I am getting closer to a solution but I am still having some difficulties with syntax. I have the following definitions in python

```
nvar = 4
with pm.Model() as surrogate_model:
ls = pm.HalfNormal('l',2./3,shape=nvar)
eta = pm.HalfNormal('sigma',1.0,shape=nvar)
cov = eta**2.*pm.gp.cov.ExpQuad(nvar,ls_inv=ls)
gp = pm.gp.Marginal(cov_func=cov)
y_ = gp.marginal_likelihood('y',X=X0.T,y=errors,noise=sigma)
mp = pm.find_MAP()
```

where X0 has a shape of (30,4). My intention is to create a gaussian process model which has a scalar output and four parameter input. This code executes without any difficulty.

Now I am trying to sample this gaussian process using random variables so I define

```
with surrogate_model:
a_ = pm.Normal('a',params[0], 0.1*abs(params[0]))
b_ = pm.Normal('b',params[1], 0.1*abs(params[1]))
c_ = pm.Normal('c',params[2], 0.1*abs(params[2]))
d_ = pm.Normal('d',params[3], 0.1*abs(params[3]))
sigma2 = pm.HalfNormal('sigma2',1.0)
params_ = tt.stack([a_,b_,c_,d_]).reshape([1,4])
f_cond = gp.conditional("f_cond",Xnew=params_)
error = pm.Normal('error',mu=f_cond,sigma=sigma2,observed=0.)
trace = pm.sample(1000,njobs=1)
```

This code fails at the definition of f_cond where I get

```
ValueError: Input dimension mis-match. (input[0].shape[1] = 1, input[1].shape[1] = 4)
```

It seems like what may be happening is that I’m not defining my gp to have four inputs at all. However when I look at

```
print(gp.cov_func.active_dims)
[0 1 2 3]
```

which seems to show that I do have four degrees of freedom. I tried the following to see if I could reproduce my training data

```
with surrogate_model:
f_pred = gp.conditional('f_pred',Xnew=X0.T)
with surrogate_model:
pred_samples = pm.sample_ppc([mp],vars=[f_pred],samples=2000)
```

and I received the error

```
ValueError: incompatible dimensions
```

Any advice would be greatly appreciated.