# Maximum Likelihood Model Parameter Estimation

\begin{equation*}
D = \text{data} \\
\theta = \text{model parameters} \\
y(\theta, x) = \text{model}
\end{equation*}

A common parameter-estimation strategy is to minimize error,

\begin{equation*}
\sum_{i}(D_i - y_i(\theta, x))^2
\end{equation*}

E.g.

Create some phony straight line data with $\theta=(a, b)$ and normally distributed noise in $y$,


```
n_points = 10
x = linspace(0,10,n_points)
def y(x, theta=(2, 4), noise=0):
    return  theta[0] * x + theta[1] + noise*randn(n_points)
D = y(x, noise=1)
scatter(x,D)
```

We are trying to discover parameters, $\theta$, so write error in terms of only $theta$,

$SSE(\theta; y(x), D)$


```
# cheat, do this with global scope of y, x and D
def ssq_err(theta):
    return np.sum((D - y(x,theta))*(D - y(x,theta)).transpose())

from scipy.optimize import minimize
fit = minimize(ssq_err,[1,1])
print fit
scatter(x, D)
plot(x,y(x, fit.x))
print fit
```

Not amazing so far $\ldots$

__Leading quesiton__:  Is normally distributed noise "unusually compatible" with least squares errors minimization?

Maybe we should do something that takes into account the expected distibution of the model errors? What if the errors are distributed by Poisson, for example, does sum-squared-error still make sense?


```
from scipy.stats import norm, poisson
rvn, rvp = norm(5,1), poisson(5)
x = np.linspace(-2,10,140)
plot(x, rvn.pdf(x))
xd = np.arange(0, np.minimum(rvp.dist.b, 15))
vlines(xd, 0, rvp.pmf(xd), lw=2, color="red")
```

__New Idea__:  The parameter fitting problem can be framed as "Can we find the mostly likely value of the parameters $\theta$ of the model given the data?"

i.e. maximize $P(\theta|D)$

### Frequentist Statistics

Significance is specified by probability that
the data was drawn from the “wrong”
model (null hypothesis): p value

* Not the probability a model is correct
* Only random variables have probabilities

That is, from a frequentist perspective, the PDF for normal noise is in terms of $P(D|\theta)$,

E.g.
\begin{equation*}
P(D|\theta) \propto exp(\sum_i -\frac{(D_i - y_i(\theta))^2}{2\sigma^2})
\end{equation*}

We need to flip this around.

### Bayesian inference

* Frequentist: Probability is a frequency
(of a random variable)
* Bayesian: Probability is a measure of
belief or certainty

Some statistics Rules

\begin{equation}
P(A,B) = P(A|B)P(B) = P(B|A)P(A) \\
P(A) = \int_\Omega P(A,B)dB = \int_\Omega P(A|B)P(B)dB \\
\end{equation}

Bayesian statistics is the application of these rules

Back to data and parameters:

\begin{equation*}
P(\theta|D) = \frac{P(D|\theta)P(\theta)}{P(D)}
\end{equation*}

\begin{equation*}
P(D|\theta) : \text{Likelihood} \\
P(\theta) : \text{Prior} \\
P(D) : \text{Evidence} \\
P(\theta|D) : \text{Posterier} \\
\end{equation*}

\begin{equation*}
P(D) = \int_\Omega P(D|\theta)P(\theta)d\theta \\
\end{equation*}

Integration over the parameter space of the model.

## Simple Discrete Model: Coin Toss

In 10 trials, we get 7 heads: What's our best estimate of the probability of flipping heads with this coin?  

The distribution of coin toss outcomes is the Binomial Distribution.

$pmf(y|t,w) = \binom{n}{k} p^k (1-p)^{(n-k)}$

Let's not forget how obvious this is and guess 0.7.


```
from scipy.stats import binom
rv = binom(10, 0.7)
x = np.arange(0, np.minimum(rv.dist.b, 10))
h = plt.vlines(x, 0, rv.pmf(x), lw=2, color="orange")
```

"Speciﬁcally, the PDF in Fig. 1 is a function of the data given a particular set of parameter values, deﬁned on the data scale. On the other hand, the likelihood function is a function of the parameter given a particular set of observed data, deﬁned on the
parameter scale. In short, Fig. 1 tells us the probability of a particular data value for a ﬁxed parameter, whereas Fig. 2 tells us the likelihood (‘‘unnormalized probability’’) of a particular parameter value for a ﬁxed data set." (.J. Myung / Journal of Mathematical Psychology 47 (2003) 90–100)


```
def f(w):
    rv = binom(10, w)
    return rv.pmf(7)
# Hold the data constant, vary the parameter
x = np.linspace(0,1,100)
plot(x, f(x))
```

It looks like the maximum likelihood on the data over the parameter is around 0.7?

Our data could have been (6,10) or (8,10) with a 0.7 coin.  But we started by asking about the coin, to make any headway, we need to flip the question around.


```
dat = [10, 7]
def L_binom(w):
    rv = binom(dat[0], w)
    x = dat[1]
    # -1 to make max problem a min problem
    return -rv.pmf(x)
```


