Make stochastic analyses reproducible

notebooks/logit_all_models.ipynb

@@ -106,7 +106,10 @@
    "source": [
     "# Two-sided 95% confidence intervals used throughout the notebook\n",
     "alpha = 0.05\n",
-    "probs = [alpha / 2, 1 - alpha / 2]\n"
+    "probs = [alpha / 2, 1 - alpha / 2]\n",
+    "\n",
+    "# Reproducibility: pass this seed to every stochastic calculation.\n",
+    "random_seed = logistic.DEFAULT_RANDOM_SEED\n"
    ]
   },
   {
@@ -170,7 +173,7 @@
     "logit = logistic.LogisticPolyRegression(degree=1, mono=False)\n",
     "\n",
     "# fits\n",
-    "fit_results = [logit.fit(x, y, method=\"local\") for x in xs]\n",
+    "fit_results = [logit.fit(x, y, method=\"local\", seed=random_seed) for x in xs]\n",
     "if not all(result[\"success\"] for result in fit_results):\n",
     "    raise RuntimeError(\"At least one linear logistic-regression fit failed.\")\n",
     "\n",
@@ -562,7 +565,8 @@
     "# goodness of fit measures\n",
     "index = pd.Index(scales, name = 'scale')\n",
     "pd.DataFrame(\n",
-    "    [logit.goodness_of_fit(x, y, theta) for x, theta in zip(xs, thetas)],\n",
+    "    [logit.goodness_of_fit(x, y, theta, bootstrap_seed=random_seed)\n",
+    "     for x, theta in zip(xs, thetas)],\n",
     "    index=index,\n",
     ")"
    ]
@@ -601,7 +605,9 @@
     "fig, axs = plt.subplots(ncols = len(scales), figsize = (6*len(scales), 5))\n",
     "\n",
     "quantile_methods = {\n",
-    "    \"normal\": logit.get_model_quantiles_normal,\n",
+    "    \"normal\": lambda x, p, theta, cov: logit.get_model_quantiles_normal(\n",
+    "        x, p, theta, cov, seed=random_seed\n",
+    "    ),\n",
     "    \"delta\": logit.get_model_quantiles_delta,\n",
     "}\n",
     "\n",
@@ -715,7 +721,7 @@
     "logit = logistic.LogisticPolyRegression(degree=3, mono=False)\n",
     "\n",
     "# fits\n",
-    "fit_results = [logit.fit(x, y, method=\"local\") for x in xs]\n",
+    "fit_results = [logit.fit(x, y, method=\"local\", seed=random_seed) for x in xs]\n",
     "if not all(result[\"success\"] for result in fit_results):\n",
     "    raise RuntimeError(\"At least one cubic logistic-regression fit failed.\")\n",
     "\n",
@@ -1174,7 +1180,8 @@
     "# goodness of fit measures\n",
     "index = pd.Index(scales, name = 'scale')\n",
     "pd.DataFrame(\n",
-    "    [logit.goodness_of_fit(x, y, theta) for x, theta in zip(xs, thetas)],\n",
+    "    [logit.goodness_of_fit(x, y, theta, bootstrap_seed=random_seed)\n",
+    "     for x, theta in zip(xs, thetas)],\n",
     "    index=index,\n",
     ")"
    ]
@@ -1211,7 +1218,9 @@
     "fig, axs = plt.subplots(ncols = len(scales), figsize = (6*len(scales), 5))\n",
     "\n",
     "quantile_methods = {\n",
-    "    \"normal\": logit.get_model_quantiles_normal,\n",
+    "    \"normal\": lambda x, p, theta, cov: logit.get_model_quantiles_normal(\n",
+    "        x, p, theta, cov, seed=random_seed\n",
+    "    ),\n",
     "    \"delta\": logit.get_model_quantiles_delta,\n",
     "}\n",
     "\n",
@@ -1447,7 +1456,7 @@
     "logit = logistic.LogisticPolyRegression(degree=3, mono=True, lam=(0.0, 1e-6))\n",
     "\n",
     "# fits\n",
