n = 1000
mu1 = 99
mu2 = 150
sd1 = sd2 = exp(3)
(phi = logit(-1, inverse=TRUE))
y = ifelse(runif(n) < phi, rnorm(n, mu1, sd1), rnorm(n, mu2, sd2))
fit = vglm(y ~ 1, mix2normal1(equalsd=TRUE))
# Compare the results
cf = coef(fit)
round(rbind('Estimated'=c(logit(cf[1], inv=TRUE),
cf[2], exp(cf[3]), cf[4]), 'Truth'=c(phi, mu1, sd1, mu2)), dig=2)
# Plot the results
xx = seq(min(y), max(y), len=200)
plot(xx, (1-phi)*dnorm(xx, mu2, sd2), type="l", xlab="y",
main="Red=estimate, blue=truth", col="blue", ylab="Density")
phi.est = logit(coef(fit)[1], inverse=TRUE)
sd.est = exp(coef(fit)[3])
lines(xx, phi*dnorm(xx, mu1, sd1), col="blue")
lines(xx, phi.est * dnorm(xx, Coef(fit)[2], sd.est), col="red")
lines(xx, (1-phi.est) * dnorm(xx, Coef(fit)[4], sd.est), col="red")
abline(v=Coef(fit)[c(2,4)], lty=2, col="red")
abline(v=c(mu1, mu2), lty=2, col="blue")
Run the code above in your browser using DataLab