<jupyter_start><jupyter_text>Chapter 2<jupyter_code>import pymc3 as pm
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>2.1<jupyter_code>ways = np.array([0, 3, 8, 9, 0])
ways / ways.sum()<jupyter_output><empty_output><jupyter_text>2.2Analytical solution for the poster$$Pr(w \mid n, p) =  \frac{n!}{w!(n − w)!} p^w (1 − p)^{n−w}$$The probability of observing six W’s in nine tosses—under a value of p=0.5<jupyter_code>probabilites = []
for i in range(10):
    probabilites.append(stats.binom.pmf(i, n=9, p=0.5))
fig, ax = plt.subplots(1, 1, figsize=(6, 4), constrained_layout=True)

ax.bar(range(10), probabilites)
plt.xticks(range(10))

ax.set_title('Binomial pmf')
ax.set_xlabel('# successes')
ax.set_ylabel('Probability')

plt.show();
sum(probabilites)<jupyter_output><empty_output><jupyter_text>2.3 and 2.5Grid approximation to compute the posterior.<jupyter_code>def posterior_grid_approx(grid_points=5, k=6, n=9, prior_type='uniform', likelihood='binomial'):
    """Generate an estimate of the posterior distribution."""
    # define grid from 0, 1 with our number of equally spaced values
    p_grid = np.linspace(0, 1, grid_points)

    # compute likelihood at each point in the grid
    likelihood = stats.binom.pmf(k, n, p_grid)
    
    # define prior probability of p 
    if prior_type == 'uniform':
        prior = np.repeat(1, repeats=grid_points)
    elif prior_type == 'truncated':
        prior = (p_grid >= 0.5).astype(int)
    elif prior_type == 'double_exp':
        prior = np.exp(- 5 * abs(p_grid - 0.5))

    # compute unstandardised posterior -  product of likelihood and prior
    unstd_posterior = likelihood * prior

    # standardize the posterior, so it sums to 1
    posterior = unstd_posterior / unstd_posterior.sum()
    
    return p_grid, posterior<jupyter_output><empty_output><jupyter_text>2.3<jupyter_code>def plot_grid_approx(points, w, n, ax, *args, ** kwargs):
    p_grid, posterior = posterior_grid_approx(points, k=w, n=n, *args, ** kwargs)

    ax.set_title(f'{points} points')
    ax.set_xlabel('probability of water', fontsize=14)
    ax.set_ylabel('posterior probability', fontsize=14)

    ax.plot(p_grid, posterior, 'o-', label=f'success = {w}\ntosses = {n}')
w, n, = 6, 9
points = [3, 5, 10, 20, 1_000, 100_000]
fig, axes, = plt.subplots(nrows=2, ncols=3, figsize=(10, 6),
                          constrained_layout=True, sharex=True)

fig.suptitle('Grid Approximations')

for i, ax in enumerate(axes.flatten()):
    plot_grid_approx(points[i], w, n, ax)
    
plt.show();
fig, axes, = plt.subplots(nrows=2, ncols=3, figsize=(10, 6),
                          constrained_layout=True, sharex=True)

fig.suptitle('Grid Approximations')

for i, ax in enumerate(axes.flatten()):
    plot_grid_approx(points[i], w, n, ax, prior_type='truncated')
    
plt.show();
fig, axes, = plt.subplots(nrows=2, ncols=3, figsize=(10, 6),
                          constrained_layout=True, sharex=True)

fig.suptitle('Grid Approximations')

for i, ax in enumerate(axes.flatten()):
    plot_grid_approx(points[i], w, n, ax, prior_type='double_exp')
    
plt.show();
fig, axes, = plt.subplots(nrows=2, ncols=3, figsize=(10, 6),
                          constrained_layout=True, sharex=True)

fig.suptitle('Grid Approximations')

for i, ax in enumerate(axes.flatten()):
    plot_grid_approx(points[i], w=8, n=15, ax=ax)
    
plt.show();<jupyter_output><empty_output><jupyter_text>2.6Computing the posterior using the quadratic aproximation<jupyter_code>data = np.repeat((0, 1), (3, 6))
with pm.Model() as normal_aproximation:
    p = pm.Uniform('p', 0, 1)
    w = pm.Binomial('w', n=len(data), p=p, observed=data.sum())
    mean_q = pm.find_MAP()
    std_q = ((1/pm.find_hessian(mean_q, vars=[p]))**0.5)[0]
mean_q['p'], std_q
norm = stats.norm(mean_q, std_q)
prob = .89
z = stats.norm.ppf([(1-prob)/2, (1+prob)/2])
pi = mean_q['p'] + std_q * z 
pi<jupyter_output><empty_output><jupyter_text>2.8<jupyter_code># analytical calculation
w, n = 6, 9
x = np.linspace(0, 1, 100)
plt.plot(x, stats.beta.pdf(x , w+1, n-w+1),
         label='True posterior')

# quadratic approximation
plt.plot(x, stats.norm.pdf(x, mean_q['p'], std_q),
         label='Quadratic approximation')
plt.legend(loc=0, fontsize=13)

plt.title('n = {}'.format(n), fontsize=14)
plt.xlabel('Proportion water', fontsize=14)
plt.ylabel('Density', fontsize=14);<jupyter_output><empty_output><jupyter_text>Versions<jupyter_code>import sys, IPython, scipy, matplotlib
print(f"This notebook was created using:\nPython {sys.version[:5]}\nIPython {IPython.__version__}\nPyMC3 {pm.__version__}\nNumPy {np.__version__}\nSciPy {scipy.__version__}\nMatplotlib {matplotlib.__version__}\n")<jupyter_output>This notebook was created using:
Python 3.7.5
IPython 7.9.0
PyMC3 3.7
NumPy 1.17.3
SciPy 1.3.1
Matplotlib 3.1.1
{'fraction_non_alphanumeric': 0.15438308355573804, 'fraction_numerical': 0.030999794703346337, 'mean_word_length': 4.536363636363636, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '601308', 'license': 'MIT', 'n_tokens_mistral': 1983, 'n_tokens_neox': 1857, 'n_words': 400, 'path': 'Rethinking/Chp_02.ipynb', 'repo_name': 'David-Loughnane/resources', 'size': '4871'}