Survival Analysis - Computing the mean residual lifetime

Questions v5
,
1

Hello, I am trying to compute the mean residual lifetime from the survival function. The trouble I am having is in regards with taking an integral. \frac{\int_{x}^{\infty}e^{(-(\frac{t}{\sigma})^a)}dt}{e^{(-(\frac{x}{\sigma})^a)}}. Here is the example I am using to try to figure this all out:

# Imports
import numpy as np
import pandas as pd
import pymc as pm

import matplotlib.pyplot as plt
import arviz as az

import seaborn as sns

from matplotlib import pyplot as plt
from matplotlib.ticker import StrMethodFormatter

# Load the data
df = pd.read_csv('https://gist.github.com/epogrebnyak/7933e16c0ad215742c4c104be4fbdeb1/raw/c932bc5b6aa6317770c4cbf43eb591511fec08f9/lamps.csv')

# Process data
event_times = np.concatenate([np.repeat(i, j) for i, j in zip(df.h, df.f)]).ravel()

# Construct the model
with pm.Model() as weibull_model:
    k = pm.Uniform('k', lower=0, upper=20)
    lambda_ = pm.Uniform('lambda_', lower=1000, upper=2000)
    lifetimes = pm.Weibull('lifetimes', alpha=k, beta=lambda_, observed=event_times)

# Sample
with weibull_model:
    idata = pm.sample(draws=8000, tune=2000, target_accept=0.99)

# Draw PPC
with weibull_model:
    ppc = pm.sample_posterior_predictive(idata, var_names=['lifetimes'])

# General survival curve
t_plot = np.linspace(500, 4000, 100)

weibull_pp_surv = np.greater_equal.outer(
    ppc.posterior_predictive['lifetimes'].mean(axis=0).to_numpy(), t_plot
)
weibull_pp_surv_mean = weibull_pp_surv.mean(axis=0)

sns.set()
blue, green, red, purple, gold, teal = sns.color_palette(n_colors=6)

pct_formatter = StrMethodFormatter("{x:.1%}")

fig, ax = plt.subplots(figsize=(8, 6))

ax.plot(t_plot, weibull_pp_surv_mean[0], c=blue)
ax.plot(t_plot, np.repeat(0.5, 100), c=blue, linestyle='--', label='Median Life')

ax.set_xlim(500, 2000)
ax.set_xlabel("Hours")

ax.set_ylim(top=1.05)
ax.yaxis.set_major_formatter(pct_formatter)
ax.set_ylabel("Survival probability")

ax.legend(loc=1)
ax.set_title("Weibull $S(t)$");

I am not sure how to approach the integration problem. Any help would really be appreciated.