7

假设我们有一个关于X的先验（例如 X ~ Gaussian）和一个前向算子y = f(x)。进一步假设我们通过一个实验观察了y，并且这个实验可以无限重复。假设输出Y为高斯 (Y ~ Gaussian) 或无噪声 (Y ~ Delta(observation))。

给定观察结果，如何持续更新我们对X的主观知识程度？我用 PyMC 尝试了以下模型，但似乎我遗漏了一些东西：

from pymc import *

xtrue = 2                        # this value is unknown in the real application
x = rnormal(0, 0.01, size=10000) # initial guess

for i in range(5):
    X = Normal('X', x.mean(), 1./x.var())
    Y = X*X                        # f(x) = x*x
    OBS = Normal('OBS', Y, 0.1, value=xtrue*xtrue+rnormal(0,1), observed=True)
    model = Model([X,Y,OBS])
    mcmc = MCMC(model)
    mcmc.sample(10000)

    x = mcmc.trace('X')[:]       # posterior samples

后验没有收敛到xtrue。

4

2 回答 2

7

@user1572508 的功能现在是 PyMC 的一部分，名称为stochastic_from_data()or Histogram()。然后该线程的解决方案变为：

from pymc import *
import matplotlib.pyplot as plt

xtrue = 2 # unknown in the real application
prior = rnormal(0,1,10000) # initial guess is inaccurate
for i in range(5):
  x = stochastic_from_data('x', prior)
  y = x*x
  obs = Normal('obs', y, 0.1, xtrue*xtrue + rnormal(0,1), observed=True)

  model = Model([x,y,obs])
  mcmc = MCMC(model)
  mcmc.sample(10000)

  Matplot.plot(mcmc.trace('x'))
  plt.show()

  prior = mcmc.trace('x')[:]
于 2013-07-29T21:55:26.590 回答
5

问题是您的函数 $y = x^2$ 不是一对一的。具体来说，因为您在平方 X 的符号时丢失了所有有关 X 符号的信息，所以无法从您的 Y 值判断您最初使用 2 还是 -2 来生成数据。如果您在第一次迭代之后为 X 创建了轨迹的直方图，您将看到：第一次迭代后迹线的直方图

此分布有 2 种模式，一种为 2（您的真实值），另一种为 -2。在下一次迭代中，x.mean() 将接近于零（对双峰分布进行平均），这显然不是您想要的。

于 2013-07-02T00:33:50.617 回答