# Fit truncated Poisson distribution in R

I am trying to fit a Poisson distribution with zero-truncated, imperfectly observed data. I have a set of samples for which I can observe 1 counts and >1 counts, although I don't know what the >1 counts are. I would like to use the proportion of 1s to fit a lambda parameter so I can estimate a Poisson distribution for the >1 counts. Any advice would be very appreciated! Thank you!

Data: 1 counts: 20% of observations, >1 80% of observations. What is lambda?

• You have Binomial data (where a $1$ is a success and anything else is a "failure"). Use your favorite estimator of $p = \lambda\exp(-\lambda)/(1 - \exp(-\lambda))$.
– whuber
Commented Aug 20, 2015 at 20:17
• It is not binomial data because all the data are non-zero counts. I can observe 1's but I cannot observe the counts of anything larger than 1. Thanks! Commented Aug 21, 2015 at 2:33
• And that, by definition, is binomial: two outcomes with fixed probabilities.
– whuber
Commented Aug 21, 2015 at 12:39

This is just an extended comment as @whuber has already given the answer in a comment.

The maximum likelihood estimator for $\lambda$ is the solution to $p=\lambda e^{-\lambda}/(1-e^{-\lambda})$ or equivalently $p=\lambda/(e^{\lambda}-1)$ where $p$ is the proportion of 1's. You essentially have a Poisson distribution truncated at 1 and censored at 2.

By having just two possible data values one can just obtain a unique maximum likelihood estimator of $\lambda$ but that estimator is highly dependent on the underlying distribution actually being Poisson and with just two data values possible, there is no way to check on that assumption. So I would call this a potentially fragile estimator.

The log likelihood function is given by

$$(n-n_1)\log(1+\lambda/(1-e^{\lambda}))-n_1\lambda-n_1 \log(1-e^{-\lambda})$$

where there are $n$ samples and $n_1$ of those samples have a count of 1. Where $n_1/n=0.2$ (i.e., 20% of the counts are 1's) I get the mle of $\lambda$ as 2.660399.

Some R code to generate such data and provide the maximum likelihood estimator and an estimate of the asymptotic standard error follows:

# MLE estimation of lambda for a Poisson left-truncated at 1 and
# right-censored at 2

# Generate some samples from this odd distribution
lambda <- 3   # Poisson mean
n <- 100       # Sample size
y <- rpois(n*5,lambda) # Oversample to be more certain of getting n samples
y <- y[y >= 1][1:n]

# Determine number of 1's
n1 = sum(y==1)

# MLE initial estimate for lambda
lambda0 <- 2  # Probably any starting value greater than zero will work

# Define log of the likelihood function
logL <- function(lambda,n,n1) {
n1*(-lambda - log(1 - exp(-lambda)) + log(lambda)) +
(n - n1)*log(1 + lambda/(1 - exp(lambda)))
}

# Find maximum likelihood estimate of lambda
trPoisson <- optim(lambda0, logL, n=n, n1=n1,
method="L-BFGS-B", lower=0.00001, upper=Inf, control=list(fnscale=-1),
hessian=TRUE)

# Show results
cat(trPoisson$par,"= MLE of lambda\n") cat((-1/trPoisson$hessian)^0.5,"= Asymptotic std. err. of MLE of lambda\n")