Translating Censored Data example from PyMC2 to PyMC3

Questions
1

Hi there,

I have found a great example of Bayesian parametric survival analysis in pymc2 here:
http://blog.yhat.com/posts/estimating-user-lifetimes-with-pymc.html

I want to reproduce it in pymc3. I am struggling with encoding the censored events. How would you translate this code to PyMC3?

N = 2500 #125 times more data than before
lifetime = mc.rweibull( 2, 5, size = N )
birth = mc.runiform(0, 10, N)
censor = ((birth + lifetime) >= 10) #an individual is right-censored if this is True 
lifetime_ = lifetime.copy()
lifetime_[censor] = 10 - birth[censor] #we only see this part of their lives.


alpha = mc.Uniform('alpha', 0, 20) #lets just use uninformative priors again
beta = mc.Uniform('beta', 0, 20)

@mc.observed
def survival(value=lifetime_, alpha = alpha, beta = beta ):
    return sum( (1-censor)(log( alpha/beta) + (alpha-1)*log(value/beta)) - (value/beta)**(alpha)

mcmc = mc.MCMC([alpha, beta, survival ] )
mcmc.sample(50000, 30000)

My understanding is that the @mc.observed decorator is not supported in PyMC3.

How would you translate this code?
What is the most elegant and scalable way to handle censored data in PyMC3?

Thanks,
Tokoro

2

Hi Tokoro,
There is a quite detail example of modelling censored data in the PyMC3/examples, it should answer most of your question.

3

Junpenglao,

Thanks for the quick reply and for the pointer. I do think the answer is in that example. I will tinker a bit with it and come back here if I have any further questions.

Thank you!
Tokoro

2 Likes