Renaming files and trying to finish figures for paper

python/logistic/logit_reg_boots.ipynb → python/logistic/logit_boots.ipynb

@@ -28,7 +28,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "3b351af7",
    "metadata": {},
    "outputs": [],
@@ -45,7 +45,7 @@
     "\n",
     "# our libs\n",
     "import data_utils\n",
-    "import logit_utils\n",
+    "import logit\n",
     "\n",
     "# testing multivariate normal\n",
     "import pingouin as pg\n",
@@ -125,7 +125,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "id": "239ae789",
    "metadata": {
     "tags": []
@@ -142,13 +142,13 @@
    "source": [
     "# fits\n",
     "logr = linear_model.LogisticRegression(penalty = None)\n",
-    "pars = logit_utils.logit_poly_fit(logr, x, y, degree = 1)\n",
+    "pars = logit.logit_poly_fit(logr, x, y, degree = 1)\n",
     "print(f\"{pars = }\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": null,
    "id": "9cf6fec1",
    "metadata": {},
    "outputs": [
@@ -163,7 +163,7 @@
    ],
    "source": [
     "# asymptotic covariance matrix of parameters\n",
-    "cov_pars = logit_utils.logit_poly_cov(x, pars)\n",
+    "cov_pars = logit.logit_poly_cov(x, pars)\n",
     "print(f\"{cov_pars = }\")"
    ]
   },
@@ -181,7 +181,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "id": "4ba4d5d4",
    "metadata": {},
    "outputs": [
@@ -196,7 +196,7 @@
    ],
    "source": [
     "# CI of params (assuming asymptotic distr of parameters)\n",
-    "pars_CI = logit_utils.logit_poly_pars_quantiles_normal(probs, pars, cov_pars)\n",
+    "pars_CI = logit.logit_poly_pars_quantiles_normal(probs, pars, cov_pars)\n",
     "print(f\"{pars_CI = }\")"
    ]
   },
