1
$\begingroup$

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:

enter image description here

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.

$\endgroup$
6
  • 7
    $\begingroup$ Since the MLE for the Poisson rate parameter is the mean of the observed data... why not just mean(a) ? $\endgroup$
    – duckmayr
    Feb 1, 2019 at 12:46
  • 3
    $\begingroup$ 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 ...) $\endgroup$
    – Ben Bolker
    Feb 1, 2019 at 12:50
  • $\begingroup$ @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 $\endgroup$
    – Feng Chen
    Feb 1, 2019 at 13:08
  • $\begingroup$ @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. $\endgroup$
    – Feng Chen
    Feb 1, 2019 at 13:09
  • 1
    $\begingroup$ @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. $\endgroup$
    – duckmayr
    Feb 1, 2019 at 13:17

1 Answer 1

11
$\begingroup$

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)))
$\endgroup$
2
  • 1
    $\begingroup$ This is what I meant! So I can use this way to deal with complex parameter estimate. Thanks $\endgroup$
    – Feng Chen
    Feb 1, 2019 at 13:25
  • 4
    $\begingroup$ No glm(a~1,family=poisson(link="identity"))$coefficients or exp(glm(a~1,family=poisson)$coefficients) ? $\endgroup$
    – Glen_b
    Feb 2, 2019 at 2:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.