{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8d31132f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from fitter import Fitter, get_common_distributions, get_distributions\n",
    "from scipy import io"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "945db1eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "suv_filename = \"../data/suv_percentilesSLOthenUWM.mat\"\n",
    "flags_filename=\"../data/flags_combined.mat\"\n",
    "normal_range_filename =\"../data/normal_range.mat\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e3ac38f3",
   "metadata": {},
   "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": "code",
   "execution_count": 5,
   "id": "b8dadd66",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(58, 17, 100)\n",
      "(58,)\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    " Extracting only data associated to lungs\n",
    "\"\"\"\n",
    "lung_suv = suv_dict['lung_SUVperc_COMBINED'][0:58,:]\n",
    "flags_lung = flags_dict['flags'][0:58,3]  # 0 == NC, 1 == AE\n",
    "\n",
    "print(lung_suv.shape)\n",
    "print(flags_lung.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f92743a3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      "NC_lung: 53\n",
      "AE_lung: 5\n"
     ]
    }
   ],
   "source": [
    "print(flags_lung)\n",
    "print(\"NC_lung:\", np.count_nonzero(flags_lung==0))\n",
    "print(\"AE_lung:\", np.count_nonzero(flags_lung==1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ee2874bd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(58, 17, 100)\n",
      "(58,)\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    " Extracting only data associated to bowel\n",
    "\"\"\"\n",
    "bowel_suv = suv_dict['bowel_SUVperc_COMBINED'][0:58,:]\n",
    "flags_bowel = flags_dict['flags'][0:58,1]  # 0 == NC, 1 == AE\n",
    "\n",
    "print(bowel_suv.shape)\n",
    "print(flags_bowel.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7f427125",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0\n",
      " 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0]\n",
      "NC_bowel: 52\n",
      "AE_bowel: 6\n"
     ]
    }
   ],
   "source": [
    "print(flags_bowel)\n",
    "print(\"NC_bowel:\", np.count_nonzero(flags_bowel==0))\n",
    "print(\"AE_bowel:\", np.count_nonzero(flags_bowel==1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "07fe6dcd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(58, 17, 100)\n",
      "(58,)\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    " Extracting only data associated to thyroid\n",
    "\"\"\"\n",
    "thyroid_suv = suv_dict['thyroid_SUVperc_COMBINED'][0:58,:]\n",
    "flags_thyroid = flags_dict['flags'][0:58,5]  # 0 == NC, 1 == AE\n",
    "\n",
    "print(thyroid_suv.shape)\n",
    "print(flags_thyroid.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9bbea04e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0]\n",
      "NC_thyroid: 49\n",
      "AE_thyroid: 9\n"
     ]
    }
   ],
   "source": [
    "print(flags_thyroid)\n",
    "print(\"NC_thyroid:\", np.count_nonzero(flags_thyroid==0))\n",
    "print(\"AE_thyroid:\", np.count_nonzero(flags_thyroid==1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8e0492f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 95  # percentile of interest\n",
    "X_lung = np.nanmax(lung_suv[:,:,p], axis=1)  # taking SUV max and ignoring NaNs\n",
    "\n",
    "X_bowel = np.nanmax(bowel_suv[:,:,p], axis=1)\n",
    "\n",
    "X_thyroid = np.nanmax(thyroid_suv[:,:,75], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "55666bc8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(53,)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_lung_NC = X_lung[flags_lung == 0]\n",
    "X_lung_AE = X_lung[flags_lung == 1]\n",
    "\n",
    "X_lung_NC.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "bfb51a60",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52,)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_bowel_NC = X_bowel[flags_bowel == 0]\n",
    "X_bowel_AE = X_bowel[flags_bowel == 1]\n",
    "\n",
    "X_bowel_NC.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "32430668",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(49,)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_thyroid_NC = X_thyroid[flags_thyroid == 0]\n",
    "X_thyroid_AE = X_thyroid[flags_thyroid == 1]\n",
    "\n",
    "X_thyroid_NC.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "dc4a8459",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(X_bowel_NC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "16fba333",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(X_lung_NC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a63a15cd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(X_thyroid_NC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "36ca3761",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 2.37582922  8.32889462  2.49695492  2.62231517  3.1779449   2.48935366\n",
      "  2.64083982  3.21895695  2.65162754  2.81365871 13.13622665  4.00094509\n",
      "  6.43130398  7.14293861  2.63764477  4.91679096  3.0133276   4.00735998\n",
      "  2.82554626  2.3756032   3.05244803  5.04009581  2.75375438  2.96612239\n",
      "  2.84024048  3.02956247  3.56154656  2.63551474  3.54600501  3.10588789\n",
      "  5.65029192  3.35723162  2.60515237  2.84188676  2.41808271  2.51086831\n",
      "  2.64900827  2.2847507   3.04761338  2.94424272  2.26703072  2.83462191\n",
      "  3.44818497  2.7778244   5.96339226  3.04745531  2.66029692  5.41591644\n",
      "  2.07756186  2.52008891  1.9738239   7.01769114]\n"
     ]
    }
   ],
   "source": [
    "print(X_bowel_NC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "412eeb4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_NC = np.concatenate((X_lung_NC,X_bowel_NC,X_thyroid_NC))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "527906e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(154,)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_NC.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "97acc5a0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.43247677, 0.29007588, 0.03164464, 0.02109643, 0.01054821,\n",
       "        0.01582232, 0.00527411, 0.        , 0.        , 0.00527411]),\n",
       " array([ 0.82417732,  2.05538226,  3.28658719,  4.51779212,  5.74899706,\n",
       "         6.98020199,  8.21140692,  9.44261186, 10.67381679, 11.90502172,\n",
       "        13.13622665]),\n",
       " <BarContainer object of 10 artists>)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 4548 characters]",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    " Plotting stuff\n",
    "\"\"\"\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.hist(X_NC,density=True, label=\"NC\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "7d5d4b62",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0.,  3., 44., 35., 23., 23., 11.,  2.,  2.,  1.,  2.,  2.,  1.,\n",
       "         0.,  3.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "         1.,  0.]),\n",
       " array([ 0. ,  0.5,  1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,  5. ,\n",
       "         5.5,  6. ,  6.5,  7. ,  7.5,  8. ,  8.5,  9. ,  9.5, 10. , 10.5,\n",
       "        11. , 11.5, 12. , 12.5, 13. , 13.5, 14. ]),\n",
       " <BarContainer object of 28 artists>)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 4500 characters]",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins =[0,0.5,1,1.5,2,2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5,7, 7.5, 8,8.5, 9,9.5, 10, 10.5, 11, 11.5, 12,12.5, 13, 13.5, 14 ]\n",
    "plt.hist(X_NC, bins, label=\"NC\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "92a3e7e2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Fitting 106 distributions:  11%|█▍           | 12/106 [00:01<00:11,  8.53it/s]WARNING:root:SKIPPED kstwo distribution (taking more than 30 seconds)\n",
      "Fitting 106 distributions:  21%|██▋          | 22/106 [00:02<00:07, 10.95it/s]WARNING:root:SKIPPED rv_continuous distribution (taking more than 30 seconds)\n",
      "Fitting 106 distributions:  23%|██▉          | 24/106 [00:02<00:09,  9.10it/s]WARNING:root:SKIPPED rv_histogram distribution (taking more than 30 seconds)\n",
      "Fitting 106 distributions:  94%|███████████▎| 100/106 [00:22<00:06,  1.08s/it]WARNING:root:SKIPPED kappa4 distribution (taking more than 30 seconds)\n",
      "Fitting 106 distributions:  95%|███████████▍| 101/106 [00:31<00:16,  3.27s/it]/home/marija/.conda/envs/pytorch/lib/python3.7/site-packages/scipy/integrate/quadpack.py:880: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
      "  If increasing the limit yields no improvement it is advised to analyze \n",
      "  the integrand in order to determine the difficulties.  If the position of a \n",
      "  local difficulty can be determined (singularity, discontinuity) one will \n",
      "  probably gain from splitting up the interval and calling the integrator \n",
      "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
      "  **opt)\n",
      "WARNING:root:SKIPPED levy_stable distribution (taking more than 30 seconds)\n",
      "Fitting 106 distributions:  96%|███████████▌| 102/106 [00:31<00:09,  2.47s/it]WARNING:root:SKIPPED recipinvgauss distribution (taking more than 30 seconds)\n",
      "Fitting 106 distributions:  97%|███████████▋| 103/106 [00:32<00:05,  1.89s/it]WARNING:root:SKIPPED studentized_range distribution (taking more than 30 seconds)\n",
      "Fitting 106 distributions:  98%|███████████▊| 104/106 [00:32<00:03,  1.56s/it]WARNING:root:SKIPPED vonmises distribution (taking more than 30 seconds)\n",
      "Fitting 106 distributions:  99%|███████████▉| 105/106 [00:33<00:01,  1.32s/it]WARNING:root:SKIPPED vonmises_line distribution (taking more than 30 seconds)\n",
      "Fitting 106 distributions: 100%|████████████| 106/106 [00:33<00:00,  3.13it/s]\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>sumsquare_error</th>\n",
       "      <th>aic</th>\n",
       "      <th>bic</th>\n",
       "      <th>kl_div</th>\n",
       "      <th>ks_statistic</th>\n",
       "      <th>ks_pvalue</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>laplace_asymmetric</th>\n",
       "      <td>0.834069</td>\n",
       "      <td>1007.300261</td>\n",
       "      <td>-788.521547</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.094202</td>\n",
       "      <td>0.121887</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>burr12</th>\n",
       "      <td>0.836352</td>\n",
       "      <td>888.541745</td>\n",
       "      <td>-783.063512</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.062840</td>\n",
       "      <td>0.555837</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>skewcauchy</th>\n",
       "      <td>0.857340</td>\n",
       "      <td>801.601692</td>\n",
       "      <td>-784.283549</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.083710</td>\n",
       "      <td>0.218161</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>johnsonsu</th>\n",
       "      <td>0.861766</td>\n",
       "      <td>884.977607</td>\n",
       "      <td>-778.453770</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.068345</td>\n",
       "      <td>0.448399</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>exponnorm</th>\n",
       "      <td>0.874751</td>\n",
       "      <td>1006.719459</td>\n",
       "      <td>-781.187455</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.095852</td>\n",
       "      <td>0.110524</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                    sumsquare_error          aic         bic  kl_div  \\\n",
       "laplace_asymmetric         0.834069  1007.300261 -788.521547     inf   \n",
       "burr12                     0.836352   888.541745 -783.063512     inf   \n",
       "skewcauchy                 0.857340   801.601692 -784.283549     inf   \n",
       "johnsonsu                  0.861766   884.977607 -778.453770     inf   \n",
       "exponnorm                  0.874751  1006.719459 -781.187455     inf   \n",
       "\n",
       "                    ks_statistic  ks_pvalue  \n",
       "laplace_asymmetric      0.094202   0.121887  \n",
       "burr12                  0.062840   0.555837  \n",
       "skewcauchy              0.083710   0.218161  \n",
       "johnsonsu               0.068345   0.448399  \n",
       "exponnorm               0.095852   0.110524  "
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/marija/.conda/envs/pytorch/lib/python3.7/site-packages/scipy/integrate/quadpack.py:880: IntegrationWarning: The integral is probably divergent, or slowly convergent.\n",
      "  **opt)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 23336 characters]",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = Fitter(X_NC)\n",
    "f.fit()\n",
    "f.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "b6db22d9",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Fitting 3 distributions: 100%|██████████████████| 3/3 [00:02<00:00,  1.32it/s]\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>sumsquare_error</th>\n",
       "      <th>aic</th>\n",
       "      <th>bic</th>\n",
       "      <th>kl_div</th>\n",
       "      <th>ks_statistic</th>\n",
       "      <th>ks_pvalue</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>lognorm</th>\n",
       "      <td>0.919659</td>\n",
       "      <td>921.257538</td>\n",
       "      <td>-773.477696</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.081209</td>\n",
       "      <td>2.480260e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gamma</th>\n",
       "      <td>1.778432</td>\n",
       "      <td>656.279200</td>\n",
       "      <td>-671.917107</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.285486</td>\n",
       "      <td>1.331359e-11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>norm</th>\n",
       "      <td>1.890541</td>\n",
       "      <td>1676.350736</td>\n",
       "      <td>-667.539883</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.190164</td>\n",
       "      <td>2.375160e-05</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         sumsquare_error          aic         bic  kl_div  ks_statistic  \\\n",
       "lognorm         0.919659   921.257538 -773.477696     inf      0.081209   \n",
       "gamma           1.778432   656.279200 -671.917107     inf      0.285486   \n",
       "norm            1.890541  1676.350736 -667.539883     inf      0.190164   \n",
       "\n",
       "            ks_pvalue  \n",
       "lognorm  2.480260e-01  \n",
       "gamma    1.331359e-11  \n",
       "norm     2.375160e-05  "
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 20080 characters]",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = Fitter(X_NC, distributions=['gamma',\n",
    "                          'lognorm',\n",
    "                          \"norm\"])\n",
    "f.fit()\n",
    "f.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "1beea73b",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Set up a helper function for checking p-values against an alpha level, and printing result\n",
    "\"\"\"\n",
    "def check_p_val(p_val, alpha):\n",
    "\n",
    "    if p_val < alpha:\n",
    "        print('We have evidence to reject the null hypothesis.')\n",
    "    else:\n",
    "        print('We do not have evidence to reject the null hypothesis.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "6d9288f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as stats\n",
    "from scipy.stats import normaltest\n",
    "from scipy.stats import kstest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "0a8f8714",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(154,)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = np.log(X_NC)\n",
    "L.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "83ee07de",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Statistic: \t0.50 \n",
      "P-Value: \t9.17e-36\n",
      "\n",
      "We have evidence to reject the null hypothesis.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/marija/.conda/envs/pytorch/lib/python3.7/site-packages/scipy/stats/_continuous_distns.py:4837: IntegrationWarning: The algorithm does not converge.  Roundoff error is detected\n",
      "  in the extrapolation table.  It is assumed that the requested tolerance\n",
      "  cannot be achieved, and that the returned result (if full_output = 1) is \n",
      "  the best which can be obtained.\n",
      "  intg = integrate.quad(f, -xi, np.pi/2, **intg_kwargs)[0]\n"
     ]
    }
   ],
   "source": [
    "stat, p_val = kstest(L, 'norm')\n",
    "print('Statistic: \\t{:1.2f} \\nP-Value: \\t{:1.2e}\\n'.format(stat, p_val))\n",
    "check_p_val(p_val, alpha=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "90b3ec11",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "distribution is not normal\n",
      "distribution is not log-normal\n"
     ]
    }
   ],
   "source": [
    "import scipy as sp\n",
    "\n",
    "p=sp.stats.mstats.normaltest(X_NC, axis=0).pvalue\n",
    "if p<0.025:\n",
    "   print ('distribution is not normal')\n",
    "p=sp.stats.mstats.normaltest(np.log(X_NC), axis=0).pvalue\n",
    "if p<0.025:\n",
    "   print ('distribution is not log-normal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "53177ff3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7529322013966726\n",
      "0.694634906696082\n",
      "1.204700345809746\n",
      "0.186230860341136\n",
      "0.7529322013966726\n"
     ]
    }
   ],
   "source": [
    "shape, location, scale = sp.stats.lognorm.fit(X_NC)\n",
    "mu, sigma = np.log(scale), shape\n",
    "print(shape)\n",
    "print(location)\n",
    "print(scale)\n",
    "print(mu)\n",
    "print(sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "adc791d9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KstestResult(statistic=0.08120944008227199, pvalue=0.24802604184218757)"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.stats.kstest(X_NC, \"lognorm\", sp.stats.lognorm.fit(X_NC))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdff3710",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60af4dca",
   "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.6"
  },
  "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
}
