{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "be915d56-9a8a-422b-9da4-d7ca37dd0b85",
   "metadata": {
    "id": "be915d56-9a8a-422b-9da4-d7ca37dd0b85"
   },
   "source": [
    "# Bayesian modelling\n",
    "\n",
    "Two separate-group discussion: \n",
    "\n",
    "    1. fitting each groups separately \n",
    "    2. obtained distributions of each group are then used in conditional probability to encounter measurement in the second group."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fa3ccb1-05c0-4ab9-a74a-a4e20d1e38f8",
   "metadata": {
    "id": "6fa3ccb1-05c0-4ab9-a74a-a4e20d1e38f8"
   },
   "source": [
    "## Data & common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ee1ce969-7f96-4ca9-882e-866805bd316c",
   "metadata": {},
   "outputs": [],
   "source": [
    "choice_setup = 1  # 0 for MD, 1 for MH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "uq21BzGM7YHZ",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 26515,
     "status": "ok",
     "timestamp": 1715067825951,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -120
    },
    "id": "uq21BzGM7YHZ",
    "outputId": "fed613c4-2be9-49ba-bdc5-0a09399dc2fd"
   },
   "outputs": [],
   "source": [
    "if choice_setup == 0:\n",
    "    from google.colab import drive\n",
    "    drive.mount('/content/drive')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0de55416",
   "metadata": {
    "executionInfo": {
     "elapsed": 1711,
     "status": "ok",
     "timestamp": 1715067832083,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -120
    },
    "id": "0de55416",
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "import scipy\n",
    "from scipy import io\n",
    "from scipy import stats\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "71b466ba",
   "metadata": {
    "executionInfo": {
     "elapsed": 305,
     "status": "ok",
     "timestamp": 1715067835069,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -120
    },
    "id": "71b466ba",
    "tags": []
   },
   "outputs": [],
   "source": [
    "wd = os.getcwd()   # working directory\n",
    "\n",
    "topic = \"bayes\"\n",
    "\n",
    "if choice_setup == 0:\n",
    "    #data_path = \"/home/marija/maki/\"\n",
    "    data_path = \"drive/MyDrive/maki_slovenija/data/\"\n",
    "    figsdir = f\"drive/MyDrive/maki_slovenija/figs_{topic}/\"  # figure directory\n",
    "    figsext = \".eps\"                                         # figure extension\n",
    "elif choice_setup == 1:\n",
    "    data_path = \"../data\"\n",
    "    figsdir = f\"./figs/{topic}\"\n",
    "    figsext = \".pdf\"  # figure extension\n",
    "    \n",
    "\n",
    "if not os.path.exists(figsdir): os.makedirs(figsdir)\n",
    "    \n",
    "suv_filename = os.path.join(data_path, \"suv_percentilesSLOthenUWM.mat\")\n",
    "flags_filename=  os.path.join(data_path, \"flags_combined.mat\")\n",
    "normal_range_filename =  os.path.join(data_path, \"normal_range.mat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c72c613c",
   "metadata": {
    "executionInfo": {
     "elapsed": 2756,
     "status": "ok",
     "timestamp": 1715067840744,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -120
    },
    "id": "c72c613c",
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Loading all data\n",
    "\"\"\"\n",
    "suv_dict = io.loadmat(suv_filename)\n",
    "flags_dict = io.loadmat(flags_filename)\n",
    "normal_range_dict = io.loadmat(normal_range_filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c4a5c84-c72a-4b48-a92c-7e25f5b752fc",
   "metadata": {
    "id": "8c4a5c84-c72a-4b48-a92c-7e25f5b752fc"
   },
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "eeeaa8fe",
   "metadata": {
    "executionInfo": {
     "elapsed": 11,
     "status": "ok",
     "timestamp": 1715067844510,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -120
    },
    "id": "eeeaa8fe",
    "tags": []
   },
   "outputs": [],
   "source": [
    "def myfitter(data, fun, cov = True):\n",
    "\n",
    "    nllf = lambda p: -fun(data, *p).sum()\n",
    "\n",
    "    p0 = (1, 0.6)\n",
    "    bounds= ((0, 10), (0.1, 4))\n",
    "    res = minimize(nllf, p0, bounds = bounds, method=\"L-BFGS-B\")\n",
    "\n",
    "    if not cov: return res.x\n",
    "\n",
    "    return res.x, res.hess_inv.todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8262d64c",
   "metadata": {
    "executionInfo": {
     "elapsed": 528,
     "status": "ok",
     "timestamp": 1715067865707,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -120
    },
    "id": "8262d64c",
    "tags": []
   },
   "outputs": [],
   "source": [
    "def gen_plots(pert, organ, dist, prefix, alpha, title = \"\",\n",
    "              bootstrapping = \"nonparam\", Nsamples = 10**4, \n",
    "              xvert = None):\n",
    "    \n",
    "    suv = suv_dict[organ + '_SUVperc_COMBINED'][0:58,:,:]\n",
    "\n",
    "    if organ == \"lung\":\n",
    "        idx = 3\n",
    "    elif organ == \"bowel\":\n",
    "        idx = 1\n",
    "    elif organ == \"thyroid\":\n",
    "        idx = 5\n",
    "    else:\n",
    "        print(\"This organ is not supported\")\n",
    "        return\n",
    "        \n",
    "    if dist == \"trunc_Gauss\":\n",
    "        # create a custom function (\"truncated\" Gaussian distribution)\n",
    "        fun = lambda x, mu, sigma: -(x-mu)**2/(2*sigma**2) + np.log(2/np.pi)/2 - np.log(sigma*(1 + scipy.special.erf(mu/(np.sqrt(2)*sigma))))\n",
    "        g = lambda x, mu, sigma: np.exp(-(x-mu)**2/(2*sigma**2))*np.sqrt(2/np.pi)/(sigma*(1 + scipy.special.erf(mu/(np.sqrt(2)*sigma))))\n",
    "    elif dist == \"Gauss\":\n",
    "        fun =  lambda x, mu, sigma: -(x-mu)**2/(2*sigma**2) - np.log(sigma*np.sqrt(2*np.pi))\n",
    "        g = lambda x, mu, sigma: np.exp(-(x-mu)**2/(2*sigma**2))/(sigma*np.sqrt(2*np.pi))\n",
    "    elif dist == \"log_normal\":\n",
    "        fun =  lambda x, mu, sigma: -(np.log(x) - mu)**2/(2*sigma**2) - np.log(x*sigma*np.sqrt(2*np.pi))\n",
    "        g = lambda x, mu, sigma: np.exp(-(np.log(x)-mu)**2/(2*sigma**2))/(x*sigma*np.sqrt(2*np.pi))\n",
    "    else:\n",
    "        print(\"This distribution is not supported\")\n",
    "        return\n",
    "    \n",
    "    flags = flags_dict['flags'][0:58, idx]\n",
    "\n",
    "    # choose percentile of interest (for lung and bowel 95, for thyroid 75),\n",
    "    # numerating here is from 0 so take the one less\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",
    "    \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",
    "    if bootstrapping == \"nonparam\":\n",
    "        # obtaining parameters via bootstrapping method\n",
    "        rng = np.random.default_rng()\n",
    "\n",
    "        # helper function to sample with replacement of same size sample\n",
    "        R = lambda X: X[rng.integers(len(X), size=len(X))]\n",
    "\n",
    "        ps = np.array([np.concatenate((myfitter(R(X_NC), fun, False), \n",
    "                                       myfitter(R(X_AE), fun, False))) \n",
    "                       for _ in range(Nsamples)])\n",
    "    elif bootstrapping == \"param\":\n",
    "        # approximationg parameters with multivatiate distribution\n",
    "        Z = np.zeros((len(popt_NC), len(popt_AE)))\n",
    "        CovMatrix = np.bmat([[pcov_NC, Z], [Z, pcov_AE]])\n",
    "\n",
    "        #generate parameters and points in x axis\n",
    "        ps = np.random.multivariate_normal(popt, CovMatrix, Nsamples)\n",
    "    else:\n",
    "        print(\"This boostrapping  is not supported\")\n",
    "        return\n",
    "\n",
    "    \n",
    "    # correct for outliers\n",
    "    mask = (ps[:,1]>0) & (ps[:,3]>0)\n",
    "    ps = ps[mask]\n",
    "\n",
    "    print(\"faulty:\", np.count_nonzero(~mask)/len(mask))\n",
    "\n",
    "    # probability for NC and AE\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",
    "    #firstly for optimal parameter values\n",
    "    f_XatNC_opt = lambda x: g(x, popt_NC[0], popt_NC[1])   # f(x | NC)\n",
    "    f_XatAE_opt = lambda x: g(x, popt_AE[0], popt_AE[1])   # f(x | AE)\n",
    "\n",
    "    # probability density for X \n",
    "    f_X_opt = lambda x: P_NC*f_XatNC_opt(x) + P_AE*f_XatAE_opt(x)\n",
    "\n",
    "    # conditional probability for AE at X\n",
    "    P_AEatX_opt = lambda x: P_AE*f_XatAE_opt(x)/f_X_opt(x)\n",
    "\n",
    "    # conditional probability for AE at X for generated parameters\n",
    "    P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n",
    "\n",
    "    #generate parameters and points in x axis\n",
    "    xs = np.linspace(1e-3, X.max()+1, 200)\n",
    "\n",
    "    # generate values in y axis\n",
    "    ys = np.array([np.fromiter(map(lambda p: P_AEatX(x, p), ps),float) for x in xs])\n",
    "    yopts = P_AEatX_opt(xs)\n",
    "\n",
    "    # doing statistics on values\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",
    "    # generate plots\n",
    "    fs = 16\n",
    "    params = {'font.size': fs, 'axes.labelsize': fs,'axes.titlesize':fs, \n",
    "              'legend.fontsize': fs, 'xtick.labelsize': fs, 'ytick.labelsize': fs}\n",
    "    plt.rcParams.update(params)\n",
    "\n",
    "    fig, ax = plt.subplots(figsize = (7, 5))\n",
    "\n",
    "    plt.xlabel(f\"$x = \\\\mathrm{{max}}(SUV(visit,p = {{{pert}}} \\\\%))$\")\n",
    "    plt.ylabel(r\"$P(Y = AE|X = x)$\")\n",
    "   \n",
    "    ax.set_xlabel(f\"$x = \\\\mathrm{{max}}(SUV(visit,p = {{{pert}}} \\\\%))$\")\n",
    "    ax.set_ylabel(f\"$P(Y = AE|X = x)$\")\n",
    "\n",
    "    ax.scatter(X[flags==0], flags[flags==0], label=\"NC\",zorder=2.5)\n",
    "    ax.scatter(X[flags==1], flags[flags==1], label=\"AE\",zorder=2.5)\n",
    "\n",
    "    ax.plot(xs, yopts, label=\"optimal\", color=\"red\",zorder=3)\n",
    "    ax.plot(xs, mys, label=\"median\", color=\"green\",zorder=3)\n",
    "\n",
    "    #plt.fill_between(xs, yerr_min, yerr_max, alpha=0.5, label=\"confidence band\")\n",
    "    ax.fill_between(xs, yerr_min, yerr_max, color='xkcd:beige', alpha=0.7, label=\"confidence interval\")\n",
    "\n",
    "    if title != \"\": ax.set_title(title+\":\" + organ)\n",
    "    if xvert is not None: \n",
    "        ax.axvline(x = xvert, linestyle = \"--\", color = \"k\", alpha = 0.7, label=\"example\")\n",
    "        I = lambda v: scipy.interpolate.interp1d(xs, v)(xvert)        \n",
    "        display(pd.DataFrame({\"type\": [\"optimal\", \"median\", \"LCL\", \"UCL\"], \n",
    "                              f\"value@{xvert}\": [I(yopts), I(mys), I(yerr_min), I(yerr_max)]}))\n",
    "\n",
    "    ax.text(0.1, 0.93, \"(\" + organ + \")\", fontsize = 16, \n",
    "            horizontalalignment = 'center', verticalalignment = 'center', \n",
    "            transform = fig.transFigure)\n",
    "\n",
    "    ax.legend(loc = 'center right', fancybox=True, fontsize='10', framealpha=0.5)\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    \n",
    "    basename =  f\"{prefix}_{dist}_{organ}{figsext}\"\n",
    "    plt.savefig(os.path.join(figsdir,basename))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a974ab00-e992-4de8-8a9f-431264fe376a",
   "metadata": {
    "id": "a974ab00-e992-4de8-8a9f-431264fe376a"
   },
   "source": [
    "### Non-parametric boostrapping: sampling with replacement (primary choice)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "00ea6e5e-5e09-45f1-a8e3-0499261ac781",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e79e7dc7-c3b3-46bc-a11b-fc44c3b2f305",
   "metadata": {},
   "source": [
    "#### Lung"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "BDoKoZ0dMekG",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "executionInfo": {
     "elapsed": 75311,
     "status": "ok",
     "timestamp": 1715067952192,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -120
    },
    "id": "BDoKoZ0dMekG",
    "outputId": "4498b0d8-e711-414b-9f12-2ea62281d491"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [1.30761056 0.25845619]\n",
      "popt_AE: [2.45119943 0.8708103 ]\n",
      "pcov_NC: [[1.18087018e-03 1.19806173e-06]\n",
      " [1.19806173e-06 6.30252819e-04]]\n",
      "pcov_AE: [[0.16908412 0.03501477]\n",
      " [0.03501477 0.07755902]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>type</th>\n",
       "      <th>value@2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>optimal</td>\n",
       "      <td>0.4699474053976193</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>median</td>\n",
       "      <td>0.6400362738480477</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>LCL</td>\n",
       "      <td>0.12335997081396466</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>UCL</td>\n",
       "      <td>0.9945968429081302</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      type              value@2\n",
       "0  optimal   0.4699474053976193\n",
       "1   median   0.6400362738480477\n",
       "2      LCL  0.12335997081396466\n",
       "3      UCL   0.9945968429081302"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 72916 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [0.25012563 0.18862822]\n",
      "popt_AE: [0.84495162 0.3054373 ]\n",
      "pcov_NC: [[ 6.83750245e-04 -8.23179537e-06]\n",
      " [-8.23179537e-06  3.41031588e-04]]\n",
      "pcov_AE: [[0.02184912 0.00118759]\n",
      " [0.00118759 0.00977597]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n",
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>type</th>\n",
       "      <th>value@2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>optimal</td>\n",
       "      <td>0.4482002564354698</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>median</td>\n",
       "      <td>0.5427501460938366</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>LCL</td>\n",
       "      <td>0.1843888655804964</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>UCL</td>\n",
       "      <td>0.9213506594067791</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      type             value@2\n",
       "0  optimal  0.4482002564354698\n",
       "1   median  0.5427501460938366\n",
       "2      LCL  0.1843888655804964\n",
       "3      UCL  0.9213506594067791"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 72552 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3min 32s, sys: 1.71 s, total: 3min 34s\n",
      "Wall time: 3min 32s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "gen_plots(95, \"lung\", \"Gauss\", \"Bayes_bs\", alpha, bootstrapping = \"nonparam\", xvert = 2)\n",
    "gen_plots(95, \"lung\", \"log_normal\", \"Bayes_bs\", alpha, bootstrapping = \"nonparam\", xvert = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2b16dac-0952-40ad-b991-12dc7b9d8745",
   "metadata": {},
   "source": [
    "#### Bowel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d3bea1f3-a899-44aa-b410-32b0e3163a7e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [3.29019865 1.61425315]\n",
      "popt_AE: [4.86590897 2.14788649]\n",
      "pcov_NC: [[0.05087085 0.00014864]\n",
      " [0.00014864 0.02513256]]\n",
      "pcov_AE: [[ 0.80285003 -0.00968201]\n",
      " [-0.00968201  0.38917906]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n",
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 68036 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [1.11554265 0.35077766]\n",
      "popt_AE: [1.49665276 0.40295539]\n",
      "pcov_NC: [[2.57517154e-03 1.39935384e-05]\n",
      " [1.39935384e-05 1.18591508e-03]]\n",
      "pcov_AE: [[2.70791363e-02 8.31548861e-05]\n",
      " [8.31548861e-05 1.40269527e-02]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n",
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 74264 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3min 45s, sys: 1.32 s, total: 3min 46s\n",
      "Wall time: 3min 45s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "gen_plots(95, \"bowel\", \"Gauss\", \"Bayes_bs\", alpha, bootstrapping = \"nonparam\")\n",
    "gen_plots(95, \"bowel\", \"log_normal\", \"Bayes_bs\", alpha, bootstrapping = \"nonparam\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9276cee-096a-40ac-8011-74b87c64b9a9",
   "metadata": {},
   "source": [
    "#### Thyroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6463d837-a8f8-406b-9403-4fccbe0a5594",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [1.93132713 0.86265893]\n",
      "popt_AE: [3.0787539  0.80272578]\n",
      "pcov_NC: [[ 0.01472579 -0.00035012]\n",
      " [-0.00035012  0.00723928]]\n",
      "pcov_AE: [[7.10424716e-02 7.45174552e-05]\n",
      " [7.45174552e-05 3.55500847e-02]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 69028 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [0.59700666 0.32538502]\n",
      "popt_AE: [1.08690741 0.28239884]\n",
      "pcov_NC: [[ 3.60331259e-03 -2.12283067e-06]\n",
      " [-2.12283067e-06  1.08060047e-03]]\n",
      "pcov_AE: [[ 8.86542388e-03 -7.03060232e-06]\n",
      " [-7.03060232e-06  4.59559687e-03]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n",
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 68876 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4min 57s, sys: 1.65 s, total: 4min 59s\n",
      "Wall time: 4min 58s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "gen_plots(75, \"thyroid\", \"Gauss\", \"Bayes_bs\", alpha, bootstrapping = \"nonparam\")\n",
    "gen_plots(75, \"thyroid\", \"log_normal\", \"Bayes_bs\", alpha, bootstrapping = \"nonparam\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bb81864-0179-4a0a-b557-c1d2ae069d64",
   "metadata": {},
   "source": [
    "### Parametric boostrapping: MVN approx of params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8efb3048-08a3-4762-b9fe-c26b96c22beb",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efb248a7-030e-4d97-b195-4cab5fd2ee11",
   "metadata": {},
   "source": [
    "#### Lung"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "50f56580",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "executionInfo": {
     "elapsed": 33129,
     "status": "ok",
     "timestamp": 1709041941629,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -60
    },
    "id": "50f56580",
    "outputId": "fdc915b1-73b5-46c2-92a4-1f428bbd2134",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [1.30761056 0.25845619]\n",
      "popt_AE: [2.45119943 0.8708103 ]\n",
      "pcov_NC: [[1.18087018e-03 1.19806173e-06]\n",
      " [1.19806173e-06 6.30252819e-04]]\n",
      "pcov_AE: [[0.16908412 0.03501477]\n",
      " [0.03501477 0.07755902]]\n",
      "faulty: 0.0012\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 62324 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [0.25012563 0.18862822]\n",
      "popt_AE: [0.84495162 0.3054373 ]\n",
      "pcov_NC: [[ 6.83750245e-04 -8.23179537e-06]\n",
      " [-8.23179537e-06  3.41031588e-04]]\n",
      "pcov_AE: [[0.02184912 0.00118759]\n",
      " [0.00118759 0.00977597]]\n",
      "faulty: 0.0011\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n",
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 68824 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min, sys: 1.54 s, total: 1min 1s\n",
      "Wall time: 1min\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "#gen_plots(95, \"lung\", \"trunc_Gauss\", \"Bayes_mvn\", alpha,  title)\n",
    "gen_plots(95, \"lung\", \"Gauss\", \"Bayes_mvn\", alpha, bootstrapping = \"param\")\n",
    "gen_plots(95, \"lung\", \"log_normal\", \"Bayes_mvn\", alpha, bootstrapping = \"param\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61f8e440-4cae-427d-8864-c208f6975d6d",
   "metadata": {},
   "source": [
    "#### Bowel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "02ac712c-1658-4120-92bb-84cb78d9bb9d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [3.29019865 1.61425315]\n",
      "popt_AE: [4.86590897 2.14788649]\n",
      "pcov_NC: [[0.05087085 0.00014864]\n",
      " [0.00014864 0.02513256]]\n",
      "pcov_AE: [[ 0.80285003 -0.00968201]\n",
      " [-0.00968201  0.38917906]]\n",
      "faulty: 0.0001\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 67700 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [1.11554265 0.35077766]\n",
      "popt_AE: [1.49665276 0.40295539]\n",
      "pcov_NC: [[2.57517154e-03 1.39935384e-05]\n",
      " [1.39935384e-05 1.18591508e-03]]\n",
      "pcov_AE: [[2.70791363e-02 8.31548861e-05]\n",
      " [8.31548861e-05 1.40269527e-02]]\n",
      "faulty: 0.0005\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n",
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 76024 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 37s, sys: 1.75 s, total: 1min 38s\n",
      "Wall time: 1min 38s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "#gen_plots(95, \"bowel\", \"trunc_Gauss\", \"Bayes_mvn\", alpha,  title, bootstrapping = \"param\")\n",
    "gen_plots(95, \"bowel\", \"Gauss\", \"Bayes_mvn\", alpha, bootstrapping = \"param\")\n",
    "gen_plots(95, \"bowel\", \"log_normal\", \"Bayes_mvn\", alpha, bootstrapping = \"param\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26314f04-8a6f-441a-880b-ddb84d862950",
   "metadata": {},
   "source": [
    "#### Thyroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d4ff1bf6-9705-49a9-b1f2-5f6cb12d4740",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [1.93132713 0.86265893]\n",
      "popt_AE: [3.0787539  0.80272578]\n",
      "pcov_NC: [[ 0.01472579 -0.00035012]\n",
      " [-0.00035012  0.00723928]]\n",
      "pcov_AE: [[7.10424716e-02 7.45174552e-05]\n",
      " [7.45174552e-05 3.55500847e-02]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 64868 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [0.59700666 0.32538502]\n",
      "popt_AE: [1.08690741 0.28239884]\n",
      "pcov_NC: [[ 3.60331259e-03 -2.12283067e-06]\n",
      " [-2.12283067e-06  1.08060047e-03]]\n",
      "pcov_AE: [[ 8.86542388e-03 -7.03060232e-06]\n",
      " [-7.03060232e-06  4.59559687e-03]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
      "  P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 68800 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 51.4 s, sys: 1.31 s, total: 52.8 s\n",
      "Wall time: 51.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "#gen_plots(75, \"thyroid\", \"trunc_Gauss\", \"Bayes_mvn\", alpha,  title, bootstrapping = \"param\")\n",
    "gen_plots(75, \"thyroid\", \"Gauss\", \"Bayes_mvn\", alpha, bootstrapping = \"param\")\n",
    "gen_plots(75, \"thyroid\", \"log_normal\", \"Bayes_mvn\", alpha, bootstrapping = \"param\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
