This brief lesson introduces the concept of discrete distribution (probability mass function) on the example of Poisson distribution. We will also do more likelihood stuff.

# 1 Poisson distribution

Poisson distribution is a simple model for discrete and positive count data.

$p\left( x \right) = \frac{{e^{ - \lambda } \lambda ^x }}{{x!}}$

Where $$\lambda$$ is the only parameter, which is also the mean of the distribution.

# 2 Simulating poisson data

Let’s explore the Poisson model a bit and simulate some data using rpois():

rpois(n=100, lambda=1)
##   [1] 0 1 0 2 0 0 0 0 0 1 0 0 3 0 1 1 1 2 2 0 2 0 4 2 1 0 0 0 0 0 0 1 1 1 2
##  [36] 0 0 2 0 1 1 3 0 0 1 2 0 1 1 0 1 3 0 2 0 1 0 1 0 1 0 1 1 0 1 1 1 2 1 2
##  [71] 0 2 0 0 2 2 1 2 1 0 0 2 1 2 1 1 0 2 1 1 1 1 0 0 2 0 1 0 0 2

Now try different values of $$\lambda$$.

We can plot the “empirical frequency histogram” for simulated data with 3 different values of parameter $$\lambda$$:

  plot(c(0,35), c(0,1), type="n",
xlab="x", ylab="Frequency")

hist(rpois(n=100, lambda=20), add=T, freq=FALSE, col="green")

Now we can plot the respective probability mass functions:

  plot(0:35, dpois(0:35, lambda=1),
col="blue", pch=19,
ylab="Probability mass",
xlab="x")
points(0:35, dpois(0:35, lambda=5), col="red", pch=19)
points(0:35, dpois(0:35, lambda=20), col="green", pch=19)

# 3 Poisson likelihood

I will generate the artificial data again:

x <- rpois(n=100, lambda=8)
x
##   [1]  7  8  7  5 18  7  7  4  6  7  9  8  7 11 15  7  7  4  6  5 12 12  9
##  [24]  3  6  9  5 12  6  8 11  7  8 10  9 11  9 11  4  7  7  7  3  7  5 14
##  [47] 11 11 11  6  2 11  9 11 10  9  7  2 11  6  8  7  7  6 10  5 10  7  5
##  [70]  9 14  8 14  9  7  7  9  3 10  9 12 11  8  5  3  5  9  8  4 11  7  8
##  [93]  9  5  9  4  8  5  6  9

Poisson distribution is didactically convenient as it only has the one parameter $$\lambda$$. We can write the deviance (negative log-likelihood) function for a given dataset $$x$$:

  deviance.function <- function(x, lambda)
{
LL <- dpois(x, lambda, log=TRUE) # the log likelihood
deviance <- -sum(LL)
return(deviance)
}

We can try different values of $$\lambda$$:

  deviance.function(x, lambda=10)
## [1] 270.6436

And we can evaluate the deviance.function for values of $$\lambda$$ between 0 and 35:

  lambdas <- 0:35
deviances <- numeric(36) # empty numeric vector

for(i in 1:36)
{
deviances[i] <- deviance.function(x, lambdas[i])
}

Let’s check what we have:

  deviances
##  [1]       Inf 1198.8962  748.5373  526.5980  398.1785  321.0025  276.2392
##  [8]  253.8435  247.8196  254.2999  270.6436  294.9673  325.8803  362.3264
## [15]  403.4847  448.7043  497.4607  549.3248  603.9410  661.0116  720.2848
## [22]  781.5454  844.6085  909.3138  975.5215 1043.1088 1111.9675 1182.0017
## [29] 1253.1258 1325.2633 1398.3455 1472.3104 1547.1019 1622.6692 1698.9659
## [36] 1775.9498

Our maximum likelihood exstimate (MLE) of $$\lambda$$ is:

MLE <- lambdas[which.min(deviances)]
MLE
## [1] 8

And this is the plot of the deviance function, together with the MLE (vertical line):

plot(lambdas, deviances, ylab="Deviance (Negative Log Likelihood)")
abline(v=MLE)