2

是否可以使用 PyMC 执行“伪实验”？

通过伪实验，我的意思是通过从先验抽样生成随机“观察”，然后，给定每个伪实验，从后验抽取样本。之后，将每个参数的轨迹与从后验抽样中使用的样本（从先验获得）进行比较。

一个更具体的例子：假设我想知道进程 X 的速率。我计算在某个时间段内发生了多少次。但是，我知道过程 Y 有时也会发生并且会污染我的计数。过程 Y 的速率是已知的，但存在一定的不确定性。所以，我建立了一个模型，包括我的观察结果，并从后验样本：

import pymc

class mymodel:
    rate_x = pymc.Uniform('rate_x', lower=0, upper=100)
    rate_y = pymc.Normal('rate_y', mu=150, tau=1./(15**2))
    total_rate = pymc.LinearCombination('total_rate', [1,1], [rate_x, rate_y])
    data = pymc.Poisson('data', mu=total_rate, value=193, observed=True)

Mod = pymc.Model(mymodel)
MCMC = pymc.MCMC(Mod)
MCMC.sample(100000, burn=5000, thin=5)
print MCMC.stats()['rate_x']['quantiles']

但是，在我进行实验之前（或在我“揭开”我的分析并查看我的数据之前），我想知道我期望的敏感度——我的 rate_x 测量的不确定性是什么？

为了回答这个问题，我可以从之前的样本中

Mod.draw_from_prior()

但这仅对rate_x,rate_y和计算进行采样total_rate。但是一旦这些值由 设置draw_from_prior()，我可以绘制一个伪实验：

Mod.data.random()

这只是返回一个数字，所以我必须将值设置Mod.data为随机样本。因为设置Mod.data了observed标志，我也必须“强制”它：

Mod.data.set_value(Mod.data.random(), force=True)

现在我可以再次从后部采样

MCMC.sample(100000, burn=500, thin=5)
print MCMC.stats()['rate_x']['quantiles']

所有这些都有效，所以我想我的问题的简单答案是“是”。但感觉非常hacky。有没有更好或更自然的方法来实现这一点？

4

0 回答 0