# Problems with finding normal approximation to binomial

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)
## 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")

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

The problem is the hessian is defined on the transformed space in versions of PyMC>=4. More context here: Numerical differences between 3.11.4 and 4.0.0b · Issue #5443 · pymc-devs/pymc · GitHub

1 Like

Hi Ricardo,