Hi there,
I’m quite new to probabilistic programming and pymc3. But what I’ve seen so far looks very cool. Thank you for providing this great toolbox!
Currently, I want to implement the Kennedy-O’Hagan framework in pymc3.
As I could see some people did the same thing before:
https://discourse.pymc.io/t/bayesian-model-calibration-with-gaussian-process/948
https://discourse.pymc.io/t/problem-with-multiprocessing-in-pymc3/1440
The setup is as follows (copied from kejzlarv’s post):
I have n observations z_i of the form
z_i = f(x_i, \theta) + g(x_i) + e_i,
where x_i are known imputs and \theta are unknown calibration parameters and e_i are iid error terms. We also have m model evaluations y_j of the form
y_j = f(x^*_j, \theta_j), where both x^*_j (different than x_i above) and \theta_j are known. Therefore, the data consists of all z_i and y_j. In the paper, Kennedy-O’Hagan model f, g using gaussian processes:
f \sim GP\{m_1(.,.), \Sigma_1[(.,.),(.,.)] \}
g \sim GP\{m_2(.), \Sigma_2[(.),(.)] \}.
Among other things, the goal is to get posterior samples for the unknow calibration iparameters \theta.
I tried to implement the framework and in the beginning it seemed easy but now I am confused. Would be great if you could give me some advise! The tricky part for me is the formulation of the marginal_likelihoods…
Here is what I did so far:
import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
from multiprocessing import freeze_support
import sys
import theano
import theano.tensor as tt
from mpl_toolkits.mplot3d import Axes3D
def physical_system(x):
return 0.65 * x / (1 + x / 5)
def observation(x):
return physical_system(x[:]) + np.random.normal(0,0.01,len(x))
def computational_system(input):
return input[:,0]*input[:,1]
if __name__ == "__main__":
freeze_support()
# observations with noise
x_obs = np.linspace(0,4,10)
y_real = physical_system(x_obs[:])
y_obs = observation(x_obs[:])
# computation model
N = 20
x_comp = np.random.uniform(0,4,N)[:,None]
t_comp = np.random.uniform(0,1,N)[:,None]
input_comp = np.hstack((x_comp,t_comp))
y_comp = computational_system(input_comp)
x_obs_shared = theano.shared(x_obs[:, None])
with pm.Model() as model:
theta = pm.Normal('theta', mu=0, sd=10)
noise = pm.HalfCauchy('noise',beta=5)
ls_1 = pm.Gamma('ls_1', alpha=1, beta=1, shape=2)
cov = pm.gp.cov.ExpQuad(2,ls=ls_1)
gp1 = pm.gp.Marginal(cov_func=cov)
gp2 = pm.gp.Marginal(cov_func=cov)
gp = gp1 + gp2
input_1 = tt.concatenate([x_obs_shared, tt.tile(theta, len(x_obs), ndim=2).T], axis=1)
f_0 = gp1.marginal_likelihood('f_0', X=input_comp, y=y_comp, noise=noise)
f_1 = gp1.marginal_likelihood('f_1', X=input_1, y=y_obs, noise=noise)
f = gp.marginal_likelihood('f', X=input_1, y=y_obs, noise=noise)
trace = pm.sample(1000,chains=1,step=pm.Metropolis())
pm.traceplot(trace)
plt.show()
Summary
This text will be hidden