3

我对概率编程和 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 1 (.,.), 西格玛1 [(.,.),(.,.)] }

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 吗？我很绝望......将不胜感激任何帮助。谢谢！

4

1 回答 1

0

您可能已经找到了解决方案，但如果没有，那就太好了（建筑能源模型的贝叶斯校准指南）

于 2021-12-06T09:20:25.770 回答