lambda1 = exp(1); lambda2 = exp(3)
(phi = logit(-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, mdata, trace = TRUE)
coef(fit, matrix = TRUE)
# Compare the results with the truth
round(rbind('Estimated' = Coef(fit),
'Truth' = c(phi, lambda1, lambda2)), dig = 2)
# Plot the results
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