The truncated Poisson

February 21, 2010

(This article was first published on Stats raving mad » R, and kindly contributed to R-bloggers)

A common model for counts data is the Poisson. There are cases however that we only record positive counts, ie there is a truncation of 0. This is the truncated Poisson model.

To study this model we only need the total counts and the sample size. This comes from the sufficient statistic principle as the likelihood is

 logL=-n-frac{e^{-lambda } n}{1-e^{-lambda }}+frac{T}{lambda },

where  T=sum _X. Let’s set  T=160 and   n=50.


loglik <- function(theta, n, sum.x) {
- n* log(exp(theta) - 1) + sum.x * log(theta)

The gradient is

 frac{e^{-2 lambda } n}{left(1-e^{-lambda }right)^2}+frac{e^{-lambda } n}{1-e^{-lambda }}-frac{T}{lambda ^2}.

gradlik <- function(theta, n, sum.x){

The second derivative (the hessian) is

 -frac{2 e^{-2 lambda } n}{left(1-e^{-lambda }right)^2}-frac{e^{-lambda } n}{1-e^{-lambda }}-e^{-lambda } left(frac{2 e^{-2 lambda }}{left(1-e^{-lambda }right)^3}+frac{e^{-lambda }}{left(1-e^{-lambda }right)^2}right) n+frac{2 T}{lambda ^3}.

hesslik <- function(theta, n, sum.x) {

Now, we simply call the s/t function of maxLik package to fit the model.

a <- maxLik(loglik, start=3.2,method="Newton-Raphson",
n=n,sum.x=sum.x, print.level=2)
# ----- Initial parameters: -----
# fcn value: 28.18494
#      parameter initial gradient free
# [1,]       3.2        -2.124718    1
# Condition number of the (active) hessian: 1
# -----Iteration 1 -----
# -----Iteration 2 -----
# -----Iteration 3 -----
# --------------
# gradient close to zero. May be a solution
# 3  iterations
# estimate: 3.048175
# Function value: 28.34853

There is an nice & handy function I found somewhere in the net (sorry I don’t remember the author;() to plot the likelihood of the truncated Poisson.

plot.logL = function(Left, Right, Step){
 grid.x = seq(Left, Right, by = Step)
 plot(grid.x, loglik(grid.x, n, sum.x), type = "l", ,col="red",
 xlab = "theta", ylab = "logL",
 main = "Log likelihood for truncated Poisson")
plot.logL(2, 4, 0.00001);grid()

To leave a comment for the author, please follow the link and comment on their blog: Stats raving mad » R. offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , , , , , ,

Comments are closed.


Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)