Байесовская калибровка с PyMC3, Кеннеди О'Хаган

Question

Я довольно новичок в вероятностном программировании c и pymc3 ... В настоящее время я хочу внедрить фреймворк Kennedy-O'Hagan в pymc3.

Установка в соответствии с работой Кеннеди и О'Хагана выглядит следующим образом:

Мы имеем n наблюдений z _i вида

z _i = f (x _i , тета) + g (x _i ) + e _i ,

, где x _i - известные импульсы, theta - неизвестные параметры калибровки и e _i - термины ошибки iid. У нас также есть m оценки модели y _j вида

y _j = f (x ^' _j , тета _j ) , где оба x ^' _j (отличные от x _i выше) и theta _j известны. Следовательно, данные состоят из всех z _i и y _j . В работе Кеннеди-О'Хагана модель f, g использует гауссовские процессы:

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

г ~ GP {м2 (.), Sigma2 [(.), (.)]}

Помимо прочего, цель состоит в том, чтобы получить последующие образцы для неизвестных параметров тета калибровки.

То, что я сделал до сих пор:

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 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)

Может кто-нибудь дать мне совет, как правильно сформулировать KOH с помощью pymc3? Я в отчаянии ... Буду признателен за любую помощь. Спасибо!