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=1), add=T, freq=FALSE, col="blue")
  hist(rpois(n=100, lambda=5), add=T, freq=FALSE, col="red")
  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)

Univariate Poisson Model

Petr Keil

January 2016

1 Poisson distribution

2 Simulating poisson data

3 Poisson likelihood