Obtención de un número aleatorio a partir de una mezcla de dos distribuciones normales

Question

Quiero muestrear de una distribución normal mixta, la primera es $N(1,2)$ El segundo es $N(5,4)$ . He utilizado rnorm(100, c(mean=c(1,5), sd=c(2,4))) . ¿Es correcto?

El problema que intento resolver es el muestreo a partir de las 2 distribuciones anteriores, la primera con un 75% y la segunda con un 25%. ¿Estoy en el camino correcto?

Edita: Reescribiré el problema para despejar, con números más fáciles. :)

Quiero tomar una muestra de $N(0,1)$ con un 70% de probabilidad, y $N(100,10)$ con un 30% de probabilidad. Por supuesto, esto es sólo por el bien de la discusión, la distribución real con la que estoy trabajando es n(21, 3,3), n(26,4).

Answer 1

Para lograr el objetivo de tomar muestras de una mezcla desigual de distribuciones, lo más sencillo es tomar muestras por separado, en proporción a la relación deseada:

 p <- 0.70 #P(from N(mu1, sd1)) 
 n.samps <- 10000
 mu1 <- 0
 sd1 <- 1
 mu2 <- 100
 sd2 <- 10

 x <- vector()
 for(i in 1:n.samps){
    b <- runif(1, 0, 1)
    if(b < p){
       x[i] <- rnorm(1, mu1, sd1)
     } else { 
       x[i] <- rnorm(1, mu2, sd2)
     }
   }

esto puede hacerse ~50 veces más rápido:

 binary <- runif(n.samps, 0, 1) > p
 x <- c(rnorm(sum(binary), 1, 2), rnorm(sum(!binary), 5, 4)

a continuación, para extraer una muestra:

sample(x, 1)

o para reorganizar:

x <- sample(x, n.samp)

Answer 2

Si desea tomar muestras desiguales (con probabilidad 0,7 y 0,3) de dos gaussianos con parámetros $(\mu_1,\sigma_1^2)$ y $(\mu_2,\sigma_2^2)$ entonces probablemente puedas intentar algo así:

n <- 100
yn <- rbinom(n, 1, .7)
# draw n units from a mixture of N(0,1) and N(100,3^2)
s <- rnorm(n, 0 + 100*yn, 1 + 2*yn)

De hecho, esta es una de las ilustraciones proporcionadas en Estadística Aplicada Moderna con S de Venables y Ripley (Springer, 2002; §5.2, pp. 110-111).

Con diferentes parámetros, puede utilizar un ifelse para seleccionar la media y la DE según la secuencia binomial dada en yn por ejemplo rnorm(n, mean=ifelse(yn, 21, 26), sd=ifelse(yn, 3.3, 4)) . (No es necesario lanzar yn a una lógica con as.logical .)