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

根据 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:]

到这里为止，一切都很好。我的 GP 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 吗？我很绝望......将不胜感激任何帮助。谢谢!

最佳答案

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

https://stackoverflow.com/questions/61065989/