Hi, I am implementing non-stationary GP regression using Gibbs kernel. The lengthscale and signal covariance are also GP.

```
with pm.Model() as model:
# Mean
mean = pm.gp.mean.Zero()
# Noise
noise_obs = pm.Normal('noise_obs', mu=0, sd=1)
# Prior for hyperparameters
l_l = pm.Normal('l_l', mu=0, sd=1)
l_sig = pm.Normal('l_sig', mu=0, sd=1)
signal_cov_l = pm.Normal('signal_cov_l', mu=0, sd=1)
signal_cov_sig = pm.Normal('signal_cov_sig', mu=0, sd=1)
# Covariance functions for Level 2
cov_l = signal_cov_l**2 * pm.gp.cov.ExpQuad(1, l_l)
cov_sig = signal_cov_sig**2 * pm.gp.cov.ExpQuad(1, l_sig)
# Level 2
gp_l = pm.gp.Marginal(mean_func=mean, cov_func=cov_l)
gp_sig = pm.gp.Marginal(mean_func=mean, cov_func=cov_sig)
# Covariance functions for Level 1
cov_obs = gp_sig**2 * pm.gp.cov.Gibbs(1,gp_l)
# Level 1
gp_obs = pm.gp.Marginal(mean_func=mean, cov_func=cov_obs)
# Level 0
y_obs = gp_obs.marginal_likelihood("y_obs", X=Xtrain[:,None], y=ytrain, noise=noise_obs)
marginal_post = pm.find_MAP()
```

I am getting the error: unsupported operand type(s) for ** or pow(): ‘Marginal’ and ‘int’ for covariance function in level 1. I do understand it.

My question is: Is there a workaround to this problem or will I have to make my own custom implementation? The model needs to be the way it is so I can’t change it for now.

Thanks and Regards!

Rohan