变量可以用作 PyMC3 模型中的 'observed' 数据吗?
Can a variable be used as 'observed' data in a PyMC3 model?
我是贝叶斯世界和 PyMC3 的新手,正在为简单的模型设置而苦苦挣扎。具体来说,如何处理 'observed' 数据本身被随机变量修改的设置?举个例子,假设我有一个二维点 [Xi, Yi] 的集合,这些点形成一个圆弧,圆的中心点是 [Xc,Yc],我不知道。但是,我希望点与圆心 Ri 之间的距离应该呈正态分布,大约已知半径 R。因此,我最初认为我可以分配 Xc 和 Yc 统一先验(在任意大的范围内)然后在模型中重新计算 Ri 并将 Ri 指定为 'observed' 数据以获得 Xc 和 Yc 的后验估计:
import pymc3 as pm
import numpy as np
points = np.array([[2.95, 4.98], [3.28, 4.88], [3.84, 4.59], [4.47, 4.09], [2.1,5.1], [5.4, 1.8]])
Xi = points[:,0]
Yi = points[:,1]
#known [Xc,Yc] = [2.1, 1.8]
R = 3.3
with pm.Model() as Cir_model:
Xc = pm.Uniform('Xc', lower=-20, upper=20)
Yc = pm.Uniform('Yc', lower=-20, upper=20)
Ri = pm.math.sqrt((Xi-Xc)**2 + (Yi-Yc)**2)
y = pm.Normal('y', mu=R, sd=1.0, observed=Ri)
samples = pm.fit(random_seed=2020).sample(1000)
pm.plot_posterior(samples, var_names=['Xc'])
pm.plot_posterior(samples, var_names=['Yc']);
虽然这段代码运行并给了我一些东西,但它显然不能正常工作,这并不奇怪,因为在 'observed' 数据中输入变量 (Ri) 似乎是不对的.然而,虽然我知道我的模型设置存在严重错误(以及我更普遍的理解),但我似乎无法识别它。非常感谢任何帮助!
此模型实际上运行良好,但您可能需要改进以下几点:
- 使用变量作为观察值并不好,因为您应该考虑它对您拟合的分布有何影响。它将符合 a 分布,但您应该考虑是否在先验和似然中重复计算变量。不过对于这个玩具模型来说,这并不重要!
- 您正在使用
pm.fit(...)
,它使用 variational inference,但 MCMC 在这里没问题,所以用 samples = pm.sample()
替换整行是可行的。
- 你提供的
points
几乎正好在一个圆上——经验标准差在0.004左右,但你提供的标准差很可能是1:真实值的250倍左右!按原样从模型中采样允许点的中心位于两个不同的位置:
如果将可能性更改为 y = pm.Normal('y', mu=R, sd=0.01, observed=Ri)
,您仍然会得到两个可能的中心,尽管在真正的中心附近有更多的质量:
最后,您可以采取一种方法,在量表上放置先验知识,并学习到最有原则的 和 给出的结果最接近真实结果.这是模型:
with pm.Model():
Xc = pm.Uniform('Xc', lower=-20, upper=20)
Yc = pm.Uniform('Yc', lower=-20, upper=20)
Ri = pm.math.sqrt((Xi-Xc)**2 + (Yi-Yc)**2)
obs_sd = pm.HalfNormal('obs_sd', 1)
y = pm.Normal('y', mu=R, sd=obs_sd, observed=Ri)
samples = pm.sample()
这是输出:
我是贝叶斯世界和 PyMC3 的新手,正在为简单的模型设置而苦苦挣扎。具体来说,如何处理 'observed' 数据本身被随机变量修改的设置?举个例子,假设我有一个二维点 [Xi, Yi] 的集合,这些点形成一个圆弧,圆的中心点是 [Xc,Yc],我不知道。但是,我希望点与圆心 Ri 之间的距离应该呈正态分布,大约已知半径 R。因此,我最初认为我可以分配 Xc 和 Yc 统一先验(在任意大的范围内)然后在模型中重新计算 Ri 并将 Ri 指定为 'observed' 数据以获得 Xc 和 Yc 的后验估计:
import pymc3 as pm
import numpy as np
points = np.array([[2.95, 4.98], [3.28, 4.88], [3.84, 4.59], [4.47, 4.09], [2.1,5.1], [5.4, 1.8]])
Xi = points[:,0]
Yi = points[:,1]
#known [Xc,Yc] = [2.1, 1.8]
R = 3.3
with pm.Model() as Cir_model:
Xc = pm.Uniform('Xc', lower=-20, upper=20)
Yc = pm.Uniform('Yc', lower=-20, upper=20)
Ri = pm.math.sqrt((Xi-Xc)**2 + (Yi-Yc)**2)
y = pm.Normal('y', mu=R, sd=1.0, observed=Ri)
samples = pm.fit(random_seed=2020).sample(1000)
pm.plot_posterior(samples, var_names=['Xc'])
pm.plot_posterior(samples, var_names=['Yc']);
虽然这段代码运行并给了我一些东西,但它显然不能正常工作,这并不奇怪,因为在 'observed' 数据中输入变量 (Ri) 似乎是不对的.然而,虽然我知道我的模型设置存在严重错误(以及我更普遍的理解),但我似乎无法识别它。非常感谢任何帮助!
此模型实际上运行良好,但您可能需要改进以下几点:
- 使用变量作为观察值并不好,因为您应该考虑它对您拟合的分布有何影响。它将符合 a 分布,但您应该考虑是否在先验和似然中重复计算变量。不过对于这个玩具模型来说,这并不重要!
- 您正在使用
pm.fit(...)
,它使用 variational inference,但 MCMC 在这里没问题,所以用samples = pm.sample()
替换整行是可行的。 - 你提供的
points
几乎正好在一个圆上——经验标准差在0.004左右,但你提供的标准差很可能是1:真实值的250倍左右!按原样从模型中采样允许点的中心位于两个不同的位置:
如果将可能性更改为 y = pm.Normal('y', mu=R, sd=0.01, observed=Ri)
,您仍然会得到两个可能的中心,尽管在真正的中心附近有更多的质量:
最后,您可以采取一种方法,在量表上放置先验知识,并学习到最有原则的 和 给出的结果最接近真实结果.这是模型:
with pm.Model():
Xc = pm.Uniform('Xc', lower=-20, upper=20)
Yc = pm.Uniform('Yc', lower=-20, upper=20)
Ri = pm.math.sqrt((Xi-Xc)**2 + (Yi-Yc)**2)
obs_sd = pm.HalfNormal('obs_sd', 1)
y = pm.Normal('y', mu=R, sd=obs_sd, observed=Ri)
samples = pm.sample()
这是输出: