使用 PyMC3 进行贝叶斯校准，Kennedy O'Hagan

Question

我对概率编程和 pymc3 还是很陌生... 目前，我想在 pymc3 中实现 Kennedy-O'Hagan 框架。

根据肯尼迪和奥哈根的论文设置如下：

我们有 n 个观察结果 z_i 的形式

z_i = f(x_i , theta) + g(x_i) + e_i,

其中 x_i 是已知输入，theta 是未知校准参数，e_i 是 iid 错误术语。我们还有 m 模型评估 y_j 的形式

y_j = f(x^'_j, theta_j) ，其中 x^'_j （不同于上面的 x_i）和 theta_j 是已知的。因此，数据由所有 z_i 和 y_j 组成。在论文中，Kennedy-O'Hagan 使用高斯过程对 f、g 建模：

f ~ GP{m₁ (.,.), Sigma₁[(.,.),(.,.)] }

g ~ GP{m2 (.), Sigma2[(.),(.)] }

除其他外，目标是获取未知校准参数 theta 的后验样本。

到目前为止我做了什么：

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
import pyDOE
from scipy.stats.distributions import uniform

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 = 60
  design = pyDOE.lhs(2, samples=N, criterion='center')

  left = [-0.2,-0.2]; right = [4.2,1.2]
  for i in range(2):
    design[:,i] = uniform(loc=left[i],scale=right[i]-left[i]).ppf(design[:,i])

  x_comp = design[:,0][:,None]; t_comp = design[:,1][:,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:

    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)
    f = pm.gp.Marginal(cov_func=cov)
    # train the gp f with data from computer model:
    f_0 = f.marginal_likelihood('f_0', X=input_comp, y=y_comp, noise=noise)
    trace = pm.sample(500, pm.Metropolis(), chains=4)
    burned_trace = trace[300:]

到这里为止，一切都很好。我的全科医生 f 是根据计算机模型进行训练的。 现在，我想测试我是否可以将这个训练有素的 GP 拟合到我观察到的数据中：

   #gp f is now trained to data from computer model
   #now I want to fit this trained gp to observed data and find posterior for theta

   with model:
    sd = pm.Gamma('eta', alpha=1, beta=1)
    theta = pm.Normal('theta', mu=0, sd=sd)
    sigma = pm.Gamma('sigma', alpha=1, beta=1)
    input_1 = tt.concatenate([x_obs_shared, tt.tile(theta, len(x_obs[:,None]), ndim=2).T], axis=1)
    f_1 = gp1.conditional('f_1', Xnew=input_1, shape=(10,))
    y_ = pm.Normal('y_', mu=f_1,sd=sigma, observed=y_obs)
    step = pm.Metropolis()
    trace_ = pm.sample(30000, step,start=pm.find_MAP(), chains=4)

这个公式正确吗？我得到非常不稳定的结果...... 根据 KOH 的完整配方应该是这样的：

    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)

有人能给我一些建议如何用 pymc3 正确地制定 KOH 吗？我很绝望......将不胜感激任何帮助。谢谢！

Answer 1

您可能已经找到了解决方案，但如果没有，that's 是一个很好的解决方案（建筑能源模型的贝叶斯校准指南）