{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1b60bfaa-5249-4ae1-890d-2ffaf49decbe",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Bayesian modelling using Box-Cox distributions with uncertainties"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23282233-0150-46be-a360-13704b955048",
   "metadata": {},
   "source": [
    "## Data & common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "939c61f5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "from scipy import io\n",
    "from scipy import stats\n",
    "from scipy.optimize import minimize\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "35de6424",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#data_path = \"/home/marija/maki/\"\n",
    "data_path = \"../data/\"\n",
    "suv_filename = data_path + \"suv_percentilesSLOthenUWM.mat\"\n",
    "flags_filename=  data_path + \"flags_combined.mat\"\n",
    "normal_range_filename =  data_path + \"normal_range.mat\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4faafd2f",
   "metadata": {
    "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": "1d5bcc0a",
   "metadata": {},
   "source": [
    "# Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dd57b318",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(58, 17, 100)\n",
      "(58,)\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Extracting only data associated to organ of interest, you need to choose data and column of flags that shoud be taken depended on data \n",
    "\"\"\"\n",
    "# in suv values can be (choose one of them): 'lung_SUVperc_COMBINED' , 'bowel_SUVperc_COMBINED' or 'thyroid_SUVperc_COMBINED'\n",
    "# in flags choose the column depended on the data of interest (1 for bowel, 3 for lung, 5 for thyroid)\n",
    "\n",
    "suv = suv_dict['lung_SUVperc_COMBINED'][0:58,:,:]\n",
    "flags = flags_dict['flags'][0:58,3]\n",
    "\n",
    "print(suv.shape)\n",
    "print(flags.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "642aa5be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(58, 1)\n"
     ]
    }
   ],
   "source": [
    "p = 94  # choose percentile of interest (for lung and bowel 95, for thyroid 75), numerating here is from 0 so take the one less\n",
    "X = np.nanmax(suv[:,:,p], axis=1).reshape(-1,1)\n",
    "\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "077147b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# log of  distribution\n",
    "fun = lambda x, mu, sigma, l: (l-1)*np.log(np.abs(x)) - np.log(sigma*np.sqrt(2*np.pi)) - (((x**l-1)/l)-mu)**2/(2*sigma**2)\n",
    "# distribution\n",
    "g = lambda x, mu, sigma, l: np.abs(x**(l-1)) * np.exp(-(((x**l-1)/l)-mu)**2/(2*sigma**2))/(np.sqrt(2*np.pi)*sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "11a7c5cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def myfitter(data, fun):\n",
    "\n",
    "    nllf = lambda p: -fun(data, *p).sum()\n",
    "\n",
    "    p0 = (1, 0.6, 0.3)\n",
    "    bounds= ((0, 10), (0.2, 2), (-0.5,0.5))\n",
    "    res = minimize(nllf, p0, bounds = bounds, method=\"L-BFGS-B\")\n",
    "    \n",
    "    return res.x, res.hess_inv.todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b84274b5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.26584345  0.21037974 -0.5       ]\n",
      "[[0.00202926 0.00056789 0.00676633]\n",
      " [0.00056789 0.00145807 0.01237667]\n",
      " [0.00676633 0.01237667 0.14698481]]\n"
     ]
    }
   ],
   "source": [
    "popt, pcov = myfitter(X,fun)\n",
    "print(popt)\n",
    "print(pcov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3dac5fb9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.5"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popt[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b195be4",
   "metadata": {},
   "source": [
    "# Final results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "dd48b114-1b24-44f9-82af-648f6de5bd38",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# distribution\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",
    "\n",
    "# log of distribution\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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "feecfa11",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def myfitter(data, fun, cov = True):\n",
    "\n",
    "    nllf = lambda p: -fun(data, *p).sum()\n",
    "\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",
    "\n",
    "    if not cov: return res.x\n",
    "        \n",
    "    return res.x, res.hess_inv.todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "04338a01",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "figsdir  = \"./figs/\"\n",
    "\n",
    "# function for generating results\n",
    "def gen_plots(pert, organ, prefix, alpha, title, bootstrapping = False):\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",
    "  \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",
    "    #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",
    "        # 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), myfitter(R(X_AE), fun, False))) for _ in range(Nsamples)])\n",
    "        \n",
    "    else: \n",
    "        # approximationg parameters with multivatiate distribution\n",
    "        Z = np.zeros((3, 3))\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",
    "\n",
    "    # correct for outliers\n",
    "    mask = (ps[:,1]>0) & (ps[:,4]>0)\n",
    "    ps = ps[mask]\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], 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",
    "    # 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",
    "    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",
    "     # conditional probability for AE at X\n",
    "    P_AEatX_opt = lambda x: P_AEatX(x, popt)\n",
    "             \n",
    "    # values on x axis\n",
    "    xs = np.linspace(1e-3, np.max(X) + 2, 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 = np.array([P_AEatX_opt(x) for x in 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",
    "    #plt.scatter(X[flags==0], flags[flags==0],label=\"NC\")  \n",
    "    #plt.scatter(X[flags==1], flags[flags==1],label=\"AE\")\n",
    "    #plt.plot(xs, yopts, label=\"fit model\",color=\"red\")\n",
    "    #plt.errorbar(xs, mys, yerr = [mys-yerr_min, yerr_max-mys], label=\"direct variance\",fmt=' ')\n",
    "    #plt.xlim([0,3])\n",
    "    #plt.legend()\n",
    "    \n",
    "    # plot of conditional distribution for NC\n",
    "    plt.clf()\n",
    "    plt.xlabel(r\"$x = max_{visit} SUV(visit, p)$\")\n",
    "    plt.ylabel(\"P(X=x| NC)\")\n",
    "    \n",
    "    plt.hist(X_NC,label=\"NC\", density=True)  \n",
    "    plt.plot(xs, f_XatNC_opt(xs), label=\"NC fit\")  \n",
    "   \n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "    # plot of conditional distribution for AE\n",
    "    plt.clf()\n",
    "    plt.xlabel(r\"$x = max_{visit} SUV(visit, p)$\")\n",
    "    plt.ylabel(\"P(X=x| AE)\")\n",
    "    \n",
    "    plt.hist(X_AE,label=\"AE\", density=True)  \n",
    "    plt.plot(xs, f_XatAE_opt(xs), label=\"AE fit\")  \n",
    "   \n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "    # generate plots of P(X|AE)\n",
    "    plt.clf()\n",
    "    plt.xlabel(r\"$x = max_{visit} SUV(visit, p)$\")\n",
    "    plt.ylabel(\"P(AE|X=x)\")\n",
    "    \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",
    "            \n",
    "    plt.savefig(figsdir + prefix + \"_\" + organ + \".pdf\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "606a133f-27e4-43e0-aca6-47fdf5bec12f",
   "metadata": {},
   "source": [
    "## Multivariate normal (MVN) approx of parameter pdf -- \"parametric BS\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "3bd68173-e600-49a9-9275-0e72a34e089d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [ 0.24241008  0.18097193 -0.16096351]\n",
      "popt_AE: [ 0.39013585  0.04654895 -2.02453434]\n",
      "pcov_NC: [[1.32581435e-02 7.27444201e-03 3.30504381e-01]\n",
      " [7.27444201e-03 1.37974708e-02 4.37917022e-01]\n",
      " [3.30504381e-01 4.37917022e-01 1.53211786e+01]]\n",
      "pcov_AE: [[0.09003669 0.0253148  0.6911939 ]\n",
      " [0.0253148  0.00746009 0.19644863]\n",
      " [0.6911939  0.19644863 5.34630204]]\n",
      "faulty: 0.3366\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 14308 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 13992 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 26004 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [ 0.43785532  0.02927488 -1.98912856]\n",
      "popt_AE: [ 0.92125498  0.14138068 -0.67747973]\n",
      "pcov_NC: [[0.0540009  0.01294842 0.3089749 ]\n",
      " [0.01294842 0.00312167 0.07412538]\n",
      " [0.3089749  0.07412538 1.76873243]]\n",
      "pcov_AE: [[0.87174071 0.33898716 1.24929408]\n",
      " [0.33898716 0.13569242 0.48849789]\n",
      " [1.24929408 0.48849789 1.80520867]]\n",
      "faulty: 0.5544\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 15072 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 13928 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 31636 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [ 0.50962594  0.24370391 -0.42937612]\n",
      "popt_AE: [2.18426726 0.86465323 1.06853018]\n",
      "pcov_NC: [[0.00548475 0.00233591 0.02114286]\n",
      " [0.00233591 0.00279503 0.01237538]\n",
      " [0.02114286 0.01237538 0.10656002]]\n",
      "pcov_AE: [[4.03191476 2.38986582 2.5599894 ]\n",
      " [2.38986582 1.51162447 1.57893243]\n",
      " [2.5599894  1.57893243 1.69508472]]\n",
      "faulty: 0.2394\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13784 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 12276 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 35944 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "title=\"Box-Cox with estimated uncertainties\"\n",
    "gen_plots(95, \"lung\", \"boxcox_mvn\", 0.05, title)\n",
    "gen_plots(95, \"bowel\", \"boxcox_mvn\", 0.05, title)\n",
    "gen_plots(75, \"thyroid\", \"boxcox_mvn\", 0.05, title)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8937c74a-425c-40b6-bf90-650d41cd9201",
   "metadata": {},
   "source": [
    "## Non-parametric Bootstrapping (BS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "2f16c726",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [ 0.24241008  0.18097193 -0.16096351]\n",
      "popt_AE: [ 0.39013585  0.04654895 -2.02453434]\n",
      "pcov_NC: [[1.32581435e-02 7.27444201e-03 3.30504381e-01]\n",
      " [7.27444201e-03 1.37974708e-02 4.37917022e-01]\n",
      " [3.30504381e-01 4.37917022e-01 1.53211786e+01]]\n",
      "pcov_AE: [[0.09003669 0.0253148  0.6911939 ]\n",
      " [0.0253148  0.00746009 0.19644863]\n",
      " [0.6911939  0.19644863 5.34630204]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 14308 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 13992 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 26360 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [ 0.43785532  0.02927488 -1.98912856]\n",
      "popt_AE: [ 0.92125498  0.14138068 -0.67747973]\n",
      "pcov_NC: [[0.0540009  0.01294842 0.3089749 ]\n",
      " [0.01294842 0.00312167 0.07412538]\n",
      " [0.3089749  0.07412538 1.76873243]]\n",
      "pcov_AE: [[0.87174071 0.33898716 1.24929408]\n",
      " [0.33898716 0.13569242 0.48849789]\n",
      " [1.24929408 0.48849789 1.80520867]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 15072 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 13928 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 31780 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [ 0.50962594  0.24370391 -0.42937612]\n",
      "popt_AE: [2.18426726 0.86465323 1.06853018]\n",
      "pcov_NC: [[0.00548475 0.00233591 0.02114286]\n",
      " [0.00233591 0.00279503 0.01237538]\n",
      " [0.02114286 0.01237538 0.10656002]]\n",
      "pcov_AE: [[4.03191476 2.38986582 2.5599894 ]\n",
      " [2.38986582 1.51162447 1.57893243]\n",
      " [2.5599894  1.57893243 1.69508472]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13784 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 12276 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 37592 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "title=\"Box-Cox with estimated uncertainties\"\n",
    "gen_plots(95, \"lung\", \"boxcox_bs\", 0.05, title, bootstrapping = True)\n",
    "gen_plots(95, \"bowel\", \"boxcox_bs\", 0.05, title, bootstrapping = True)\n",
    "gen_plots(75, \"thyroid\", \"boxcox_bs\", 0.05, title, bootstrapping = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95998d6e-a741-417f-a491-9e7d64794124",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "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.10.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
