I was trying to follow this post, but github gists miss few variables, correct code shall be:

import pymc as mc import numpy as np N=20 #create some artificial data, like the data in the figure above. 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_ = np.ma.masked_array( lifetime, censor ) #create the censorship event. lifetime_.set_fill_value( 10 ) #good for computations later. #this begins the model alpha = mc.Uniform("alpha", 0,20) #lets just use uninformative priors beta = mc.Uniform("beta", 0,20) obs = mc.Weibull( 'obs', alpha, beta, value = lifetime_, observed = True )

Posteriory

@mc.potential def censorfactor(obs=obs): if np.any((obs + birth < 10)[lifetime_.mask] ): return -100000 else: return 0 #perform Markov Chain Monte Carlo - see chapter 3 of BMH mcmc = mc.MCMC([alpha, beta, obs, censorfactor ] ) mcmc.sample(50000, 30000)

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)

## Comments

comments powered by Disqus