{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "850e7703",
   "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[1], 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 with weak priors\n",
    "with pm.Model() as logistic_model:\n",
    "    # Weak priors (wide normal distribution for intercept and slope)\n",
    "    intercept = pm.Normal('intercept', mu=0, sigma=1e6)  # Weak prior for intercept\n",
    "    slope = pm.Normal('slope', mu=0, sigma=1e6)  # Weak prior for slope\n",
    "    \n",
    "    # Logistic regression model: logit link\n",
    "    linear_combination = intercept + slope * X_stdized.flatten()\n",
    "    likelihood = pm.Bernoulli('likelihood', logit_p=linear_combination, observed=flags)\n",
    "    \n",
    "    # Sample from the posterior using MCMC\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"
   ]
  }
 ],
 "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
}