@@ -310,7 +310,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": null,
    "id": "6ef6f46b",
    "metadata": {},
    "outputs": [
@@ -451,12 +451,12 @@
    ],
    "source": [
     "# goodness of fit measures\n",
-    "pd.DataFrame([logit_utils.logit_poly_goodness_of_fit(x, y, pars)])"
+    "pd.DataFrame([logit.logit_poly_goodness_of_fit(x, y, pars)])"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": null,
    "id": "54f463d0",
    "metadata": {
     "tags": []
@@ -487,13 +487,13 @@
     "\n",
     "# plotting fitted model \n",
     "xp = np.linspace(0, 6, 100)\n",
-    "yp = logit_utils.logit_poly_model(xp, pars)\n",
+    "yp = logit.logit_poly_model(xp, pars)\n",
     "ax.plot(xp, yp, label = \"logistic reg\")\n",
     "\n",
     "for lab, c in zip([\"normal\", \"delta\"], [\"red\", \"green\"]):\n",
     "    \n",
     "    # define quantile model\n",
-    "    qm = eval(f\"logit_utils.logit_poly_model_quantiles_{lab}\")\n",
+    "    qm = eval(f\"logit.logit_poly_model_quantiles_{lab}\")\n",
     "    \n",
     "    # get result of model quantiles\n",
     "    res = qm(xp, probs, pars, cov_pars)\n",
@@ -598,7 +598,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": null,
    "id": "2a860703",
    "metadata": {},
    "outputs": [
@@ -614,7 +614,7 @@
    "source": [
     "# performing bootstrapping\n",
     "m = n_boots\n",
-    "bpars_nonpar = logit_utils.get_nonparam_boots_pars(logr, x, y, m)\n",
+    "bpars_nonpar = logit.get_nonparam_boots_pars(logr, x, y, m)\n",
     "\n",
     "print(f\"{bpars_nonpar.shape = }\")\n",
     "\n",
@@ -776,7 +776,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": null,
    "id": "b0e015ea",
    "metadata": {},
    "outputs": [
@@ -791,7 +791,7 @@
    "source": [
     "# performing stratified bootstrapping\n",
     "m = n_boots\n",
-    "bpars_nonpar_strat = logit_utils.get_nonparam_stratified_boots_pars(logr, x, y, m)\n",
+    "bpars_nonpar_strat = logit.get_nonparam_stratified_boots_pars(logr, x, y, m)\n",
     "\n",
     "print(f\"{bpars_nonpar_strat.shape= }\")"
    ]
@@ -903,7 +903,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": null,
    "id": "88415c37",
    "metadata": {},
    "outputs": [
@@ -918,7 +918,7 @@
    "source": [
     "# performing stratified bootstrapping\n",
     "m = n_boots\n",
-    "bpars_param = logit_utils.get_parametric_boots_pars(logr, x, y, m)\n",
+    "bpars_param = logit.get_parametric_boots_pars(logr, x, y, m)\n",
     "\n",
     "print(f\"{bpars_param.shape= }\")"
    ]
@@ -990,7 +990,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": null,
    "id": "924d6463",
    "metadata": {},
    "outputs": [
@@ -1027,7 +1027,7 @@
     "\n",
     "# plotting fitted model \n",
     "xp = np.linspace(min(x), max(x), 100)\n",
-    "yp = logit_utils.logit_poly_model(xp, pars)\n",
+    "yp = logit.logit_poly_model(xp, pars)\n",
     "ax.plot(xp, yp, label = \"logistic reg\")\n",
     "\n",
     "# CI based on asymptotic MLE distribution\n",
@@ -1035,7 +1035,7 @@
     "                  [\"red\", \"green\"]):\n",
     "    \n",
     "    # define quantile model\n",
-    "    qm = eval(f\"logit_utils.logit_poly_model_quantiles_{lab}\")\n",
+    "    qm = eval(f\"logit.logit_poly_model_quantiles_{lab}\")\n",
     "    \n",
     "    # get result of model quantiles\n",
     "    res = qm(xp, probs, pars, cov_pars)\n",
@@ -1049,7 +1049,7 @@
     "                         [\"green\", \"orange\", \"purple\"] ):\n",
     "\n",
     "    # calculate model values of for all bootstrapped parameters\n",
-    "    model_values = [logit_utils.logit_poly_model(xp, p) for p in bpars]\n",
+    "    model_values = [logit.logit_poly_model(xp, p) for p in bpars]\n",
     "    res = np.quantile(model_values, probs, axis = 0)\n",
     "\n",
     "    # make plot of quantiles\n",

python/logistic/logit_reg_fit.ipynb → python/logistic/logit_fit.ipynb

@@ -20,7 +20,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "id": "3b351af7",
    "metadata": {},
    "outputs": [],
@@ -35,7 +35,7 @@
     "\n",
     "# our libs\n",
     "import data_utils\n",
-    "import logit_utils"
+    "import logit"
    ]
   },
   {
@@ -157,7 +157,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "id": "239ae789",
    "metadata": {
     "tags": []
@@ -174,13 +174,13 @@
    ],
    "source": [
     "# fits\n",
-    "parss = np.array([logit_utils.logit_poly_fit(x, y, degree = 1) for x in xs])\n",
+    "parss = np.array([logit.logit_poly_fit(x, y, degree = 1) for x in xs])\n",
     "print(f\"{parss = }\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "id": "9cf6fec1",
    "metadata": {},
    "outputs": [
@@ -198,7 +198,7 @@
    ],
    "source": [
     "# asymptotic covariance matrix of parameters\n",
-    "cov_parss = np.array([logit_utils.logit_poly_cov(x, pars) for x, pars in zip(xs, parss)])\n",
+    "cov_parss = np.array([logit.logit_poly_cov(x, pars) for x, pars in zip(xs, parss)])\n",
     "print(f\"{cov_parss = }\")"
    ]
   },
@@ -216,7 +216,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "id": "4ba4d5d4",
    "metadata": {},
    "outputs": [
@@ -234,7 +234,7 @@
    ],
    "source": [
     "# CI of params (assuming asymptotic distr of parameters)\n",
-    "pars_CIs = np.array([logit_utils.logit_poly_pars_quantiles_normal(probs, pars, cov_pars) for pars, cov_pars in zip(parss, cov_parss)])\n",
+    "pars_CIs = np.array([logit.logit_poly_pars_quantiles_normal(probs, pars, cov_pars) for pars, cov_pars in zip(parss, cov_parss)])\n",
     "print(f\"{pars_CIs = }\")"
    ]
   },
@@ -387,7 +387,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": null,
    "id": "6ef6f46b",
    "metadata": {},
    "outputs": [
@@ -569,12 +569,12 @@
    "source": [
     "# goodness of fit measures\n",
     "index = pd.Index(labs, name = 'scale')\n",
-    "pd.DataFrame([logit_utils.logit_poly_goodness_of_fit(x, y, pars) for x, pars in zip(xs, parss)], index = index)"
+    "pd.DataFrame([logit.logit_poly_goodness_of_fit(x, y, pars) for x, pars in zip(xs, parss)], index = index)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": null,
    "id": "54f463d0",
    "metadata": {
     "tags": []
@@ -612,14 +612,14 @@
     "\n",
     "    # plotting fitted model \n",
     "    xp = np.linspace(min(x), max(x), 100)\n",
-    "    yp = logit_utils.logit_poly_model(xp, pars)\n",
+    "    yp = logit_sim.logit_poly_model(xp, pars)\n",
     "\n",
     "    ax.plot(xp, yp, label = \"logistic reg\")\n",
     "\n",
     "    for lab, c in zip([\"normal\", \"delta\"], [\"red\", \"green\"]):\n",
     "        \n",
     "        # define quantile model\n",
-    "        fname = f\"logit_utils.logit_poly_model_quantiles_{lab}\"\n",
+    "        fname = f\"logit.logit_poly_model_quantiles_{lab}\"\n",
     "        \n",
     "        # get result of model quantiles at probs\n",
     "        # pars and cov_pars are obtained via MLE method\n",
@@ -680,7 +680,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": null,
    "id": "896ca504",
    "metadata": {},
    "outputs": [
@@ -695,13 +695,13 @@
    ],
    "source": [
     "# fits\n",
-    "parss = np.array([logit_utils.logit_poly_fit(x, y, degree = 3) for x in xs])\n",
+    "parss = np.array([logit.logit_poly_fit(x, y, degree = 3) for x in xs])\n",
     "print(f\"{parss = }\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": null,
    "id": "45404e51",
    "metadata": {},
    "outputs": [
@@ -731,7 +731,7 @@
    ],
    "source": [
     "# asymptotic covariance matrix of parameters\n",
-    "cov_parss = np.array([logit_utils.logit_poly_cov(x, pars) for x, pars in zip(xs, parss)])\n",
+    "cov_parss = np.array([logit.logit_poly_cov(x, pars) for x, pars in zip(xs, parss)])\n",
     "print(f\"{cov_parss = }\")"
    ]
   },
@@ -749,7 +749,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": null,
    "id": "f3f7a0cc",
    "metadata": {},
    "outputs": [
@@ -767,7 +767,7 @@
    ],
    "source": [
     "# CI of params (assuming asymptotic distr of parameters)\n",
-    "pars_CIs = np.array([logit_utils.logit_poly_pars_quantiles_normal(probs, pars, cov_pars) for pars, cov_pars in zip(parss, cov_parss)])\n",
+    "pars_CIs = np.array([logit.logit_poly_pars_quantiles_normal(probs, pars, cov_pars) for pars, cov_pars in zip(parss, cov_parss)])\n",
     "print(f\"{pars_CIs = }\")"
    ]
   },
@@ -972,7 +972,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": null,
    "id": "d88f2c91",
    "metadata": {},
    "outputs": [
@@ -1154,12 +1154,12 @@
    "source": [
     "# goodness of fit measures\n",
     "index = pd.Index(labs, name = 'scale')\n",
-    "pd.DataFrame([logit_utils.logit_poly_goodness_of_fit(x, y, pars) for x, pars in zip(xs, parss)], index = index)"
+    "pd.DataFrame([logit.logit_poly_goodness_of_fit(x, y, pars) for x, pars in zip(xs, parss)], index = index)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": null,
    "id": "c2f9bbbd",
    "metadata": {},
    "outputs": [
@@ -1195,14 +1195,14 @@
     "\n",
     "    # plotting fitted model \n",
     "    xp = np.linspace(min(x), max(x), 100)\n",
-    "    yp = logit_utils.logit_poly_model(xp, pars)\n",
+    "    yp = logit.logit_poly_model(xp, pars)\n",
     "\n",
     "    ax.plot(xp, yp, label = \"logistic reg\")\n",
     "\n",
     "    for lab, c in zip([\"normal\", \"delta\"], [\"red\", \"green\"]):\n",
     "        \n",
     "        # define quantile model\n",
-    "        fname = f\"logit_utils.logit_poly_model_quantiles_{lab}\"\n",
+    "        fname = f\"logit.logit_poly_model_quantiles_{lab}\"\n",
     "        \n",
     "        # get result of model quantiles at probs\n",
     "        # pars and cov_pars are obtained via MLE method\n",

python/logistic/logit_reg_fit_gen.ipynb → python/logistic/logit_gen.ipynb

@@ -22,7 +22,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "3b351af7",
    "metadata": {},
    "outputs": [],
@@ -37,7 +37,7 @@
     "\n",
     "# our libs\n",
     "import data_utils\n",
-    "import logit_utils_gen\n",
+    "import logit_gen\n",
     "\n",
     "np.set_printoptions(precision=16)"
    ]
@@ -116,7 +116,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "id": "239ae789",
    "metadata": {
     "tags": []
@@ -136,7 +136,7 @@
    ],
    "source": [
     "# fits\n",
-    "lg = logit_utils_gen.LogisticPolyRegression(3, mono = True, lam = (0, 1e-3))\n",
+    "lg = logit_gen.LogisticPolyRegression(3, mono = True, lam = (0, 1e-3))\n",
     "ress = [lg.fit(x, y, method=\"diff_evol\") for x in xs]\n",
     "\n",
     "pd.DataFrame(ress, index = pd.Index(scales, name = 'scale'))"

python/logistic/logit_utils_gen.py → python/logistic/logit_gen.py

@@ -63,6 +63,8 @@ import scipy.optimize
 
 import mono_cubic2 as mc
 
+
+
 def resize_with_const(v, n, val=0):
     """
     Resize a 1D vector to a specified length `n`.
@@ -87,25 +89,6 @@ def resize_with_const(v, n, val=0):
     
     return np.concatenate([v, np.full(n - len(v), val)])
 
-def safe_exp(x, max_exp = 700):
-    """
-    A numerically robust version of np.exp that avoids overflow by clipping the input.
-    
-    Parameters:
-        x : array_like
-            Input value or array.
-        max_exp : float
-            Maximum allowed exponent value. np.exp(709) ≈ 8.2e307 (close to float64 max).
-    
-    Returns:
-        array_like
-            The exponential of the input with overflow protection.
-    """
-
-    return np.exp(np.clip(x, -max_exp, max_exp))
-
-def safe_expit(x, max_exp = 700): return 1/(1 + safe_exp(-x, max_exp))
-
 """
     Fitting data
 
@@ -241,7 +224,7 @@ class LogisticPolyRegression:
         X = np.column_stack([x**i for i in range(len(beta))])
         F = X @ beta              # decision function, X beta
 
-        return safe_expit(F)
+        return scipy.special.expit(F)
 
     """
         Calculate negative log-likelihood function
@@ -265,26 +248,24 @@ class LogisticPolyRegression:
             nllf           : if jac is false
             (nllf, grad)   : if jac is true
     """
-    def get_nllf(self, x, y, pars, jac = False):
+    def get_nllf(self, x, y, pars, jac=False):
+        beta = self.get_beta(pars)              # shape (p,)
+        J = self.get_jac_beta(pars) if jac else None   # shape (p, q)
 
-        beta = self.get_beta(pars)
-        X = np.column_stack([x**i for i in range(len(beta))])
-    
-        # signs
-        s = 2.0*y - 1
-        
-        # decision function for conditional probability Prob(Y = y| x)
-        F = s*(X @ beta)
+        X = np.column_stack([x**i for i in range(len(beta))])   # shape (n, p)
+        s = 2.0 * y - 1.0
 
-        nllf = np.sum(np.log(1 + safe_exp(-F))) 
+        eta = X @ beta
+        z = s * eta
 
-        if not jac: return nllf
-        
-        J = self.get_jac_beta(pars)
+        # stable negative log-likelihood
+        nllf = -np.sum(scipy.special.log_expit(z))
 
-        grad = -(s*safe_expit(-F)) @ (X @ J)
+        if not jac:
+            return nllf
 
-        return (nllf, grad)
+        grad = -(s * (1.0 - scipy.special.expit(z))) @ (X @ J)
+        return nllf, grad
 
     """
         Penalty function
@@ -485,7 +466,7 @@ class LogisticPolyRegression:
         F = V @ beta
 
         # probabilities p_i = P(Y=y_i | x_i)
-        p = safe_expit(F)
+        p = scipy.special.expit(F)
         q = 1 - p
 
         # calculate hessian
@@ -611,7 +592,7 @@ class LogisticPolyRegression:
             # quantiles of decision function
             Q = np.quantile(X@beta.T, probs, axis = 1)
             
-            return np.apply_along_axis(safe_expit, 1, Q)
+            return np.apply_along_axis(scipy.special.expit, 1, Q)
 
         # J = d(beta)/d(pars)
         J = self.get_jac_beta(mean_pars)
@@ -627,7 +608,7 @@ class LogisticPolyRegression:
         Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)
         
         # convert logit to expit
-        return safe_expit(Q)
+        return scipy.special.expit(Q)
 
     """
         Calculating quantiles using delta method of the model values 
@@ -672,8 +653,10 @@ class LogisticPolyRegression:
         S = X@J
 
         # attributes of normal distribution of model values
-        locs = safe_expit(F)
-        scales = np.sqrt(np.diag(S@cov_pars@S.T))/(4*np.cosh(F/2)**2)
+        locs = scipy.special.expit(F)
+
+        derivative = locs * (1 - locs)
+        scales = np.sqrt(np.diag(S @ cov_pars @ S.T)) * derivative
 
         # computing quantiles of logit
         Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)

python/logistic/report_utils.py

@@ -0,0 +1,619 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import os
+
+
+def build_plot_data(df_data, df_fit_index, alpha, lg, n_grid=100, verbose=False):
+    """
+    Compute all quantities needed for plotting and return them as dataframes.
+
+    Parameters
+    ----------
+    df_data : pd.DataFrame
+        Input data with columns including scale, dataset, X, Y.
+    df_fit_index : pd.DataFrame
+        Fit results indexed by (scale, dataset), with columns 'pars' and 'cov'.
+    alpha : float
+        Two-sided significance level for confidence intervals.
+        Example: alpha = 0.05 gives a 95% CI.
+    lg : object
+        Model object with methods:
+            - model(x, pars)
+            - get_model_quantiles_normal(x, probs, pars, cov)
+            - get_model_quantiles_delta(x, probs, pars, cov)
+    n_grid : int, default 100
+        Number of x-grid points for plotting.
+    verbose : bool, default False
+        Print progress information.
+
+    Returns
+    -------
+    df_points : pd.DataFrame
+        Original observed data, one row per observation.
+    df_fit_plot : pd.DataFrame
+        Fitted curve on a grid, one row per x-grid value.
+    df_ci_plot : pd.DataFrame
+        Confidence interval bands on a grid, one row per x-grid value and method.
+    """
+    probs = [alpha / 2, 1 - alpha / 2]
+    confidence = 100 * (1 - alpha)
+
+    points_list = []
+    fit_list = []
+    ci_list = []
+
+    for (scale, dataset), g in df_data.groupby(["scale", "dataset"]):
+        if verbose:
+            print(f"Calculating: {scale}, {dataset}")
+
+        X = g["X"].to_numpy()
+        Y = g["Y"].to_numpy()
+
+        # store raw points
+        tmp_points = g[["X", "Y"]].copy()
+        tmp_points["scale"] = scale
+        tmp_points["dataset"] = dataset
+        tmp_points["class_label"] = tmp_points["Y"].map({0: "NC", 1: "AE"})
+        points_list.append(tmp_points[["scale", "dataset", "X", "Y", "class_label"]])
+
+        # get fit objects
+        fit_pars = df_fit_index.loc[(scale, dataset), "pars"]
+        fit_cov = df_fit_index.loc[(scale, dataset), "cov"]
+
+        # grid for smooth curve / bands
+        x_fit = np.linspace(X.min(), X.max(), n_grid)
+
+        # fitted curve
+        y_fit = lg.model(x_fit, fit_pars)
+        fit_list.append(
+            pd.DataFrame({
+                "scale": scale,
+                "dataset": dataset,
+                "X": x_fit,
+                "y_fit": y_fit,
+            })
+        )
+
+        # confidence bands
+        for method in ["normal", "delta"]:
+            qfun = getattr(lg, f"get_model_quantiles_{method}")
+            y_low, y_high = qfun(x_fit, probs, fit_pars, fit_cov)
+
+            ci_list.append(
+                pd.DataFrame({
+                    "scale": scale,
+                    "dataset": dataset,
+                    "method": method,
+                    "confidence": confidence,
+                    "X": x_fit,
+                    "y_low": y_low,
+                    "y_high": y_high,
+                })
+            )
+
+    df_points = pd.concat(points_list, ignore_index=True)
+    df_fit_plot = pd.concat(fit_list, ignore_index=True)
+    df_ci_plot = pd.concat(ci_list, ignore_index=True)
+
+    return df_points, df_fit_plot, df_ci_plot
+
+
+def plot_from_dataframes(
+    df_points,
+    df_fit_plot,
+    df_ci_plot,
+    xlabs,
+    scales=None,
+    datasets=None,
+    results_path=None,
+    filename="logit_fit_CI_simple_paper.pdf",
+    verbose=False,
+    color_map={
+        "data_NC": "blue",
+        "data_AE": "orange",
+        "fit": "black",
+        "normal": "red",
+        "delta": "green",
+    },
+    ls_map={
+        "fit": "-",
+        "normal": "--",
+        "delta": "-.",
+    },
+    ci_alpha=0.15,
+    ci_linewidth=1.5,
+    fit_linewidth=2.0,
+    width=5,
+    height=4,
+):
+    
+
+    if scales is None:
+        scales = sorted(df_points["scale"].unique())
+    if datasets is None:
+        datasets = sorted(df_points["dataset"].unique())
+
+    fig, axs = plt.subplots(
+        ncols=len(scales),
+        nrows=len(datasets),
+        figsize=(width * len(scales), height * len(datasets)),
+        layout="constrained",
+        squeeze=False,
+    )
+
+    plot_labs = dict(
+        zip(
+            [(scale, dataset) for dataset in datasets for scale in scales],
+            ["A", "B", "C", "D"][: len(scales) * len(datasets)],
+        )
+    )
+    if verbose:
+        print(plot_labs)
+
+    for dataset_i, dataset in enumerate(datasets):
+        for scale_i, scale in enumerate(scales):
+            ax = axs[dataset_i, scale_i]
+
+            ax.set_xlabel(xlabs[scale])
+            ax.set_ylabel("P(AE|X = x)")
+            ax.text(
+                0.05,
+                0.9,
+                plot_labs[(scale, dataset)],
+                fontsize=14,
+                transform=ax.transAxes,
+            )
+
+            # raw points
+            g_points = df_points[
+                (df_points["scale"] == scale) & (df_points["dataset"] == dataset)
+            ]
+
+            for lab in ["NC", "AE"]:
+                sub = g_points[g_points["class_label"] == lab]
+                ax.scatter(
+                    sub["X"],
+                    sub["Y"],
+                    label=f"data: {lab}",
+                    alpha=0.5,
+                    color=color_map[f"data_{lab}"],
+                )
+
+            # fitted curve
+            g_fit = df_fit_plot[
+                (df_fit_plot["scale"] == scale) & (df_fit_plot["dataset"] == dataset)
+            ].sort_values("X")
+
+            ax.plot(
+                g_fit["X"],
+                g_fit["y_fit"],
+                label="fit",
+                color=color_map["fit"],
+                linestyle=ls_map["fit"],
+                linewidth=fit_linewidth,
+            )
+
+            # CI bands + border lines
+            g_ci = df_ci_plot[
+                (df_ci_plot["scale"] == scale) & (df_ci_plot["dataset"] == dataset)
+            ]
+
+            for method, sub in g_ci.groupby("method"):
+                sub = sub.sort_values("X")
+
+                color = color_map[method]
+                ls = ls_map[method]
+                confidence = sub["confidence"].iloc[0]
+
+                # format nicely: 95 instead of 95.0 when possible
+                conf_str = f"{confidence:.0f}" if float(confidence).is_integer() else f"{confidence:.1f}"
+
+                # light filled band
+                ax.fill_between(
+                    sub["X"],
+                    sub["y_low"],
+                    sub["y_high"],
+                    color=color,
+                    alpha=ci_alpha,
+                    linewidth=0,
+                )
+
+                # lower border
+                ax.plot(
+                    sub["X"],
+                    sub["y_low"],
+                    color=color,
+                    linestyle=ls,
+                    linewidth=ci_linewidth,
+                    label=f"CI: {method} {conf_str}%",
+                )
+
+                # upper border
+                ax.plot(
+                    sub["X"],
+                    sub["y_high"],
+                    color=color,
+                    linestyle=ls,
+                    linewidth=ci_linewidth,
+                )
+
+    for ax in axs.flat:
+        ax.label_outer()
+
+    handles, labels = axs[0, 0].get_legend_handles_labels()
+    by_label = dict(zip(labels, handles))
+    fig.legend(
+        by_label.values(),
+        by_label.keys(),
+        bbox_to_anchor=(0.975, 0.2),
+        framealpha=0.5,
+        loc="center right",
+    )
+
+    if results_path is not None:
+        plt.savefig(
+            os.path.join(results_path, filename),
+            bbox_inches="tight",
+        )
+
+    plt.show()
+
+
+import os
+import pickle
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+
+def load_or_build_bootstrap_pars(
+    df_data,
+    lg,
+    results_path,
+    filename="boots_pars_results.pkl",
+    n_boots=10000,
+    methods=("normal", "nonparam_boots", "nonparam_stratified_boots", "parametric_boots"),
+    verbose=False,
+):
+    """
+    Load bootstrap parameter results from disk if available, otherwise compute and save.
+
+    Returns
+    -------
+    boots_pars_results : dict
+        Dictionary with keys (scale, dataset, method) and values bootstrap parameter arrays.
+    """
+    result_file = os.path.join(results_path, filename)
+
+    if os.path.isfile(result_file):
+        with open(result_file, "rb") as f:
+            boots_pars_results = pickle.load(f)
+        return boots_pars_results
+
+    boots_pars_results = {}
+
+    for (scale, dataset), g in df_data.groupby(["scale", "dataset"]):
+        if verbose:
+            print(f"Bootstrap parameters: {scale}, {dataset}")
+
+        X = g["X"].to_numpy()
+        Y = g["Y"].to_numpy()
+
+        for method in methods:
+            if verbose:
+                print(f"\tmethod: {method}")
+
+            par_fun = getattr(lg, f"get_{method}_pars")
+            boots_pars_results[(scale, dataset, method)] = par_fun(X, Y, m=n_boots)
+
+    with open(result_file, "wb") as f:
+        pickle.dump(boots_pars_results, f)
+
+    return boots_pars_results
+
+
+def build_plot_data_full_ci(
+    df_data,
+    df_fit_index,
+    lg,
+    alpha,
+    boots_pars_results,
+    n_grid=100,
+    analytic_methods=("normal", "delta"),
+    bootstrap_methods=("nonparam_boots", "parametric_boots"),
+    verbose=False,
+):
+    """
+    Compute all quantities needed for plotting and return them as dataframes.
+
+    Parameters
+    ----------
+    df_data : pd.DataFrame
+        Input data with columns including scale, dataset, X, Y.
+    df_fit_index : pd.DataFrame
+        Fit results indexed by (scale, dataset), with columns 'pars' and 'cov'.
+    lg : object
+        Model object with methods:
+            - model(x, pars)
+            - get_model_quantiles_<method>(x, probs, pars, cov)
+    alpha : float
+        Two-sided significance level. Example alpha=0.05 gives 95% CI.
+    boots_pars_results : dict
+        Dictionary with keys (scale, dataset, method) and bootstrap parameter samples as values.
+    n_grid : int
+        Number of x-grid points.
+    analytic_methods : tuple
+        CI methods based on fitted covariance.
+    bootstrap_methods : tuple
+        CI methods based on bootstrap parameter samples.
+    verbose : bool
+        Print progress.
+
+    Returns
+    -------
+    df_points : pd.DataFrame
+    df_fit_plot : pd.DataFrame
+    df_ci_plot : pd.DataFrame
+        Long dataframe with both analytic and bootstrap confidence intervals.
+    """
+    probs = [alpha / 2, 1 - alpha / 2]
+    confidence = 100 * (1 - alpha)
+
+    points_list = []
+    fit_list = []
+    ci_list = []
+
+    for (scale, dataset), g in df_data.groupby(["scale", "dataset"]):
+        if verbose:
+            print(f"Calculating plot data: {scale}, {dataset}")
+
+        X = g["X"].to_numpy()
+        Y = g["Y"].to_numpy()
+
+        # raw points
+        tmp_points = g[["X", "Y"]].copy()
+        tmp_points["scale"] = scale
+        tmp_points["dataset"] = dataset
+        tmp_points["class_label"] = tmp_points["Y"].map({0: "NC", 1: "AE"})
+        points_list.append(tmp_points[["scale", "dataset", "X", "Y", "class_label"]])
+
+        # fitted objects
+        fit_pars = df_fit_index.loc[(scale, dataset), "pars"]
+        fit_cov = df_fit_index.loc[(scale, dataset), "cov"]
+
+        # x-grid
+        x_fit = np.linspace(X.min(), X.max(), n_grid)
+
+        # fitted curve
+        y_fit = lg.model(x_fit, fit_pars)
+        fit_list.append(
+            pd.DataFrame({
+                "scale": scale,
+                "dataset": dataset,
+                "X": x_fit,
+                "y_fit": y_fit,
+            })
+        )
+
+        # analytic CIs
+        for method in analytic_methods:
+            if verbose:
+                print(f"\tanalytic CI: {method}")
+
+            qfun = getattr(lg, f"get_model_quantiles_{method}")
+            y_low, y_high = qfun(x_fit, probs, fit_pars, fit_cov)
+
+            ci_list.append(
+                pd.DataFrame({
+                    "scale": scale,
+                    "dataset": dataset,
+                    "method": method,
+                    "ci_source": "analytic",
+                    "confidence": confidence,
+                    "X": x_fit,
+                    "y_low": y_low,
+                    "y_high": y_high,
+                })
+            )
+
+        # bootstrap CIs
+        for method in bootstrap_methods:
+            if verbose:
+                print(f"\tbootstrap CI: {method}")
+
+            bpars = boots_pars_results[(scale, dataset, method)]
+            quant = np.quantile([lg.model(x_fit, p) for p in bpars], probs, axis=0)
+            y_low, y_high = quant
+
+            ci_list.append(
+                pd.DataFrame({
+                    "scale": scale,
+                    "dataset": dataset,
+                    "method": method,
+                    "ci_source": "bootstrap",
+                    "confidence": confidence,
+                    "X": x_fit,
+                    "y_low": y_low,
+                    "y_high": y_high,
+                })
+            )
+
+    df_points = pd.concat(points_list, ignore_index=True)
+    df_fit_plot = pd.concat(fit_list, ignore_index=True)
+    df_ci_plot = pd.concat(ci_list, ignore_index=True)
+
+    return df_points, df_fit_plot, df_ci_plot
+
+
+def plot_from_dataframes_full_ci(
+    df_points,
+    df_fit_plot,
+    df_ci_plot,
+    xlabs,
+    scales=None,
+    datasets=None,
+    results_path=None,
+    filename="logit_fit_CI_paper.pdf",
+    verbose=False,
+    color_map=None,
+    ls_map=None,
+    ci_alpha=0.12,
+    ci_linewidth=1.5,
+    fit_linewidth=2.0,
+    width=5,
+    height=4,
+):
+    """
+    Plot data, fitted curve, and confidence intervals from prepared dataframes.
+    """
+    if scales is None:
+        scales = list(df_points["scale"].drop_duplicates())
+    if datasets is None:
+        datasets = list(df_points["dataset"].drop_duplicates())
+
+    if color_map is None:
+        color_map = {
+            "data_NC": "blue",
+            "data_AE": "orange",
+            "fit": "black",
+            "normal": "red",
+            "delta": "green",
+            "nonparam_boots": "blue",
+            "nonparam_stratified_boots": "purple",
+            "parametric_boots": "cyan",
+        }
+
+    if ls_map is None:
+        ls_map = {
+            "fit": "-",
+            "normal": "--",
+            "delta": "-.",
+            "nonparam_boots": ":",
+            "nonparam_stratified_boots": (0, (3, 1, 1, 1)),
+            "parametric_boots": (0, (5, 2)),
+        }
+
+    fig, axs = plt.subplots(
+        ncols=len(scales),
+        nrows=len(datasets),
+        figsize=(width * len(scales), height * len(datasets)),
+        layout="constrained",
+        squeeze=False,
+    )
+
+    plot_labs = dict(
+        zip(
+            [(scale, dataset) for dataset in datasets for scale in scales],
+            list("ABCDEFGHIJKLMNOPQRSTUVWXYZ")[: len(scales) * len(datasets)],
+        )
+    )
+    if verbose:
+        print(plot_labs)
+
+    for dataset_i, dataset in enumerate(datasets):
+        for scale_i, scale in enumerate(scales):
+            ax = axs[dataset_i, scale_i]
+
+            ax.set_xlabel(xlabs[scale])
+            ax.set_ylabel("P(AE|X = x)")
+            ax.text(
+                0.05,
+                0.9,
+                plot_labs[(scale, dataset)],
+                fontsize=14,
+                transform=ax.transAxes,
+            )
+
+            # raw points
+            g_points = df_points[
+                (df_points["scale"] == scale) & (df_points["dataset"] == dataset)
+            ]
+
+            for lab in ["NC", "AE"]:
+                sub = g_points[g_points["class_label"] == lab]
+                ax.scatter(
+                    sub["X"],
+                    sub["Y"],
+                    label=f"data: {lab}",
+                    alpha=0.5,
+                    color=color_map[f"data_{lab}"],
+                )
+
+            # fitted curve
+            g_fit = df_fit_plot[
+                (df_fit_plot["scale"] == scale) & (df_fit_plot["dataset"] == dataset)
+            ].sort_values("X")
+
+            ax.plot(
+                g_fit["X"],
+                g_fit["y_fit"],
+                label="fit",
+                color=color_map["fit"],
+                linestyle=ls_map["fit"],
+                linewidth=fit_linewidth,
+            )
+
+            # CI bands + border lines
+            g_ci = df_ci_plot[
+                (df_ci_plot["scale"] == scale) & (df_ci_plot["dataset"] == dataset)
+            ]
+
+            for method, sub in g_ci.groupby("method", sort=False):
+                sub = sub.sort_values("X")
+
+                color = color_map[method]
+                ls = ls_map[method]
+                confidence = sub["confidence"].iloc[0]
+                conf_str = (
+                    f"{confidence:.0f}"
+                    if float(confidence).is_integer()
+                    else f"{confidence:.1f}"
+                )
+
+                ax.fill_between(
+                    sub["X"],
+                    sub["y_low"],
+                    sub["y_high"],
+                    color=color,
+                    alpha=ci_alpha,
+                    linewidth=0,
+                )
+
+                ax.plot(
+                    sub["X"],
+                    sub["y_low"],
+                    color=color,
+                    linestyle=ls,
+                    linewidth=ci_linewidth,
+                    label=f"CI: {method} {conf_str}%",
+                )
+
+                ax.plot(
+                    sub["X"],
+                    sub["y_high"],
+                    color=color,
+                    linestyle=ls,
+                    linewidth=ci_linewidth,
+                )
+
+    for ax in axs.flat:
+        ax.label_outer()
+
+    handles, labels = axs[0, 0].get_legend_handles_labels()
+    by_label = dict(zip(labels, handles))
+    fig.legend(
+        by_label.values(),
+        by_label.keys(),
+        bbox_to_anchor=(0.975, 0.2),
+        framealpha=0.5,
+        loc="center right",
+    )
+
+    if results_path is not None:
+        plt.savefig(
+            os.path.join(results_path, filename),
+            bbox_inches="tight",
+        )
+
+    plt.show()