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

9. Type-1 Censoring*#

Censoring occurs when we have incomplete information about the value of an observation. That is, instead of observing the exact value of a variable, we only know that it falls above, below, or between certain thresholds. This is common in many fields:

  • In survival analysis, if a patient leaves a study before the event (e.g., death) occurs, we only know their survival time exceeds the time they were last seen.

  • In industrial testing, sometimes instruments can only measure up to a certain maximum (or minimum).

    • A scale might return a value of “too heavy”, over a certain weight threshold. This is known as right censoring.

    • A water testing kit that unable to detect a chemical substance if the concentration is below \(x\) parts per million (ppm) might return “not detected”. This is an example of left censoring.

    • Notice that, in both cases, though we cannot observe the value exactly, we are resonably certain it is above (right) or below (left) a certain threshold.

Bartholomew (1957)#

This is our first introduction to censoring. It is adapted from Unit 6: exponential1.odc.

Results of a life test on ten pieces of equipment#

The observations above are the result of a test on the lifespan of some piece of equipment. The test had to end on August 31, but items 2, 4, and 10 did not fail by that date (the values in parentheses are the eventual lifespans of those pieces of equipment, but they were not known on August 31). So all we know for the purposes of the experiment is that they worked for 72, 60, and 21 days, respectively, and that they will continue working for some unknown additional time period.

Dataset and table from: A Problem in Life Testing (Bartholemew 1957).

Item No.

1

2

3

4

5

6

7

8

9

10

Date of Installation

11 June

21 June

22 June

2 July

21 July

31 July

31 July

1 Aug.

2 Aug.

10 Aug.

Date of Failure

13 June

—

12 Aug.

—

23 Aug.

27 Aug.

14 Aug.

25 Aug.

6 Aug.

—

Life in Days (\(t_i\))

2

(119)

51

(77)

33

27

14

24

4

(37)

No. Days to End of Period (\(T_i\))

81

72

70

60

41

31

31

30

29

21

The goal of the experiment is to provide a lifespan estimate. We could:

  1. Take the censored data as observed:

    • Divide the total observed working time (308 days) and divide by the equipment count (10) to get an average lifetime of 308 for an estimate of 30.8 days.

    • Problem: underestimates actual average lifespan.

  2. Ignore the equipment that didn’t fail:

    • Take mean lifetime of the pieces of equipment that broke within the experiment period for an estimate of 22.1 days.

    • Problems: selection bias, underestimates even more.

  3. Use the classical method:

    • Maximum Likelihood Estimation (MLE) under an exponential model. Total observed lifetime divided by 7 (number of observed failures) for an estimate of 44.0 days.

    • Problem: small sample size.

The actual mean lifespan of all pieces of equipment was 38.8 days.

A Bayesian Approach : Two Ways#

We will right-censor the three machines that still hadn’t failed by August 31, following this example.

We specify our priors as above, and use pm.Deterministic for \(\frac{1}{\lambda}\). Notice we create introduce an exponential distribution without supplying an observed argument. Then, wrapping that in pm.Censored with an observed argument tells PyMC that this is a likelihood with censoring. Specifying combinations of lower and upper arguments will give the desired censoring:

  • lower=None and upper specified : right censoring

  • lower specified and upper=None : left censoring

  • both lower and upper specified: interval censoring

Imputing Censored Points#

We can impute the 3 censored points by using a left-truncated version of the latent distribution. This makes sense because we know the lifetime reached the censored value, and then we sample the latent distribution on the right of each censored point.

# gamma dist parameters
α = 0.01
β = 0.1

# max possible observed life for each piece of equipment (before end of experiment)
censored = (81, 72, 70, 60, 41, 31, 31, 30, 29, 21)

# observed life within experiment dates
y = (2, 72, 51, 60, 33, 27, 14, 24, 4, 21)


