PyMC 3中的生存分析

Question

我试图将一个简单的生存模型从这里（介绍中的第一个）移植到PyMC 3。然而，我没有找到“observed”装饰器的等效物，我的尝试编写一个新的分布也失败了。有人能提供在PyMC 3中如何完成这个任务的示例吗？

- Marigold

Answer 1

这是一个棘手的端口，需要三个新概念：

1. 使用theano张量
2. 使用DensityDist
3. 将dict作为observed传递

此代码提供了与您上面链接的PyMC2版本相同的模型:

import pymc3 as pm
from pymc.examples import melanoma_data as data
import theano.tensor as t

times = data.t # not to be confused with the theano tensor t!
failure = (data.censored==0).astype(int)

with pm.Model() as model:

    beta0 = pm.Normal('beta0', mu=0.0, tau=0.0001)
    beta1 = pm.Normal('beta1', mu=0.0, tau=0.0001)
    lam = t.exp(beta0 + beta1*data.treat)

    def survival_like(failure, value):
        return t.sum(failure * t.log(lam) - lam * value)

    survive = pm.DensityDist('survive', survival_like,
                        observed={'failure': failure, 'value': times})

with model:

    start = pm.find_MAP()
    step = pm.NUTS(scaling=start)
    trace = pm.sample(10000, step=step, start=start)

pm.traceplot(trace);

以下是输出结果：

- inversion

---

当传递 observed={...} 时，参数如何传递到 survival_like？这些参数必须按字母顺序排列吗？谢谢！ - Stefan Novak

我不认为这有什么关系。但是您可以进行简单的测试来确认一下。 - inversion

谢谢！我已经深入研究了 PyMC 代码库，发现 logp 将使用 **data 调用，因此无论顺序如何，字典中的值都将传递到正确的参数。https://github.com/pymc-devs/pymc3/blob/master/pymc3/model.py#L535 - Stefan Novak

我认为这行代码： return t.sum(failure * t.log(lam) - lam * value) 应该改为： return t.sum(failure * (t.log(lam) - lam * value)) - Yetti