y <- c(rnorm(1000), rnorm(2000, 1, 0.5))
es una mezcla de dos normales (N(0, 1) y N(1, 0.5)) con pesos 1/3 y 2/3 respectivamente. Pero, ¿cómo podríamos estimar los parámetros a partir de esos datos?
Se puede usar, p.e., flexmix, que implementa eso del EM. Pero en el librillo de este maestrillo dice
library(rstan)
y <- c(rnorm(1000), rnorm(2000, 1, 0.5))
codigo <- "
data {
int<lower=1> K; // number of mixture components
int<lower=1> N; // number of data points
real y[N]; // observations
}
parameters {
simplex[K] theta; // mixing proportions
real mu[K]; // locations of mixture components
real<lower=0> sigma[K]; // scales of mixture components
}
model {
real ps[K]; // temp for log component densities
sigma ~ cauchy(0,2.5);
mu ~ normal(0,10);
for (n in 1:N) {
for (k in 1:K) {
ps[k] <- log(theta[k]) + normal_log(y[n],mu[k], sigma[k]);
}
increment_log_prob(log_sum_exp(ps));
}
}"
fit <- stan(model_code = codigo,
data = list(K = 2, N = length(y), y = y),
iter=48000, warmup=2000,
chains=1, thin=10)
En el código anterior no sé si queda claro cómo cada punto $y_i$ sigue una distribución (condicionada a los parámetros) con densidad $\theta_1 \phi(y_i, \mu_1, \sigma_1) + \theta_2 \phi(y_i, \mu_2, \sigma_2)$.
