{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "850e7703",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [intercept, slope]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "570fd1615d6640639a2795169dea6d2a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">\n",
       "</pre>\n"
      ],
      "text/plain": [
       "\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 162 seconds.\n",
      "There were 15 divergences after tuning. Increase `target_accept` or reparameterize.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Frequency')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 19452 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Visualization of the distribution of X (Max SUV at P94)\n",
    "\n",
    "\n",
    "plt.hist(X, bins=20)\n",
    "plt.title(\"Distribution of X (Max SUV at P94)\")\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "757919c7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [intercept, slope]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "87bf72eb6aa542f08070386bbedf1a28",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">\n",
       "</pre>\n"
      ],
      "text/plain": [
       "\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 200 seconds.\n",
      "There were 8 divergences after tuning. Increase `target_accept` or reparameterize.\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 72528 characters]",
      "text/plain": [
       "<Figure size 800x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pymc as pm\n",
    "import arviz as az\n",
    "from scipy.io import loadmat\n",
    "\n",
    "\n",
    "# Load .mat files\n",
    "suv_file = \"../data/suv_percentilesSLOthenUWM.mat\"\n",
    "flags_file = \"../data/flags_combined.mat\"\n",
    "\n",
    "suv_data = loadmat(suv_file)\n",
    "flags_data = loadmat(flags_file)\n",
    "\n",
    "# Extract data\n",
    "suv = suv_data['lung_SUVperc_COMBINED'][0:58, :, :]  # SUV shape: (58, patients, percentiles)\n",
    "flags = flags_data['flags'][0:58, 3]  # Binary adverse effect (0 or 1)\n",
    "\n",
    "# Choose SUV percentile (94th) and reshape\n",
    "p = 94  # Use 94th percentile (0-indexed)\n",
    "X = np.nanmax(suv[:, :, p], axis=1).reshape(-1, 1)  # Shape (58, 1)\n",
    "y = flags.flatten().astype(int)                    # Shape (58,)\n",
    "\n",
    "X_mean = X.mean()\n",
    "X_std = X.std()\n",
    "X_stdized = (X - X_mean) / X_std\n",
    "\n",
    "\n",
    "\n",
    "with pm.Model() as logistic_model:\n",
    "    # Priors on coefficients (robust heavy-tailed)\n",
    "    intercept = pm.StudentT(\"intercept\", mu=0, sigma=5, nu=3)\n",
    "    slope = pm.StudentT(\"slope\", mu=0, sigma=5, nu=3)\n",
    "\n",
    "\n",
    "    logits = intercept + slope * X_stdized.flatten()\n",
    "\n",
    "\n",
    "    theta = pm.Deterministic(\"theta\", pm.math.sigmoid(logits))\n",
    "\n",
    "\n",
    "    y_obs = pm.Bernoulli(\"y_obs\", p=theta, observed=y)\n",
    "\n",
    "\n",
    "    trace = pm.sample(2000, tune=1000, target_accept=0.95, return_inferencedata=True)\n",
    "\n",
    "\n",
    "\n",
    "# Grid of x values for prediction\n",
    "x_grid = np.linspace(X_stdized.min(), X_stdized.max(), 200)\n",
    "\n",
    "# Extract posterior samples\n",
    "intercept_samples = trace.posterior['intercept'].stack(samples=(\"chain\", \"draw\")).values\n",
    "slope_samples = trace.posterior['slope'].stack(samples=(\"chain\", \"draw\")).values\n",
    "\n",
    "# Compute predicted probabilities\n",
    "logits_samples = intercept_samples[:, None] + slope_samples[:, None] * x_grid[None, :]\n",
    "probs_samples = 1 / (1 + np.exp(-logits_samples))\n",
    "\n",
    "# Mean and 95% credible interval\n",
    "mean_probs = np.mean(probs_samples, axis=0)\n",
    "lower = np.percentile(probs_samples, 2.5, axis=0)\n",
    "upper = np.percentile(probs_samples, 97.5, axis=0)\n",
    "\n",
    "# De-standardize x-axis\n",
    "x_plot = x_grid * X_std + X_mean\n",
    "\n",
    "# Plot\n",
    "plt.figure(figsize=(8, 5))\n",
    "plt.plot(x_plot, mean_probs, label='Mean prediction (sigmoid)', color='blue')\n",
    "plt.fill_between(x_plot, lower, upper, color='blue', alpha=0.3, label='95% CI')\n",
    "plt.scatter(X, y, color='red', alpha=0.6, label='Observed data')\n",
    "plt.xlabel(\"Max SUV at 94th Percentile\")\n",
    "plt.ylabel(\"Probability of Adverse Effect\")\n",
    "plt.title(\"Bayesian Logistic Regression\\n(Sigmoid Curve with 95% Credible Interval)\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "112beb58",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 36204 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.hist(X_NC, bins=20, density=True, alpha=0.5, label=\"NC histogram\")\n",
    "plt.hist(X_AE, bins=20, density=True, alpha=0.5, label=\"AE histogram\")\n",
    "plt.plot(xs, f_XatNC_opt(xs, popt_NC), label=\"NC fit\", color='blue')\n",
    "plt.plot(xs, f_XatAE_opt(xs, popt_AE), label=\"AE fit\", color='red')\n",
    "plt.legend()\n",
    "plt.title(\"Fit of Custom Distribution to NC and AE Groups\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1f40976",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NC params: [ 0.24241006  0.18097191 -0.16096395]\n",
      "AE params: [ 0.39013584  0.04654894 -2.02453446]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 54912 characters]",
      "text/plain": [
       "<Figure size 800x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import loadmat\n",
    "from scipy.optimize import minimize\n",
    "import os\n",
    "\n",
    "\n",
    "suv_file = \"../data/suv_percentilesSLOthenUWM.mat\"\n",
    "flags_file = \"../data/flags_combined.mat\"\n",
    "\n",
    "suv_data = loadmat(suv_file)\n",
    "flags_data = loadmat(flags_file)\n",
    "\n",
    "suv = suv_data['lung_SUVperc_COMBINED'][0:58,:,:]\n",
    "flags = flags_data['flags'][0:58, 3]\n",
    "\n",
    "\n",
    "p = 94  # percentile index\n",
    "X = np.nanmax(suv[:,:,p], axis=1).reshape(-1, 1)\n",
    "flags = flags.flatten()\n",
    " \n",
    "\n",
    "g = lambda x, mu, sigma, l: np.abs(x)**(l-1) * np.exp(-(((np.abs(x)**l - 1)/l)-mu)**2/(2*sigma**2)) / (np.sqrt(2*np.pi)*sigma)\n",
    "fun = lambda x, mu, sigma, l: (l-1)*np.log(np.abs(x)) - np.log(sigma*np.sqrt(2*np.pi)) - (((np.abs(x)**l - 1)/l)-mu)**2/(2*sigma**2)\n",
    "\n",
    "def safe_divide(numerator, denominator, eps=1e-10):\n",
    "    return numerator / np.clip(denominator, eps, np.inf)\n",
    "\n",
    "def myfitter(data, fun, cov=True):\n",
    "    nllf = lambda p: -fun(data, *p).sum()\n",
    "    p0 = (1, 0.6, 0.3)\n",
    "    bounds = ((0, 10), (0.01, 10), (-5, 5))\n",
    "    res = minimize(nllf, p0, bounds=bounds, method=\"L-BFGS-B\")\n",
    "    if not res.success:\n",
    "        raise RuntimeError(\"Fitting failed\")\n",
    "    if not cov:\n",
    "        return res.x, None\n",
    "    return res.x, res.hess_inv.todense()\n",
    "\n",
    "\n",
    "X_NC = X[flags == 0].flatten()\n",
    "X_AE = X[flags == 1].flatten()\n",
    "flags_NC = flags[flags == 0]\n",
    "flags_AE = flags[flags == 1]\n",
    "\n",
    "\n",
    "popt_NC, _ = myfitter(X_NC, fun)\n",
    "popt_AE, _ = myfitter(X_AE, fun)\n",
    "\n",
    "print(\"NC params:\", popt_NC)\n",
    "print(\"AE params:\", popt_AE)\n",
    "\n",
    "P_NC = np.count_nonzero(flags == 0) / len(flags)\n",
    "P_AE = np.count_nonzero(flags == 1) / len(flags)\n",
    "\n",
    "f_XatNC_opt = lambda x: g(x, *popt_NC)\n",
    "f_XatAE_opt = lambda x: g(x, *popt_AE)\n",
    "f_X_opt = lambda x: P_NC * f_XatNC_opt(x) + P_AE * f_XatAE_opt(x)\n",
    "P_AEatX_opt = lambda x: safe_divide(P_AE * f_XatAE_opt(x), f_X_opt(x))\n",
    "\n",
    "xs = np.linspace(1e-3, X.max()+1, 200)\n",
    "yopts = P_AEatX_opt(xs)\n",
    "\n",
    "\n",
    "n_boot = 500\n",
    "ysamples = []\n",
    "\n",
    "for i in range(n_boot):\n",
    "    try:\n",
    "        X_NC_boot = np.random.choice(X_NC, size=len(X_NC), replace=True)\n",
    "        popt_NC_boot, _ = myfitter(X_NC_boot, fun, cov=False)\n",
    "\n",
    "        X_AE_boot = np.random.choice(X_AE, size=len(X_AE), replace=True)\n",
    "        popt_AE_boot, _ = myfitter(X_AE_boot, fun, cov=False)\n",
    "\n",
    "        f_XatNC_boot = lambda x: g(x, *popt_NC_boot)\n",
    "        f_XatAE_boot = lambda x: g(x, *popt_AE_boot)\n",
    "        f_X_boot = lambda x: P_NC * f_XatNC_boot(x) + P_AE * f_XatAE_boot(x)\n",
    "        P_AEatX_boot = lambda x: safe_divide(P_AE * f_XatAE_boot(x), f_X_boot(x))\n",
    "\n",
    "        ysamples.append(P_AEatX_boot(xs))\n",
    "    except Exception as e:\n",
    "        continue  # skip failed fits\n",
    "\n",
    "ysamples = np.array(ysamples)\n",
    "\n",
    "if ysamples.shape[0] < 10:\n",
    "    print(\"Too many failed fits. Try reducing complexity or increasing data.\")\n",
    "    lower = upper = np.full_like(xs, np.nan)\n",
    "else:\n",
    "    lower = np.percentile(ysamples, 2.5, axis=0)\n",
    "    upper = np.percentile(ysamples, 97.5, axis=0)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.scatter(X_NC, flags_NC, label=\"NC\", alpha=0.5)\n",
    "plt.scatter(X_AE, flags_AE, label=\"AE\", alpha=0.5)\n",
    "plt.plot(xs, yopts, label=\"Fit model\", color=\"red\")\n",
    "if not np.all(np.isnan(lower)):\n",
    "    plt.fill_between(xs, lower, upper, color='red', alpha=0.3, label=\"95% CI\")\n",
    "plt.xlabel(\"Max SUV at 94th Percentile\")\n",
    "plt.ylabel(\"Probability of Adverse Effect\")\n",
    "plt.title(\"Bayesian-like model with Bootstrap Confidence Interval\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "cbf195fe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting pymc3\n",
      "  Using cached pymc3-3.11.6-py3-none-any.whl.metadata (15 kB)\n",
      "Requirement already satisfied: arviz>=0.11.0 in c:\\users\\zahra\\anaconda3\\lib\\site-packages (from pymc3) (0.21.0)\n",
      "Requirement already satisfied: cachetools>=4.2.1 in c:\\users\\zahra\\anaconda3\\lib\\site-packages (from pymc3) (5.3.3)\n",
      "Collecting deprecat (from pymc3)\n",
      "  Using cached deprecat-2.1.3-py2.py3-none-any.whl.metadata (1.6 kB)\n",
      "Requirement already satisfied: dill in c:\\users\\zahra\\anaconda3\\lib\\site-packages (from pymc3) (0.3.8)\n",
      "Collecting fastprogress>=0.2.0 (from pymc3)\n",
      "  Using cached fastprogress-1.0.3-py3-none-any.whl.metadata (5.6 kB)\n",
      "Collecting numpy<1.22.2,>=1.15.0 (from pymc3)\n",
      "  Using cached numpy-1.22.1.zip (11.4 MB)\n",
      "  Installing build dependencies: started\n",
      "  Installing build dependencies: still running...\n",
      "  Installing build dependencies: finished with status 'done'\n",
      "  Getting requirements to build wheel: started\n",
      "  Getting requirements to build wheel: finished with status 'error'\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  error: subprocess-exited-with-error\n",
      "  \n",
      "  × Getting requirements to build wheel did not run successfully.\n",
      "  │ exit code: 1\n",
      "  ╰─> [33 lines of output]\n",
      "      Traceback (most recent call last):\n",
      "        File \"C:\\Users\\zahra\\anaconda3\\Lib\\site-packages\\pip\\_vendor\\pyproject_hooks\\_in_process\\_in_process.py\", line 353, in <module>\n",
      "          main()\n",
      "        File \"C:\\Users\\zahra\\anaconda3\\Lib\\site-packages\\pip\\_vendor\\pyproject_hooks\\_in_process\\_in_process.py\", line 335, in main\n",
      "          json_out['return_val'] = hook(**hook_input['kwargs'])\n",
      "                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n",
      "        File \"C:\\Users\\zahra\\anaconda3\\Lib\\site-packages\\pip\\_vendor\\pyproject_hooks\\_in_process\\_in_process.py\", line 112, in get_requires_for_build_wheel\n",
      "          backend = _build_backend()\n",
      "                    ^^^^^^^^^^^^^^^^\n",
      "        File \"C:\\Users\\zahra\\anaconda3\\Lib\\site-packages\\pip\\_vendor\\pyproject_hooks\\_in_process\\_in_process.py\", line 77, in _build_backend\n",
      "          obj = import_module(mod_path)\n",
      "                ^^^^^^^^^^^^^^^^^^^^^^^\n",
      "        File \"C:\\Users\\zahra\\anaconda3\\Lib\\importlib\\__init__.py\", line 90, in import_module\n",
      "          return _bootstrap._gcd_import(name[level:], package, level)\n",
      "                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n",
      "        File \"<frozen importlib._bootstrap>\", line 1387, in _gcd_import\n",
      "        File \"<frozen importlib._bootstrap>\", line 1360, in _find_and_load\n",
      "        File \"<frozen importlib._bootstrap>\", line 1310, in _find_and_load_unlocked\n",
      "        File \"<frozen importlib._bootstrap>\", line 488, in _call_with_frames_removed\n",
      "        File \"<frozen importlib._bootstrap>\", line 1387, in _gcd_import\n",
      "        File \"<frozen importlib._bootstrap>\", line 1360, in _find_and_load\n",
      "        File \"<frozen importlib._bootstrap>\", line 1331, in _find_and_load_unlocked\n",
      "        File \"<frozen importlib._bootstrap>\", line 935, in _load_unlocked\n",
      "        File \"<frozen importlib._bootstrap_external>\", line 995, in exec_module\n",
      "        File \"<frozen importlib._bootstrap>\", line 488, in _call_with_frames_removed\n",
      "        File \"C:\\Users\\zahra\\AppData\\Local\\Temp\\pip-build-env-pniegcuk\\overlay\\Lib\\site-packages\\setuptools\\__init__.py\", line 16, in <module>\n",
      "          import setuptools.version\n",
      "        File \"C:\\Users\\zahra\\AppData\\Local\\Temp\\pip-build-env-pniegcuk\\overlay\\Lib\\site-packages\\setuptools\\version.py\", line 1, in <module>\n",
      "          import pkg_resources\n",
      "        File \"C:\\Users\\zahra\\AppData\\Local\\Temp\\pip-build-env-pniegcuk\\overlay\\Lib\\site-packages\\pkg_resources\\__init__.py\", line 2172, in <module>\n",
      "          register_finder(pkgutil.ImpImporter, find_on_path)\n",
      "                          ^^^^^^^^^^^^^^^^^^^\n",
      "      AttributeError: module 'pkgutil' has no attribute 'ImpImporter'. Did you mean: 'zipimporter'?\n",
      "      [end of output]\n",
      "  \n",
      "  note: This error originates from a subprocess, and is likely not a problem with pip.\n",
      "error: subprocess-exited-with-error\n",
      "\n",
      "× Getting requirements to build wheel did not run successfully.\n",
      "│ exit code: 1\n",
      "╰─> See above for output.\n",
      "\n",
      "note: This error originates from a subprocess, and is likely not a problem with pip.\n"
     ]
    }
   ],
   "source": [
    "!pip install pymc3\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79080223",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [5.28692634 2.06426605 1.        ]\n",
      "popt_AE: [5.45406231 1.27809753 1.        ]\n",
      "pcov_NC: [[0.1 0.  0. ]\n",
      " [0.  0.1 0. ]\n",
      " [0.  0.  0.1]]\n",
      "pcov_AE: [[0.1 0.  0. ]\n",
      " [0.  0.1 0. ]\n",
      " [0.  0.  0.1]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 34736 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 32420 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 57000 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "figsdir = \"./figs/\"\n",
    "\n",
    "\n",
    "def myfitter(X, fun, return_cov=True):\n",
    "    \n",
    "    p = np.array([np.mean(X), np.std(X), 1.0])\n",
    "    pcov = np.eye(3) * 0.1  #covariance matrix\n",
    "    return p, pcov\n",
    "\n",
    "def fun(x, mu, sigma, l):\n",
    "\n",
    "    return np.exp(-((x - mu)**2) / (2*sigma**2))\n",
    "\n",
    "def g(x, mu, sigma, l):\n",
    "    return fun(x, mu, sigma, l)\n",
    "\n",
    "def load_example_data():\n",
    "    suv_dict = {\n",
    "        'lung_SUVperc_COMBINED': np.abs(np.random.randn(58, 10, 3)) * 3  # some positive random SUV values\n",
    "    }\n",
    "    flags_dict = {\n",
    "        'flags': np.random.choice([0, 1], size=(58, 6))  # binary flags for 6 organs\n",
    "    }\n",
    "    return suv_dict, flags_dict\n",
    "\n",
    "def gen_plots(pert, prefix, alpha, title, suv_dict, flags_dict, bootstrapping=False):\n",
    "\n",
    "    organ = \"lung\"   # fixed organ\n",
    "\n",
    "    suv = suv_dict[organ + '_SUVperc_COMBINED'][0:58,:,:]\n",
    "    idx = 3  # lung index\n",
    "    \n",
    "    flags = flags_dict['flags'][0:58, idx]\n",
    "\n",
    "    # choose percentile of interest (for lung and bowel 95, for thyroid 75),\n",
    "    X = np.nanmax(suv[:, :, pert - 1], axis=1).reshape(-1,1)\n",
    "\n",
    "    # extracting X associated to NC and AE\n",
    "    X_NC = X[flags == 0]\n",
    "    X_AE = X[flags == 1]\n",
    "    # and flags associated to NC and AE\n",
    "    flags_NC = flags[flags==0]\n",
    "    flags_AE = flags[flags==1]\n",
    "\n",
    "    popt_NC, pcov_NC = myfitter(X_NC, fun)\n",
    "    popt_AE, pcov_AE = myfitter(X_AE, fun)\n",
    "\n",
    "    # print optimization results\n",
    "    print(\"popt_NC:\", popt_NC)\n",
    "    print(\"popt_AE:\", popt_AE)\n",
    "\n",
    "    print(\"pcov_NC:\", pcov_NC)\n",
    "    print(\"pcov_AE:\", pcov_AE)\n",
    "\n",
    "    popt = np.concatenate((popt_NC, popt_AE))\n",
    "\n",
    "    Nsamples = 10**4\n",
    "\n",
    "    if bootstrapping:\n",
    "        rng = np.random.default_rng()\n",
    "        R = lambda X: X[rng.integers(len(X), size=len(X))]\n",
    "        ps = np.array([np.concatenate((myfitter(R(X_NC), fun, False), myfitter(R(X_AE), fun, False))) for _ in range(Nsamples)])\n",
    "    else:\n",
    "\n",
    "\n",
    "        Z = np.zeros((3, 3))\n",
    "        CovMatrix = np.bmat([[pcov_NC, Z], [Z, pcov_AE]])\n",
    "\n",
    "        ps = np.random.multivariate_normal(popt, CovMatrix, Nsamples)\n",
    "\n",
    "\n",
    "\n",
    "    mask = (ps[:,1]>0) & (ps[:,4]>0)\n",
    "    ps = ps[mask]\n",
    "    print(\"faulty:\", np.count_nonzero(~mask)/len(mask))\n",
    "\n",
    "\n",
    "    P_NC = np.count_nonzero(flags == 0)/len(flags)  # Prob(Y = NC)\n",
    "    P_AE = np.count_nonzero(flags == 1)/len(flags)  # Prob(Y = AE)\n",
    "\n",
    "\n",
    "    f_XatNC_opt = lambda x: g(x, popt_NC[0], popt_NC[1], popt_NC[2])   # f(x | NC)\n",
    "    f_XatAE_opt = lambda x: g(x, popt_AE[0], popt_AE[1], popt_AE[2])   # f(x | AE)\n",
    "\n",
    "\n",
    "    f_X_opt = lambda x: P_NC*f_XatNC_opt(x) + P_AE*f_XatAE_opt(x)\n",
    "\n",
    "\n",
    "    P_AEatX_opt = lambda x: P_AE*f_XatAE_opt(x)/f_X_opt(x)\n",
    "\n",
    "\n",
    "    def P_AEatX(x, p): \n",
    "        mu_NC, s_NC, l_NC, mu_AE, s_AE, l_AE = p\n",
    "        tmp = ((x**l_AE - 1)/l_AE - mu_AE)**2/s_AE**2 - ((x**l_NC - 1)/l_NC - mu_NC)**2/s_NC**2\n",
    "\n",
    "        if tmp > 300: \n",
    "            return 0\n",
    "        elif tmp < -300:\n",
    "            return 1\n",
    "\n",
    "        return 1/(1 + P_NC/P_AE*s_AE/s_NC*np.abs(x)**(l_NC-l_AE)*np.exp(0.5*tmp))\n",
    "\n",
    "\n",
    "    P_AEatX_opt = lambda x: P_AEatX(x, popt)\n",
    "\n",
    "    xs = np.linspace(1e-3, np.max(X) + 2, 200)\n",
    "\n",
    "    ys = np.array([np.fromiter(map(lambda p: P_AEatX(x, p), ps), float) for x in xs])\n",
    "    yopts = np.array([P_AEatX_opt(x) for x in xs])\n",
    "\n",
    "    mys = np.median(ys, axis=1)\n",
    "    yerr_min = np.quantile(ys, alpha/2, axis=1)\n",
    "    yerr_max = np.quantile(ys, 1-alpha/2, axis=1)\n",
    "\n",
    "    plt.clf()\n",
    "    plt.xlabel(r\"$x = max_{visit} SUV(visit, p)$\")\n",
    "    plt.ylabel(\"P(X=x| NC)\")\n",
    "    plt.hist(X_NC,label=\"NC\", density=True)  \n",
    "    plt.plot(xs, f_XatNC_opt(xs), label=\"NC fit\")  \n",
    "    plt.legend()\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "    plt.clf()\n",
    "    plt.xlabel(r\"$x = max_{visit} SUV(visit, p)$\")\n",
    "    plt.ylabel(\"P(X=x| AE)\")\n",
    "    plt.hist(X_AE,label=\"AE\", density=True)  \n",
    "    plt.plot(xs, f_XatAE_opt(xs), label=\"AE fit\")  \n",
    "    plt.legend()\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "    plt.clf()\n",
    "    plt.xlabel(r\"$x = max_{visit} SUV(visit, p)$\")\n",
    "    plt.ylabel(\"P(AE|X=x)\")\n",
    "    plt.scatter(X_NC, flags_NC, label=\"NC\")  \n",
    "    plt.scatter(X_AE, flags_AE, label=\"AE\")\n",
    "\n",
    "    plt.plot(xs, yopts, label=\"fit model\", color=\"red\")\n",
    "    plt.plot(xs, mys, label=\"median\", color=\"blue\")\n",
    "\n",
    "    plt.fill_between(xs, yerr_min, yerr_max, alpha=0.5, label=\"prediction interval\")\n",
    "\n",
    "    plt.title(title+\": \" + organ)\n",
    "    plt.legend(loc = 'center right')\n",
    "    plt.tight_layout()\n",
    "\n",
    "    if not os.path.exists(figsdir): os.makedirs(figsdir)\n",
    "    plt.savefig(os.path.join(figsdir, prefix + \"_\" + organ + \".pdf\"))\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    suv_dict, flags_dict = load_example_data()\n",
    "    gen_plots(\n",
    "        pert=1,\n",
    "        prefix=\"lung_test\",\n",
    "        alpha=0.05,\n",
    "        title=\"Lung SUV Analysis\",\n",
    "        suv_dict=suv_dict,\n",
    "        flags_dict=flags_dict,\n",
    "        bootstrapping=False\n",
    "    )\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
