Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Time-to-Event Data, Censoring, and Hazard Functions › Implementation in NumPy and SciPy

Michael BrenndoerferOctober 25, 20267 min read

Machine Learning from ScratchAll chapters

Part 6 of 15Survival Analysis

Chapter 48Time-to-Event Data, Censoring, and Hazard Functions › Implementation in NumPy and SciPy

Reading section 7 of 10 · Implementation in NumPy and SciPy · 7 minSection 7 of 10 · 7 min

Section 7 of 10 · 7 min read

Time-to-Event Data, Censoring, and Hazard Functions › Implementation in NumPy and SciPySection view · Sections in this chapter7 / 10 ↓

Sections in this chapter10

Reading tools Definitions, links & notes

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Notes & comments

Notes & commentsLeave notes and join the discussionSelection toolbarShow comment actions when selecting text. Turn off for the browser menu only.Comment indicatorsShow small counts beside annotated passagesPassage highlightsGently shade annotated text selections

For larger datasets, code can accumulate the event count and exposure while retaining the same closed-form exponential estimator. We now build the log-likelihood from primitives, fit it with SciPy, and compare with the analytic estimator. This is a reusable pattern for continuous parametric survival models fitted by full right-censoring likelihood; other estimators use different objectives or risk-set calculations. We encode the hazard, encode the cumulative hazard, sum the log-likelihood contributions, and minimize the negative. The hand calculation gives us a known answer to check against. The implementation gives us a procedure that scales.

Our simulation returns $(Y,\Delta)$ with exponential latent lifetimes and independent uniform follow-up windows. As in the swimlane example, varying administrative follow-up can represent uniformly staggered installations with a common calendar cutoff. Latent event times are retained only as a teaching diagnostic; fitting uses the observed pair.

PythonIn [14]
CollapseExpanddef simulate_exponential_censored(n, lam_true, censoring_window, rng):
    """Simulate right-censored exponential event times and a binary event indicator."""
    T_true = rng.exponential(scale=1.0 / lam_true, size=n)
    C_time = rng.uniform(0.0, censoring_window, size=n)
    Y_obs = np.minimum(T_true, C_time)
    delta = (T_true <= C_time).astype(int)
    return Y_obs, delta, T_true


lam_true = 0.15
censoring_window = 12.0
rng = np.random.default_rng(7)

Y_train, d_train, T_train = simulate_exponential_censored(
    500, lam_true, censoring_window, rng
)
Y_test, d_test, _ = simulate_exponential_censored(
    200, lam_true, censoring_window, rng
)

We now derive the log-likelihood from scratch. For the exponential model with rate $\lambda$ , $h(y) = \lambda$ and $H(y) = \lambda y$ , so the negative log-likelihood is

-\ell(\lambda) = \lambda \sum_i Y_i - \sum_i \Delta_i \log \lambda.

where:

$\lambda$ : the exponential rate (constant hazard), with $\lambda > 0$
$Y_i$ : the observed time for record $i$
$\Delta_i$ : the event indicator for record $i$
$\sum_i Y_i$ : total time under observation, including censored records
$\sum_i \Delta_i$ : total number of observed events

We minimize this with SciPy's minimize. A small guard rejects any non-positive $\lambda$ , because $\log \lambda$ is undefined there.

PythonIn [15]
CollapseExpanddef negative_log_likelihood_exponential(lam, Y, delta):
    lam = float(np.asarray(lam).reshape(-1)[0])
    if not np.isfinite(lam) or lam <= 0:
        return np.inf
    return -(delta.sum() * np.log(lam) - lam * Y.sum())


fit_result = minimize(
    fun=negative_log_likelihood_exponential,
    x0=np.array([0.10]),
    args=(Y_train, d_train),
    method="Nelder-Mead",
    options=dict(xatol=1e-10, fatol=1e-10),
)

if not fit_result.success:
    raise RuntimeError(f"Exponential optimizer failed: {fit_result.message}")

lam_mle_numeric = float(fit_result.x[0])
lam_mle_analytic = float(d_train.sum() / Y_train.sum())
np.testing.assert_allclose(
    lam_mle_numeric, lam_mle_analytic, rtol=1e-7, atol=1e-10
)

ConsoleOut [16]
CollapseExpandevents in training set: 250 / 500
total observed time:    1743.3
analytic MLE lambda:    0.143406
numeric  MLE lambda:    0.143406
true     lambda:        0.150000

The numerical estimate agrees with $D/\sum_iY_i$ to the asserted relative tolerance of $10^{-7}$ ; this is a numerical check, not a machine-precision guarantee. The optimizer also reports successful termination. Agreement with the analytic solution checks this objective and fitting call. The estimate is near the generating rate $0.15$ in this sample, but sampling variability prevents exact recovery of a population parameter.

To complete the workflow, we compute the negative log-likelihood on the held-out test set at the fitted $\lambda$ . This tells us how the model generalizes to new records, not just how well it fit the training sample.

PythonIn [17]
CollapseExpandnll_train = float(
    negative_log_likelihood_exponential(lam_mle_analytic, Y_train, d_train)
)
nll_test = float(
    negative_log_likelihood_exponential(lam_mle_analytic, Y_test, d_test)
)
per_record_train = nll_train / len(Y_train)
per_record_test = nll_test / len(Y_test)

Out[18]:

	metric	train	test
0	negative log-likelihood	735.5183	310.7176
1	per-record NLL	1.4710	1.5536

