Based on the PyMC port of the book “Rethinking Statistics - e2”, I have been trying to reproduce the normal approximation to a binomial. I get the correct result for the mean, but the SD is off by a factor of 4 - and I cannot find the problem. This is also relevant for the
https://github.com/pymc-devs/pymc-resources/blob/main/Rethinking_2/Chp_02.ipynb
import matplotlib.pyplot as plt
import numpy as np
import pymc as pm
import scipy.stats as stats
## Generate some data
throws = 9
successes = 6
failures = throws - successes
data = np.repeat((0, 1), (failures, successes))
print(data)
## The model -------------------------------------------------------------------------
with pm.Model() as normal_approximation:
p = pm.Uniform("p", 0, 1) # uniform priors
w = pm.Binomial("w", n=len(data), p=p, observed=data.sum()) # binomial likelihood
mean_q = pm.find_MAP()
std_q = (1 / ( pm.find_hessian(mean_q, vars=[p])) ** 0.5)[0][0]
## End of model ---------------------------------------------------------------------
## display summary of quadratic approximation
print("Mean, Standard deviation\np {:.2}, {:.2}".format(mean_q["p"], std_q))
# analytical calculation
w, n = successes, throws
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)
plt.title(f"n = {n}")
plt.xlabel("p");
results in
