Shape of posterior predictive samples of GP

I have the following to model my univariate time series data:

with pm.Model() as model:
    ℓ = pm.Gamma("ℓ", alpha=2, beta=1)
    η = pm.HalfCauchy("η", beta=5)

    cov = η**2 * pm.gp.cov.Matern52(1, ℓ)
    gp = pm.gp.Latent(cov_func=cov)
    f = gp.prior("f", X = train_x)

    σ = pm.HalfCauchy("σ", beta=5)
    ν = pm.Gamma("ν", alpha=2, beta=0.1)
    y_ = pm.StudentT("y", mu=f, lam=1.0/σ, nu=ν, observed=train_y)

    trace = pm.sample(50, tune=10, cores=4)

However, when I sample from the posterior by using the block below, the shape I get for pred_samples[‘y’] is (50, 48, 48):

with model:
    f_pred = gp.conditional("f_pred", train_x)
    # Sample from the GP conditional distribution
    pred_samples = pm.sample_posterior_predictive(trace, vars=[f_pred, y_], samples=50)

So I don’t quite understand why each sample is 48 x 48. Shouldn’t it be simple 50x48 (which is the shape of pred_samples['f_pred']). I have included a simple test script here in this colab notebook which runs the above.

Any thoughts/ comments would be highly appreciated.

@lucianopaz could you take a look at this example? Seems like there is a broadcasting error?

@sachinruk, looks like a broadcasting error. Does it work if you try train_y = np.ones(48) instead of train_y = np.ones(48)[:, None]? I’ll be able to take a close look during the week

Perfect! That fixed the problem.