```python
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
from typing import List, Tuple
```


```python
data = np.array([0.5, -0.32, -0.55, -0.76, -0.07, 0.44, -0.48])
theta = np.linspace(-10, 10, 10000)
```

(a) Let $x_1,...,x_n \sim DExp(\theta)$. The contribution of $x_i$ to the likelihood is:
$$L_i(\theta) = \frac{1}{2} e^{|x_i - \theta|} $$

and the total log-likelihood is:
\begin{align}
\log L(\theta) & = \sum_{i=1}^n \log L_i(\theta) \\
& = -n \log 2 + \sum_{i=1}^n |x_i - \theta|
\end{align}


```python
def log_likelihood_double_exp(theta: List[float], data: List[float]) -> List[float]:
    log_like = []
    for th in theta:
        ll = 0
        for dt in data:
            ll += np.abs(dt - th)
        ll *= -np.log(2)
        log_like.append(ll)
    return log_like
```


```python
def plot_log_likelihood(theta: List[float], data: List[float]):
    log_like = log_likelihood_double_exp(theta, data)
    plt.plot(theta, log_like)
    plt.title('Double exponential likelihood')
    plt.xlabel(r'$\theta$')
    plt.ylabel('Log likelihood');
```


```python
plot_log_likelihood(theta, data)
```

(b) Then $\hat{\theta}$ the MLE of $\theta$ is given by the solution to the score equation:

\begin{align}
S(\theta) = 0 & = \frac{\partial}{\partial \theta} \log L(\theta) \\
& = \sum_{i=1}^n sgn |x_i - \theta| \\
& = median(x)
\end{align}

where the last statement is due to the following argument: https://math.stackexchange.com/questions/1678740/mle-of-double-exponential.


```python
print(f"MLE = {np.median(data)}")
```

or alternatively we can use the log-likelihood directly:


```python
log_like = log_like = log_likelihood_double_exp(theta, data)
print(f"MLE = {np.round(theta[np.argmax(log_like)], 2)}") 
```

The likelihood based interval:


```python
def likelihood_interval(theta: List[float],
                        likelihood: List[float],
                        cutoff: float) -> Tuple[float, float]:
    # intersection points occur below and above the maximum likelihood estimate
    mle_index = np.argmax(likelihood)
    interp_below_max = interp1d(likelihood[:mle_index], theta[:mle_index])
    interp_above_max = interp1d(likelihood[mle_index:], theta[mle_index :])
    lower_int = np.round(interp_below_max(cutoff).flatten()[0], 2)
    upper_int = np.round(interp_above_max(cutoff).flatten()[0], 2)
    return (lower_int, upper_int)
```

at a cut-off of:


```python
c = 0.15 # 95% confidence interval
print(f'Likelihood interval for c = {c} is {likelihood_interval(theta, np.exp(log_like), c)}') 
```

(c) If the largest value in the data is now 2.5 instead of 0.5:


```python
data[0] = 2.5
data
```

then the log-loikelihood function is now:


```python
plot_log_likelihood(theta, data)
```

Despite the change in the largest value, the inference of the MLE $\hat{\theta}$ does not change. This is expected as it is the median of the data which remains unchanged. The median is known to be robust to outliers.


```python
print(f"MLE = {np.median(data)}")
```

7986fb8c679b29d88c2cbf81e551561296881c5f
ipynb
Jupyter Notebook
python/chapter-2/exercises/EX2-13.ipynb
covuworie/in-all-likelihood
6638bec8bb4dde7271adb5941d1c66e7fbe12526
MIT
python/chapter-2/exercises/EX2-13.ipynb
covuworie/in-all-likelihood
6638bec8bb4dde7271adb5941d1c66e7fbe12526
MIT
2020-03-24T17:53:04.000Z
2021-08-23T20:16:17.000Z
python/chapter-2/exercises/EX2-13.ipynb
covuworie/in-all-likelihood
6638bec8bb4dde7271adb5941d1c66e7fbe12526
MIT
2021-11-21T10:24:59.000Z
2021-11-21T10:24:59.000Z
Qwen/Qwen-72B
1. YES
2. YES