Rethinking_2/Chp_02.ipynb: Mistake in quadratic approximation ("Code 2.6")?

[thomas-haslwanter]

When I run the code for the quadratic approximation (attached below), I get a SD that is 4 times larger than the one in the book (0.64, instead of 0.16 as expected). Similarly, the "approximations" of the binomial with a gaussian have a much-too-large SD. Since I don't have an equation to compare, I assume that there is a mistake in the current code:

data = np.repeat((0, 1), (3, 6))
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()

    p_value = normal_approximation.rvs_to_values[p]
    p_value.tag.transform = None
    p_value.name = p.name

    std_q = ((1 / pm.find_hessian(mean_q, vars=[p])) ** 0.5)[0]

# display summary of quadratic approximation
print("Mean, Standard deviation\np {:.2}, {:.2}".format(mean_q["p"], std_q[0]))

produces on my computer

Mean, Standard deviation
p 0.67, 0.64

instead of the expected SD=0.16

I have been running the code on Win64, Python 3.10, PyMC 4

[LukeBenjaminT]

I think I hit the same problem. I'm running PyMC 5.0.2. I dug online and got to this. Worked for me but doesn't look efficient (two models).

with pm.Model() as m:
    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()
    
#calculating std_dev without trasforming... I think it was in the log space. Unbounded???
    
with pm.Model() as untransformed_m:
    p = pm.Uniform("p", 0, 1,transform=None)  # uniform priors
    w = pm.Binomial("w", n=len(data), p=p, observed=data.sum(),transform=None)  # this seems to be the difference
    std_q = ((1 / pm.find_hessian(mean_q, vars=[p])) ** 0.5)[0]
    
print("Mean, Standard deviation\np {:.2}, {:.2}".format(mean_q["p"], std_q[0]))

[jakobpete]

Had the same issue! Did you find a solution for that?
Very Best!

[LukeBenjaminT]

Something like the below works for me. Still no fundamental solutions. This is trying to work with actual samples from the posterior (there is some variability because they are just samples). The original method I think is faster and fancier but running into difficulties with how its handled behind the scenes.

data = np.repeat((0, 1), (3, 6))
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]
    idata = pm.sample() #taking actual samples here with default values
# display summary of quadratic approximation

print("Mean, Standard deviation\np {:.2}, {:.2}".format(mean_q["p"], std_q[0]))
print("SAMPLE:\nMean, Standard deviation\np {:.2}, {:.2}".format(idata.posterior['p'].mean(),idata.posterior['p'].std()))

I'm basing this on a discussion in the pymc forum.

[peterdizo]

This issue seems to describe what is the problem: pymc-devs/pymc#5443
And this code snippet gave me the correct answers (using PyMC 5.1.1): pymc-devs/pymc#5443 (comment)

It seems like we need to use .rvs_to_tranforms to access and remove the log transform, instead of accessing it through .tag.transform
Just be careful, as this will probably also change the underlying model

Hope it helps

[ivaquero]

Hi, @LukeBenjaminT. Your solution works. Could you please provide the link of your post on the forum?

Something like the below works for me. Still no fundamental solutions. This is trying to work with actual samples from the posterior (there is some variability because they are just samples). The original method I think is faster and fancier but running into difficulties with how its handled behind the scenes.

data = np.repeat((0, 1), (3, 6))
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]
    idata = pm.sample() #taking actual samples here with default values
# display summary of quadratic approximation

print("Mean, Standard deviation\np {:.2}, {:.2}".format(mean_q["p"], std_q[0]))
print("SAMPLE:\nMean, Standard deviation\np {:.2}, {:.2}".format(idata.posterior['p'].mean(),idata.posterior['p'].std()))

I'm basing this on a discussion in the pymc forum.