Bayesian GP PyMC3 PPC Problem

gp_ppc_broken.ipynb

I was able to get the PPC to start sampling with this (hacky) change:

with pm.Model() as model:
    tau = pm.StudentTpos('tau', nu=4)
    l = pm.StudentTpos('l', nu=4)
    K = pm.Deterministic('K', tau**2 * tt.exp(-x_diffs_shared**2 / l**2))

    log_sigma = pm.Uniform('log_sigma', -2., 2.)
    sigma = pm.Deterministic('sigma', tt.power(10., log_sigma))

    cov = pm.Deterministic('cov', K + sigma**2 * tt.eye(x_diffs_shared.shape[0]))
    y_obs = pm.MvNormal('y_obs', tt.zeros_like(cov)[0, :], cov, observed=y)

However, the output does not seem correct:

Actually seems to work somewhat:

Interesting! I'll see if I can take this any further.

I am looking for an example of Gaussian Process classification in PyMC3 with PPC sampling. The only example of Gaussian Process for classification that is provided in the documentation[1] does not provide an example of PPC sampling. When I use that example with PPC I receive the following error: TypeError: object of type 'NoneType' has no len(). Here is the code that I use at the moment:

# ******************** Modelling *******************************
print("Modelling ...")
model_input = theano.shared(np.array(X_train))
model_output = theano.shared(y_train.values)
n, p = X_train.shape
with pm.Model() as gp_model:
    # prior over alpha
    alpha = pm.Normal('alpha', mu=0.5, sd=5 ** -5, shape=(1))
    # prior over betas, notice that each beta has its own distribution
    betas = pm.Normal('betas', mu=0, sd=10 ** -5, shape=(1, p))
    # Logistic Regression model, you have alpha and beta * x^i
    f_transform = pm.invlogit(alpha + tt.sum(betas * model_input, 1))

    # Binomial likelihood
    y = pm.Binomial('y', observed=model_output, n=np.ones(n), p=f_transform, 
    shape=n)
    # Interpolate function values using noiseless covariance matrix
    L = tt.slinalg.cholesky(K_stable)
    f_pred = pm.Deterministic('f_pred', tt.dot(K_s.T, tt.slinalg.solve(L.T, 
    tt.slinalg.solve(L, f_transform))))

with gp_model:
    if sampler == "1":
        print("Mini_batch does not work ATM!")

    if sampler == "2":
        step = pm.NUTS()
        trace = pm.sample(advi_iter, step)

    # save the trace plot
    pm.traceplot(trace[n_burn_in:])
    plt.savefig(output_dir + "trace_" + str(trace_draws) + "_draw.png")
    print("Test the model ...")
    # Replace shared variables with testing set
    model_input.set_value(X_test)
    model_output.set_value(y_test)
    n = X_test.shape[0]
    print("OK ....")
    # Create posterior predictive samples from ADVI
    ppc = pm.sample_ppc(trace, model=gp_model, samples=ppc_samples)
    # Use probability of > 0.5 to assume prediction of class 1
    pred = ppc['f_pred'].mean(axis=0) > 0.5

[1]https://pymc-devs.github.io/pymc3/notebooks/GP-slice-sampling.html#Gaussian-Process-Classification