Maximum Likelihood Estimation

1. Definition of a Likelihood Function

A likelihood function \(L(\theta|X)\) as a function of a parameter \(\theta\) given a sample \(X = (x_1,x_2, \cdots, x_n)^T, x_i\in\mathbb{R}\) for \(i\in[n]\), which is independently and identically distributed () is the product of probabilities of observations \(x_i\) given \(\theta\): \[L(\theta|X) = \prod_{i=1}^n p(x_i|\theta)\].

2. Example

2.1 Samples from a Bernoulli distribution with parameter \(\theta\)

Suppose we have an sample \(X = (x_1,x_2, \cdots, x_n)^T, x_i\in\{0,1\}\) for \(i\in[n]\). We assume that the sample is generated from a Bernoulli distribution with parameter \(\theta \in (0,1)\). The probability of observing \(x_i\) given \(\theta\) is \[p(x_i|\theta) = \theta^{x_i}(1-\theta)^{(1-x_i)}\].

Then the likelihood function \(L\) of parameter \(\theta\) given \(X\) is \[L(\theta|X) = \prod_{i=1}^n p(x_i|\theta) = \theta^{\sum_{i=1}^n x_i}(1-\theta)^{\sum_{i=1}^n (1-x_i)} = \theta^{n_1}(1-\theta)^{n_0}\] where \(n_1 = \sum_{i=1}^n x_i\) and \(n_0 = \sum_{i=1}^n (1-x_i)\), \(n = n_1 + n_0\).

3. Estimation

3.1 Maximum Likelihood Estimation

We want to find \(\theta\) that maximizes the likelihood \(L(\theta|X)\). In other words, \[\widehat{\theta}_{\text{MLE}} = \arg\max_{\theta \in \Theta}L(\theta|X)\].

This is equivalent to \[\widehat{\theta}_{\text{MLE}} = \arg\max_{\theta \in \Theta}\text{log} L(\theta|X)\], which is a monotonic transformation of the likelihood function \(L(\theta|X)\).

teta= seq(0,1, by =1e-3)
sample1 = rbinom(100,1,1/3) # interesting observation: if p becomes bigger or smaller, likelihood useless
bernlike= function(y,x){
  temp =  y^(sum(x))*(1-y)^(length(x)-sum(x))
  #temp = -(2*log(y) + log(1-y))
  return(temp)
}
sum(sample1)

## [1] 37

plot(teta, bernlike(teta, sample1), xlab = "theta", ylab = "L(theta|X)", main = "Likelihood function of Bernoulli distribution with sample")

Let’s find the \(\widehat{\theta}_{\text{MLE}}\) by three approaches:

1. Analytic approach

By solving the first order condition of objective function with respect to \(\theta\), we can find the \(\widehat{\theta}_{\text{MLE}}\):

F.O.C. [\(\theta\)]

\(\begin{align*} \frac{d}{d\theta} \text{log}L(\theta|X) =& \frac{d}{d\theta}\sum_{i=1}^n\big( x_i \text{log}\theta + (1-x_i) \text{log}(1-\theta)\big)\\ =& \frac{1}{\theta}\sum_{i=1}^n x_i - \frac{1}{1-\theta}\sum_{i=1}^n(1- x_i) = 0\\ \widehat{\theta}_{\text{MLE}}=& \frac{1}{n}\sum_{i=1}^n x_i = \overline{x} \end{align*}\)

theta.hat = sum(sample1)/length(sample1)
theta.hat

## [1] 0.37

2. Built-in solver optimize

Also we can use a built-in solver optimize to find the \(\widehat{\theta}_{\text{MLE}}\):

#install.packages("stats")
#library(stats)
bernloglike= function(y,x){
  #temp =  y^2*(1-y)
  temp = -(sum(x)*log(y) + (length(x)-sum(x))*log(1-y))
  return(temp)
}

stats::optimize(bernloglike, c(0,1), x = sample1)$minimum # is this true? why different value from above? need to construct gradient

## [1] 0.3699958

3. Gradient Descent

Also we can apply famous optimization method Gradient Descent (GD).

grad = function(y, x){
  temp = - (1/(y*(1-y))*sum(x) - length(x)/(1-y))
  return(temp)
}

tol = 1e-6
diff = 1
theta0 = 0.1 # if the initial value is too extreme, then it does not work well.
alpha = 1e-3
theta.collect = list()
while (TRUE){
  theta1 = theta0 - alpha*grad(theta0, sample1)
  theta.collect[length(theta.collect)+1] = theta1
  #if (abs(bernloglike(theta1, sample)-bernloglike(theta0, sample)) <=tol){
  if (abs(theta1 - theta0)<=tol){
    break
    print(theta1)
  }else 
  {theta0 = theta1
  }
}
plot(seq(from =1, to = length(theta.collect)), theta.collect, xlab = "Iteration", ylab = expression(theta), main = expression(paste("Convergence of ", theta, " by gradient descent")))

theta.collect[length(theta.collect)]

## [[1]]
## [1] 0.3700013