[Reliability Bayesian pymc-devs#474] fixed empirical coverage estimate

Signed-off-by: Nathaniel <[email protected]>

Author: NathanielF

examples/case_studies/reliability_and_calibrated_prediction.ipynb

@@ -33,7 +33,6 @@
     "import pandas as pd\n",
     "import pymc as pm\n",
     "\n",
-    "from joblib import Parallel, delayed\n",
     "from lifelines import KaplanMeierFitter, LogNormalFitter, WeibullFitter\n",
     "from lifelines.utils import survival_table_from_events\n",
     "from scipy.stats import binom, lognorm, norm, weibull_min"
@@ -114,15 +113,13 @@
     }
    ],
    "source": [
-    "from scipy.stats import lognorm\n",
-    "\n",
     "mu, sigma = 6, 0.3\n",
     "\n",
     "\n",
     "def plot_ln_pi(mu, sigma, xy=(700, 75), title=\"Exact Prediction Interval for Known Lognormal\"):\n",
     "    failure_dist = lognorm(s=sigma, scale=np.exp(mu))\n",
     "    samples = failure_dist.rvs(size=1000, random_state=100)\n",
-    "    fig, axs = plt.subplots(1, 3, figsize=(20, 10))\n",
+    "    fig, axs = plt.subplots(1, 3, figsize=(20, 8))\n",
     "    axs = axs.flatten()\n",
     "    axs[0].hist(samples, ec=\"black\", color=\"slateblue\", bins=30)\n",
     "    axs[0].set_title(f\"Failure Time Distribution: LN({mu}, {sigma})\")\n",
@@ -1970,15 +1967,15 @@
     "def bayes_boot(df, lb, ub, seed=100):\n",
     "    w = np.random.dirichlet(np.ones(len(df)), 1)[0]\n",
     "    lnf = LogNormalFitter().fit(df[\"t\"] + 1e-25, df[\"failed\"], weights=w)\n",
-    "    ## Sample random choice from 95% percentile interval of bootstrapped dist\n",
-    "    # choices = draws['t'].values\n",
-    "    choices = np.linspace(df[\"t\"].min(), df[\"t\"].max(), 1000)\n",
+    "    rv = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_))\n",
+    "    ## Sample random choice from implied bootstrapped distribution\n",
+    "    choices = rv.rvs(1000)\n",
     "    future = random.choice(choices)\n",
     "    ## Check if choice is contained within the MLE 95% PI\n",
     "    contained = (future >= lb) & (future <= ub)\n",
     "    ## Record 95% interval of bootstrapped dist\n",
-    "    lb = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_)).ppf(0.025)\n",
-    "    ub = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_)).ppf(0.975)\n",
+    "    lb = rv.ppf(0.025)\n",
+    "    ub = rv.ppf(0.975)\n",
     "    return lb, ub, contained, future, lnf.sigma_, lnf.mu_"
    ]
   },
@@ -2162,16 +2159,15 @@
     "    draws.sort_values(\"t\", inplace=True)\n",
     "    ## Fit Lognormal Dist to\n",
     "    lnf = LogNormalFitter().fit(draws[\"t\"] + 1e-25, draws[\"failed\"])\n",
-    "    ## Sample random choice from 95% percentile interval of bootstrapped dist\n",
-    "    # choices = draws['t'].values\n",
-    "    ## Essentially sampling from a uniform interval\n",
-    "    choices = np.linspace(draws[\"t\"].min(), draws[\"t\"].max(), 1000)\n",
+    "    rv = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_))\n",
+    "    ## Sample random choice from implied distribution\n",
+    "    choices = rv.rvs(1000)\n",
     "    future = random.choice(choices)\n",
     "    ## Check if choice is contained within the MLE 95% PI\n",
     "    contained = (future >= lb) & (future <= ub)\n",
     "    ## Record 95% interval of bootstrapped dist\n",
-    "    lb = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_)).ppf(0.025)\n",
-    "    ub = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_)).ppf(0.975)\n",
+    "    lb = rv.ppf(0.025)\n",
+    "    ub = rv.ppf(0.975)\n",
     "    return lb, ub, contained, future, lnf.sigma_, lnf.mu_\n",
     "\n",
     "\n",