with pm.Model() as m:
    λ = pm.Gamma("λ", α, β)
    μ = pm.Deterministic("μ", 1 / λ)
    obs_latent = pm.Exponential.dist(λ)

    observed = pm.Censored(
        "observed",
        obs_latent,
        lower=None,
        upper=censored,
        observed=y,
        shape=len(y),
    )

    # impute censored points
    censored_imputed = pm.Truncated(
        "censored_imputed", obs_latent, lower=[72, 60, 21]
    )

    trace = pm.sample(5000)
Hide code cell output 
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [λ, censored_imputed]



Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 1 seconds.















az.summary(trace, hdi_prob=0.95)













  
    

      
      mean
      sd
      hdi_2.5%
      hdi_97.5%
      mcse_mean
      mcse_sd
      ess_bulk
      ess_tail
      r_hat
    

  
  
    

      censored_imputed[0]
      123.111
      61.181
      72.000
      234.111
      0.580
      1.443
      9004.0
      6434.0
      1.0
    

    

      censored_imputed[1]
      111.495
      61.776
      60.008
      223.221
      0.541
      1.360
      10197.0
      7827.0
      1.0
    

    

      censored_imputed[2]
      71.592
      58.517
      21.001
      181.962
      0.508
      1.152
      11514.0
      8874.0
      1.0
    

    

      λ
      0.023
      0.009
      0.008
      0.040
      0.000
      0.000
      11541.0
      9435.0
      1.0
    

    

      μ
      50.986
      22.636
      18.702
      93.074
      0.240
      0.452
      11541.0
      9435.0
      1.0
    

  











The Conjugate Prior Approach:#


The priors chosen above also allow us to use analytical methods (conjugate priors).  This will give us some intuition for what happens to the likelihood function for censored observations.  We will use an \(\text{Exponential}(\lambda)\) likelihood, where \(\lambda\) has a \(\text{Gamma}(\alpha = 0.01,\beta = 0.1)\) prior distribution.  First, consider the likelihood function for the case that all values were fully observed (no censoring):



\[
f(Y=y\vert \lambda ) = \prod_{i=1}^{n}\lambda \exp\big\{-\lambda y_{i}\big\}
\]


We evaluate the exponential PDF at the observation value and multiply them all together.  Consider now the cases where a failure was not observed.  Note that for these, \(T_{i}\) in the table above denotes the point after which a failure would be observed had it occurred.


In the censored case, points above the censoring threshold \(T_{i}\) will contribute to the likelihood a bit differently:



\[\begin{align*}
f(Y \geq T_i \mid \lambda)
&= \int_{T_i}^{\infty} \lambda \exp\big\{ -\lambda y \big\} \, dy = \exp\big\{ -\lambda T_i \big\}
\end{align*}\]


The censoring-adjusted likelihood function becomes:



\[
f(y\vert \lambda) = \prod_{\{i\vert y_{i} < T_{i}\}}\lambda \exp\big\{-\lambda y_{i}\big\}\times \prod_{\{i\vert y_{i}\geq T_{i}\}}\exp\big\{-\lambda T_i\big\}
\]


We introduce an indicator variable \(\delta_{i} = 1\) for \(y_{i}< T_{i}\), otherwise \(\delta_{i} = 0\).  The posterior distribution is:



\[\begin{split}
\begin{align*}
\pi(\lambda \vert \mathbf{y})&\propto \prod_{\{i\vert y_{i} < T_{i}\}}\lambda \exp\big\{-\lambda y_{i}\big\}\times \prod_{\{i\vert y_{i}\geq T_{i}\}}\exp\big\{-\lambda T_i\big\}\times \lambda^{\alpha - 1}\exp\big\{ - \lambda \beta \big\}\\
&\phantom{;}\\ 
&= \lambda^{(\alpha + \sum_{i}\delta_{i}) - 1}\exp\big\{-\lambda (\beta + \sum_{\{i\vert y_{i} < T_{i}\}}y_{i} +  \sum_{\{i\vert y_{i} \geq T_{i}\}}T_{i})\big\}\\
&\sim \text{Gamma}(\alpha + \sum_{i=1}^{n}\delta_{i}\text{ , } \beta + \sum_{\{i\vert y_{i} < T_{i}\}}y_{i} +  \sum_{\{i\vert y_{i} \geq T_{i}\}}T_{i})
\end{align*}
\end{split}\]


