{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8f50efc2",
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'pymc3'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[6], line 4\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpyplot\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mplt\u001b[39;00m\n\u001b[0;32m      3\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mscipy\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m io\n\u001b[1;32m----> 4\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mpymc3\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpm\u001b[39;00m\n\u001b[0;32m      5\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01marviz\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01maz\u001b[39;00m\n\u001b[0;32m      7\u001b[0m \u001b[38;5;66;03m# Load the data\u001b[39;00m\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'pymc3'"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import io\n",
    "import pymc3 as pm\n",
    "import arviz as az\n",
    "\n",
    "# Load the data\n",
    "data_path = \"../data/\"\n",
    "suv_file = data_path + \"suv_percentilesSLOthenUWM.mat\"\n",
    "flags_file = data_path + \"flags_combined.mat\"\n",
    "\n",
    "# Load .mat files\n",
    "suv_dict = io.loadmat(suv_file)\n",
    "flags_dict = io.loadmat(flags_file)\n",
    "\n",
    "# Extract relevant data\n",
    "suv = suv_dict['lung_SUVperc_COMBINED'][0:58, :, :]\n",
    "flags = flags_dict['flags'][0:58, 3]\n",
    "\n",
    "# Choose percentile of interest\n",
    "p = 94  # For lung and bowel, use 94 (0-indexed)\n",
    "X = np.nanmax(suv[:, :, p], axis=1).reshape(-1, 1)\n",
    "\n",
    "# Handle missing values\n",
    "valid_indices = ~np.isnan(X).flatten()\n",
    "X = X[valid_indices]\n",
    "flags = flags[valid_indices]\n",
    "\n",
    "# Standardize the predictor\n",
    "X_mean = X.mean()\n",
    "X_std = X.std()\n",
    "X_stdized = (X - X_mean) / X_std\n",
    "\n",
    "# Define the Bayesian logistic regression model\n",
    "with pm.Model() as logistic_model:\n",
    "    # Priors for the intercept and slope\n",
    "    intercept = pm.Normal('intercept', mu=0, sigma=10)\n",
    "    slope = pm.Normal('slope', mu=0, sigma=10)\n",
    "    \n",
    "    # Linear combination\n",
    "    linear_combination = intercept + slope * X_stdized.flatten()\n",
    "    \n",
    "    # Likelihood\n",
    "    likelihood = pm.Bernoulli('likelihood', logit_p=linear_combination, observed=flags)\n",
    "    \n",
    "    # Sampling\n",
    "    trace = pm.sample(2000, tune=1000, target_accept=0.95, return_inferencedata=True)\n",
    "\n",
    "# Summary of the posterior distributions\n",
    "az.summary(trace, hdi_prob=0.95)\n",
    "\n",
    "# Plotting the posterior distributions\n",
    "az.plot_trace(trace)\n",
    "plt.show()\n",
    "\n",
    "# Predictive probabilities\n",
    "x_vals = np.linspace(X_stdized.min(), X_stdized.max(), 100)\n",
    "with logistic_model:\n",
    "    pm.set_data({\"X\": x_vals})\n",
    "    posterior_predictive = pm.sample_posterior_predictive(trace, var_names=[\"likelihood\"])\n",
    "\n",
    "# Compute the mean predicted probability\n",
    "mean_predicted_probs = posterior_predictive[\"likelihood\"].mean(axis=0)\n",
    "\n",
    "# Plotting\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.title(\"Bayesian Logistic Regression: Prediction Probability of AE Depending on maxSUV Value\")\n",
    "plt.xlabel(r\"$x = \\max_{visit} SUV(visit, p)$ (standardized)\")\n",
    "plt.ylabel(\"P(AE | X = x)\")\n",
    "\n",
    "# Scatter data points\n",
    "plt.scatter(X_stdized[flags == 0], flags[flags == 0], label=\"NC\", color='blue', alpha=0.6)\n",
    "plt.scatter(X_stdized[flags == 1], flags[flags == 1], label=\"AE\", color='red', alpha=0.6)\n",
    "\n",
    "# Plot the mean logistic curve (S-shape)\n",
    "plt.plot(x_vals, mean_predicted_probs, color='green', label=\"Bayesian Logistic Regression\")\n",
    "\n",
    "plt.legend(loc='center right')\n",
    "plt.grid(True)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e66b489",
   "metadata": {},
   "source": [
    "We should not do regularization"
   ]
  }
 ],
 "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
}