{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "be915d56-9a8a-422b-9da4-d7ca37dd0b85",
   "metadata": {},
   "source": [
    "# Bayesian modelling using normal distributions with uncertainties"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fa3ccb1-05c0-4ab9-a74a-a4e20d1e38f8",
   "metadata": {},
   "source": [
    "## Data & common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0de55416",
   "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": 3,
   "id": "71b466ba",
   "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": 4,
   "id": "c72c613c",
   "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": "8c4a5c84-c72a-4b48-a92c-7e25f5b752fc",
   "metadata": {},
   "source": [
    "## Final results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "eeeaa8fe",
   "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)\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": 13,
   "id": "8262d64c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "figsdir = \"./figs/\"\n",
    "\n",
    "def gen_plots(pert, organ, dist, 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",
    "    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",
    "        \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**3\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",
    "    else: \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",
    "        \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: P_AE*g(x, p[2], p[3])/(P_NC*g(x, p[0], p[1]) + P_AE*g(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",
    "    #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",
    "    # generate plots\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",
    "    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 + \"_\"+  dist + \"_\" + organ + \".pdf\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a974ab00-e992-4de8-8a9f-431264fe376a",
   "metadata": {},
   "source": [
    "### Multivaiate normal approx of parameter pdf == \"parametric BS\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "50f56580",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [1.30761035 0.25845691]\n",
      "popt_AE: [2.44364741 0.88137533]\n",
      "pcov_NC: [[1.25913825e-03 5.06705885e-06]\n",
      " [5.06705885e-06 6.06774810e-04]]\n",
      "pcov_AE: [[ 0.17534308 -0.003193  ]\n",
      " [-0.003193    0.08975444]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 30296 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: [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.002\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 29796 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.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.002\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_27868/1449691649.py:89: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  P_AEatX = lambda x, p: P_AE*g(x, p[2], p[3])/(P_NC*g(x, p[0], p[1]) + P_AE*g(x, p[2], p[3]))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 32160 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: [3.14290066 1.75796851]\n",
      "popt_AE: [4.76760397 2.25648474]\n",
      "pcov_NC: [[14.70862832 12.92525704]\n",
      " [12.92525704 11.49269931]]\n",
      "pcov_AE: [[ 0.99567237 -0.09317099]\n",
      " [-0.09317099  0.64015636]]\n",
      "faulty: 0.311\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 32280 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: [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"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 32452 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: [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.001\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 37092 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: [1.88855494 0.90927873]\n",
      "popt_AE: [3.07854855 0.80312247]\n",
      "pcov_NC: [[ 0.02367611 -0.00551307]\n",
      " [-0.00551307  0.01323124]]\n",
      "pcov_AE: [[0.29051356 0.57028087]\n",
      " [0.57028087 2.82295922]]\n",
      "faulty: 0.32\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 33224 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: [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, 30928 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.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"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 33536 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "title = \"Bayesian approach with estimated uncertainties\"\n",
    "alpha = 0.05\n",
    "gen_plots(95, \"lung\", \"trunc_Gauss\", \"Bayes_mvn\", alpha,  title)\n",
    "gen_plots(95, \"lung\", \"Gauss\", \"Bayes_mvn\", alpha, title)\n",
    "gen_plots(95, \"lung\", \"log_normal\", \"Bayes_mvn\", alpha, title)\n",
    "\n",
    "gen_plots(95, \"bowel\", \"trunc_Gauss\", \"Bayes_mvn\", alpha,  title)\n",
    "gen_plots(95, \"bowel\", \"Gauss\", \"Bayes_mvn\", alpha, title)\n",
    "gen_plots(95, \"bowel\", \"log_normal\", \"Bayes_mvn\", alpha, title)\n",
    "\n",
    "gen_plots(75, \"thyroid\", \"trunc_Gauss\", \"Bayes_mvn\", alpha,  title)\n",
    "gen_plots(75, \"thyroid\", \"Gauss\", \"Bayes_mvn\", alpha, title)\n",
    "gen_plots(75, \"thyroid\", \"log_normal\", \"Bayes_mvn\", alpha, title)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1dec924-7c34-4c4e-97a0-d6cdc4acf699",
   "metadata": {},
   "source": [
    "### Non-parametric bootstrapping (BS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "079cc43f-19b4-4abc-ad4d-70e098b4e553",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [1.30761035 0.25845691]\n",
      "popt_AE: [2.44364741 0.88137533]\n",
      "pcov_NC: [[1.25913825e-03 5.06705885e-06]\n",
      " [5.06705885e-06 6.06774810e-04]]\n",
      "pcov_AE: [[ 0.17534308 -0.003193  ]\n",
      " [-0.003193    0.08975444]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 33732 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: [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": {
      "image/png": "[encoded data, 33564 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.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_27868/1449691649.py:89: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  P_AEatX = lambda x, p: P_AE*g(x, p[2], p[3])/(P_NC*g(x, p[0], p[1]) + P_AE*g(x, p[2], p[3]))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 33460 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "alpha = 0.05\n",
    "gen_plots(95, \"lung\", \"trunc_Gauss\", \"Bayes_bs\", alpha,  title, bootstrapping = True)\n",
    "gen_plots(95, \"lung\", \"Gauss\", \"Bayes_bs\", alpha, title, bootstrapping = True)\n",
    "gen_plots(95, \"lung\", \"log_normal\", \"Bayes_bs\", alpha, title, bootstrapping = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "71521a47-ce8c-4f17-a3d7-dd34df9b4169",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [3.14290066 1.75796851]\n",
      "popt_AE: [4.76760397 2.25648474]\n",
      "pcov_NC: [[14.70862832 12.92525704]\n",
      " [12.92525704 11.49269931]]\n",
      "pcov_AE: [[ 0.99567237 -0.09317099]\n",
      " [-0.09317099  0.64015636]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 33816 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: [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"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 33304 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: [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_27868/1449691649.py:89: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  P_AEatX = lambda x, p: P_AE*g(x, p[2], p[3])/(P_NC*g(x, p[0], p[1]) + P_AE*g(x, p[2], p[3]))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 35764 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gen_plots(95, \"bowel\", \"trunc_Gauss\", \"Bayes_bs\", alpha,  title, bootstrapping = True)\n",
    "gen_plots(95, \"bowel\", \"Gauss\", \"Bayes_bs\", alpha, title, bootstrapping = True)\n",
    "gen_plots(95, \"bowel\", \"log_normal\", \"Bayes_bs\", alpha, title, bootstrapping = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "31d40eb5-0d21-4feb-b21f-5788680cc52b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt_NC: [1.88855494 0.90927873]\n",
      "popt_AE: [3.07854855 0.80312247]\n",
      "pcov_NC: [[ 0.02367611 -0.00551307]\n",
      " [-0.00551307  0.01323124]]\n",
      "pcov_AE: [[0.29051356 0.57028087]\n",
      " [0.57028087 2.82295922]]\n",
      "faulty: 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 32960 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: [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, 33632 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.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_27868/1449691649.py:89: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  P_AEatX = lambda x, p: P_AE*g(x, p[2], p[3])/(P_NC*g(x, p[0], p[1]) + P_AE*g(x, p[2], p[3]))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 33716 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gen_plots(75, \"thyroid\", \"trunc_Gauss\", \"Bayes_bs\", alpha,  title, bootstrapping = True)\n",
    "gen_plots(75, \"thyroid\", \"Gauss\", \"Bayes_bs\", alpha, title, bootstrapping = True)\n",
    "gen_plots(75, \"thyroid\", \"log_normal\", \"Bayes_bs\", alpha, title, bootstrapping = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dedc1382-7e16-4e15-9986-a8c345672221",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
