Rethinking_2/Chp_02.ipynb: Mistake in quadratic approximation ("Code 2.6")? about pymc-resources HOT 5 OPEN

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

from pymc-resources.

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]))

from pymc-resources.

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

from pymc-resources.

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.

from pymc-resources.

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.

from pymc-resources.