采用不同测量方法建立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

多个观测值

PyMC3支持多个observates，也就是说，可以将多个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)

可能会有一些额外的语法来让矩阵正确相乘，但这至少包括正确的形状。你知道吗

---

前科

一旦你得到设置工作，这将是在你的利益重新考虑前科。我怀疑你有更多的关于这些的典型值的信息，而不是现在包括的，特别是因为这看起来像是一个化学或物理模型。你知道吗

例如，现在模型说

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。你知道吗