Thanks for the reply, here’s the basic of it. The goal is to come up with a Gaussian approx to the marginal posterior P(θ|x) = ∫dz P(θ,z|x). The following is the mean (the covariance of the estimate is different but uses the identical pieces):
θ = # initial guess for θ
H = # some guess for Hessian of θ -> logP(θ|x)
while norm(θ - θlast) < θtol:
# a bunch of simulated x's generated from P(x,z|θ)
x_sims = [sample_prior_predictive(θ).x for i in 1:nsims]
# zMAP maximizes the function z -> logP(x,z|θ)
zMAP_data = zMAP(x, θ)
zMAP_sims = [zMAP(x_sim, θ) for x_sim in x_sims]
# ∇θ_logLike is gradient of θ -> logP(x,z|θ)
g_data = ∇θ_logLike(zMAP_data, x, θ)
g_sims = [∇θ_logLike(zMAP_sim, x_sim, θ) for (zMAP_sim, x_sim) in zip(zMAP_sims,x_sims)]
# gradient of θ -> logP(θ)
g_prior = ∇θ_logPrior(θ)
θlast = θ
θ -= H \ (g_data - mean(g_sims) + g_prior)
The one thing I’ve already figured out what its called in PyMC is sample_prior_predictive, although I haven’t found the right way to condition it on θ. The gradient of the likelihood with respect to the latent z would be used in zMAP(), then I also need ∇θ_logLike() and ∇θ_logPrior() to evaluate the gradient of the likelihood and prior w.r.t. θ. That’s everything.
Yes sorry, I guess my question really boils down to, would I expect that the code I write here is compatible with both the current latest release and the future PyMC v4, or is the relevant API changing such that I’d need to update / support both versions?
And thanks for the SMC example, am taking a look now.