Take note of the subset indexing that appears in the rate parameter.  We only want to add the \(y_{i}\)’s that are fully observed.  For observations that are right censored, we add the \(T_{i}\)’s instead.


The posterior distribution for \(\lambda\) (the failure rate) summarizes our updated beliefs about the typical rate at which equipment failures occur, given both the observed and censored data. Lower values of \(\lambda\) indicate a lower failure rate, and longer equipment lifetimes.  Next, for interpretability, let’s define \(\mu = \frac{1}{\lambda}\).  Then \(\pi(\mu\vert \mathbf{y})\) is a posterior distribution of expected lifetimes until equipment failure.  Next we will:


    

  • Calculate the mean and standard deviation of \(\pi(\lambda \vert \mathbf{y})\) using their closed form solutions
    • 

  • Calculate an estimate of the 95% HDI credible interval using np.random.gamma samples and az.hdi
    • 

  • Calculate the mean, standard deviation, and 95% HDI credible interval for \(\pi(\mu\vert \mathbf{y})\) using the reciprocals of our samples above.
    • 



An estimate of the mean, standard deviation, and 95% HDI credible interval for \(\pi(\frac{1}{\lambda}\vert \mathbf{y})\).  This is a posterior distribution for the mean lifetime of the equipment.






# prior parameters:
α = 0.01
β = 0.1

# data:

# max possible observed life for each piece of equipment (before end of experiment)
T = (81, 72, 70, 60, 41, 31, 31, 30, 29, 21)

# observed life within experiment dates
Y = (2, 72, 51, 60, 33, 27, 14, 24, 4, 21)

ty_paired = list(zip(T, Y))

# prior alpha    +   number of fully observed
α_posterior = α + sum(y < t for t, y in zip(T, Y))


# prior beta     +    #sum of fully observed values      +    #sum of T's for right censored values ONLY
β_posterior = (
    β
    + sum((y < t) * (y) for t, y in ty_paired)
    + sum((y >= t) * (t) for t, y in ty_paired)
)


# Signature: np.random.gamma(shape, scale=1.0, size=None)
# Docstring:
# gamma(shape, scale=1.0, size=None)
# TLDR: need to use "scale= 1/β_posterior"

samples = np.random.gamma(α_posterior, 1 / β_posterior, 100000)


print("λ Posterior:")
hdi = az.hdi(samples, hdi_prob=0.95)
summary = f"Posterior Mean: {α_posterior / β_posterior:<.4f}   Posterior SD: {(α_posterior / ((β_posterior)**(2)))**(1/2):<.4f}   HDI Bounds[{hdi[0]:<.4f},{hdi[1]:<.4f}]"
print("-" * len(summary))
print(summary)
print("")
print("μ = 1/λ Posterior:")
hdi = az.hdi(1 / samples, hdi_prob=0.95)
summary = f"Posterior Mean: {(1/samples).mean():<.4f}   Posterior SD: {(1/samples).std():<.4f}   HDI Bounds[{hdi[0]:<.4f},{hdi[1]:<.4f}]"
print("-" * len(summary))
print(summary)









λ Posterior:
-------------------------------------------------------------------------
Posterior Mean: 0.0228   Posterior SD: 0.0086   HDI Bounds[0.0078,0.0400]

μ = 1/λ Posterior:
-----------------------------------------------------------------------------
Posterior Mean: 51.2348   Posterior SD: 22.9443   HDI Bounds[19.2470,95.1005]









The results closely match what we found using PyMC.






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









Last updated: Wed Jun 11 2025

Python implementation: CPython
Python version       : 3.13.3
IPython version      : 9.2.0

pytensor: 2.30.3

numpy: 2.2.6
pymc : 5.22.0
arviz: 0.21.0