intro-MLEs.qmd

# Introduction to Maximum Likelihood Inference {#sec-intro-MLEs}

---

These notes are derived primarily from @dobson4e (mostly chapters 1-5).

Some material was also taken from @mclachlan2007em and @CaseBerg01.

---

{{< include shared-config.qmd >}}

## Overview of maximum likelihood estimation

### The likelihood function

:::{#def-lik-obs}
#### Likelihood of a single observation

Let $X$ be a random variable and let $x$ be $X$'s observed data value.
Let $\p_{\Th}(X=x)$ be a probability model for the distribution of $X$, with parameter vector $\Th$.

Then the **likelihood** of parameter value $\th$, for model $\p_{\Th}(X=x)$ and data $X = x$, is simply the probability of the event $X=x$ given $\Th = \th$:
 
$$
\ba
\Lik(\theta) &\eqdef \P_{\theta}(X = x)
\ea
$$
:::

---

:::{#def-lik}
#### Likelihood of a dataset

Let $\vec x$ be a dataset with corresponding random variable $\vec X$.
Let $\p_{\Th}(\vec X)$ be a probability model for the distribution of $\vX$ with unknown parameter vector $\Th$.

Then the **likelihood** of parameter value $\th$, for model $\p_{\Th}(X)$ and data $\vX = \vx$, is the *joint probability* of $\vX = \vx$ given $\Th = \th$:
 
$$
\ba
\Lik(\theta) &\eqdef p_{\theta}(\vX = \vx)
\\&=p_{\theta}(X_1=x_1, ..., X_n = x_n)
\ea
$$
:::

::: callout-note

#### Notation for the likelihood function

The likelihood function can be written as:

- $\Lik(\theta)$
- $\Lik(\vec x;\theta)$
- $\Lik(\theta; \vec x)$
- $\Lik_{\vec x}(\theta)$
- $\Lik_{\theta}(\vec x)$
- $\Lik(\vec x | \theta)$

All of these notations mean the same thing.
:::

::: notes
The likelihood is a function that takes $\theta$ (and implicitly, $\vec X$) as inputs and outputs a single number, the joint probability of $\vec x$ for model $p_\Theta(\vX=\vx)$ with $\Theta = \theta$.
:::

---

:::{#thm-lik-iid}

#### Likelihood of an independent sample

For [mutually independent](probability.qmd#def-indpt) data $X_1, ..., X_n$:

$$\Lik(\vec x|\theta) = \prod_{i=1}^n \p(X_i=x_i|\theta)$$ {#eq-Lik}

:::

:::{.proof}

$$
\ba
\Lik(\vec x|\theta) 
&\eqdef \p(X_1 = x_1, …,X_n =x_n|\theta) 
\\&= \prod_{i=1}^n \p(X_i=x_i|\theta)
\ea
$$
The second equality is by the definition of statistical independence.

:::

---

:::{#def-lik-factor}
#### Likelihood components

Given an $\iid$ dataset $\vec x$, the **likelihood component** or **likelihood factor** of observation $X_i=x_i$ is the marginal likelihood of $X_i=x_i$: 

$$\Lik_i(\theta) = \P(X_i=x_i)$$

:::

---

:::{#thm-ds-lik-obs-lik}

For $\iid$ data $\vx \eqdef \x1n$, 
the likelihood of the dataset is equal to the product of the observation-specific likelihood factors:

$$\Lik(\theta) = \prodin \Lik_i(\theta)$$
:::

---

### The maximum likelihood estimate

:::{#def-mle}

#### Maximum likelihood estimate

The **maximum likelihood estimate** of a parameter vector $\Theta$, denoted $\hthml$, is the value of $\Theta$ that maximizes the likelihood:

$$
\hthml \eqdef \arg \max_\Th \Lik(\Th)
$$ {#eq-mle}
:::

### Finding the maximum of a function

Recall from calculus: the maxima of a continuous function $f(x)$ over a range of input values $\rangef{x}$ can be found either:

- at the edges of the range of input values, *OR*:
- where the function is flat (i.e. where the gradient function $f'(x) = 0$) *AND* the second derivative is negative definite ($f''(x) < 0$).

### Directly maximizing the likelihood function for *iid* data

To find the maximizer(s) of the likelihood function, we need to solve $\Lik'(\th) = 0$ for $\theta$. However, even for mutually independent data, we quickly run into a problem:

$$
\ba
\Lik'(\th) 
&= \deriv{\th} \Lik(\th)
\\ &= \deriv{\th} \prod_{i=1}^n p(X_i=x_i|\theta)
\ea
$$ {#eq-deriv-Lik}

The derivative of the likelihood of independent data is the derivative of a product. 
We will have to perform a massive application of the product rule to evaluate this derivative.

### The log-likelihood function

It is typically easier to work with the log of the likelihood function:

:::{#def-loglik}

#### Log-likelihood

The **log-likelihood** of parameter value $\theta$, for model $\p_{\Theta}(\vX)$ and data $\vX = \vx$, is the natural logarithm of the likelihood^[ <https://en.wikipedia.org/wiki/Does_exactly_what_it_says_on_the_tin>]:

$$\lik(\th) \eqdef \logf{\Lik(\th)}$$
:::

---

:::{#def-loglik}

#### Log-likelihood components

Given a dataset $\vX = \vx$, the **log-likelihood component of observation $X_i=x_i$** is the natural logarithm of the likelihood component:

$$\lik_i(\th) \eqdef \logf{\Lik_i(\th)}$$
:::

---

:::{#thm-mle-use-log}
####

The likelihood and log-likelihood have the same maximizer:

$$
\am_\th \Lik(\th) = \am_\th \lik(\th)
$$
::: 

::: proof
Left to the reader.
:::

---

:::{#thm-llik-iid}

#### Log-likelihood of an $\iid$ sample

For $\iid$ data $X_1, ..., X_n$ with shared distribution $\p(X=x)$:

$$\ell(x|\theta) = \sum_{i=1}^n \log{p(X=x_i|\theta)}$$ {#eq-loglik}
:::

:::{.proof}
$$
\ba
\ell(x|\theta) 
&\eqdef \log{\Lik(\vec x|\theta)}
\\&= \log{\prod_{i=1}^n \p(X_i=x_i|\theta)}
\\&= \sum_{i=1}^n \log{p(X=x_i|\theta)}
\ea
$$
:::

---

::: notes

For $\iid$ data, we will have a much easier time taking the derivative of the log-likelihood:

:::

:::{#thm-deriv-llik-iid}
#### Derivative of the log-likelihood function for $\iid$ data

For $\iid$ data:

$$\ell'(\theta) = \sumin \deriv{\theta} \log{\p(X=x_i|\theta)}$$ {#eq-deriv-llik}

:::

:::{.proof}
$$
\ba
\lik'(\th) 
&= \deriv{\th} \lik(\th)
\\ &= \deriv{\th} \sum_{i=1}^n \log{\p(X=x_i|\theta)}
\\ &= \sum_{i=1}^n \deriv{\th} \log{\p(X=x_i|\theta)}
\ea
$$
:::

---

### The score function

The first derivative^[a.k.a. the [gradient](https://en.wikipedia.org/wiki/Gradient)] of the log-likelihood, $\lik'(\th)$, is important enough to have its own name: the *score function*.

:::{#def-score}
#### Score function

The **score function** of a statistical model $\pr(\vec X=\vec x)$ is the gradient (i.e., first derivative) of the log-likelihood of that model:

$$\lik'(\th) \eqdef \deriv{\th} \lik(\th)$$
:::

::: notes
We often
skip writing the arguments $x$ and/or $\theta)$, so
$\ell' \eqdef \ell'(\vec x;\theta) \eqdef \ell'(\theta)$.[^1] Some statisticians
use $U$ or $S$ instead of $\ell'$. I prefer $\ell'$.
Why use up extra letters?
:::

### Asymptotic distribution of the maximum likelihood estimate

::: notes
We learned how to quantify our uncertainty about these maximum likelihood estimates; with sufficient sample size, $\hthml$ has an approximately Gaussian distribution [@newey1994large]:
:::
$$
\hat\theta_{ML} \dot \sim N(\theta,\mathcal I(\theta)^{-1})
$$

Recall:

- $\einf(\theta) \eqdef \E{\oinf(\vX;\theta)}$
- $\oinf(\vX,\theta) \eqdef -\ell''(\vX;\theta)$

We can estimate $\einf(\th)$ using either $\einf(\hthml)$ or $\oinf(\vec x; \hthml)$.

So we can estimate the standard error of $\hth_k$ as:

$$
\HSE{\hth_k} = \sqrt{\sb{\inv{\heinff{\hthml}}}_{kk}}
$$

### The (Fisher) (expected) information matrix

The variance of $\ell'(x,\theta)$,
${Cov}\left\{ \ell'(x,\theta) \right\}$, is also very
important; we call it the "expected information matrix", "Fisher
information matrix", or just "information matrix", and we represent it
using the symbol $\einff{I}$ (`\frakturI` in Unicode, `\mathfrak{I}` in LaTeX).

$$
\ba
\einf 
\eqdef \einf(\theta) 
\\ &\eqdef \Covf{\ell'|\theta} 
\\ &= \Expp[ \ell'{\ell'}\' ] - \Expp[ \ell' ] \ \Expp[ \ell' ]\'
\ea
$$

The elements of $\mathfrak{I}$ are:

$$
\ba
\mathfrak{I}_{ij} 
&\eqdef \Covf{{\ell'}_{i},{\ell'}_{j}}
\\ &= \Expp[ \ell_{i}'\ell_{j}' ] - \Expp[ {\ell'}_{i} ] \Expp[ {\ell'}_{j} ]
\ea
$$

Here, 

$$
\ba
\E{\ell'}
&\eqdef \int_{x \in \rangef{x}}
{
\ell'(x,\th) \p(X = x | \th) dx
}
\\ &= \int_{x \in \rangef{X}}
{
\paren
{
\deriv{\th}
\log{\p(X = x | \th)}
}
\p(X = x | \theta) dx
}
\\ &= 
\int_{x \in \rangef{X}}
{
\frac
{\deriv{\theta} \p(X = x | \th)}
{\p(X = x | \theta)}
\p(X = x | \theta) dx
}
\\ &= 
\int_{x \in \rangef{X}}
{
\deriv{\theta} \p(X = x | \th) dx
}
\ea
$$



And similarly

$$
\Exp{\ell' \ell'^{\top}} 
\eqdef 
\int_{x \in R(x)}
{\ell'(x,\theta)\ell'(x,\theta)^{\top}\ 
\pf{X = x | \th}\ dx}
$$

Note that $\Exp{\ell'}$ and
$\Exp{\ell'{\ell'}^{\top}}$
are functions of $\theta$ but not of $x$; 
the expectation operator removed $x$.

Also note that for most of the distributions you are familiar with
(including Gaussian, binomial, Poisson, exponential):

$$\Exp{\ell'} = 0$$

So

$$\einff{\theta} = \Exp{\ell'{\ell'}^{\top} }$$

Moreover, for those distributions (called the "exponential family"), we
have:

$$
\mathfrak{I} = -\Exp{\ell''}
= \Exp{- \ell''}
$$

(see @dobson4e, §3.17), where

$$\ell'' \eqdef \deriv{\theta}\ell^{'(x,\theta)^{\top}} = \deriv{\theta}\deriv{\theta^{\top}}\ell(x,\theta)$$

is the $p \times p$ matrix whose elements are:

$$\ell_{ij}'' \eqdef \deriv{\theta_{i}}\deriv{\theta_{j}}\log{ p\left( X = x \mid \theta \right)}$$

$\ell''$ is called the "Hessian"^[named after mathematician [Otto Hesse](https://en.wikipedia.org/wiki/Otto_Hesse)] of the log-likelihood
function.

Sometimes, we use $I(\theta;x) \eqdef - \ell''$ (note the
standard-font "I" here). $I(\theta;x)$ is the observed information, precision, or concentration
matrix (Negative Hessian).

:::{.callout-important}
#### Key point
 The asymptotics of MLEs gives us
${\widehat{\theta}}_{ML} \sim N\left( \theta,\mathfrak{I}^{- 1}(\theta) \right)$,
approximately, for large sample sizes.
:::

We can estimate $\einf^{- 1}(\theta)$ by working out
$\Ef{-\ell''}$ or
$\Ef{\ell'{\ell'}^{\top}}$
and plugging in $\hthml$, but sometimes we instead use
$\oinf(\hthml,\vx)$ for convenience; there are
some cases where it’s provably better according to some criteria (@efron1978assessing).

### Iterative maximization {#sec-newton-raphson}

(c.f., @dobson4e, Chapter 4)

::: notes
Later, 
when we are trying to find MLEs for likelihoods which we can’t easily differentiate, 
we will "hill-climb" using the Newton-Raphson algorithm:
:::

$$
\begin{aligned}
\esttmp{\theta} 
&\leftarrow \esttmp{\theta} + \inv{\oinff{\vec y;\esttmp{\theta}}}
\scoref{\vec y;\esttmp{\theta}}
\\ &= \esttmp{\theta} - \inv{\hessf{\vec y;\esttmp{\theta}}} 
\scoref{\vec y;\esttmp{\theta}}
\end{aligned}
$$

---

::: notes
The reasoning for this algorithm is that we can approximate the the score function using the first-order [Taylor polynomial](https://en.wikipedia.org/wiki/Taylor%27s_theorem):
:::

$$
\ba
\score(\th) 
&\approx \score^*(\th)
\\ &\eqdef \score(\esttmp{\th}) + \hessian(\esttmp{\th})(\th - \esttmp{\th})
\ea
$$

---

::: notes
The approximate score function, $\score^*(\th)$, is a linear function of $\th$, so it is easy to solve the corresponding approximate score equation, $\score^*(\th) = 0$, for $\th$:

:::

$$
\ba
\th 
&= \esttmp{\th} - \score(\esttmp{\th}) \cd \inv{\hessian(\esttmp{\th})}
\ea
$$

---

For computational simplicity, we will sometimes use
$\mathfrak{I}^{- 1}(\theta)$ in place of
$I\left( \widehat{\theta},y \right)$; 
doing so is called "Fisher scoring" or the "method of scoring". 
Note that this is the opposite of the substitution that we are making for estimating the variance of the MLE; 
this time we should technically use the observed information but we use the expected information instead.


---

There’s also an "empirical information matrix" (see @mclachlan2007em):

$$I_{e}(\theta,y) \eqdef \sum_{i = 1}^{n}{\ell_{i}'\ {\ell_{i}'}^{\top}} - \frac{1}{n}\ell'{\ell'}^{\top}$$

where $\ell_{i}$ is the log-likelihood of the ith observation.
Note that $\ell' = \sum_{i = 1}^{n}\ell_{i}'$.

$\frac{1}{n}I_{e}(\theta,y)$ is the sample equivalent of

$$\mathfrak{I \eqdef I(}\theta) \eqdef {Cov}\left( \ell'|\theta \right) = E[ \ell'{\ell'}^{\top} ] - E[ \ell' ]\ E[ \ell' ]^{\top}$$

$$\left\{ \mathfrak{I}_{jk} \eqdef {Cov}\left( {\ell'}_{j},{\ell'}_{k} \right) = E[ \ell_{j}'\ell_{k}' ] - E[ {\ell'}_{j} ] E[ {\ell'}_{k} ] \right\}$$

$I_{e}(\theta,y)$ is sometimes computationally easier to compute for
Newton-Raphson-type maximization algorithms.

c.f. <https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization>

### Quantifying (un)certainty of MLEs

#### Confidence intervals for MLEs

An asymptotic approximation of a 95% confidence interval for $\theta_k$ is

$$
\hthml \pm z_{0.975} \times \HSE{\hth_k}
$$

where $z_\beta$ the $\beta$ quantile of the standard Gaussian distribution.

#### p-values and hypothesis tests for MLEs

(to add)

#### Likelihood ratio tests for MLEs

log(likelihood ratio) tests [c.f. @dobson4e §5.7]:

$$-2\ell_{0} \sim \chi^{2}(p - q)$$

See also <https://online.stat.psu.edu/stat504/book/export/html/657>

#### Prediction intervals for MLEs

$$\overline{X} \in [ \widehat{\mu} \pm z_{1 - \alpha\text{/}2}\frac{\sigma}{m} ]$$

Where $m$ is the sample size of the new data to be predicted (typically
1, except for binary outcomes, where it needs to be bigger for
prediction intervals to make sense)

[^1]: I might sometimes switch the order of $x,$ $\theta$; this is
    unintentional and not meaningful.

## Example: Maximum likelihood for Tropical Cyclones in Australia

{{< include dobson-cyclone-example.qmd >}}

## Maximum likelihood inference for univariate Gaussian models

Suppose $X_{1}, ..., X_{n} \siid N(\mu, \sigma^{2})$.
Let $X = (X_{1},\ldots,X_{n})^{\top}$ be these random
variables in vector format. Let $x_{i}$ and $x$ denote the corresponding
observed data. Then $\theta = (\mu,\sigma^{2})$ is
the vector of true parameters, and $\Theta = (\Mu, \Sigma^2)$ is the vector of parameters as a random
vector.

Then the log-likelihood
is:

$$
\begin{aligned}
\ell 
&\propto - \frac{n}{2}\log{\sigma^{2}} - \frac{1}{2}\sum_{i = 1}^{n}\frac{( x_{i} - \mu)^{2}}{\sigma^{2}}\\
&= - \frac{n}{2}\log{\sigma^{2}} - \frac{1}{2\sigma^{2}}\sum_{i = 1}^{n}{x_{i}^{2} - 2x_{i}\mu + \mu^{2}}
\end{aligned}
$$


### The score function

$$\ell'(x,\theta) \eqdef \deriv{\theta}\ell(x,\theta) = \left( \begin{array}{r}
\deriv{\mu}\ell(\theta;x) \\
\deriv{\sigma^{2}}\ell(\theta;x)
\end{array} \right) = \left( \begin{array}{r}
\ell_{\mu}'(\theta;x) \\
\ell_{\sigma^{2}}'(\theta;x)
\end{array} \right)$$.

$\ell'(x,\theta)$ is the function we set equal to 0 and solve
to find the MLE:

$${\widehat{\theta}}_{ML} = \left\{ \theta:\ell'(x,\theta) = 0 \right\}$$

### MLE of $\mu$

$$
\ba
\frac{d\ell}{d\mu} 
&= - \frac{1}{2}\sum_{i = 1}^{n}
\frac{- 2(x_{i} - \mu)}{\sigma^{2}}
\\ &= \frac{1}{\sigma^{2}}
\sb{
    \paren{
        \sum_{i = 1}^{n}x_{i}
    } 
    - n\mu
}
\ea
$$

If $\frac{d\ell}{d\mu} = 0$, then
$\mu = \overline{x} \eqdef \frac{1}{n}\sum_{i = 1}^{n}x_{i}$.

$$\frac{d^{2}\ell}{(d\mu)^{2}} = \frac{- n}{\sigma^{2}} < 0$$

So ${\widehat{\mu}}_{ML} = \overline{x}$.

### MLE of $\sigma^{2}$

:::{.callout-tip}
#### Reparametrizing the Gaussian distribution

When solving for ${\widehat{\sigma}}_{ML}$, you can treat
$\sigma^{2}$ as an atomic variable (don’t differentiate with respect to
$\sigma$ or things get messy). In fact, you can replace $\sigma^{2}$
with $1/\tau$ and differentiate with respect to $\tau$ instead, and the
process might be even easier.
:::

$$\frac{d\ell}{d\sigma^{2}} = \deriv{\sigma^{2}}\left( - \frac{n}{2}\log{\sigma^{2}} - \frac{1}{2}\sum_{i = 1}^{n}\frac{\left( x_{i} - \mu \right)^{2}}{\sigma^{2}} \right)\ $$

$$= - \frac{n}{2}\left( \sigma^{2} \right)^{- 1} + \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2}$$

If $\frac{d\ell}{d\sigma^{2}} = 0$, then:

$$\frac{n}{2}\left( \sigma^{2} \right)^{- 1} = \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2}$$

$$\sigma^{2} = \frac{1}{n}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2}$$

We plug in ${\widehat{\mu}}_{ML} = \overline{x}$ to maximize globally (a
technique called profiling):

$$\sigma_{ML}^{2} = \frac{1}{n}\sum_{i = 1}^{n}\left( x_{i} - \overline{x} \right)^{2}$$

Now:

$$
\begin{aligned}
\frac{d^{2}\ell}{\left( d\sigma^{2} \right)^{2}} 
&= \deriv{\sigma^{2}}\left\{ - \frac{n}{2}\left( \sigma^{2} \right)^{- 1} + \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2} \right\}\\
&= \left\{ - \frac{n}{2}\deriv{\sigma^{2}}\left( \sigma^{2} \right)^{- 1} + \frac{1}{2}\deriv{\sigma^{2}}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2} \right\}\\
&= \left\{ \frac{n}{2}\left( \sigma^{2} \right)^{- 2} - \left( \sigma^{2} \right)^{- 3}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2} \right\}\\
&= \left( \sigma^{2} \right)^{- 2}\left\{ \frac{n}{2} - \left( \sigma^{2} \right)^{- 1}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2} \right\}
\end{aligned}
$$

Evaluated at
$\mu = \overline{x},\sigma^{2} = \frac{1}{n}\sum_{i = 1}^{n}\left( x_{i} - \overline{x} \right)^{2}$,
we have:

$$
\begin{aligned}
\frac{d^{2}\ell}{\left( d\sigma^{2} \right)^{2}} 
&= \left( {\widehat{\sigma}}^{2} \right)^{- 2}\left\{ \frac{n}{2} - \left( {\widehat{\sigma}}^{2} \right)^{- 1}\sum_{i = 1}^{n}\left( x_{i} - \overline{x} \right)^{2} \right\}\\
&= \left( {\widehat{\sigma}}^{2} \right)^{- 2}\left\{ \frac{n}{2} - \left( {\widehat{\sigma}}^{2} \right)^{- 1}n{\widehat{\sigma}}^{2} \right\}\\
&= \left( {\widehat{\sigma}}^{2} \right)^{- 2}\left\{ \frac{n}{2} - n \right\}\\
&= \left( {\widehat{\sigma}}^{2} \right)^{- 2}n\left\{ \frac{1}{2} - 1 \right\}\\
&= \left( {\widehat{\sigma}}^{2} \right)^{- 2}n\left( - \frac{1}{2} \right) < 0
\end{aligned}
$$

Finally, we have:

$$
\begin{aligned}
\frac{d^{2}\ell}{d\mu\ d\sigma^{2}} 
&= \deriv{\mu}\left\{ - \frac{n}{2}\left( \sigma^{2} \right)^{- 1} + \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2} \right\}\\
&= \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\deriv{\mu}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2}\\
&= \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}{- 2(x_{i} - \mu)}\\
&= - \left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}{(x_{i} - \mu)}
\end{aligned}
$$

Evaluated at
$\mu = \widehat{\mu} = \overline{x},\sigma^{2} = {\widehat{\sigma}}^{2} = \frac{1}{n}\sum_{i = 1}^{n}\left( x_{i} - \overline{x} \right)^{2}$,
we have:

$$\frac{d^{2}\ell}{d\mu\ d\sigma^{2}} = - \left( {\widehat{\sigma}}^{2} \right)^{- 2}\left( n\overline{x} - n\overline{x} \right) = 0$$

### Covariance matrix

$$I = \begin{bmatrix}
\frac{n}{\sigma^{2}} & 0 \\
0 & \left( {\widehat{\sigma}}^{2} \right)^{- 2}n\left( - \frac{1}{2} \right)
\end{bmatrix} = \begin{bmatrix}
a & 0 \\
0 & d
\end{bmatrix}$$

So:

$$I^{- 1} = \frac{1}{ad}\begin{bmatrix}
d & 0 \\
0 & a
\end{bmatrix} = \begin{bmatrix}
\frac{1}{a} & 0 \\
0 & \frac{1}{d}
\end{bmatrix}$$

$$I^{- 1} = \begin{bmatrix}
\frac{{\widehat{\sigma}}^{2}}{n} & 0 \\
0 & \frac{{2\left( {\widehat{\sigma}}^{2} \right)}^{2}}{n}
\end{bmatrix}$$

See @CaseBerg01 p322, example 7.2.12.

To prove it’s a maximum, we need:

- $\ell' = 0$

- At least one diagonal element of $\ell''$ is negative.

- Determinant of $\ell''$ is positive.


{{< include HERS-example.qmd >}}