Generando variables aleatorias a partir de una mezcla de distribuciones normales

20

¿Cómo puedo tomar muestras de una distribución de mezcla, y en particular una mezcla de distribuciones normales en R? Por ejemplo, si quisiera probar de:

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

¿Cómo podría hacer eso?

gung - Restablece a Monica
fuente
3
Realmente no me gusta esta forma de denotar una mezcla. Sé que se hace convencionalmente de esta manera, pero me resulta engañoso. La notación sugiere que para muestrear, debe muestrear las tres normales y sopesar los resultados por esos coeficientes que obviamente no serían correctos. Alguien sabe una notación mejor?
StijnDeVuyst
Nunca tuve esa impresión. Pienso en las distribuciones (en este caso las tres distribuciones normales) como funciones y luego el resultado es otra función.
roundsquare
@StijnDeVuyst es posible que desee visitar esta pregunta originada en su comentario: stats.stackexchange.com/questions/431171/…
ankii
@ankii: ¡gracias por señalarlo!
StijnDeVuyst

Respuestas:

32

Es una buena práctica evitar forbucles Rpor razones de rendimiento. Una solución alternativa que explota el hecho rnormes vectorizada:

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])
M. Berk
fuente
3
Alternativamente, puede usar las propiedades de la distribución normal para reemplazar la última línea por samples <- rnorm(N)*sds[components]+mus[components]. Me resulta más fácil de leer :)
Elvis
Muy elegante (cc @ Elvis)!
Itamar
18

En general, una de las formas más fáciles de tomar muestras de una distribución de mezcla es la siguiente:

Pasos de algoritmo

1) Generar una variable aleatoria UUniforme(0 0,1)

U[yo=1kpagk,yo=1k+1pagk+1)pagkkthkth

3) Repita los pasos 1) y 2) hasta obtener la cantidad deseada de muestras de la distribución de la mezcla.

Ahora, usando el algoritmo general dado anteriormente, puede tomar muestras de su mezcla de ejemplo de normales usando el siguiente Rcódigo:

#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)

Lo que genera:

ingrese la descripción de la imagen aquí

y como control de cordura:

ingrese la descripción de la imagen aquí


fuente
¡Hola! ¡Muchas gracias! Esta respuesta me ayudó mucho. Estoy usando esto en un proyecto de investigación. Deseo citar una referencia para lo anterior. ¿Puede sugerir una cita de artículo de investigación?
Abhishek Bhatia
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))

ingrese la descripción de la imagen aquí

gung - Restablece a Monica
fuente
Buena respuesta, me
1
Gracias por el consejo, @BabakP. No estoy seguro de qué era. Era algo en la ifelse()declaración, pero tendré que resolverlo más tarde. Reemplacé ese código con un bucle.
gung - Restablece a Monica
66
(cc @BabakP) Estas son buenas respuestas y obviamente son correctas (+ 1s). Solo un Rtruco de programación: también puede usar los comandos findInterval()y cumsum()para simplificar el código y, lo que es más importante, facilitar la generalización a un número diferente de dimensiones. Por ejemplo, para un vector de entrada de mediasμ( mu) y variacionesσ2( s), y las probabilidades de la mezcla ( p), una función simple para generar n muestras a partir de esta mezcla seríamix <- function(n,mu,s,p) { ii <- findInterval(runif(n),cumsum(p))+1; x <- rnorm(n,mean=mu[ii],sd=sqrt(s[ii])); return(x); }
Macro
1
@Macro, muy cierto y muy lindo código! No he visto el findInterval()comando antes, sin embargo, me gusta escribir código aquí de la manera más simple posible porque quiero que sea una herramienta de comprensión en lugar de eficiencia.
1
Dije que estas eran buenas respuestas. Mi propósito no era criticarte, sino ofrecerte un enfoque que se generalice fácilmente a más de tres dimensiones cambiando solo un argumento, no ningún código. No me queda claro por qué lo que has escrito es más transparente que lo que he escrito, pero seguramente no quiero discutir sobre eso. Salud.
Macro
0

Ya tengo respuestas perfectas, así que para aquellos que quieren lograr esto en Python, aquí está mi solución:

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)

ingrese la descripción de la imagen aquí

UNA RATA
fuente