Генерация случайных величин из смеси нормальных распределений

20

Как я могу сделать выборку из распределения смеси, и в частности из смеси нормальных распределений в R? Например, если я хотел сделать выборку из:

0.3×N(0,1)+0.5×N(10,1)+0.2×N(3,.1)

как я мог это сделать?

Gung - Восстановить Монику
источник
3
Мне действительно не нравится этот способ обозначения смеси. Я знаю, что это обычно делается так, но я нахожу это вводящим в заблуждение. Нотация предполагает, что для выборки вам необходимо выбрать все три нормали и взвесить результаты по тем коэффициентам, которые, очевидно, были бы неверными. Кто-нибудь знает лучшую запись?
StijnDeVuyst
У меня никогда не было такого впечатления. Я думаю о распределениях (в данном случае о трех нормальных распределениях) как о функциях, а затем результатом является другая функция.
roundsquare
@StijnDeVuyst вы можете захотеть посетить этот вопрос возник из вашего комментария: stats.stackexchange.com/questions/431171/…
ankii
@ankii: спасибо за указание на это!
StijnDeVuyst

Ответы:

32

Рекомендуется избегать forциклов Rпо причинам производительности. Альтернативное решение, которое использует этот факт rnorm, векторизовано:

N <- 100000

components <- sample(1:3,prob=c(0.3,0.5,0.2),size=N,replace=TRUE)
mus <- c(0,10,3)
sds <- sqrt(c(1,1,0.1))

samples <- rnorm(n=N,mean=mus[components],sd=sds[components])
М. Берк
источник
3
Кроме того, вы можете использовать свойства нормального распределения, чтобы заменить последнюю строку на samples <- rnorm(N)*sds[components]+mus[components]. Мне легче читать :)
Элвис
Очень элегантно (см @ Элвис)!
Итамар
18

В общем, один из самых простых способов отбора проб из распределения смеси заключается в следующем:

Шаги алгоритма

UUniform(0,1)

U[i=1kpk,i=1k+1pk+1)pkkthkth

3) Повторите шаги 1) и 2), пока не получите желаемое количество образцов из распределения смеси

Теперь, используя общий алгоритм, приведенный выше, вы можете взять из вашего примера смесь нормалей, используя следующий Rкод:

#The number of samples from the mixture distribution
N = 100000                 

#Sample N random uniforms U
U =runif(N)

#Variable to store the samples from the mixture distribution                                             
rand.samples = rep(NA,N)

#Sampling from the mixture
for(i in 1:N){
    if(U[i]<.3){
        rand.samples[i] = rnorm(1,0,1)
    }else if(U[i]<.8){
        rand.samples[i] = rnorm(1,10,1)
    }else{
        rand.samples[i] = rnorm(1,3,.1)
    }
}

#Density plot of the random samples
plot(density(rand.samples),main="Density Estimate of the Mixture Model")

#Plotting the true density as a sanity check
x = seq(-20,20,.1)
truth = .3*dnorm(x,0,1) + .5*dnorm(x,10,1) + .2*dnorm(x,3,.1)
plot(density(rand.samples),main="Density Estimate of the Mixture Model",ylim=c(0,.2),lwd=2)
lines(x,truth,col="red",lwd=2)

legend("topleft",c("True Density","Estimated Density"),col=c("red","black"),lwd=2)

Который генерирует:

введите описание изображения здесь

и как проверка работоспособности:

введите описание изображения здесь


источник
Здравствуй! Спасибо! Этот ответ мне очень помог. Я использую это в исследовательском проекте. Я хочу процитировать ссылку на выше. Можете ли вы предложить цитату из научной статьи.
Абхишек Бхатия
7

kR

set.seed(8)               # this makes the example reproducible
N     = 1000              # this is how many data you want
probs = c(.3,.8)          # these are *cumulative* probabilities; since they 
                          #   necessarily sum to 1, the last would be redundant
dists = runif(N)          # here I'm generating random variates from a uniform
                          #   to select the relevant distribution

# this is where the actual data are generated, it's just some if->then
#   statements, followed by the normal distributions you were interested in
data = vector(length=N)
for(i in 1:N){
  if(dists[i]<probs[1]){
    data[i] = rnorm(1, mean=0, sd=1)
  } else if(dists[i]<probs[2]){
    data[i] = rnorm(1, mean=10, sd=1)
  } else {
    data[i] = rnorm(1, mean=3, sd=.1)
  }
}

# here are a couple of ways of looking at the results
summary(data)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# -3.2820  0.8443  3.1910  5.5350 10.0700 13.1600 

plot(density(data))

введите описание изображения здесь

Gung - Восстановить Монику
источник
Хороший ответ, вы побили меня до публикации: P
1
Спасибо за совет, @BabakP. Я не уверен, что это было. Это было что-то в ifelse()заявлении, но я должен это выяснить позже. Я заменил этот код с циклом.
gung - Восстановить Монику
6
RfindInterval()cumsum()μmuσ2spmix <- function(n,mu,s,p) { ii <- findInterval(runif(n),cumsum(p))+1; x <- rnorm(n,mean=mu[ii],sd=sqrt(s[ii])); return(x); }
1
@Macro, очень верный и очень хороший код! Я не видел эту findInterval()команду раньше, однако мне нравится писать код здесь максимально упрощенно, потому что я хочу, чтобы она была инструментом для понимания, а не эффективности.
1
Я сказал, что это были хорошие ответы. Моя цель состояла не в том, чтобы критиковать вас, а в том, чтобы предложить подход, который легко обобщает более трех измерений, изменяя только один аргумент, а не какой-либо код. Мне не ясно, почему то, что вы написали, более прозрачно, чем то, что я написал, но я, конечно, не хочу спорить об этом. Приветствия.
Макрос
0

Уже даны отличные ответы, поэтому для тех, кто хочет достичь этого в Python, вот мое решение:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

mu = [0, 10, 3]
sigma = [1, 1, 1]
p_i = [0.3, 0.5, 0.2]
n = 10000

x = []
for i in range(n):
    z_i = np.argmax(np.random.multinomial(1, p_i))
    x_i = np.random.normal(mu[z_i], sigma[z_i])
    x.append(x_i)

def univariate_normal(x, mean, variance):
    """pdf of the univariate normal distribution."""
    return ((1. / np.sqrt(2 * np.pi * variance)) * 
            np.exp(-(x - mean)**2 / (2 * variance)))

a = np.arange(-7, 18, 0.01)
y = p_i[0] * univariate_normal(a, mean=mu[0], variance=sigma[0]**2) + p_i[1] * univariate_normal(a, mean=mu[1], variance=sigma[0]**2)+ p_i[2] * univariate_normal(a, mean=mu[2], variance=sigma[0]**2)

fig, ax = plt.subplots(figsize=(8, 4))

ax.hist(x, bins=100, density=True)
ax.plot(a, y)

введите описание изображения здесь

КРЫСА
источник