mu1 <- 99; mu2 <- 150; nn <- 1000
sd1 <- sd2 <- exp(3)
(phi <- logit(-1, inverse = TRUE))
mdata <- data.frame(y = ifelse(runif(nn) < phi, rnorm(nn, mu1, sd1),
rnorm(nn, mu2, sd2)))
fit <- vglm(y ~ 1, mix2normal1(equalsd = TRUE), mdata)
# Compare the results
cfit <- coef(fit)
round(rbind('Estimated' = c(logit(cfit[1], inverse = TRUE),
cfit[2], exp(cfit[3]), cfit[4]),
'Truth' = c(phi, mu1, sd1, mu2)), dig = 2)
# Plot the results
xx <- with(mdata, seq(min(y), max(y), len = 200))
plot(xx, (1-phi)*dnorm(xx, mu2, sd2), type = "l", xlab = "y",
main = "Orange = 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 = "orange")
lines(xx, (1-phi.est) * dnorm(xx, Coef(fit)[4], sd.est), col = "orange")
abline(v = Coef(fit)[c(2,4)], lty = 2, col = "orange")
abline(v = c(mu1, mu2), lty = 2, col = "blue")
Run the code above in your browser using DataLab