```
x = np.linspace(0.1, 0.9, 9)
scatter(x, L_binom(x))
from scipy.optimize import minimize_scalar
minimize_scalar(L_binom, bounds=[0,1], method="bounded")
```

## Slightly More Complicated Model: Exponential with Normal Error


```
n_points = 10
x = linspace(0,100,n_points)
def y(x, theta=(10, 0.05), noise=0):
    return  theta[0] * np.exp(-x * theta[1]) + noise*randn(n_points)
D = y(x, noise=1)
scatter(x,D)
```

Likelihood:
    
\begin{equation}
P(D|\theta) \propto \prod_i {\exp(-\frac{(D_i - y_i(\theta))^2}{2\sigma^2})} \propto \prod_i {\exp(-\frac{(D_i - (10\exp(-wx)))^2}{2})}
\end{equation}




```
# use x, D from global scope
def L_cost(w):
    return np.prod(np.exp(((D - y(x,(10, w)))*(D - y(x,(10, w))).transpose())/2.))
w = np.linspace(0,0.1,10)
L_c = [L_cost(i) for i in w]
scatter(w,L_c)
```

Numerically troubling!  Turn the product into a sum with log!


```
# use x, D from global scope
def L_cost(w):
    return np.sum(((D - y(x,(10, w)))*(D - y(x,(10, w))).transpose())/2.)
w = np.linspace(0,0.1,10)
L_c = [L_cost(i) for i in w]
scatter(w,L_c)
from scipy.optimize import minimize_scalar
minimize_scalar(L_cost, bounds=[0,0.1], method="bounded")
```

So it might be useful to use the $\log$ likelihood:

\begin{equation}
\log(P(D|\theta)) \propto \sum_i \log(exp(-\frac{(D_i - (10\exp(-wx)))^2}{2})) = - \sum_i \frac{(D_i - (10\exp(-wx)))^2}{2}
\end{equation}

__Normally distributed errors + Minimum-Log-Likelihood is the same solution as Lease Squares.__

## Exponential Decay and Poisson Likelihood

\begin{equation}
\log(P) \propto \sum_i \log(\frac{\lambda^{x_i}e^{-\lambda}}{\Gamma(x_i)}) = \sum_i (x_i) \log\lambda + (-\lambda) - \log\Gamma(x_i)
\end{equation}

Sterling $\ldots$ or realize that $\Gamma(x_i)$ is a constant w.r.t. MLE estimation of $\lambda$.


\begin{equation}
\log(P(D|\theta)) \propto \sum_i \log(\frac{y_i(\theta))^{D_i}e^{-y_i(\theta))}}{\Gamma(D_i)}) = \sum_i D_i \log y_i(\theta) -y_i(\theta)
\end{equation}


```
from numpy.random import poisson
n_points = 20
x = linspace(0,300,n_points)
def y(x, theta=(110, 0.03), noise=0):
    ty = theta[0] * np.exp(-x * theta[1])
    if noise != 0:
        ty = [poisson(k)for k in ty] 
    return ty
D = y(x, noise=5)
scatter(x,D)
plot(x,y(x))
```


```
# use x, D from global scope
def P_cost(*w):
    return -np.sum(D*np.log(y(x,*w)) - y(x,*w))
fit_ml = minimize(P_cost, (105,0.05))
print fit_ml
```

      status: 2
     success: False
        njev: 40
        nfev: 172
         fun: -901.2355508250871
           x: array([  1.09423618e+02,   2.74569336e-02])
     message: 'Desired error not necessarily achieved due to precision loss.'
        hess: array([[  1.84852921e-05,  -1.78246740e-05],
           [ -1.78246740e-05,   1.78203684e-05]])
         jac: array([  2.28881836e-05,   2.34222412e-03])



```
# cheat, do this with global scope of y, x and D
def ssq_err(*w):
    return np.sum((D - y(x,*w))*(D - y(x,*w)).transpose())
fit_ssq = minimize(ssq_err,[105,0.05])
print fit_ssq
scatter(x, D)
plot(x,y(x,fit_ml.x), color="green")
plot(x,y(x,fit_ssq.x), color="red")
```

Need more?

http://people.physics.anu.edu.au/~tas110/Teaching/Lectures/L3/Material/Myung03.pdf

http://sciencehouse.files.wordpress.com/2013/07/gordonconf13.pdf

http://www.math.utah.edu/~levin/M5080/mle.pdf



```

```

c2e2b9341da211d9d244355e60c7dbf14df12f34
ipynb
Jupyter Notebook
max-likelihood/Maximum-Likelihood.ipynb
NoonienSoong/Data-Science-45min-Intros
8de3538945f7bad3e0b73554c64b8c38c1bdfdc2
Unlicense
2020-06-17T16:12:24.000Z
2020-06-17T16:12:24.000Z
max-likelihood/Maximum-Likelihood.ipynb
NoonienSoong/Data-Science-45min-Intros
8de3538945f7bad3e0b73554c64b8c38c1bdfdc2
Unlicense
max-likelihood/Maximum-Likelihood.ipynb
NoonienSoong/Data-Science-45min-Intros
8de3538945f7bad3e0b73554c64b8c38c1bdfdc2
Unlicense
Qwen/Qwen-72B
1. YES
2. YES