Software 2.5.1 (The Weibull functions). \(S\), \(h\) and \(\exp (-H)\) at \(t = 5\) for Weibull distributions with scale \(10\) and shapes \(0.5\), \(1\), \(2\); a check that \(f = hS\) integrates to \(1\); and the mean of the \(k = 2\) distribution as the area under \(S\) and from the formula \(\lambda \Gamma (1 + 1/k)\):

> S <- function(t, k, lam) exp(-(t / lam)^k)          # Weibull survival
> h <- function(t, k, lam) (k / lam) * (t / lam)^(k - 1)  # and its hazard
> H <- function(t, k, lam) (t / lam)^k                 # cumulative hazard
> for (k in c(0.5, 1, 2)) {
+   cat(sprintf("k = %3.1f  S(5) = %.4f  h(5) = %.4f  exp(-H(5)) = %.4f\n",
+               k, S(5, k, 10), h(5, k, 10), exp(-H(5, k, 10))))
+ }
k = 0.5  S(5) = 0.4931  h(5) = 0.0707  exp(-H(5)) = 0.4931
k = 1.0  S(5) = 0.6065  h(5) = 0.1000  exp(-H(5)) = 0.6065
k = 2.0  S(5) = 0.7788  h(5) = 0.1000  exp(-H(5)) = 0.7788
> f <- function(t) h(t, 2, 10) * S(t, 2, 10)           # density = h S
> integrate(f, 0, Inf)$value                           # integrates to 1
[1] 1
> c(mean_by_integral = integrate(function(t) S(t, 2, 10), 0, Inf)$value,
+   mean_by_formula = 10 * gamma(1 + 1/2))
mean_by_integral  mean_by_formula
        8.862269         8.862269

Reading the output. In every case \(S(5) = \exp (-H(5))\), as Theorem 2.3.2 requires. At \(t = 5\) the decreasing hazard (\(k = 0.5\)) has already removed half the population, while the increasing one (\(k = 2\)) has removed only \(22\%\). The area under \(S\) gives the mean \(8.862\), equal to \(10\,\Gamma (1.5)\).

Python: scipy.stats.weibull_min(c=k, scale=10) provides sf, pdf and mean.

Software 2.5.2 (Medians, memory and residual life). For the exponential with mean \(10\): the median, the conditional probability of surviving \(5\) more units having survived \(10\), and the unconditional probability of surviving \(5\). Then the same comparison, the median and the mean residual life at ages \(0, 5, 10, 20\) for the Weibull with \(k = 2\):

> lam <- 0.1  # exponential, mean 10
> c(median = log(2) / lam, mean = 1 / lam)
   median      mean
 6.931472 10.000000
> c(P_survive_5_more_given_10 = exp(-lam * 15) / exp(-lam * 10),
+   P_survive_5_from_start = exp(-lam * 5))            # memoryless
P_survive_5_more_given_10    P_survive_5_from_start
                0.6065307                 0.6065307
> Sw <- function(t) S(t, 2, 10)  # Weibull k = 2: ageing
> c(Sw(15) / Sw(10), Sw(5))  # conditional is lower
[1] 0.2865048 0.7788008
> med <- uniroot(function(t) Sw(t) - 0.5, c(0, 100), tol = 1e-10)$root
> c(median = med, formula = 10 * log(2)^(1/2))
  median  formula
8.325546 8.325546
> mrl <- function(t) integrate(Sw, t, Inf)$value / Sw(t)  # mean residual life
> round(sapply(c(0, 5, 10, 20), mrl), 3)
[1] 8.862 5.456 3.789 2.263

Reading the output. The exponential’s median, \(6.93\), is below its mean \(10\), and its conditional and unconditional probabilities are both \(0.607\): it is memoryless. For the Weibull with increasing hazard, surviving \(5\) more units after \(10\) has probability \(0.287\), much less than \(0.779\) from new: it ages. Its median is \(8.33 = 10\sqrt {\log 2}\), and its mean residual life falls from \(8.86\) when new to \(2.26\) at age \(20\).

Python: scipy.integrate.quad for the residual life integral.

Software 2.5.3 (A Gompertz life table). Example 2.4.2: the probabilities of surviving from \(30\) to \(60\), \(80\) and \(90\), the median age at death, the probability of reaching \(80\) from \(60\), and the doubling time of the hazard:

> a <- 5e-5; b <- 0.09  # Gompertz hazard a e^(b x)
> Hg <- function(x, x0) (a / b) * (exp(b * x) - exp(b * x0))
> Sg <- function(x, x0 = 30) exp(-Hg(x, x0))           # from age 30
> round(c(to_60 = Sg(60), to_80 = Sg(80), to_90 = Sg(90)), 4)
 to_60  to_80  to_90
0.8916 0.4791 0.1617
> uniroot(function(x) Sg(x) - 0.5, c(30, 120), tol = 1e-8)$root  # median
[1] 79.34316
> round(Sg(80) / Sg(60), 4)  # P(reach 80 | alive at 60)
[1] 0.5373
> round(log(2) / b, 2)  # hazard doubling time
[1] 7.7

Reading the output. \(89\%\) of \(30\)-year-olds reach \(60\) and \(48\%\) reach \(80\); the median age at death is \(79.3\); of those alive at \(60\), \(54\%\) reach \(80\). The hazard doubles every \(7.7\) years, close to the eight years observed in many populations.

Python: the same functions with numpy.exp and scipy.optimize.brentq for the median.