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.

# HMC Model

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