用PyMC解决逆问题

小编典典

用PyMC解决逆问题

python

假设我们给定了 X 的先验值（例如X〜高斯）和正向运算符 y = f（x） 。进一步假设我们通过实验观察了 y
，并且该实验可以无限期地重复。假定输出 Y 为高斯（Y〜高斯）或无噪声（Y〜Delta（观测值））。

给定观察结果，如何持续更新我们对 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 。

阅读 139

2020-12-20

共1个答案

小编典典

@
user1572508提供的功能现在是PyMC的一部分，名称为stochastic_from_data()或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')[:]

2020-12-20