-    "ress = [logit.fit(x, y, method=\"diff_evol\") for x in xs]\n",
+    "ress = [logit.fit(x, y, method=\"diff_evol\", seed=random_seed) for x in xs]\n",
     "if not all(result[\"success\"] for result in ress):\n",
     "    raise RuntimeError(\"At least one monotonic cubic logistic-regression fit failed.\")\n",
     "\n",
@@ -2143,7 +2152,8 @@
    "source": [
     "# goodness of fit measures\n",
     "pd.DataFrame(\n",
-    "    [logit.goodness_of_fit(x, y, theta) for x, theta in zip(xs, thetas)], \n",
+    "    [logit.goodness_of_fit(x, y, theta, bootstrap_seed=random_seed)\n",
+    "     for x, theta in zip(xs, thetas)], \n",
     "    index = pd.Index(scales, name = 'scale')\n",
     ")"
    ]
@@ -2182,7 +2192,9 @@
     "fig, axs = plt.subplots(ncols = len(scales), figsize = (6*len(scales), 5))\n",
     "\n",
     "quantile_methods = {\n",
-    "    \"normal\": logit.get_model_quantiles_normal,\n",
+    "    \"normal\": lambda x, p, theta, cov: logit.get_model_quantiles_normal(\n",
+    "        x, p, theta, cov, seed=random_seed\n",
+    "    ),\n",
     "    \"delta\": logit.get_model_quantiles_delta,\n",
     "}\n",
     "\n",
@@ -2275,7 +2287,9 @@
     "\n",
     "boots_theta_results = {}\n",
     "quantile_methods = {\n",
-    "    \"normal\": logit.get_model_quantiles_normal,\n",
+    "    \"normal\": lambda x, p, theta, cov: logit.get_model_quantiles_normal(\n",
+    "        x, p, theta, cov, seed=random_seed\n",
+    "    ),\n",
     "    \"delta\": logit.get_model_quantiles_delta,\n",
     "}\n",
     "parameter_generators = {\n",
@@ -2328,7 +2342,9 @@
     "        print(f\"model CI:scale:{scale},boots-method:{method}\")\n",
     "\n",
     "        # generate sampled or bootstrapped parameters\n",
-    "        btheta = parameter_generators[method](x, y, n_parameter_samples)\n",
+    "        btheta = parameter_generators[method](\n",
+    "            x, y, n_parameter_samples, seed=random_seed\n",
+    "        )\n",
     "    \n",
     "        boots_theta_results[(scale, method)] = btheta\n",
     "\n",

notebooks/logit_mono-cubic4paper.ipynb

@@ -57,6 +57,8 @@
     "\n",
     "# general setting\n",
     "np.set_printoptions(precision=16)\n",
+    "# Reproducibility: pass this seed to every stochastic calculation.\n",
+    "random_seed = logistic.DEFAULT_RANDOM_SEED\n",
     "width=5\n",
     "\n",
     "height=3.6\n",
@@ -1157,7 +1159,7 @@
     "    X, Y = g['X'].to_numpy(), g['Y'].to_numpy()\n",
     "\n",
     "    # fit\n",
-    "    res_fit = logit.fit(X, Y, method=\"diff_evol\")\n",
+    "    res_fit = logit.fit(X, Y, method=\"diff_evol\", seed=random_seed)\n",
     "    \n",
     "    # calc nllf and beta for fitted theta\n",
     "    res_nllf = logit.get_nllf(X, Y, res_fit['theta'], jac = True)\n",
@@ -1178,7 +1180,9 @@
     "    theta_CI = logit.get_theta_quantiles_normal(probs, res_fit['theta'], cov)\n",
     "    \n",
     "    # goodness of fit\n",
-    "    res_gof = logit.goodness_of_fit(X, Y, res_fit['theta'])\n",
+    "    res_gof = logit.goodness_of_fit(\n",
+    "        X, Y, res_fit['theta'], bootstrap_seed=random_seed\n",
+    "    )\n",
     "\n",
     "    lst.append(\n",
     "        {\"scale\": scale, \"dataset\": dataset, \"beta\": beta} | \n",
@@ -6711,6 +6715,7 @@
     "    width=width,\n",
     "    height=height,\n",
     "    jitter_y=0.02,\n",
+    "    random_state=random_seed,\n",
     "    alpha_full=0.50,\n",
     "    alpha_trim=0.50,\n",
     "    plot_labs={\"log\": \"C\", \"plain\": \"D\"}\n",
@@ -7042,6 +7047,7 @@
     "    lg=logit,\n",
     "    n_grid=100,\n",
     "    verbose=verbose,\n",
+    "    seed=random_seed,\n",
     ")\n",
     "\n",
     "reporting.plot_from_dataframes(\n",
@@ -7114,6 +7120,7 @@
     "    n_boots=10000,\n",
     "    methods=(\"normal\", \"nonparam_boots\", \"nonparam_stratified_boots\", \"parametric_boots\"),\n",
     "    verbose=True,\n",
+    "    seed=random_seed,\n",
     ")"
    ]
   },
@@ -7138,6 +7145,7 @@
     "    analytic_methods=(\"normal\", \"delta\"),\n",
     "    bootstrap_methods=(\"nonparam_boots\", \"parametric_boots\"),\n",
     "    verbose=verbose,\n",
+    "    seed=random_seed,\n",
     ")"
    ]
   },

notebooks/logit_review_boots.ipynb

@@ -868,6 +868,8 @@
     "y = df_xy['y'].to_numpy(dtype=int)\n",
     "\n",
     "n_boots = 10_000\n",
+    "# Reproducibility: derive all stochastic calculations from one fixed seed.\n",
+    "random_seed = logistic.DEFAULT_RANDOM_SEED\n",
     "df_xy"
    ]
   },
@@ -1038,7 +1040,7 @@
    ],
    "source": [
     "logit = logistic.LogisticPolyRegression(degree=1)\n",
-    "res_fit = logit.fit(x, y, method='diff_evol')\n",
+    "res_fit = logit.fit(x, y, method='diff_evol', seed=random_seed)\n",
     "if not res_fit['success']:\n",
     "    raise RuntimeError(f\"Logistic fit failed: {res_fit}\")\n",
     "\n",
@@ -1197,7 +1199,9 @@
     }
    ],
    "source": [
-    "pd.DataFrame([logit.goodness_of_fit(x, y, theta)])"
+    "pd.DataFrame([\n",
+    "    logit.goodness_of_fit(x, y, theta, bootstrap_seed=random_seed)\n",
+    "])"
    ]
   },
   {
@@ -1225,9 +1229,16 @@
     "\n",
     "ax.plot(xp, logit.model(xp, theta), color='black', label='logistic fit', zorder=4)\n",
     "\n",
-    "for label, method, color in [('normal', logit.get_model_quantiles_normal, 'tab:red'),\n",
-    "                             ('delta', logit.get_model_quantiles_delta, 'tab:green')]:\n",
-    "    lo, hi = method(xp, probs, theta, cov_theta)\n",
+    "quantile_methods = {\n",
+    "    'normal': lambda xgrid: logit.get_model_quantiles_normal(\n",
+    "        xgrid, probs, theta, cov_theta, seed=random_seed\n",
+    "    ),\n",
+    "    'delta': lambda xgrid: logit.get_model_quantiles_delta(\n",
+    "        xgrid, probs, theta, cov_theta\n",
+    "    ),\n",
+    "}\n",
+    "for label, color in [('normal', 'tab:red'), ('delta', 'tab:green')]:\n",
+    "    lo, hi = quantile_methods[label](xp)\n",
     "    ax.fill_between(xp, lo, hi, color=color, alpha=0.18, label=f'95% CI ({label})')\n",
     "\n",
     "ax.set(xlabel='maximum SUV percentile across visits', \n",
@@ -1264,9 +1275,15 @@
     "# Each method returns exactly n_boots bootstrap estimates. The original-data\n",
     "# MLE is used as the initial point for fitting each replicate but is not\n",
     "# included as a row in the returned bootstrap array.\n",
-    "btheta_nonpar = logit.get_nonparam_boots_theta(x, y, n_boots, seed=1977)\n",
-    "btheta_nonpar_strat = logit.get_nonparam_stratified_boots_theta(x, y, n_boots, seed=1978)\n",
-    "btheta_param = logit.get_parametric_boots_theta(x, y, n_boots, seed=1979)\n",
+    "btheta_nonpar = logit.get_nonparam_boots_theta(\n",
+    "    x, y, n_boots, seed=random_seed\n",
+    ")\n",
+    "btheta_nonpar_strat = logit.get_nonparam_stratified_boots_theta(\n",
+    "    x, y, n_boots, seed=random_seed + 1\n",
+    ")\n",
+    "btheta_param = logit.get_parametric_boots_theta(\n",
+    "    x, y, n_boots, seed=random_seed + 2\n",
+    ")\n",
     "\n",
     "print(btheta_nonpar.shape, btheta_nonpar_strat.shape, btheta_param.shape)"
    ]
@@ -1501,7 +1518,9 @@
     "ax.plot(xp, logit.model(xp, theta), color='black', lw=2, label='logistic fit', zorder=5)\n",
     "\n",
     "bands = [\n",
-    "    ('normal', logit.get_model_quantiles_normal(xp, probs, theta, cov_theta), 'tab:red'),\n",
+    "    ('normal', logit.get_model_quantiles_normal(\n",
+    "        xp, probs, theta, cov_theta, seed=random_seed\n",
+    "    ), 'tab:red'),\n",
     "    ('delta', logit.get_model_quantiles_delta(xp, probs, theta, cov_theta), 'tab:green'),\n",
     "    ('non-parametric bootstrap', bootstrap_band(xp, btheta_nonpar, probs), 'tab:blue'),\n",
     "    ('stratified bootstrap', bootstrap_band(xp, btheta_nonpar_strat, probs), 'tab:purple'),\n",

src/irae_risk/logistic.py

@@ -63,6 +63,13 @@ import scipy.optimize
 
 from . import monotonic as mc
 
+
+# Shared deterministic seed for every stochastic calculation in this module.
+# Callers can still supply a different seed explicitly when independent runs
+# are required.
+DEFAULT_RANDOM_SEED = 1977
+
+
 def resize_with_const(v, n, val=0):
     """
     Resize a 1D vector to a specified length `n`.
@@ -577,12 +584,21 @@ class LogisticPolyRegression:
         Return:
             dict {"theta", "cost", "success"}
     """
-    def fit(self, x, y, theta0 = None, method = "local", seed = None):
+    def fit(
+        self,
+        x,
+        y,
+        theta0=None,
+        method="local",
+        seed=DEFAULT_RANDOM_SEED,
+    ):
         """Fit the model using a local or global optimization method.
 
         ``theta0`` is used as the initial point for local optimization.  If it is
         omitted, an initial estimate is calculated from the data.  ``seed`` is
-        passed to stochastic global optimizers for reproducible fits.
+        passed to stochastic global optimizers for reproducible fits and
+        defaults to ``DEFAULT_RANDOM_SEED``.  Pass ``None`` explicitly only
+        when a nondeterministic run is desired.
         """
         x, y = self._validate_data(x, y)
         bnds = self._get_bounds()
@@ -715,7 +731,7 @@ class LogisticPolyRegression:
         thresh = 0.5,
         regularization = False,
         bootstrap_samples = 1000,
-        bootstrap_seed = 1977,
+        bootstrap_seed = DEFAULT_RANDOM_SEED,
         bootstrap_max_attempts = None,
     ):
         """Calculate fit summaries and a bootstrap goodness-of-fit test.
@@ -977,8 +993,16 @@ class LogisticPolyRegression:
             array of mxn floats
     """
 
-    def get_model_quantiles_normal(self, x, probs, mean_theta, cov_theta,
-                                   exact = True, seed = 1977, m = 10**5):
+    def get_model_quantiles_normal(
+        self,
+        x,
+        probs,
+        mean_theta,
+        cov_theta,
+        exact=True,
+        seed=DEFAULT_RANDOM_SEED,
+        m=10**5,
+    ):
 
         mean_beta = self.get_beta(mean_theta)
         X = np.column_stack([x**i for i in range(len(mean_beta))])
