import arviz as az
import numpy as np
import pymc as pm

rng = np.random.default_rng(6420)

5. Time-to-event Models: Example 2*

Consider the model

\[\begin{split} \begin{align*} Y_i & \sim Wei(\nu,\lambda) \\ \nu & = \frac{3}{2} \\ \lambda & \sim Ga(2,3) \\ \end{align*}\end{split}\]

with data

\[ Y_1 = 2, Y_2 = 3, Y_3 = 1^*, Y_4 = \frac{5}{2}, Y_5 = 3^*\]

where *’s represent censored points. The goal is to estimate the parameter \(\lambda\) in Bayesian fashion. The likelihood is the Weibull PDF and Survival function for non-censored and censored points, respectively

\[\begin{split} \begin{align*} L(\lambda | \nu, y_1, ...,y_n) & = \prod_{i=1}^k \nu \lambda y_i ^ {\nu-1} e^{-\lambda y_{i} ^\nu}\times \prod_{i=k+1}^n e^{-\lambda y_i ^nu} \\ & = \nu^k \lambda^k \left ( \prod_{i=1}^k y_i\right )^{\nu-1} e^{-\lambda \sum_{i=1}^n y_i^\nu} \end{align*} \end{split}\]

We can elimate constants, and see that the likelihood is proportional to

\[ L(\lambda | \nu, y_1, ...,y_n) \propto \lambda^k e^{-\lambda \sum_{i=1}^n y_i^\nu} \]

Prior of \(\lambda\) is specified as \(Ga(\alpha,\beta)\), which is also simplified

\[ \pi(\lambda | \alpha, \beta) \propto \lambda^{\alpha-1} e^{-\beta \lambda}\]

With the joint likelihood and prior, we can find the posterior distribution

\[\pi(\lambda | \nu, y_1,...,y_n,\alpha,\beta) \propto \lambda^{k+\alpha-1} e^{-\lambda (\beta + \sum_{i=1}^n y_i^\nu)}\]

And by substituting in the known values, we arrive at

\[\pi(\lambda | \nu, y_1,...,y_n,\alpha,\beta) \propto Ga(3 + 2, 18.1736 +3)\]

The estimator \(\hat{\lambda}_{bayes}\) is the mean of this Gamma distribution

\[ \hat{\lambda}_{bayes} = \frac{5}{21.1736} = 0.2361\]

PyMC implementation

Now, let’s use PyMC to find the \( \hat{\lambda}_{bayes} \). Note, we are converting between Weibull parameterizations to be compatible with the courses historical solution in WINBUGS.

censored = np.array([np.inf, np.inf, 1, np.inf, 3])
y = np.array([2, 3, 1, 2.5, 3])

with pm.Model() as m:
    λ = pm.Gamma("λ", 2, 3)
    α = 1.5
    β = λ ** (-1 / α)
    obs_latent = pm.Weibull.dist(α, β)
    likelihood = pm.Censored(
        "likelihood", obs_latent, lower=None, upper=censored, observed=y
    )

    trace = pm.sample(1000, random_seed=rng)

az.summary(trace)

NUTS[nutpie]: [λ]

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
λ	0.237	0.106	0.095	0.43	1706	1704	1.00	0.0024	0.0022

%load_ext watermark
%watermark -n -u -v -iv -p pytensor

Last updated: Mon, 19 Oct 2026

Python implementation: CPython
Python version       : 3.14.7
IPython version      : 9.16.1

pytensor: 3.3.0

arviz: 1.3.0
numpy: 2.4.6
pymc : 6.3.1