构建包含不同测量值的 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 支持多个 observable，也就是说，您可以使用 observed 参数集将多个 RandomVariable 对象添加到图中。

单人审判

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

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)

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

先验

一旦您使该设置正常工作，重新考虑先验知识将符合您的利益。我怀疑您对这些典型值的信息比目前包含的更多，尤其是因为这看起来像是化学或物理模型。

例如，现在模型说，

We believe the true value of B remains fixed for the duration of a trial, but across trials is a completely arbitrary value between 0 and 1000, and measuring it repeatedly within a trial would be Poisson distributed.

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

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

building PyMC3 model incorporating different measurements

python

bayesian

pymc3

多个 Observables

单人审判

多次试验

先验