@@ -1078,7 +1102,7 @@ class LogisticPolyRegression:
         Return:
             array of mx(degree + 1)
     """
-    def get_normal_theta(self, x, y, m, seed = 1977):
+    def get_normal_theta(self, x, y, m, seed=DEFAULT_RANDOM_SEED):
         x, y = self._validate_data(x, y)
         m = self._validate_sample_count(m)
 
@@ -1119,7 +1143,7 @@ class LogisticPolyRegression:
         x,
         y,
         m,
-        seed = 1977,
+        seed=DEFAULT_RANDOM_SEED,
         max_attempts = None,
     ):
         x, y = self._validate_data(x, y)
@@ -1163,7 +1187,7 @@ class LogisticPolyRegression:
         x,
         y,
         m,
-        seed = 1977,
+        seed=DEFAULT_RANDOM_SEED,
         max_attempts = None,
     ):
         x, y = self._validate_data(x, y)
@@ -1217,7 +1241,7 @@ class LogisticPolyRegression:
         x,
         y,
         m,
-        seed = 1977,
+        seed=DEFAULT_RANDOM_SEED,
         max_attempts = None,
     ):
         x, y = self._validate_data(x, y)

src/irae_risk/reporting.py

@@ -4,6 +4,8 @@ import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
 
+from .logistic import DEFAULT_RANDOM_SEED
+
 # =========================
 # small generic helpers
 # =========================
@@ -159,9 +161,20 @@ def build_cost_perturbation_data(
     return pd.concat(rows, ignore_index=True)
 
 
-def build_plot_data(df_data, df_fit_index, alpha, lg, n_grid=100, verbose=False):
+def build_plot_data(
+    df_data,
+    df_fit_index,
+    alpha,
+    lg,
+    n_grid=100,
+    verbose=False,
+    seed=DEFAULT_RANDOM_SEED,
+):
     """
     Build raw points, fitted curves, and analytic confidence intervals.
+
+    ``seed`` is forwarded to the sampling-based normal interval used by a
+    monotonic cubic model.
     """
     probs = [alpha / 2, 1 - alpha / 2]
     confidence = 100 * (1 - alpha)
@@ -184,7 +197,12 @@ def build_plot_data(df_data, df_fit_index, alpha, lg, n_grid=100, verbose=False)
 
         for method in ("normal", "delta"):
             qfun = getattr(lg, f"get_model_quantiles_{method}")
-            y_low, y_high = qfun(x_fit, probs, theta, cov)
+            if method == "normal":
+                y_low, y_high = qfun(
+                    x_fit, probs, theta, cov, seed=seed
+                )
+            else:
+                y_low, y_high = qfun(x_fit, probs, theta, cov)
 
             cis.append(pd.DataFrame({
                 "scale": scale,
@@ -211,9 +229,11 @@ def load_or_build_bootstrap_theta(
     n_boots=10000,
     methods=("normal", "nonparam_boots", "nonparam_stratified_boots", "parametric_boots"),
     verbose=False,
+    seed=DEFAULT_RANDOM_SEED,
 ):
     """
-    Load bootstrap parameter draws from disk or build them.
+    Load bootstrap parameter draws from disk or build them.  When the cache is
+    absent, ``seed`` is forwarded to every parameter generator.
     """
     result_file = os.path.join(results_path, filename)
 
@@ -236,7 +256,9 @@ def load_or_build_bootstrap_theta(
 
             theta_fun = getattr(lg, f"get_{method}_theta")
             try:
-                out[(scale, dataset, method)] = theta_fun(x, y, m=n_boots)
+                out[(scale, dataset, method)] = theta_fun(
+                    x, y, m=n_boots, seed=seed
+                )
             except Exception:
                 if verbose:
                     print("\tfailed")
@@ -258,9 +280,12 @@ def build_plot_data_full_ci(
     analytic_methods=("normal", "delta"),
     bootstrap_methods=("nonparam_boots", "parametric_boots"),
     verbose=False,
+    seed=DEFAULT_RANDOM_SEED,
 ):
     """
