Maximum likelihood

In [3]: 
import numpy
step_size = 0.0001
x= numpy.arange(0,1, step_size)
len(x)
Out[3]:
10000
In [6]: 
likelihood = [ p * p * (1-p) for p in x]
In [7]: 
from matplotlib import pyplot
%matplotlib inline
pyplot.plot(x,likelihood)
Out[7]:
[<matplotlib.lines.Line2D at 0x10a08b650>]
In [8]: 
index = numpy.argmax(likelihood)
mle_p = x[index]
mle_p
Out[8]:
0.66670000000000007

Bayesian Estimation

We use Bayes' Rule to infer the probability distribution over the parameter $\theta$. Now we use the posterior probability to tell us how good the $\theta$. $$P(\theta|D) = \frac{P(D|\theta)P(\theta)}{P(D)}$$

How do we compute $P(D)$? We can approximate it with $\sum_{\theta'} P(D|\theta') P(\theta')$

In [52]: 
from scipy.stats import beta
prior_y = beta.pdf(x, 2, 2)
pyplot.plot(x, prior_y)
Out[52]:
[<matplotlib.lines.Line2D at 0x10c203d10>]
In [53]: 
normalization = 0.0
for l, p in zip(likelihood, prior_y):
    normalization += l * p
posterior = [l * p / normalization for l, p in zip(likelihood, prior_y)]
pyplot.plot(x, posterior)
Out[53]:
[<matplotlib.lines.Line2D at 0x10c27d750>]
In [54]: 
index = numpy.argmax(posterior)
x[index]
Out[54]:
0.5999999999999502
In [55]: 
float((2 + 1)) / (2+1 + 1 +1)
Out[55]:
0.6
In [58]: 
index = numpy.argmax(likelihood)
x[index]
Out[58]:
0.6666999999999429
In [59]: 
2.0/3
Out[59]:
0.6666666666666666