Maximum likelihood

Slide 1

• Nonlinear Models
• Maximum Likelihood Primer

Slide 2

What is Maximum Likelihood-I?

• Maximum likelihood is a method for fitting a mathematical model to some data. “Fitting” involves estimating the values for the parameters of the model. This is what is happening under the hood when you fit a glm, lmer, glmer, etc.
• For some problems, you may not have canned functions, but this has become increasingly uncommon.

Slide 3: What is Maximum Likelihood-II?

• - Given a set of data, and a probability model, maximum
• likelihood chooses the values for the parameters that were most likely to have produced the data

• Probability: knowing parameters -> predicting data.
• Likelihood: knowing data -> parameter estimation.

L(p|x) = Prob(x|p)

Slide 4: What is Maximum Likelihood-II?

Parameters, data and error distributions:

Data

Parameter

Error

What are the data, parameters, and incorporation of error?

La = ℓ∞(1 − e−κa) + ε

ε ~ N(0, σ2)

Slide 5: What is Maximum Likelihood-III?

• Which of these curves
• are ”most likely” to have
• generated the data?

• These two seem pretty
• similar – is one better?

5

Slide 6: Basics of Likelihood

The likelihood of the data is the product of the likelihoods of each data point

Notation often inconsistent, but better to put parameters conditional on data

The set of data Xi

Product of the PDF/Likelihood across all data points

L˜(p | {Xi}) = ∏i=1n f(Xi | p)

Slide 7: Further Likelihood Considerations-I

• The deterministic (model) relationship between the input
• variables (covariates) and the (expected) data:
• The probability model for how the data pertain to the model prediction – in this case normal

To compute the likelihood for a given set of parameters we need:

Model prediction

Data

f(x) = 1 / (σ √(2π)) · exp( −½ ((x − μ)/σ)2 )

L(a) = ℓ∞(1 − e−κa)

L = 1 / (√(2π) σ) · exp( −(L(a) − La)2 / (2σ2) )

Slide 8: Further Likelihood Considerations-II

dnorm(observed, expected, sigma)

In R, the term:

is evaluated using the expression

or in R:

The full likelihood is the product over all the data points:

prod(dnorm(observed_vector, expected_vector, sigma))

Note: The function dnorm returns the value of the normal distribution PDF given x, μ, and σ.

L = 1 / (√(2π) σ) · exp( −(L(a) − La)2 / (2σ2) )

L = ∏a 1 / (√(2π) σ) · exp( −(L(a) − La)2 / (2σ2) )

Slide 9: Back to our original problem

• pred_vec <- Linfs[II]*(1-exp(-Kappa[II]*age_vec))
• likelihood <- prod(dnorm(len_vec, pred_vec, sigma))
• print(likelihood)

We can compare the four parameters choices using the following code (open R code):

• Results:
• 4.82e-149
• 2.74e-87
• 2.744-80
• 3.64e-82

• Best (so far)
• How do we know?

Slide 10: Overview

• Given a data set:
• Select a model for the underlying process.
• Select a model for the uncertainty of the data.
• Write the likelihood for one data point.
• Multiply the likelihoods across all data points.
• Find the values for the parameters that maximize the likelihood.
• Do this by minimizing negative log likelihood

Slide 11

• You monitor an intersection for cars that run the
• red light. Based on 50 days of observation, you detect
• 30 days on which “light running” took place. You repeat the “experiment” after a red light camera is installed and now observe 10 days on which “running” occurred.

• Questions
• What is the probability of “light running” before the camera was installed?
• 2. Has the probability of “light running” changed significantly?

• We’ll explore numerical methods later
• For now, let’s apply this to a simpler problem

Slide 12

• Let the probability of “running” be p (and hence the
• probability of not “running” be (1-p)). The probability of
• 30 “runs” and 20 “non-runs” is then p30(1-p)20
• (assuming “runs” and “non-runs” are independent).
• This leads to the following likelihood function from the Binomial:

We refer to the value for p that maximizes the likelihood function as the maximum likelihood estimate and denote it as p̂.

Is something missing? Why is this ok?

L(p) = p30 (1 − p)20

Slide 13

• In R:
• Prob <- seq(from=0.01,to=0.99,by=0.01)
• Like <- Probs^30*(1-Probs)^20
• plot(Prob,Like)

• We can find the MLE in R:
• Ibest <- which(Like==max(Like),arr.ind=T)
• print(Probs[Ibest])

See why normalizing constant doesn’t matter here?

Slide 14: Negative log-likelihoods

Rather than working with likelihoods we often work with negative log-likelihoods. Therefore, rather than maximizing the likelihood, we minimize the negative log-likelihood – this leads to same point estimates.

The estimate of p, p̂, is again 0.6.

Logs help with the tiny number problem, and negative works for numerical optimizers

−ln L = −30 ln p − 20 ln(1 − p)

Slide 15: A Poisson mixture problem

The blow rate for a given whale species is Poisson-distributed and it is “known” that the average number of blows in an interval for males and females is different. We monitor a population and find the following distribution of blow times.

15

• Questions
• What is the fraction of females in the population (females have
• higher blow rate)?
• Can we conclude that females and males really do differ in mean blow rate?

Slide 16

• Now, an animal is either male or female and letting the
• mean blow rate for males be λ1 and that for females be λ2,
• and assuming that the proportion of males is p:

The Poisson distribution is:

16

How would you do this with linear models?

P(x | λ) = λx exp(−λ) / x!

P(x | p, λ1, λ2) = p λ1x exp(−λ1) / x! + (1 − p) λ2x exp(−λ2) / x!

Slide 17

• In English:
• An animal is selected at random; the probability it is male is p.
• A blow rate is generated given the appropriate Poisson distribution.
• This type of model is called a mixture model. In this case there are two mixture components.

17

Slide 18

The likelihood function is then:

• This looks pretty complicated but R is there to help us:
• dpois(x,lambda) =

18

P(x | p, λ1, λ2) = ∏i ( p λ1xᵢ exp(−λ1) / xi! + (1 − p) λ2xᵢ exp(−λ2) / xi! )

P(x | λ) = λx exp(−λ) / x!

Slide 19: Start by assuming we know Propn We now have a three dimensional likelihood surface

• TheData <- matrix(scan(”Whale.dat",skip=1),byrow=T,ncol=2)
• TheData <- TheData[,2]
• Blow1 <- seq(from=1,to=40,by=1)
• Blow2 <- seq(from=1,to=60,by=1)
• Likes <- matrix(0,nrow=length(Blow1),ncol=length(Blow2))
• for (ii in 1:length(Blow1))
• for (jj in 1:length(Blow2))
• {
• LikeOut <- Propn*dpois(TheData,Blow1[ii])+(1-Propn)*dpois(TheData,Blow2[jj])
• Likes[ii,jj] <- -1*sum(log(LikeOut))
• }
• par(mfrow=c(1,2)); persp(x=Blow1,y=Blow2,z=Likes);contour(x=Blow1,y=Blow2,z=Likes,120)
• Ibest <- which(Likes==min(Likes),arr.ind=T)
• cat(Propn,Blow1[Ibest[1]],Blow2[Ibest[2]],min(Likes),"\n")

• GIVEN
• a value for Propn

19

Slide 20

• Notes on the code:
• dpois is taking a vector and returning a vector.
• persp is producing a 3d plot, contour is similar.
• Ibest is now a vector and not a scalar.
• We are computing the negative log-likelihood.
• takes “propn” and returns the lowest value of the negative log-likelihood.

20

Slide 21

This is a plot of the negative log-likelihood surface

• The minimum negative
• log-likelihood

21

Slide 22

22

Notice anything interesting about this likelihood surface?

Slide 23: The Optim/MLE routines

• The optim function is in the base package and the mle
• function is in the stats4 package – We will focus on mle

• ll1 <- function(x) {
• obj <- (x-0.3)^2
• return(obj)
• }
• ret1 <- mle(ll1,start=list(x=1))
• print(summary(ret1))

Function to be minimized

• Initial values for the
• parameters

23

Slide 24: Some options for MLE-I

• mle(Func, start=list(), method, lower, upper)
• Func – function to be minimized
• start – initial values for the parameters
• method – method for minimization (e.g. “L-BFGS-B”,
• “Nelder-Mead”, “BFGS”).
• lower – lower bounds for the variables (including –Inf)
• upper - upper bounds for the variables (including Inf)

24

Note, “L-BFGS-B” allows you to bound the parameters

Slide 25: Question-1

• You monitor an intersection for cars that run the
• red light. Based on 50 days of observation, you detect
• 30 days on which “light running” took place. You repeat the “experiment” after a red light camera is installed and now observe 10 days on which “running” occurred.

• Questions
• What is the probability of “light running” before the camera
• was installed?
• Has the probability of “light running” changed significantly?
• We’ll get to this question later

Slide 26: Example 1 - revisited

• LL2 <- function(p) {
• obj <- -1*30*log(p)-20*log(1-p)
• return(obj)
• }
• ret2 <- mle(LL2,start=list(p=0.1),
• lower=c(0.00001),upper=c(0.999999),
• method="L-BFGS-B")
• print(summary(ret2))

Negative log-likelihood

• “Method” must be L-BFGS-B if you
• want to bound the parameters

26

Slide 27: Example 2 - revisited

The blow rate for a given whale species is Poisson-distributed and it is “known” that the average number of blows in an interval for males and females is different. We monitor a population and find the following distribution of blow times.

• Questions
• What is the fraction of females in the population (females have
• longer blow times)?
• Can we conclude that females and males really do differ in mean
• length of blow?

27

Slide 28

• Now, an animal is either male or female and letting the
• mean blow rate for males be λ1 and that for females be λ2,
• and assuming that the proportion of males is p:

The Poisson distribution is:

28

Example 2 - revisited

The likelihood function is then:

P(x | λ) = λx exp(−λ) / x!

P(x | p, λ1, λ2) = p λ1x exp(−λ1) / x! + (1 − p) λ2x exp(−λ2) / x!

P(x | p, λ1, λ2) = ∏i ( p λ1xᵢ exp(−λ1) / xi! + (1 − p) λ2xᵢ exp(−λ2) / xi! )

Slide 29

29

Example 2 - revisited

The likelihood function is then:

Question: What are we solving for here with optimization?

P(x | p, λ1, λ2) = ∏i ( p λ1xᵢ exp(−λ1) / xi! + (1 − p) λ2xᵢ exp(−λ2) / xi! )

Slide 30: Example 2 - revisited

• LL3 <- function(lam1,lam2,p,Blows)
• {
• LikeOut <- p*dpois(Blows,lam1)+(1-p)*dpois(Blows,lam2)
• NegLogLike <- -1*sum(log(LikeOut))
• return(NegLogLike)
• }
• TheData <- matrix(scan(”Whale.dat",skip=1),byrow=T,ncol=2)
• Blows <- TheData[,2]
• ret3 <- mle(LL3,start=list(lam1=10,lam2=10,p=0.4),fixed=list(Blows=Blows),lower=c(0.01,0.01,0.00001),upper=c(Inf,Inf,0.99999),method="L-BFGS-B")
• print(summary(ret3))

30

Slide 31: Dissecting Example 2-I

• LL3 <- function(lam1,lam2,p,Blows)
• {
• LikeOut <- p*dpois(Blows,lam1)+(1-p)*dpois(Blows,lam2)
• NegLogLike <- -1*sum(log(LikeOut))
• return(NegLogLike)
• }

Parameters

The data

• This is a vector of
• likelihoods for each
• data point

• This is the negative
• log-likelihood (sum over negative
• log-likelihoods for each datum)

31

Slide 32: Dissecting Example 2-II

• ret3 <- mle(LL3,
• start=list(lam1=10,lam2=10,p=0.4),
• fixed=list(Blows=Blows),
• lower=c(0.01,0.01,0.00001),
• upper=c(Inf,Inf,0.99999),
• method="L-BFGS-B")

• Starting values
• (must be a list)

• Data
• (must be a list)

• Upper and
• lower bounds for
• each parameter

32

Slide 33: Hints for optimization

• Consider different parameterizations (e.g. log(x) rather than x).
• Methods may work better if the variables are all scaled
• Never trust a non-linear minimizer!
• Try different starting values.
• Combine different methods
• Test your code with data for which you know the correct answer (simulate!).

33

Slide 34: Question 3

The length of a fish is related to its age according to the von Bertalanffy growth curve:

• Questions
• Estimate the parameters of this equation (what are they?) for
• males and females, and compute 95% confidence intervals for
• asymptotic length.
• Can we conclude that the asymptotic sizes for females and males
• are the same?

La = ℓ∞(1 − e−κa) + ε

ε ~ N(0, σ2)

Slide 35: Example 3-I

• Estimate the values for ℓ∞, κ, and σ from the age and growth data.
• Create vectors with the age, sex, and length data.
• Create a function that takes values for the
• parameters as well as the data, and computes the
• negative log-likelihood.
• LL4 <- function(Linf1,Kappa1,Linf2,Kappa2,Sigma,Age,Length,Sex)
• use mle to find the maximum likelihood estimates
• Bound the region of parameter space to reasonable values
• E.g. Can the parameters be negative?

35

Slide 36: Likelihood Profile

Parameter of Interest

• -log
• likelihood

MLE

•

•

•

•

•

•

•

• Alternative
• values for the
• parameter

• Negative log-likelihood, fixing the
• parameter to each value in turn
• (minimize on the rest of the parameters)

36

Slide 37: Confidence Interval

• 95% confidence
• interval

1.92

37

The mle() function in R does this for you, but you can pretty easily do it by hand.