# Maximum likelihood Estimation (MLE)
based on http://python-for-signal-processing.blogspot.com/2012/10/maximum-likelihood-estimation-maximum.html
## Simulate coin flipping
- [Bernoulli distribution](https://en.wikipedia.org/wiki/Bernoulli_distribution) 
is the probability distribution of a random variable which takes the value 1 with probability $p$ and the value 0 with probability $q = 1 - p$
- [scipy.stats.bernoulli](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.bernoulli.html)


```python
import numpy as np
from scipy.stats import bernoulli 

np.random.seed(123456789)

p_true = 1/2 # this is the value we will try to estimate from the observed data
fp = bernoulli(p_true)

def sample(n=10):
    """
    simulate coin flipping
    """
    return fp.rvs(n) # flip it n times

xs = sample(100) # generate some samples
```

## Find maximum of Bernoulli distribution
Single experiment
$$\phi(x) = p ^ {x} * (1 - p) ^ { 1 - x }$$
Series of experiments
$$\mathcal{L}(p|x) = \prod_{i=1}^{n} p^{x_{i}}*(p-1)^{1-x_{i}}$$
### Hints
- [sympy.diff()](http://docs.sympy.org/dev/modules/core.html#sympy.core.function.diff)
- [sympy.expand()](http://docs.sympy.org/dev/modules/core.html#sympy.core.function.expand)
- [sympy.expand_log()](http://docs.sympy.org/dev/modules/core.html#sympy.core.function.expand_log)
- [sympy.solve()](http://docs.sympy.org/dev/modules/core.html#sympy.core.function.solve)
- [sympy.symbols()](http://docs.sympy.org/dev/modules/core.html#symbols)
- [sympy gotchas](http://docs.sympy.org/dev/tutorial/gotchas.html)


```python
import sympy
from sympy.abc import x

p = sympy.symbols('p', positive=True)
phi = p ** x * (1 - p) ** (1 - x)
L = np.prod([phi.subs(x, i) for i in xs]) # objective function to maximize
log_L = sympy.expand_log(sympy.log(L))
sol = sympy.solve(sympy.diff(log_L, p), p)[0]
```


```python
import matplotlib.pyplot as plt

x_space = np.linspace(1/100, 1, 100, endpoint=False)

plt.plot(x_space,
         list(map(sympy.lambdify(p, log_L, 'numpy'), x_space)),
         sol,
         log_L.subs(p, sol),
         'o',
         p_true,
         log_L.subs(p, p_true),
         's',
        )
plt.xlabel('$p$', fontsize=18)
plt.ylabel('Likelihood', fontsize=18)
plt.title('Estimate not equal to true value', fontsize=18)
plt.grid(True)
plt.show()
```

## Empirically examine the behavior of the maximum likelihood estimator 
- [evalf()](http://docs.sympy.org/dev/modules/core.html#module-sympy.core.evalf)


```python
def estimator_gen(niter=10, ns=100):
    """
    generate data to estimate distribution of maximum likelihood estimator'
    """
    x = sympy.symbols('x', real=True)
    phi = p**x*(1-p)**(1-x)
    for i in range(niter):
        xs = sample(ns) # generate some samples from the experiment
        L = np.prod([phi.subs(x,i) for i in xs]) # objective function to maximize
        log_L = sympy.expand_log(sympy.log(L)) 
        sol = sympy.solve(sympy.diff(log_L, p), p)[0]
        yield float(sol.evalf())
    
entries = list(estimator_gen(100)) # this may take awhile, depending on how much data you want to generate
plt.hist(entries) # histogram of maximum likelihood estimator
plt.title('$\mu={:3.3f},\sigma={:3.3f}$'.format(np.mean(entries), np.std(entries)), fontsize=18)
plt.show()
```

## Dynamic of MLE by length sample sequence


```python
def estimator_dynamics(ns_space, num_tries = 20):
    for ns in ns_space:
        estimations = list(estimator_gen(num_tries, ns))
        yield np.mean(estimations), np.std(estimations)
    
ns_space = list(range(10, 100, 5))
entries = list(estimator_dynamics(ns_space))
entries_mean = list(map(lambda e: e[0], entries))
entries_std = list(map(lambda e: e[1], entries))

plt.errorbar(ns_space, entries_mean, entries_std, fmt='-o')
plt.show()
```


```python

```

4940682496d40fe493ef3f22d66b75e23e264263
ipynb
Jupyter Notebook
mle.ipynb
hyzhak/mle
257d8046a950b7381052cc56d9931cf98aeb0a5c
MIT
2017-10-22T09:29:36.000Z
2017-10-22T09:29:36.000Z
mle.ipynb
hyzhak/mle
257d8046a950b7381052cc56d9931cf98aeb0a5c
MIT
mle.ipynb
hyzhak/mle
257d8046a950b7381052cc56d9931cf98aeb0a5c
MIT
2019-01-23T04:46:01.000Z
2020-04-21T18:38:49.000Z
Qwen/Qwen-72B
1. YES
2. YES