I’m not 100% sure yet, but it seems like it used to compute the hessian (or actually negative of the hessian, as I just noticed…) on the untransformed space, but now it computes it on the transformed space. There is also the question if the jacobian terms should be included here, they are not included in the find_MAP call, but are included in the find_hessian call…
It does seem however that the derivatives are correct on the transformed space at least:
data = np.repeat([0, 1], (10, 3))
with pm.Model() as normal_aproximation:
p = pm.Beta('p', 1., 1.)
w = pm.Binomial('w',n=1, p=p, observed=data)
mean_q = pm.find_MAP()
from scipy import optimize
jac = True
logp = normal_aproximation.logp(jacobian=jac)
func = aesara.function(normal_aproximation.value_vars, logp)
d2logp = normal_aproximation.compile_d2logp(jacobian=jac)
hesse = d2logp({"p_logodds__" : mean_q["p_logodds__"]})
x = np.linspace(mean_q["p_logodds__"] - 0.5, mean_q["p_logodds__"] + 0.5)
y = np.array([func(val) for val in x])
p = mean_q["p_logodds__"]
dlogp = normal_aproximation.compile_dlogp(jacobian=jac)
dy = dlogp({"p_logodds__": p})
ddy = d2logp({"p_logodds__": p}).ravel()
plt.plot(x, y)
plt.axvline(mean_q["p_logodds__"])
plt.plot(x, func(p) + (x - p) * dy)
plt.plot(x, func(p) + (x - p) * dy - (x - p) ** 2 * ddy / 2)
I’m teaching a Bayesian course based partly on @aloctavodia’s Bayesian Analysis in Python, but I’m updating the notebooks to use PyMC instead of PyMC3.
In Chapter 8, he uses the MAP and the Hessian in his code to demonstrate the Gaussian approximation and I needed to update it for the students.
More generally, I can imagine that there is use in getting a Gaussian approximation of the untransfomed mode. Gaussian approximations are quite generally useful.
That’s good to know. I could just tell them a couple of words about bambi (not a bad idea in any case) and then use this function.
On the other hand, looking at the source, it seems like this code is calling find_hessian on the untransformed parameter. Why can’t this work just the same outside this function (just changing the name to the result of get_untransformed_name.)
@aloctavodia I think that snippet suffers the same “problem” (assuming you don’t want the hessian on the transformed space) as in this issue with PyMC v4.
Thanks for the example. I think we need to improve our API to allow to easily define expressions on the transformed and untransformed space. I don’t think the old API was that great though…
I think the easiest solution for now is to recreate the PyMC model, while passing transform=None for every variable, and then call find_hessian on that model with the results from find_MAP that you got in the original model.
with pm.Model() as m:
x = pm.Normal("x")
y = pm.HalfNormal("y")
pm.Normal("llike", x, y, observed=[1, 2, 3])
_map = find_MAP()
with pm.Model() as untransformed_m:
x = pm.Normal("x", transform=None) # Not needed, but doesn't hurt
y = pm.HalfNormal("y", transform=None)
pm.Normal("llike", x, y, observed=[1, 2, 3], transform=None) # Not needed, but doesn't hurt
find_hessian(_map)`
Just to clarify that Bambi performs the computations on the unbounded space and then returns the results on the original scale. I agree that showing the internals to students could be too much.
From the code, yes, they are from a gaussian fit to the transformed (unconstrained space). But when they later call idata = _posterior_samples_to_idata(samples, self.model), the sampled values are transformed back to the constrained space.