构建包含不同测量值的 PyMC3 模型

Question

我正在尝试将不同类型和重复的测量合并到 PyMC3 中的一个模型中。

考虑以下模型：P(t)=P0*exp(-kBt) 其中 P(t)、P0 和 B 是浓度。 k 是一个比率。我们通过粒子计数在不同时间测量 P(t) 和 B 一次。 k 是我们试图推断的感兴趣的参数。

我的问题有两个部分： (1) 如何将 P(t) 和 B 的测量值合并到一个模型中？ (2) 如何使用可变数量的重复实验来告知k的值？

我想我可以回答第 (1) 部分，但不确定它是否正确或以正确的方式完成。我未能概括代码以包含可变数量的复制。

对于一个实验（一个重复）：

ts=np.asarray([time0,time1,...])
counts=np.asarray([countforB,countforPattime0,countforPattime1,...])
basic_model = pm.Model()
with basic_model:
    k=pm.Uniform('k',0,20)
    B=pm.Uniform('B',0,1000)
    P=pm.Uniform('P',0,1000)
    exprate=pm.Deterministic('exprate',k*B)
    modelmu=pm.math.concatenate(B*(np.asarray([1.0]),P*pm.math.exp(-exprate*ts)))
    Y_obs=pm.Poisson('Y_obs',mu=modelmu,observed=counts))

我尝试按照上述方法包含不同的复制，但无济于事：

...
    k=pm.Uniform('k',0,20) # same within replicates
    B=pm.Uniform('B',0,1000,shape=numrepl) # can vary between expts.
    P=pm.Uniform('P',0,1000,shape=numrepl) # can vary between expts.
    exprate=???
    modelmu=???

Answer 1

多个 Observables

PyMC3 支持多个 observables，也就是说，您可以将多个 RandomVariable 对象添加到带有 observed 参数集的图中。

单试

在您的第一种情况下，这将使模型更加清晰：

counts=[countforPattime0, countforPattime1, ...]

with pm.Model() as single_trial:
    # priors
    k = pm.Uniform('k', 0, 20)
    B = pm.Uniform('B', 0, 1000)
    P = pm.Uniform('P', 0, 1000)

    # transformed RVs
    rate = pm.Deterministic('exprate', k*B)
    mu = P*pm.math.exp(-rate*ts)

    # observations
    B_obs = pm.Poisson('B_obs', mu=B, observed=countforB)
    Y_obs = pm.Poisson('Y_obs', mu=mu, observed=counts)

多次试验

有了这种额外的灵活性，希望它可以更明显地过渡到多个试验。它应该是这样的：

B_cts = np.array(...) # shape (N, 1)
Y_cts = np.array(...) # shape (N, M)
ts = np.array(...)    # shape (1, M)

with pm.Model() as multi_trial:
    # priors
    k = pm.Uniform('k', 0, 20)
    B = pm.Uniform('B', 0, 1000, shape=B_cts.shape)
    P = pm.Uniform('P', 0, 1000, shape=B_cts.shape)

    # transformed RVs
    rate = pm.Deterministic('exprate', k*B)
    mu = P*pm.math.exp(-rate*ts)

    # observations
    B_obs = pm.Poisson('B_obs', mu=B, observed=B_cts)
    Y_obs = pm.Poisson('Y_obs', mu=mu, observed=Y_cts)

可能需要一些额外的语法来使矩阵正确相乘，但这至少包括正确的形状。

---

先验

一旦您的设置正常工作，重新考虑先验条件将符合您的利益。我怀疑您对这些典型值的了解比目前包含的更多，特别是因为这看起来像是一个化学或物理模型。

例如，现在模型说，

我们相信B 的真实值在试验期间保持不变，但在试验期间是一个介于 0 和 1000 之间的完全任意值，在试验中重复测量它会是泊松分布。

通常，应该避免截断，除非它们排除了无意义的值。因此，0 的下限很好，但上限是任意的。我建议看看the Stan Wiki on choosing priors。