library(MASS)
N = 2000
set.seed(1)
# Works well when the variance of the normal term is not overly large
# When the variance is very large, it tends to be underestimated
x = rbinom(N, 1, 0.5)
z = rnorm(N)
y = rpois(N, exp(-1 + x + z + 0.5 * rnorm(N)))
est = pln(y~x+z)
print(est$estimates, digits=3)
Run the code above in your browser using DataLab