if (FALSE) lambda1 <- exp(1); lambda2 <- exp(3)
(phi <- logitlink(-1, inverse = TRUE))
mdata <- data.frame(y1 = rexp(nn <- 1000, lambda1))
mdata <- transform(mdata, y2 = rexp(nn, lambda2))
mdata <- transform(mdata, Y = ifelse(runif(nn) < phi, y1, y2))
fit <- vglm(Y ~ 1, mix2exp, data = mdata, trace = TRUE)
coef(fit, matrix = TRUE)
# Compare the results with the truth
round(rbind('Estimated' = Coef(fit),
'Truth' = c(phi, lambda1, lambda2)), digits = 2)
with(mdata, hist(Y, prob = TRUE, main = "Orange=estimate, blue=truth"))
abline(v = 1 / Coef(fit)[c(2, 3)], lty = 2, col = "orange", lwd = 2)
abline(v = 1 / c(lambda1, lambda2), lty = 2, col = "blue", lwd = 2)
Run the code above in your browser using DataLab