@@ -2186,7 +2182,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We can use these bootstrapped statistics to further calculate quantities of the predictive distribution."
+    "We can use these bootstrapped statistics to further calculate quantities of the predictive distribution. In our case we could use the parametric CDF for our simple parametric model, but we'll adopt the empirical cdf here to show how this technique can be used when we have more complicated models too."
    ]
   },
   {
@@ -2234,7 +2230,7 @@
     "for i in range(1000):\n",
     "    samples = lognorm(s=draws.iloc[i][\"Sigma\"], scale=np.exp(draws.iloc[i][\"Mu\"])).rvs(1000)\n",
     "    qe, pe = ecdf(samples)\n",
-    "    ax.plot(qe, pe, color=\"grey\", alpha=0.2)\n",
+    "    ax.plot(qe, pe, color=\"skyblue\", alpha=0.2)\n",
     "    lkup = dict(zip(pe, qe))\n",
     "    hist_data.append([lkup[0.05]])\n",
     "hist_data = pd.DataFrame(hist_data, columns=[\"p05\"])\n",
@@ -2246,10 +2242,10 @@
     "ax1.hist(hist_data[\"p05\"], color=\"slateblue\", ec=\"black\", alpha=0.4, bins=30)\n",
     "ax1.set_title(\"Estimate of Uncertainty in the 5% Failure Time\", fontsize=20)\n",
     "ax1.axvline(\n",
-    "    hist_data[\"p05\"].quantile(0.025), color=\"cyan\", label=\"Lower Bound PI for 5% failure time\"\n",
+    "    hist_data[\"p05\"].quantile(0.025), color=\"cyan\", label=\"Lower Bound CI for 5% failure time\"\n",
     ")\n",
     "ax1.axvline(\n",
-    "    hist_data[\"p05\"].quantile(0.975), color=\"cyan\", label=\"Upper Bound PI for 5% failure time\"\n",
+    "    hist_data[\"p05\"].quantile(0.975), color=\"cyan\", label=\"Upper Bound CI for 5% failure time\"\n",
     ")\n",
     "ax1.legend()\n",
     "ax.legend();"
@@ -2365,7 +2361,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "These simulations should be repeated a far larger number of times than we do here. We can also vary the interval size to achieve the desired coverage level."
+    "These simulations should be repeated a far larger number of times than we do here. It should be clear to see how we can also vary the MLE interval size to achieve the desired coverage level."
    ]
   },
   {
@@ -7333,12 +7329,12 @@
     "hist_data_info = pd.DataFrame(hist_data_info, columns=[\"p10\", \"p05\"])\n",
     "draws = pm.draw(pm.Weibull.dist(alpha=np.mean(alphas), beta=np.mean(betas)), 1000)\n",
     "qe, pe = ecdf(draws)\n",
-    "ax.plot(qe, pe, color=\"purple\", label=\"Expected CDF Uninformative\")\n",
+    "ax.plot(qe, pe, color=\"purple\", label=\"Expected CDF Uninformative Prior\")\n",
     "draws = pm.draw(\n",
     "    pm.Weibull.dist(alpha=np.mean(alphas_informative), beta=np.mean(betas_informative)), 1000\n",
     ")\n",
     "qe, pe = ecdf(draws)\n",
-    "ax.plot(qe, pe, color=\"magenta\", label=\"Expected CDF Informative\")\n",
+    "ax.plot(qe, pe, color=\"magenta\", label=\"Expected CDF Informative Prior\")\n",
     "ax.plot(\n",
     "    actuarial_table_bearings[\"t\"],\n",
     "    actuarial_table_bearings[\"logit_CI_95_ub\"],\n",

examples/references.bib

@@ -431,7 +431,7 @@ @book{mcelreath2018statistical
   publisher     = {Chapman and Hall/CRC}
 }
 @book{Meeker2021,
-  author        = {Meeker, William},
+  author        = {Escobar, L.A. and Meeker, W.Q. and Pascual, Francis},
   publisher     = {Wiley},
   title         = {Statistical Methods for Reliability Data},
   year          = {2021}

myst_nbs/case_studies/reliability_and_calibrated_prediction.myst.md

@@ -34,7 +34,6 @@ import numpy as np
 import pandas as pd
 import pymc as pm
 
-from joblib import Parallel, delayed
 from lifelines import KaplanMeierFitter, LogNormalFitter, WeibullFitter
 from lifelines.utils import survival_table_from_events
 from scipy.stats import binom, lognorm, norm, weibull_min
@@ -81,15 +80,13 @@ Throughout the focus will be how the understanding of the CDF can help us unders
 In the study of reliability statistics there is a focus on location-scale based distributions with long tails. In an ideal world we'd know exactly which distribution described our failure process and the prediction interval for the next failure could be defined exactly.
 
 ```{code-cell} ipython3
-from scipy.stats import lognorm
-
 mu, sigma = 6, 0.3
 
 
 def plot_ln_pi(mu, sigma, xy=(700, 75), title="Exact Prediction Interval for Known Lognormal"):
     failure_dist = lognorm(s=sigma, scale=np.exp(mu))
     samples = failure_dist.rvs(size=1000, random_state=100)
-    fig, axs = plt.subplots(1, 3, figsize=(20, 10))
+    fig, axs = plt.subplots(1, 3, figsize=(20, 8))
     axs = axs.flatten()
     axs[0].hist(samples, ec="black", color="slateblue", bins=30)
     axs[0].set_title(f"Failure Time Distribution: LN({mu}, {sigma})")
@@ -422,15 +419,15 @@ The second method we'll use to assess coverage is to bootstrap estimates of a 95
 def bayes_boot(df, lb, ub, seed=100):
     w = np.random.dirichlet(np.ones(len(df)), 1)[0]
     lnf = LogNormalFitter().fit(df["t"] + 1e-25, df["failed"], weights=w)
-    ## Sample random choice from 95% percentile interval of bootstrapped dist
-    # choices = draws['t'].values
-    choices = np.linspace(df["t"].min(), df["t"].max(), 1000)
+    rv = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_))
+    ## Sample random choice from implied bootstrapped distribution
+    choices = rv.rvs(1000)
     future = random.choice(choices)
     ## Check if choice is contained within the MLE 95% PI
     contained = (future >= lb) & (future <= ub)
     ## Record 95% interval of bootstrapped dist
-    lb = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_)).ppf(0.025)
-    ub = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_)).ppf(0.975)
+    lb = rv.ppf(0.025)
+    ub = rv.ppf(0.975)
     return lb, ub, contained, future, lnf.sigma_, lnf.mu_
 ```
 
@@ -440,16 +437,15 @@ def bootstrap(df, lb, ub, seed=100):
     draws.sort_values("t", inplace=True)
     ## Fit Lognormal Dist to
     lnf = LogNormalFitter().fit(draws["t"] + 1e-25, draws["failed"])
-    ## Sample random choice from 95% percentile interval of bootstrapped dist
-    # choices = draws['t'].values
-    ## Essentially sampling from a uniform interval
-    choices = np.linspace(draws["t"].min(), draws["t"].max(), 1000)
+    rv = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_))
+    ## Sample random choice from implied distribution
+    choices = rv.rvs(1000)
     future = random.choice(choices)
     ## Check if choice is contained within the MLE 95% PI
     contained = (future >= lb) & (future <= ub)
     ## Record 95% interval of bootstrapped dist
