将不确定性纳入 pymc3 模型

Question

我有一组数据，其中有每个点的平均值、标准差和观察次数（即，我了解有关测量准确性的知识）。在传统的 pymc3 模型中，我只关注方法，我可能会做以下事情：

x = data['mean']

with pm.Model() as m:
    a = pm.Normal('a', mu=0, sd=1)
    b = pm.Normal('b', mu=1, sd=1)
    y = a + b*x

    eps= pm.HalfNormal('eps', sd=1)
    likelihood = pm.Normal('likelihood', mu=y, sd=eps, observed=x)

将有关观测方差的信息纳入模型的最佳方法是什么？显然，结果应该对低方差观测值赋予比高方差（不太确定）观测值更大的权重。

统计学家建议的一种方法是执行以下操作：

x = data['mean']   # mean of observation
x_sd = data['sd']  # sd of observation
x_n = data['n']    # of measures for observation
x_sem = x_sd/np.sqrt(x_n)

with pm.Model() as m:
    a = pm.Normal('a', mu=0, sd=1)
    b = pm.Normal('b', mu=1, sd=1)
    y = a + b*x

    eps = pm.HalfNormal('eps', sd=1)
    obs = mc.Normal('obs', mu=x, sd=x_sem, shape=len(x))
    likelihood = pm.Normal('likelihood', mu=y, eps=eps, observed=obs)

但是，当我运行这个时，我得到：

TypeError: observed needs to be data but got: <class 'pymc3.model.FreeRV'>

我正在运行 pymc3 的 master 分支（3.0 有一些性能问题，导致采样时间非常慢）。

Answer 1

你已经很接近了，你只需要做一些小的改变。主要原因是 PyMC3 的数据始终是恒定的。检查以下代码：

with pm.Model() as m:
    a = pm.Normal('a', mu=0, sd=1)
    b = pm.Normal('b', mu=1, sd=1)
    mu = a + b*x

    mu_est = pm.Normal('mu_est', mu, x_sem, shape=len(x))

    likelihood = pm.Normal('likelihood', mu=mu_est, sd=x_sd, observed=x)

请注意，我保持数据固定，并在两点引入观察到的不确定性：估计mu_est和可能性。当然，您可以自由地不使用x_semorx_sd而是估计它们，就像您在代码中使用变量 所做的那样eps。

从历史记录来看，带有“随机数据”的代码曾经在 PyMC3 上运行（至少对于某些模型而言），但考虑到它并不是真正设计成以这种方式运行的，开发人员决定阻止用户使用随机数据，并且解释您收到的消息。