The per-record negative log-likelihood is similar across these training and test samples. This is one held-out check under the same simulated data-generating and censoring processes. It does not prove calibration or absence of overfitting. Likelihood comparisons between event-time models should use the same observed records and observation mechanism. The chapter on evaluating survival predictions adds censoring-aware discrimination and probability-error measures.

We also fit the two naive estimators to the same observed training records. This measures their error in this simulated realization; repeated samples would be needed to estimate their sampling bias empirically.

PythonIn [19]
CollapseExpandlam_naive_events = len(Y_train) / Y_train.sum()
lam_naive_drop = d_train.sum() / Y_train[d_train == 1].sum()

Out[20]:

	estimator	lambda_hat	expected_lifetime
0	correct MLE	0.1434	6.9732
1	treat censored as events	0.2868	3.4866
2	drop censored	0.3237	3.0889

Both naive rates exceed the generating rate in this run: about $0.29$ when censorings are counted as events and $0.32$ after dropping them, versus $0.15$ . Their ordering follows from the retained exposure in this sample; closeness to the generating parameter is not guaranteed by that ordering in every sample. The figure compares the resulting survival curves.

PythonIn [21]
CollapseExpandt_grid_fit = np.linspace(0.01, 25.0, 400)
survival_true = np.exp(-lam_true * t_grid_fit)
survival_mle = np.exp(-lam_mle_analytic * t_grid_fit)
survival_naive_events = np.exp(-lam_naive_events * t_grid_fit)
survival_naive_drop = np.exp(-lam_naive_drop * t_grid_fit)

PythonIn [22]
CollapseExpandimport matplotlib.pyplot as plt
from book_plot_style import PALETTE, polish_axes, use_book_style

use_book_style()
fig, ax = plt.subplots(figsize=(6.0, 4.0))
ax.plot(
    t_grid_fit,
    survival_true,
    color=pal(0),
    linewidth=1.1,
    label=f"true  $\\lambda = {lam_true:.3f}$",
)
ax.plot(
    t_grid_fit,
    survival_mle,
    color=pal(1),
    linewidth=1.1,
    linestyle="--",
    label=f"correct MLE  $\\hat\\lambda = {lam_mle_analytic:.3f}$",
)
ax.plot(
    t_grid_fit,
    survival_naive_events,
    color=pal(2),
    linewidth=1.1,
    linestyle=":",
    label=f"treat censored as events  $\\hat\\lambda = {lam_naive_events:.3f}$",
)
ax.plot(
    t_grid_fit,
    survival_naive_drop,
    color=pal(3),
    linewidth=1.1,
    linestyle="-.",
    label=f"drop censored  $\\hat\\lambda = {lam_naive_drop:.3f}$",
)
ax.set_xlabel("Time t")
ax.set_ylabel("Survival probability S(t)")
ax.set_title("Survival fits with and without censoring")
ax.set_ylim(0, 1.02)
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2))
polish_axes(ax)
plt.tight_layout()
plt.show()

VisualizationOut [22]
CollapseExpand

Four monotonically decreasing curves. The true and MLE curves are close together and lie above the two naive curves. — Fitted exponential survival curves compared with the generating truth. The censoring-aware MLE closely tracks the true curve, while counting censorings as events or dropping them both depress the fitted survival curve at every positive plotted time in this realization.

The two naive curves lie below the generating survival curve on the positive plotted times, while the censoring-aware estimate is close in this realization. The simulation makes the consequence of the event indicator visible. It does not establish exact recovery or validate the exponential assumption for a real equipment dataset.

Key parametersLink Copied

The exponential has a single parameter, the rate $\lambda > 0$ , and its interpretation is worth spelling out because the same parameterization recurs throughout this part of the book. Its survival function is $S(t) = e^{-\lambda t}$ , its density is $f(t) = \lambda e^{-\lambda t}$ , and its hazard is constant at $h(t) = \lambda$ .

$\mathbb{E}[T] = 1/\lambda$ is the mean event time.
$\text{Median}(T) = \log 2 / \lambda$ is the median event time.
$\mathrm{Var}[T] = 1/\lambda^2$ , and the coefficient of variation (the standard deviation divided by the mean, $(1/\lambda)/(1/\lambda)$ ) is 1.
The cumulative hazard is $H(t) = \lambda t$ , linear in time.

These identities make the exponential a convenient baseline but a restrictive one. A constant hazard is restrictive when risk changes with age. The Weibull family we used in the visualization section has two parameters (shape and scale) and lets the hazard rise or fall. A later chapter will fit those parameters with the same style of likelihood.

Key methodsLink Copied

For independent right censoring without delayed entry, the same likelihood construction supports other absolutely continuous parametric event-time models:

Write the hazard $h_\theta(t)$ and the cumulative hazard $H_\theta(t)$ in terms of the model's parameters $\theta$ .
Form the negative log-likelihood $-\ell(\theta) = \sum_i H_\theta(Y_i) - \sum_i \Delta_i \log h_\theta(Y_i)$ .
Optimize with SciPy, guarding against invalid parameter regions (negative rates, positive-shape violations, etc.).
Validate on a held-out sample with the same per-record likelihood or a suitable survival metric.
Compare against the naive estimators to see how much censoring matters in this particular dataset.

Other families require their own parameter constraints, stable density or survival calculations, diagnostics, and checks that the fitted optimum is meaningful. Delayed entry, other censoring schemes, clustered observations, and competing events change the likelihood or data representation; they are not handled by substituting a new hazard alone.

Comments

Add to the conversation

A question, a correction, or an example from your own work.

No comments yet. Be the first to share your thoughts!

Key parametersLink Copied

Key methodsLink Copied

Comments

Add to the conversation

Stay up to date