-    lb = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_)).ppf(0.025)
-    ub = lognorm(s=lnf.sigma_, scale=np.exp(lnf.mu_)).ppf(0.975)
+    lb = rv.ppf(0.025)
+    ub = rv.ppf(0.975)
     return lb, ub, contained, future, lnf.sigma_, lnf.mu_
 
 
@@ -460,7 +456,7 @@ draws = pd.DataFrame(
 draws
 ```
 
-We can use these bootstrapped statistics to further calculate quantities of the predictive distribution.
+We can use these bootstrapped statistics to further calculate quantities of the predictive distribution. In our case we could use the parametric CDF for our simple parametric model, but we'll adopt the empirical cdf here to show how this technique can be used when we have more complicated models too.
 
 ```{code-cell} ipython3
 def ecdf(sample):
@@ -486,7 +482,7 @@ hist_data = []
 for i in range(1000):
     samples = lognorm(s=draws.iloc[i]["Sigma"], scale=np.exp(draws.iloc[i]["Mu"])).rvs(1000)
     qe, pe = ecdf(samples)
-    ax.plot(qe, pe, color="grey", alpha=0.2)
+    ax.plot(qe, pe, color="skyblue", alpha=0.2)
     lkup = dict(zip(pe, qe))
     hist_data.append([lkup[0.05]])
 hist_data = pd.DataFrame(hist_data, columns=["p05"])
@@ -498,10 +494,10 @@ ax.set_title("Bootstrapped CDF functions for the Shock Absorbers Data", fontsize
 ax1.hist(hist_data["p05"], color="slateblue", ec="black", alpha=0.4, bins=30)
 ax1.set_title("Estimate of Uncertainty in the 5% Failure Time", fontsize=20)
 ax1.axvline(
-    hist_data["p05"].quantile(0.025), color="cyan", label="Lower Bound PI for 5% failure time"
+    hist_data["p05"].quantile(0.025), color="cyan", label="Lower Bound CI for 5% failure time"
 )
 ax1.axvline(
-    hist_data["p05"].quantile(0.975), color="cyan", label="Upper Bound PI for 5% failure time"
+    hist_data["p05"].quantile(0.975), color="cyan", label="Upper Bound CI for 5% failure time"
 )
 ax1.legend()
 ax.legend();
@@ -587,7 +583,7 @@ axs[2].annotate(
 );
 ```
 
-These simulations should be repeated a far larger number of times than we do here. We can also vary the interval size to achieve the desired coverage level.
+These simulations should be repeated a far larger number of times than we do here. It should be clear to see how we can also vary the MLE interval size to achieve the desired coverage level.
 
 +++
 
@@ -962,12 +958,12 @@ hist_data = pd.DataFrame(hist_data, columns=["p10", "p05"])
 hist_data_info = pd.DataFrame(hist_data_info, columns=["p10", "p05"])
 draws = pm.draw(pm.Weibull.dist(alpha=np.mean(alphas), beta=np.mean(betas)), 1000)
 qe, pe = ecdf(draws)
-ax.plot(qe, pe, color="purple", label="Expected CDF Uninformative")
+ax.plot(qe, pe, color="purple", label="Expected CDF Uninformative Prior")
 draws = pm.draw(
     pm.Weibull.dist(alpha=np.mean(alphas_informative), beta=np.mean(betas_informative)), 1000
 )
 qe, pe = ecdf(draws)
-ax.plot(qe, pe, color="magenta", label="Expected CDF Informative")
+ax.plot(qe, pe, color="magenta", label="Expected CDF Informative Prior")
 ax.plot(
     actuarial_table_bearings["t"],
     actuarial_table_bearings["logit_CI_95_ub"],