Type-1 Censoring*

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

rng = np.random.default_rng(6420)

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.

Typically in this course, we’ll see this in survival analysis or time-to-event modeling. If a patient leaves a study before the event (e.g. disease recurrence, death) occurs, we only know their survival time exceeds the time they were last seen. We don’t know when the event occurred. This is where right-censoring comes in.

Another example is water-quality testing. Say we’re testing for arsenic in groundwater. The lab’s test procedure will have a method detection limit (MDL); concentrations below the limit may exist, but aren’t detectable by the test. So we can’t say there’s absolutely no arsenic even if it’s not detected by the test. For a non-detection at an MDL of 0.05 μg/L, we can only say the arsenic concentration isn’t greater than that, not whether there’s arsenic at a concentration between 0 and 0.05 μg/L. This is when we might use left-censoring.

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: censoring at both ends of the distribution

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, random_seed=rng)

Hide code cell output

 
NUTS[nutpie]: [λ, censored_imputed]













az.summary(trace, ci_prob=0.95)










  
    

      
      mean
      sd
      eti95_lb
      eti95_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      λ
      0.0226
      0.0085
      0.0091
      0.042
      11199
      9692
      1.00
      7.6e-05
      6.5e-05
    

    

      censored_imputed[0]
      123
      59
      73
      280
      9635
      8118
      1.00
      0.54
      1.1
    

    

      censored_imputed[1]
      112
      61
      61
      270
      9910
      8021
      1.00
      0.56
      1.4
    

    

      censored_imputed[2]
      73
      60
      22
      240
      9634
      8065
      1.00
      0.56
      1.3
    

    

      μ
      51.5
      23
      24
      110
      11199
      9692
      1.00
      0.25
      0.44
    

  









The Conjugate Prior Approach#


We can also solve this problem analytically using conjugate priors. We’ll use the same exponential likelihood and Gamma prior as in the PyMC model:



\[
\lambda \sim \text{Gamma}(\alpha=0.01,\beta=0.1)
\]


Here, \(\beta\) is the rate parameter of the Gamma distribution.


First, suppose we observe the exact failure time \(y_i\) for every machine. The likelihood is simply the product of the exponential PDFs:



\[
L(\lambda \mid \mathbf{y}) = \prod_{i=1}^{n}\lambda e^{-\lambda y_i}
\]


Censored observations are different. If a machine was still working after \(T_i\) days, we don’t know its failure time. We only know that it exceeds \(T_i\). Instead of the PDF, we use the survival probability:



\[
P(Y_i > T_i \mid \lambda)
= \int_{T_i}^{\infty}\lambda e^{-\lambda y}\,dy
= e^{-\lambda T_i}
\]


So our likelihood has two parts: the exponential PDF for observed failures, and the survival probability for censored observations.



\[
L(\lambda \mid \mathbf{y})
=
\prod_{\text{observed}}\lambda e^{-\lambda y_i}
\prod_{\text{censored}}e^{-\lambda T_i}
\]


Let \(\delta_i=1\) if we observed a failure and \(\delta_i=0\) if the observation was censored. Combining the likelihood with our Gamma prior gives:



\[\begin{split}
\begin{align*}
\pi(\lambda \mid \mathbf{y})
&\propto L(\lambda \mid \mathbf{y})\,\pi(\lambda)\\
&\propto
\lambda^{\alpha+\sum_i\delta_i-1}
\exp\left[
-\lambda\left(
\beta+
\sum_{\text{observed}}y_i+
\sum_{\text{censored}}T_i
\right)
\right]
\end{align*}
\end{split}\]


Recognizing the Gamma distribution, we get:



\[
\lambda \mid \mathbf{y}
\sim
\text{Gamma}\left(
\alpha+\sum_i\delta_i,\;
\beta+\sum_{\text{observed}}y_i+
\sum_{\text{censored}}T_i
\right)
\]


The posterior shape parameter increases by the number of observed failures, while the rate parameter increases by the total observed operating time, including the time contributed by censored machines.


In our data, there are 7 observed failures and 308 total days of operating time:



\[
\lambda \mid \mathbf{y}\sim\text{Gamma}(7.01,308.1)
\]


The code below calculates the posterior mean and standard deviation of \(\lambda\) directly. We’ll also draw samples from this Gamma distribution to estimate a 95% HDI.


Finally, recall that the expected lifetime of a machine under the exponential model is \(\mu=1/\lambda\). We can transform our Gamma samples by taking their reciprocals to obtain posterior estimates and a credible interval for the mean lifetime.


The results should be close to the ones from the PyMC model.






# 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)
)


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


print("λ Posterior:")
hdi = az.hdi(samples, 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, 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.0080,0.0400]

μ = 1/λ Posterior:
-----------------------------------------------------------------------------
Posterior Mean: 51.1710   Posterior SD: 22.7014   HDI Bounds[19.6744,95.5234]













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









Last updated: Sun, 18 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