# how to calculate the parameter lambda in Poisson distribution?

Let's say there is a sequence:

a <- c(1,2,3,1,2,1,1,3,1,2,3,5)


This conforms to a Poisson distribution, the formula of which is shown as:

Now I want to calculate the parameter lambda of Poisson. Usually we need to use maximum likelihood estimation to do this. But how to do it in R?

This is my first time to do statistical analysis in R, so please provide as many as details as possible.

• Since the MLE for the Poisson rate parameter is the mean of the observed data... why not just mean(a) ? Feb 1, 2019 at 12:46
• Can you give us a little more context? Is this a homework question? (Not disallowed, but in any case you should tell us what you've tried so far ...) Feb 1, 2019 at 12:50
• @Ben Bolker: I am trying to study statistics and R. This is my first try. You know, knowing how to use mle is very important for parameter estimation Feb 1, 2019 at 13:08
• @duckmayr: Thanks for your advice. I know this is work. But I want to know how to use mle to do parameter estimation, So I can use it some where else in the future. Poisson distribution is just my starting point. Feb 1, 2019 at 13:09
• @FengChen I believe then I do not understand your objective. What I described is "using mle to do parameter estimation" -- MLE stands for maximum likelihood estimate. The mean of the observed data is an estimate of the rate parameter -- specifically the estimate that maximizes the likelihood of the data. Feb 1, 2019 at 13:17

## the easy way

Since we know the mean is the MLE for $$\lambda$$:

mean(a)  ## 2.0833


## fitting distributions

MASS::fitdistr is a built-in method for ML estimation of the parameters of a variety of distributions.

MASS::fitdistr(a,"Poisson")


## brute force: optim

Define a function that returns the negative log-likelihood for a given value of $$\lambda$$:

f <- function(lambda) {
-sum(dpois(a,lambda=lambda,log=TRUE))
}
optim(par=1, ## starting value
fn=f,
method="Brent",   ## need to specify for 1-D optimization
lower=0.001, upper=10)


## mle2

The bbmle::mle2() function (a variant of stats4::mle()) does the same optimization, but has more features for doing things with the results (e.g. computing likelihood profiles, comparing models via Likelihood Ratio Test). (mle2 uses BFGS instead of Nelder-Mead optimization by default, which works in 1-D, so we don't need the method="Brent" from above [we could use it if we wanted].)

library(bbmle)
mle2(minuslogl=f,start=list(lambda=1))


## mle2 with formula

mle2 also allows some shortcuts:

mle2(a ~ dpois(lambda),
data=data.frame(a),
start=list(lambda=1))


## glm

As @glen_b points out in comments, this is also a special case of a generalized linear model. Since the Poisson model uses a log link by default, we have to be a little bit careful.

coef(glm(a~1,family=poisson(link="identity")))
exp(coef(glm(a~1,family=poisson)))

• This is what I meant! So I can use this way to deal with complex parameter estimate. Thanks Feb 1, 2019 at 13:25
• No glm(a~1,family=poisson(link="identity"))$coefficients or exp(glm(a~1,family=poisson)$coefficients) ? Feb 2, 2019 at 2:39