MLE for the Poisson distribution

Suppose that we have the following independent observations and we know that they come from the same probability density function

k<-c(32,38,42,33,27) #our observations

library('ggplot2')
dat<-data.frame(k,y=0) #we plotted our observations in the x-axis
p<-ggplot(data=dat,aes(x=k,y=y))+geom_point(size=4)
p

We don’t know the exact probability density function, but we know that the shape is the shape of the poisson distribution \[f(x|\theta)=f(k|\lambda)=\frac{\lambda^k}{k!}e^{-\lambda}\] Let’s plot two arbitrary probability density functions fixing the parameter of the parametric family

kseq<-seq(20,45,1)
f1<-dpois(kseq,33) #lambda=33
f2<-dpois(kseq,27) #lambda=27
curve1<-data.frame(kseq,f1)
curve2<-data.frame(kseq,f2)

p<-p+geom_segment(data=curve1,aes(x=kseq,xend=kseq,y=0,yend=f1))+
     geom_segment(data=curve2,aes(x=kseq+.2,xend=kseq+.2,y=0,yend=f2)) #we introduce a small jitter in the distribution to avoid overplotting
p

It is clear that the probability density function \(f\) for \(\lambda=33\) describes the observations better than the other one. We are interested in the \(f\) that maximizes \(L\).

The joint density function is \[f(k|\lambda)=f(k_1|\lambda)f(k_2|\lambda)...f(k_5|\lambda)=\] \[ =\frac{\lambda^{k_1}}{k1!}e^{-\lambda}\frac{\lambda^{k_2}}{k2!}e^{-\lambda}...\frac{\lambda^{k_5}}{k5!}e^{-\lambda}=\]

that when considered as a function of the parameter is \[L(\lambda|k)=L(\lambda|32,38,42,33,27)=f(k|\lambda)=f(32,38,42,33,27|\lambda)=\] \[=f(32|\lambda)f(38|\lambda)...f(27|\lambda)=\] \[=\frac{\lambda^{32}}{32!}e^{-\lambda}\frac{\lambda^{38}}{38!}e^{-\lambda}...\frac{\lambda^{27}}{27!}e^{-\lambda}\]

logL<-function(p) sum(dpois(k,p,log=T)) # p=lambda

We can calculate the \(log(L)\) for the two previous examples to verify that \(log(L)\) is larger for \(\lambda=33\)

logL(33)

## [1] -15.51074

logL(27)

## [1] -20.0261

logL(33)>logL(27)

## [1] TRUE

Let’s plot the \(log(L)\) for some values of \(\lambda\)

library('plyr')
lambdaSeq<-seq(20,45,.1)
logLike<-laply(lambdaSeq,function(z) logL(z))
dLogL<-data.frame(lambdaSeq,logLike)
ggplot(data=dLogL,aes(x=lambdaSeq,y=logLike))+geom_line()

So it seems that values of \(\lambda\) around 34 maximizes log(L). Let’s maximize it properly using optimize.

estPar<-optimize(logL,c(10,50),maximum=T) #with optimize we can specify that we search for a max  
MLEparameter<-estPar$maximum
MLEparameter

## [1] 34.4

So, we found that from the parametric family, the probability density function that better characterizes the observations according to MLE is the one described by the parameter \(\lambda=\) 34.3999994. Let’s plot the distribution in red in the previous graph

fMLE<-dpois(kseq,MLEparameter) 
curveMLE<-data.frame(kseq,fMLE)

p<-p+geom_segment(data=curveMLE,aes(x=kseq-.2,xend=kseq-.2,y=0,yend=fMLE),color='red') #we introduce a small jitter in the distribution to avoid overplotting
p

It could be demonstrated that, actually when the parametric family is the poisson distribution, then the MLE of \(\lambda\) is the mean of the observations

mean(k)

## [1] 34.4