Hi all,

I feel stupid because I can’t seem to find how to do such a basic thing, but I’d like to **evaluate the log density log(Pr(D, theta)) at a given theta** for the model at hand (in which D is the observed data).

FWIU this computation is at the heart of what PyMC does for running HMC, this is why I’m surprised I can’t find how to compute it.

**Motivation:** the general context for this is implementing the Laplace Method for approximating the posterior, in particular as a way of estimating Model Evidence for model comparison, in “complex-observations” situations where Information Criteria are not directly applicable. I’m already able to compute the theta_MAP and the corresponding Hessian, all I’m lacking now is a robust way of resolving `log(Pr(D, theta_MAP))`

; I guess I could use the `optres.fun`

returned by `scipy.optimize.minimize`

, but that seems like relying on a fragile implementation detail of PyMC.