meanlog = 0.5; sdlog = 0.5; y = 0:19
proby = dpolono(y, m=meanlog, sd=sdlog)
sum(proby) # Should be 1
opar = par(no.readonly = TRUE)
par(mfrow=c(2,2))
plot(y, proby, type="h", col="blue", ylab="P[Y=y]", log="",
main=paste("Poisson lognormal(meanlog=",meanlog,", sdlog=",sdlog,")",
sep=""))
# More extreme values; use the approximation and plot on a log scale
# Notice the kink at bigx.
y = 0:190
proby = dpolono(y, m=meanlog, sd=sdlog, bigx=100)
sum(proby) # Should be 1
plot(y, proby, type="h", col="blue", ylab="P[Y=y]", log="y",
main=paste("Poisson lognormal(meanlog=",meanlog,", sdlog=",sdlog,")"))
# Random number generation
table(y <- rpolono(n=1000, m=meanlog, sd=sdlog))
hist(y, breaks=((-1):max(y))+0.5, prob=TRUE, border="blue")
par(opar)
Run the code above in your browser using DataLab