-    Build raw points, fitted curves, analytic CI, and bootstrap CI.
+    Build raw points, fitted curves, analytic CI, and bootstrap CI.  ``seed``
+    is forwarded to the sampling-based normal interval used by a monotonic
+    cubic model.
     """
     probs = [alpha / 2, 1 - alpha / 2]
     confidence = 100 * (1 - alpha)
@@ -286,7 +311,12 @@ def build_plot_data_full_ci(
                 print(f"\tanalytic CI: {method}")
 
             qfun = getattr(lg, f"get_model_quantiles_{method}")
-            y_low, y_high = qfun(x_fit, probs, theta, cov)
+            if method == "normal":
+                y_low, y_high = qfun(
+                    x_fit, probs, theta, cov, seed=seed
+                )
+            else:
+                y_low, y_high = qfun(x_fit, probs, theta, cov)
 
             cis.append(pd.DataFrame({
                 "scale": scale,
@@ -354,7 +384,7 @@ def plot_from_dataframes_basic(
     jitter_y=0.0,
     jitter_y_full=None,
     jitter_y_trim=None,
-    random_state=123,
+    random_state=DEFAULT_RANDOM_SEED,
     plot_labs = None,
 ):
     """

tests/test_logistic.py

@@ -1,7 +1,9 @@
+import inspect
+
 import numpy as np
 import pytest
 
-from irae_risk.logistic import LogisticPolyRegression
+from irae_risk.logistic import DEFAULT_RANDOM_SEED, LogisticPolyRegression
 
 
 @pytest.fixture
@@ -12,6 +14,35 @@ def regression_data():
     return x, y, theta
 
 
+def test_stochastic_methods_share_a_fixed_default_seed():
+    methods_and_seed_parameters = (
+        (LogisticPolyRegression.fit, "seed"),
+        (LogisticPolyRegression.goodness_of_fit, "bootstrap_seed"),
+        (LogisticPolyRegression.get_model_quantiles_normal, "seed"),
+        (LogisticPolyRegression.get_normal_theta, "seed"),
+        (LogisticPolyRegression.get_nonparam_boots_theta, "seed"),
+        (LogisticPolyRegression.get_nonparam_stratified_boots_theta, "seed"),
+        (LogisticPolyRegression.get_parametric_boots_theta, "seed"),
+    )
+
+    for method, parameter in methods_and_seed_parameters:
+        default = inspect.signature(method).parameters[parameter].default
+        assert default == DEFAULT_RANDOM_SEED
+
+
+def test_differential_evolution_fit_is_reproducible_by_default(
+    regression_data,
+):
+    x, y, _ = regression_data
+    model = LogisticPolyRegression(degree=1)
+
+    first = model.fit(x, y, method="diff_evol")
+    second = model.fit(x, y, method="diff_evol")
+
+    np.testing.assert_array_equal(first["theta"], second["theta"])
+    assert first["cost"] == second["cost"]
+
+
 def test_x50_and_s50_for_linear_logistic_model():
     model = LogisticPolyRegression(degree=1)
     theta = np.array([-2.0, 0.5])

tests/test_reporting.py

@@ -2,7 +2,7 @@ import numpy as np
 import pandas as pd
 
 from irae_risk import reporting
-from irae_risk.logistic import LogisticPolyRegression
+from irae_risk.logistic import DEFAULT_RANDOM_SEED, LogisticPolyRegression
 
 
 def _reporting_frames():
@@ -78,8 +78,9 @@ def test_bootstrap_loader_calls_theta_api(tmp_path):
     df_data, _ = _reporting_frames()
 
     class ThetaSampler:
-        def get_normal_theta(self, x, y, m):
+        def get_normal_theta(self, x, y, m, seed):
             assert len(x) == len(y) == 4
+            assert seed == DEFAULT_RANDOM_SEED
             return np.full((m, 2), 1.5)
 
     result = reporting.load_or_build_bootstrap_theta(