zero.inflated.Poisson.R

Jim — Oct 1, 2013, 10:46 AM

# zero-truncated Poisson modeling

y <- c(3, 4, 1, 2, 2, 1, 3, 3, 3, 4, 3, 2, 1, 1, 3, 2, 1, 2, 
       2, 1, 4, 2)

# assign lambda a prior proportional to 1/lambda

zero.poisson <- function(theta, y){
  lambda <- exp(theta)
  n <- length(y)
  sum(y * log(lambda)) - n * lambda - n * log(1 - exp(-lambda))
}

# plot posterior on a grid

lam <- seq(0, 1.4, by=0.01)
g.lam <- exp(sapply(lam, zero.poisson, y))
g.lam <- g.lam / sum(g.lam) / .01
plot(lam, g.lam, type="l")

# laplace fitting

library(LearnBayes)
fit <- laplace(zero.poisson, 1, y)

# approximation:
paste("log Lambda is distributed approximately normal with mean", 
      round(fit$mode[1], 2), "and standard deviation", 
      round(sqrt(fit$var), 2))

[1] "log Lambda is distributed approximately normal with mean 0.67 and standard deviation 0.17"

curve(dnorm(x, fit$mode, sqrt(fit$var)), add=TRUE, col="red")

legend(1, 2, 
       legend=c("Exact", "Approximation"),
       col=c("black", "red"), lty=1)