Hi! I’m new to implementing Gaussian Processes in PyMC, and I have a follow-up to a question asked earlier this year on implementing a GP as an emulator for simulated data.
Instead of creating an emulator, @lilianschuster and I would like to use a GP, fitted to simulated data, as a latent model for modeling observed inputs when we only have one response observation.
To do this, I write up a latent GP model that fits the simulated inputs to the simulated response, following the example in the gallery:
with pm.Model() as model:
#define the covariance parameters and function
input_dim = 2 #number of fit variables
cov = pm.gp.cov.ExpQuad(input_dim, ls=0.1)
#instantiate latent GP
gp = pm.gp.Latent(cov_func=cov)
# np.shape(X_sim) : (663, 2)
f = gp.prior('latent_gp', X=X_train,
jitter=1e-7)
gp_sigma = pm.HalfNormal('gp_sigma', sigma=2.0)
nu = 1 + pm.Gamma('nu', alpha=2, beta=0.1)
# np.shape(X_sim) : (663, 2)
y_ = pm.StudentT('y', mu=f, lam=1.0 / gp_sigma, nu=nu,
observed=y_train)
latent_trace = pm.sample(1000, chains=2, target_accept=0.9)
The above code runs smoothly, but I get stuck on the rest of the implementation. After this I try to run gp.conditional
on two observationally informed distributions for each input dimension.
with model:
#define the "observational" distributions for inputs to gp.conditional
dim1 = pm.TruncatedNormal('dim1_fit', mu=0, sigma=2,
lower=-8, upper=8)
dim2 = pm.TruncatedNormal('dim2_fit',
mu=250,
sigma=200,
lower=1, upper=1000)
Xnew = np.array([dim1, dim2])
#instantiate the conditional fit of GP to new "data"
gp_mu = gp.conditional('conditional_fit', Xnew=Xnew)
Trying to fit gp.conditional to Xnew composed of two distribution tensors throws errors, and I’m not sure how to fix this. Even if I try to reshape or add and index of [:,None]
to the numpy array, I get the following error message:
IndexError: index 1 is out of bounds for axis 1 with size 1
Additionally, once I get gp.conditional
working, I would like to pass it as the center function to sample from a Normal distribution centered on the output of gp.conditional
. I think the following code would be able to accomplish this, but I’m not sure how to sample from it while incorporating the previous trace?
with model:
#sigma_obs and y_obs are both lists with len=1
sigma = pm.Data('sigma', sigma_obs)
observed = pm.Data('observed', y_obs)
log_l = pm.Normal('log_mb', mu=gp_mu, sigma=sigma,
observed=observed)
Thank you so much in advance for any and all help on this!