MT5762 Lecture 15

A simple start

• start with just a discrete random variable \( Y \), so the probability mass function gives probabilities directly for a value of \( Y \)

• Consider 2 independent identically distributed random variables \( Y_1 \) and \( Y_2 \) which are Poisson distributed, and which we can observe

• If we know the mean \( \theta \) of the Poisson from which they have come, then we can write down the probability of getting observed values \( y_1 \) and \( y_2 \):

\[ P ( y_1, y_2 | \theta ) = \frac{\theta ^ {y_1} e^{- \theta} }{y_1 !} \times \frac{\theta ^ {y_2} e^{- \theta} }{y_2 !} \]

• suppose that we have the two observations, \( y_1 \) and \( y_2 \), but don't know the parameter \( \theta \).

• In some sense the most 'likely' value for \( \theta \) is the one which maximises the previous expression, i.e. the value of \( \theta \) which makes our observed data look most probable

• When expressions like previous contain known data, but unknown parameters (like \( \theta \)), they are known as likelihoods for the unknown parameters

• Estimating a parameter by maximising its likelihood (given some observed data) is known as maximum likelihood estimation

• set up the likelihood for your observed data assuming a general form of Probability Function \( f(X|\boldsymbol{\theta}) \) (e.g Poisson, Binomial etc).

• Note \( \boldsymbol{\theta} \) can be a vector of several parameters \( \boldsymbol{\theta}=[ \theta_1, \theta_2, ...] \). Let us ignore \( X \) for the time being, so use \( f(X|\boldsymbol{\theta})=f(\boldsymbol{\theta}) \)

• multiply together the functions for calculating probabilities/densities for each potential data values i.e. \( \prod f(\boldsymbol{\theta}) \)

• take logs of this to simplify (i.e. now have \( g(\boldsymbol{\theta})=\sum \log(f(\boldsymbol{\theta})) \) then differentiate with respect to the parameter of interest (\( \theta_1 \) say)

• solve for the value of \( \theta_1 \) \[ \frac{d(g(\boldsymbol{\theta}))}{d\theta_1}=0 \] e.g. we find the parameter value of \( \theta_1 \) where the log-likelihood function peaks

• To actually get a log-likelihood value for some value of \( \theta \) - insert the observed values of \( X=x_1,...x_n \) into step 2

• we are usually interested in the log likelihood at the MLE e.g. this is where we might compare competing models.

• If we are dealing with iid variables, then we can often derive the form of MLE for a general value of \( X \).

Note that the theory applies to continuous RVs just as well as discrete ones - just use PDFs to create the likelihood function.

• A very (in)famous dataset relating to count data is that of Bortkiewicz, published in 1898.
• Bortkiewicz was an Economist/Statistician who took an interest in the numbers of Prussian Cavalry men dying after being kicked in the head by their horses.
• The data was collected across several cavalry corps over a 20 year period (some 90+ fatal kicks to the head!).

• Suppose this data represents the reports from one cavalry corp over this period: 10, 6, 2, 1, 1 where the first number represents the number of years with no deaths, the second is the number of years with 1 reported death etc.
• Assuming a Poisson distribution is appropriate, calculate the Maximum Likelihood Estimate for the yearly death rate parameter \( \lambda \).

• Our data comprises of 20 points i.e. \( {0,0,\cdots,1,1,\cdots,2,2,3,4} \).
• There are several ways to achieve this that differ in the details, but the process remains the same.
• First construct the likelihood function as a product of the PDF (the standard Poisson) for each of the \( n \) data points:

\[ \prod_{i=1}^n \frac{\lambda^{x_i}e^{-\lambda}}{x_i!} \]

Taking logs here will simplify matters

\[ \begin{align*} \sum_{i=1}^n \log \left( \frac{\lambda^{x_i}e^{-\lambda}}{x_i!}\right)&\quad = \quad \sum_{i=1}^n\left(\log(\lambda^{x_i})-\lambda-\log(x_i!)\right)\\ &\quad = \quad \sum_{i=1}^nx_i\log(\lambda)-n\lambda-\sum_{i=1}^n\log(x_i!)\\ \end{align*} \]

differentiate wrt parameter of interest then solve for said parameter

\[ \begin{align*} \sum_{i=1}^n\frac{x_i}{\lambda}-n=0&\quad = \quad \sum_{i=1}^n\frac{x_i}{\lambda}=n\\ &\quad = \quad\frac{\sum_{i=1}^n x_i}{n}=\lambda \end{align*} \]

which is (as you perhaps expected) simply the mean of the data.

So, summing the data gives \( 0\times 10+1\times 6 + 2\times 2 + 3\times 1 + 4\times 1=17 \) and thus \( \lambda=17/20=0.85 \). Slightly less than one fatal kick to the head per year. So given the assumption of a Poisson process, the Poisson distribution that best fits our data is one with a mean of 0.85; it is the most likely value of \( \lambda \).

From our grid search, the answer is about the mean of the data. The sample mean is the theoretical MLE for \( \theta \) in this problem.

  mleindicator<- which(likelihood==max(likelihood))

  mle<-theta[mleindicator]

  mean(data)

[1] 0.85

  mle

[1] 0.85

Likelihood has a wide range of uses in statistics:

Likelihood has a wide range of uses in statistics:

• it is often used as the basis of parameter estimation in models e.g. for a regression \( Y=\beta_0+\beta_1 X \), what are the most likely values of \( \beta_0 \) and \( \beta_1 \) given a set of data \( X=x_1, x_2, ... , x_n \) assuming Normal errors? i.e. \( Y\sim N(\mu, \sigma) \) where \( \mu=E[Y|X]=\beta_0+\beta_1 X \).
• comparing competing models, which model is best supported by the data, which of the models is more likely to have generated the data we observed e.g. likelihood ratios.

Let's try looking at simple regression. We'll

• Optimise using sums-of-squared residuals i.e. minimse these
• Re-pose as (log)-likelihood problem and optimise i.e. maximise this

Here we predict RunTime as a function of Oxygen - data that relates to running and metabolic rates.

Set up a grid to search over

  #~ grid search over parameters
  slopeGrid<- seq(-1,1, length=200)

  intGrid<- seq(10,50, length=200)

  gridSearch<- expand.grid(intGrid, slopeGrid)

Then set up the search - this is the thing I minimise (\( y-\hat{y} \) = \( y - (\beta_0 + \beta_1 x) \))

  errorSSfunction<- function(x, y, intercept, slope){sum((y-(intercept+slope*x))^2)}

This thing just passes all our intercept and slope parameters through that function:

  surfaceVals<- apply(gridSearch, 1, function(q){
    errorSSfunction(fitnessData$Oxygen, fitnessData$RunTime, q[1], q[2])})

So, which of the parameter pairs had the smallest residual sum of squares (RSS)?

  minPoint<- gridSearch[which(surfaceVals==min(surfaceVals)),]

  minPoint

          Var1       Var2
15457 21.25628 -0.2261307