I am trying to experiment with fully Bayesian inference in GP classification. The use case is to draw samples of the hypers (theta) and latent f jointly. I want to implement a scheme where the sampling is interleaved so the pseudo code looks something like this:
Step 1-> Sample some hypers, freeze those
Step 2-> Sample from the conditional p( f | theta, y).
The idea is to navigate in the constrained f space (high dimensional) using different rules than the low dimensional theta space.
I wrote some trial code using the compound step but not sure if it is giving me f samples from the conditional.
with pm.Model() as model: # covariance function log_l = pm.Normal('log_l', mu=0, sigma=3) log_s = pm.Normal('log_s', mu=0, sigma=3) l = pm.Deterministic('l', pm.math.exp(log_l)) s = pm.Deterministic('s', pm.math.exp(log_s)) cov = s**2 * pm.gp.cov.ExpQuad(1, l) gp = pm.gp.Latent(cov_func=cov) # make gp prior f = gp.prior("f", X=x[:,None], reparameterize=True) # logit link and Bernoulli likelihood p = pm.Deterministic("p", pm.math.invlogit(f)) y_ = pm.Bernoulli("y", p=p, observed=y) step_theta = pm.Metropolis([log_l, log_s]) step_latent = pm.NUTS([f]) trace_hmc = pm.sample(draws=100, tune=200, chains=1, step=[step_theta, step_latent])