Survival analysis in PyMC 3
-
21-12-2019 - |
Question
I tried to port simple survival model from here (the first one in introduction) form PyMC 2 to PyMC 3. However, I didn't find any equivalent to "observed" decorator and my attempt to write a new distribution failed. Could someone provide an example how is this done in PyMC 3?
Solution
This is a tricky port, and requires three new concepts:
- Use of the
theano
tensor - Use of the
DensityDist
- Passing a
dict
asobserved
This code provides the equivalent model as the PyMC2 version you linked to above:
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);
Output as follows:
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow