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 = pm.gp.cov.ExpQuad(1, ls=0.1)
    gp = pm.gp.Marginal(cov_func=cov_func)
    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!