EuroCoinProblem¶
Import¶
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sns
import empiricaldist
from empiricaldist import Pmf, Distribution
from ipywidgets import interact, interactive, fixed
euro = Pmf.from_seq(range(101))
euro
| probs | |
|---|---|
| 0 | 0.009901 |
| 1 | 0.009901 |
| 2 | 0.009901 |
| 3 | 0.009901 |
| 4 | 0.009901 |
| ... | ... |
| 96 | 0.009901 |
| 97 | 0.009901 |
| 98 | 0.009901 |
| 99 | 0.009901 |
| 100 | 0.009901 |
101 rows × 1 columns
Custom Likelihood¶
def likelihood_euro(data, hypo):
# print(data, hypo)
x = hypo/100
# outcome = data
if data == 'H':
return x
else:
return 1-x
euro = Pmf.from_seq(range(101))
euro.update(likelihood_euro, 'H')
euro.plot()
euro.update(likelihood_euro, 'H')
euro.plot()
euro.update(likelihood_euro, 'T')
euro.plot()
euro.update(likelihood_euro, 'H')
euro.plot()
outcome = 'HHHHHHTTT'
euro = Pmf.from_seq(range(101))
for data in outcome:
euro.update(likelihood_euro, data)
euro.plot()
outcome = 'H'*140+'T'*110
euro = Pmf.from_seq(range(101))
for data in outcome:
euro.update(likelihood_euro, data)
euro.plot()
euro.mean()
55.95238095238095
np.percentile(euro, [10,90])
array([1.62410199e-96, 3.48751559e-02])
#MAP - Maximum aposteori probability
euro.max_prob()
56
euro.quantile(0.5)
array(56.)
euro.credible_interval(0.9)
array([51., 61.])
Note
Assuming prior and credible is true
euro.mean??
np.sum(euro.ps*euro.qs)
55.95238095238095
euro.max_prob??
euro.idxmax()
56
euro.credible_interval??
Effect of Different Priors¶
def TriangularPrior():
"""Makes a Suite with a triangular prior
"""
suite = Pmf(name='triangle')
for x in range(0,51):
suite[x] = x
for x in range(51, 101):
suite[x] = 100-x
suite.normalize()
return suite
euro1 = Pmf.from_seq(range(101), name='uniform')
euro1.plot()
euro2 = TriangularPrior()
euro2.plot()
def decorate_euro(title):
"""Labels the axes.
title: string
"""
plt.xlabel('Probability of heads')
plt.ylabel('PMF')
plt.title(title)
euro1 = Pmf.from_seq(range(101), name='uniform prior')
euro1.plot(color='steelblue')
euro2 = TriangularPrior()
euro2.plot(color='red', label='triangular prior')
plt.legend()
decorate_euro("Prior Distribution")
outcome = 'H'*140+'T'*110
euro1 = Pmf.from_seq(range(101), name='uniform')
euro1.plot()
euro2 = TriangularPrior()
euro2.plot()
for data in outcome:
euro1.update(likelihood_euro, data)
euro2.update(likelihood_euro, data)
euro1.plot()
euro2.plot()
plt.legend()
decorate_euro("Posterior Distribution")
outcome = 'H'*14+'T'*11
euro1 = Pmf.from_seq(range(101), name='uniform')
euro1.plot()
euro2 = TriangularPrior()
euro2.plot()
for data in outcome:
euro1.update(likelihood_euro, data)
euro2.update(likelihood_euro, data)
euro1.plot()
euro2.plot()
plt.legend()
decorate_euro("Posterior Distribution")
outcome = 'H'*4+'T'*1
euro1 = Pmf.from_seq(range(101), name='uniform')
euro1.plot()
euro2 = TriangularPrior()
euro2.plot()
for data in outcome:
euro1.update(likelihood_euro, data)
euro2.update(likelihood_euro, data)
euro1.plot()
euro2.plot()
plt.legend()
decorate_euro("Posterior Distribution")
outcome = 'H'*140+'T'*110
euro1 = Pmf.from_seq(range(101), name='uniform')
euro1.plot()
euro2 = TriangularPrior()
euro2.plot()
for data in outcome:
euro1.update(likelihood_euro, data)
euro2.update(likelihood_euro, data)
euro1.plot()
euro2.plot()
plt.legend()
decorate_euro("Posterior Distribution")
euro1.mean() -euro2.mean()
0.20888151378589725
euro1.credible_interval(0.9)
array([51., 61.])
euro2.credible_interval(0.9)
array([51., 61.])
Swammping the priors
More or less agree for posterior
if you start with probability zero , you will always get zero
Avoid assigning probability to zero on everything
Estimated proportions
euro2.plot(color='red', label='uniform posterior')
plt.legend()
<matplotlib.legend.Legend at 0x7f98ad26e820>
list(range(51))
[0,
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12,
13,
14,
15,
16,
17,
18,
19,
20,
21,
22,
23,
24,
25,
26,
27,
28,
29,
30,
31,
32,
33,
34,
35,
36,
37,
38,
39,
40,
41,
42,
43,
44,
45,
46,
47,
48,
49,
50]
Interactive¶
def likelihood_euro(data, hypo):
# print(data, hypo)
x = hypo/100
# outcome = data
if data == 'H':
return x
else:
return 1-x
def TriangularPrior(t_peak, n_range):
"""Makes a Suite with a triangular prior
"""
suite = Pmf(name='triangle')
for x in range(0, t_peak+1):
suite[x] = x
for x in range(t_peak+1, n_range):
suite[x] = (n_range-1)-x
suite.normalize()
return suite
def compare(n_h = 140, n_t=110, t_peak=50, n_range=101):
outcomes = 'H'*n_h+'T'*n_t
euro1 = Pmf.from_seq(range(n_range), name='uniform')
euro2 = TriangularPrior(t_peak, n_range)
euro1.plot(color='steelblue', label='uniform prior', linestyle='--')
euro2.plot(color='red', label='triangular prior', linestyle='--')
for data in outcomes:
euro1.update(likelihood_euro, data)
euro2.update(likelihood_euro, data),
euro1.plot(color='steelblue', label='uniform posterior')
euro2.plot(color='red', label='triangular posterior')
plt.legend()
interactive(compare, n_range=fixed(101), t_peak=fixed(50), n_h=(0,140), n_t=(0,140))