Sample prior predictive for Gaussian process model

Hi there!
Is it possible to get prior predictive samples for Gaussian process models? I’d like something like this (MWE taken from the docs):

x = np.linspace(0,1,50)
y = 0.5 *x + np.random.randn(50) * 0.01

with pm.Model() as gp_model:
    cov_func =, ls=0.1)
    gp =
    sigma = pm.HalfCauchy("sigma", beta=5)
    y_ = gp.marginal_likelihood("y", X=x[:,None], y=y, noise=sigma)
    y_star = gp.conditional("y_star", x[:,None])

with gp_model:
    prior = pm.sample_prior_predictive()

to create prior predictive samples for y_star? However atm this always returns an array of NaNs.
Thank you!

It seems to be a bug? @bwengals?

@junpenglao @bwengals So is the way I specify the prior predictive correct? Happy to attempt to fix this if this is a bug!

To sample the prior predictive correctly I think you’d need to use gp.Latent. The gp.Marginal implementation marginalizes out the actual GP, which makes it more efficient, but the actual GP f isn’t explicitly included in the model.

1 Like