- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
我对概率编程和 pymc3 很陌生...目前,我想在 pymc3 中实现 Kennedy-O'Hagan 框架。
根据 Kennedy 和 O'Hagan 的论文设置如下:
我们有 n 个观测 zi 形式
zi = f(xi , theta) + g(xi) + ei,
其中 xi 是已知输入,theta 是未知校准参数,ei 是 iid 错误项。我们还有 m 模型评估 yj 的形式
yj = f(x'j, thetaj) ,其中 x'j(不同于上面的 xi)和 < em>thetaj 是已知的。因此,数据由所有 zi 和 yj 组成。在论文中,Kennedy-O'Hagan 模型 f,g 使用高斯过程:
f ~ GP{m1 (.,.), Sigma1[(.,.),(.,.)] }
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一个很好的(建筑能量模型贝叶斯校准指南)
关于python - 使用 PyMC3 进行贝叶斯校准,Kennedy O'Hagan,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/61065989/
使用 Kinect for Windows SDK,校准相机的最简单方法是什么? 我找到了这篇可爱的博文 http://nicolas.burrus.name/index.php/Research/K
我最初在 OpenCV 论坛上发布了这个,但不幸的是,我没有得到太多的意见/回复,所以我在这里发布是希望有人可能有方向建议? 我正在使用 Bumblebee XB3 立体相机,它有 3 个镜头。我花了
当移动很远的距离时,比如去另一个城市,关闭 GPS,打开它需要很长时间才能找到第一个正确的点。 因此,如果我通过示例注册到 gps 提供程序,则可能需要 10 分钟才能收到第一个 onLocation
在我的应用中,我使用加速度计来控制游戏中的角色。现在我只允许纵向,所以用户必须向右或向左倾斜设备才能移动角色。到目前为止一切正常。我现在想要完成的是“校准”加速度计以考虑用户正在玩的当前倾斜度。假设用
我使用 2 个 CANON EOS60d 作为立体设置来进行摄影测量。我正在使用 OpenCV 使用高精度 Circlegrid 校准模式校准两个相机。我正在使用信号发生器同时触发两个相机,早些时候我
我花了很长时间才让函数在 OpenCV 中运行,所以我想知道我的总体计划是否有意义,然后再深入研究尝试实现它的细节。 (2.3.1、Windows 7、C++)如果有任何建议,我将不胜感激。 问题:
我注意到 h2o.ai 套件的一个相对较新的添加,能够执行补充 Platt Scaling 以改进输出概率的校准。 (请参阅 calibrate_model in h2o manual 。)不过,在线
我正在尝试在 STM32F042 微 Controller 上读取 VDDA。我在 VDD 为 3.29V 时得到了意想不到的结果。我一定缺少一些基本的东西。 输出: VREFINT=1917; VR
请原谅我对编码完全陌生。首先,出于该项目的目的,我正在使用 Python 绑定(bind)到 OpenCV 的库。 我的相机已针对显示鱼眼失真进行了校准。我分别获得了 K 和 D 的以下值,即固有相机
如何使用 Netbeans 8.1 进行 JDK 校准。我用谷歌搜索了一些并找到了这个链接 here . It says "choose Profile > Advanced Commands > R
我目前正在开发增强现实应用程序。目标设备是光学透视 HMD,我需要校准其显示器以实现虚拟对象的正确注册。我用过那个implementation of SPAAM对于 android 来说,结果对于我的
我是一名优秀的程序员,十分优秀!