3-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "87632a5a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: arch in /srv/conda/envs/notebook/lib/python3.7/site-packages (4.19)\n",
      "Requirement already satisfied: statsmodels>=0.10 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from arch) (0.12.2)\n",
      "Requirement already satisfied: numpy>=1.14 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from arch) (1.20.2)\n",
      "Requirement already satisfied: pandas>=0.23 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from arch) (1.2.4)\n",
      "Requirement already satisfied: cython>=0.29.14 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from arch) (0.29.23)\n",
      "Requirement already satisfied: property-cached>=1.6.4 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from arch) (1.6.4)\n",
      "Requirement already satisfied: scipy>=1.2.3 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from arch) (1.6.2)\n",
      "Requirement already satisfied: python-dateutil>=2.7.3 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pandas>=0.23->arch) (2.8.1)\n",
      "Requirement already satisfied: pytz>=2017.3 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pandas>=0.23->arch) (2021.1)\n",
      "Requirement already satisfied: six>=1.5 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from python-dateutil>=2.7.3->pandas>=0.23->arch) (1.15.0)\n",
      "Requirement already satisfied: patsy>=0.5 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from statsmodels>=0.10->arch) (0.5.1)\n",
      "Requirement already satisfied: numpy in /srv/conda/envs/notebook/lib/python3.7/site-packages (1.20.2)\n",
      "Requirement already satisfied: pandas in /srv/conda/envs/notebook/lib/python3.7/site-packages (1.2.4)\n",
      "Requirement already satisfied: python-dateutil>=2.7.3 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pandas) (2.8.1)\n",
      "Requirement already satisfied: pytz>=2017.3 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pandas) (2021.1)\n",
      "Requirement already satisfied: numpy>=1.16.5 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pandas) (1.20.2)\n",
      "Requirement already satisfied: six>=1.5 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from python-dateutil>=2.7.3->pandas) (1.15.0)\n",
      "Requirement already satisfied: seaborn in /srv/conda/envs/notebook/lib/python3.7/site-packages (0.11.1)\n",
      "Requirement already satisfied: matplotlib>=2.2 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from seaborn) (3.4.1)\n",
      "Requirement already satisfied: scipy>=1.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from seaborn) (1.6.2)\n",
      "Requirement already satisfied: pandas>=0.23 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from seaborn) (1.2.4)\n",
      "Requirement already satisfied: numpy>=1.15 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from seaborn) (1.20.2)\n",
      "Requirement already satisfied: cycler>=0.10 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib>=2.2->seaborn) (0.10.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib>=2.2->seaborn) (1.3.1)\n",
      "Requirement already satisfied: python-dateutil>=2.7 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib>=2.2->seaborn) (2.8.1)\n",
      "Requirement already satisfied: pyparsing>=2.2.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib>=2.2->seaborn) (2.4.7)\n",
      "Requirement already satisfied: pillow>=6.2.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib>=2.2->seaborn) (8.2.0)\n",
      "Requirement already satisfied: six in /srv/conda/envs/notebook/lib/python3.7/site-packages (from cycler>=0.10->matplotlib>=2.2->seaborn) (1.15.0)\n",
      "Requirement already satisfied: pytz>=2017.3 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pandas>=0.23->seaborn) (2021.1)\n",
      "Requirement already satisfied: mlxtend in /srv/conda/envs/notebook/lib/python3.7/site-packages (0.18.0)\n",
      "Requirement already satisfied: scipy>=1.2.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from mlxtend) (1.6.2)\n",
      "Requirement already satisfied: numpy>=1.16.2 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from mlxtend) (1.20.2)\n",
      "Requirement already satisfied: matplotlib>=3.0.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from mlxtend) (3.4.1)\n",
      "Requirement already satisfied: joblib>=0.13.2 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from mlxtend) (1.0.1)\n",
      "Requirement already satisfied: setuptools in /srv/conda/envs/notebook/lib/python3.7/site-packages (from mlxtend) (49.6.0.post20210108)\n",
      "Requirement already satisfied: scikit-learn>=0.20.3 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from mlxtend) (0.24.1)\n",
      "Requirement already satisfied: pandas>=0.24.2 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from mlxtend) (1.2.4)\n",
      "Requirement already satisfied: python-dateutil>=2.7 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib>=3.0.0->mlxtend) (2.8.1)\n",
      "Requirement already satisfied: pyparsing>=2.2.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib>=3.0.0->mlxtend) (2.4.7)\n",
      "Requirement already satisfied: cycler>=0.10 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib>=3.0.0->mlxtend) (0.10.0)\n",
      "Requirement already satisfied: pillow>=6.2.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib>=3.0.0->mlxtend) (8.2.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib>=3.0.0->mlxtend) (1.3.1)\n",
      "Requirement already satisfied: six in /srv/conda/envs/notebook/lib/python3.7/site-packages (from cycler>=0.10->matplotlib>=3.0.0->mlxtend) (1.15.0)\n",
      "Requirement already satisfied: pytz>=2017.3 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pandas>=0.24.2->mlxtend) (2021.1)\n",
      "Requirement already satisfied: threadpoolctl>=2.0.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from scikit-learn>=0.20.3->mlxtend) (2.1.0)\n"
     ]
    }
   ],
   "source": [
    "%%bash \n",
    "pip install arch\n",
    "pip install numpy\n",
    "pip install pandas\n",
    "pip install seaborn\n",
    "pip install mlxtend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0d92282e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker\n",
    "import os\n",
    "import seaborn as sns\n",
    "from scipy import stats\n",
    "from mlxtend.evaluate import permutation_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "50db6a45",
   "metadata": {},
   "outputs": [],
   "source": [
    "#A=pd.read_pickle(\"data/allyears-avg_plate.pkl\")\n",
    "A=pd.read_pickle(\"data/allyears.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "9adc3479",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.999000099990002e-05\n"
     ]
    }
   ],
   "source": [
    "# Permutation test for combined off-festival vs festival period\n",
    "p_value = permutation_test( A.loc[A['festival'] == 'during', 'Bioindicator'] , \n",
    "                            A.loc[A['festival'] != 'during', 'Bioindicator'] , \n",
    "                           method='approximate',\n",
    "                           num_rounds=10000,\n",
    "                           seed=0)\n",
    "print(p_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "finished-willow",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rooms               5636.720645\n",
       "year                  17.900000\n",
       "rain                   5.025000\n",
       "Bioindicator          58.055556\n",
       "Coliform             289.583333\n",
       "total               1490.138889\n",
       "week                   4.600000\n",
       "Bioindicator UV+       9.722222\n",
       "dtype: float64"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.query('festival != \"during\"').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "777c9d6c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16vs17 = 9.999000099990002e-05, 16vs20 = 0.007099290070992901, 17vs20 = 9.999000099990002e-05 on Bioindicator\n",
      "16vs17 = 0.7629237076292371, 16vs20 = 9.999000099990002e-05, 17vs20 = 9.999000099990002e-05 on Coliform\n"
     ]
    }
   ],
   "source": [
    "for c in ['Bioindicator','Coliform']:\n",
    "    p1 = permutation_test( A.loc[A['year'] == 16 , c] , \n",
    "                                A.loc[A['year'] == 17, c] , \n",
    "                               method='approximate',\n",
    "                               num_rounds=10000,\n",
    "                               seed=0)\n",
    "    p2 = permutation_test( A.loc[A['year'] == 16 , c] , \n",
    "                                A.loc[A['year'] == 20, c] , \n",
    "                               method='approximate',\n",
    "                               num_rounds=10000,\n",
    "                               seed=0)\n",
    "    p3 = permutation_test( A.loc[A['year'] == 17 , c] , \n",
    "                                A.loc[A['year'] == 20, c] , \n",
    "                               method='approximate',\n",
    "                               num_rounds=10000,\n",
    "                               seed=0)\n",
    "    print(f\"16vs17 = {p1}, 16vs20 = {p2}, 17vs20 = {p3} on {c}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "virtual-bracelet",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(124.30555555555556, 108.68055555555556, 644.4444444444445)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.query('year == 16')['Coliform'].mean(), A.query('year == 17')['Coliform'].mean(), A.query('year == 20')['Coliform'].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "thrown-canyon",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18      0.0\n",
       "19      0.0\n",
       "20      0.0\n",
       "21    500.0\n",
       "22      0.0\n",
       "23    250.0\n",
       "24      0.0\n",
       "25      0.0\n",
       "26    250.0\n",
       "27    550.0\n",
       "28    425.0\n",
       "29    300.0\n",
       "30    800.0\n",
       "31    900.0\n",
       "32    825.0\n",
       "33    125.0\n",
       "34     50.0\n",
       "35    175.0\n",
       "Name: Bioindicator, dtype: float64"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.loc[(A['year']==16) & (A['festival'] == 'during'), 'Bioindicator']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "lyric-convergence",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "27    250.0\n",
       "28    250.0\n",
       "29    275.0\n",
       "30     25.0\n",
       "31    175.0\n",
       "32     75.0\n",
       "33    325.0\n",
       "34    200.0\n",
       "35    275.0\n",
       "36    300.0\n",
       "37    200.0\n",
       "38    200.0\n",
       "39      0.0\n",
       "40      0.0\n",
       "41      0.0\n",
       "42      0.0\n",
       "43      0.0\n",
       "44      0.0\n",
       "Name: Bioindicator, dtype: float64"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.loc[(A['year']==17) & (A['festival'] == 'during'), 'Bioindicator']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "portuguese-enemy",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "27    250.0\n",
       "28    250.0\n",
       "29    275.0\n",
       "30     25.0\n",
       "31    175.0\n",
       "32     75.0\n",
       "33    325.0\n",
       "34    200.0\n",
       "35    275.0\n",
       "36    300.0\n",
       "37    200.0\n",
       "38    200.0\n",
       "39      0.0\n",
       "40      0.0\n",
       "41      0.0\n",
       "42      0.0\n",
       "43      0.0\n",
       "44      0.0\n",
       "Name: Bioindicator UV+, dtype: float64"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.loc[(A['year']==17) & (A['festival'] == 'during'), 'Bioindicator UV+']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "3cdb3d82",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.35326467353264673\n"
     ]
    }
   ],
   "source": [
    "p_value = permutation_test( A.loc[(A['year']==16) & (A['festival'] == 'during'), 'Bioindicator'] , \n",
    "                            A.loc[(A['year']==17) & (A['festival'] == 'during'), 'Bioindicator'], \n",
    "                           method='approximate',\n",
    "                           num_rounds=10000,\n",
    "                           seed=0)\n",
    "print(p_value)\n",
    "#0.35326467353264673 BUV+"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "46204c24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0387961203879612\n"
     ]
    }
   ],
   "source": [
    "p_value = permutation_test( A.loc[(A['year']==16) & (A['festival'] == 'during'), 'Bioindicator'] , \n",
    "                            A.loc[(A['year']==16) & (A['festival'] != 'during'), 'Bioindicator'], \n",
    "                           method='approximate',\n",
    "                           num_rounds=10000,\n",
    "                           seed=0)\n",
    "print(p_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "35243d23",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0014998500149985001\n"
     ]
    }
   ],
   "source": [
    "p_value = permutation_test( A.loc[(A['year']==17) & (A['festival'] == 'during'), 'Bioindicator'] , \n",
    "                            A.loc[(A['year']==17) & (A['festival'] != 'during'), 'Bioindicator'], \n",
    "                           method='approximate',\n",
    "                           num_rounds=10000,\n",
    "                           seed=0)\n",
    "print(p_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "071991f6",
   "metadata": {},
   "outputs": [
    {
     "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>Bioindicator</th>\n",
       "      <th>Coliform</th>\n",
       "      <th>Bioindicator UV+</th>\n",
       "      <th>total</th>\n",
       "      <th>rain</th>\n",
       "      <th>rooms</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bioindicator</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.300159</td>\n",
       "      <td>0.834559</td>\n",
       "      <td>0.360940</td>\n",
       "      <td>0.123900</td>\n",
       "      <td>0.287612</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Coliform</th>\n",
       "      <td>0.300159</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.463321</td>\n",
       "      <td>0.785548</td>\n",
       "      <td>-0.029285</td>\n",
       "      <td>0.205608</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bioindicator UV+</th>\n",
       "      <td>0.834559</td>\n",
       "      <td>0.463321</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.582578</td>\n",
       "      <td>-0.183682</td>\n",
       "      <td>0.226850</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>0.360940</td>\n",
       "      <td>0.785548</td>\n",
       "      <td>0.582578</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.053015</td>\n",
       "      <td>0.205720</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rain</th>\n",
       "      <td>0.123900</td>\n",
       "      <td>-0.029285</td>\n",
       "      <td>-0.183682</td>\n",
       "      <td>-0.053015</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.399218</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rooms</th>\n",
       "      <td>0.287612</td>\n",
       "      <td>0.205608</td>\n",
       "      <td>0.226850</td>\n",
       "      <td>0.205720</td>\n",
       "      <td>0.399218</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  Bioindicator  Coliform  Bioindicator UV+     total  \\\n",
       "Bioindicator          1.000000  0.300159          0.834559  0.360940   \n",
       "Coliform              0.300159  1.000000          0.463321  0.785548   \n",
       "Bioindicator UV+      0.834559  0.463321          1.000000  0.582578   \n",
       "total                 0.360940  0.785548          0.582578  1.000000   \n",
       "rain                  0.123900 -0.029285         -0.183682 -0.053015   \n",
       "rooms                 0.287612  0.205608          0.226850  0.205720   \n",
       "\n",
       "                      rain     rooms  \n",
       "Bioindicator      0.123900  0.287612  \n",
       "Coliform         -0.029285  0.205608  \n",
       "Bioindicator UV+ -0.183682  0.226850  \n",
       "total            -0.053015  0.205720  \n",
       "rain              1.000000  0.399218  \n",
       "rooms             0.399218  1.000000  "
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "main_vars=['Bioindicator', 'Coliform', 'Bioindicator UV+',  'total' ,'rain', 'rooms']\n",
    "A[main_vars].corr(method ='spearman')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "0ac01a68",
   "metadata": {},
   "outputs": [],
   "source": [
    "def append_vars(ds, cols):\n",
    "    v=ds[cols]\n",
    "    X=ds.drop(columns=cols)\n",
    "    return pd.concat([X.assign(cfu=ds[col], name=ds['year'].apply(lambda y: f\"{col} {y}\")) for col in cols]).dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "42b959ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "A2=append_vars(A, ['Bioindicator', 'Bioindicator UV+', 'Coliform', 'total'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "6fcb603e",
   "metadata": {},
   "outputs": [
    {
     "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>Bioindicator</th>\n",
       "      <th>Coliform</th>\n",
       "      <th>Bioindicator UV+</th>\n",
       "      <th>total</th>\n",
       "      <th>rain</th>\n",
       "      <th>rooms</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bioindicator</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.463458</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.687153</td>\n",
       "      <td>0.319278</td>\n",
       "      <td>0.425990</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Coliform</th>\n",
       "      <td>0.463458</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.586674</td>\n",
       "      <td>0.000310</td>\n",
       "      <td>0.018802</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bioindicator UV+</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>0.687153</td>\n",
       "      <td>0.586674</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.083664</td>\n",
       "      <td>0.152319</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rain</th>\n",
       "      <td>0.319278</td>\n",
       "      <td>0.000310</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.083664</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.893917</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rooms</th>\n",
       "      <td>0.425990</td>\n",
       "      <td>0.018802</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.152319</td>\n",
       "      <td>0.893917</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  Bioindicator  Coliform  Bioindicator UV+     total  \\\n",
       "Bioindicator          1.000000  0.463458               NaN  0.687153   \n",
       "Coliform              0.463458  1.000000               NaN  0.586674   \n",
       "Bioindicator UV+           NaN       NaN               NaN       NaN   \n",
       "total                 0.687153  0.586674               NaN  1.000000   \n",
       "rain                  0.319278  0.000310               NaN  0.083664   \n",
       "rooms                 0.425990  0.018802               NaN  0.152319   \n",
       "\n",
       "                      rain     rooms  \n",
       "Bioindicator      0.319278  0.425990  \n",
       "Coliform          0.000310  0.018802  \n",
       "Bioindicator UV+       NaN       NaN  \n",
       "total             0.083664  0.152319  \n",
       "rain              1.000000  0.893917  \n",
       "rooms             0.893917  1.000000  "
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.query('year==16')[main_vars].corr(method ='spearman')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "196f2f38",
   "metadata": {},
   "outputs": [
    {
     "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>Bioindicator</th>\n",
       "      <th>Coliform</th>\n",
       "      <th>Bioindicator UV+</th>\n",
       "      <th>total</th>\n",
       "      <th>rain</th>\n",
       "      <th>rooms</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bioindicator</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.630557</td>\n",
       "      <td>0.834559</td>\n",
       "      <td>0.609393</td>\n",
       "      <td>-0.026901</td>\n",
       "      <td>0.250758</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Coliform</th>\n",
       "      <td>0.630557</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.463321</td>\n",
       "      <td>0.782111</td>\n",
       "      <td>-0.219374</td>\n",
       "      <td>0.439946</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bioindicator UV+</th>\n",
       "      <td>0.834559</td>\n",
       "      <td>0.463321</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.582578</td>\n",
       "      <td>-0.183682</td>\n",
       "      <td>0.226850</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>0.609393</td>\n",
       "      <td>0.782111</td>\n",
       "      <td>0.582578</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.294689</td>\n",
       "      <td>0.562313</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rain</th>\n",
       "      <td>-0.026901</td>\n",
       "      <td>-0.219374</td>\n",
       "      <td>-0.183682</td>\n",
       "      <td>-0.294689</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.013363</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rooms</th>\n",
       "      <td>0.250758</td>\n",
       "      <td>0.439946</td>\n",
       "      <td>0.226850</td>\n",
       "      <td>0.562313</td>\n",
       "      <td>-0.013363</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  Bioindicator  Coliform  Bioindicator UV+     total  \\\n",
       "Bioindicator          1.000000  0.630557          0.834559  0.609393   \n",
       "Coliform              0.630557  1.000000          0.463321  0.782111   \n",
       "Bioindicator UV+      0.834559  0.463321          1.000000  0.582578   \n",
       "total                 0.609393  0.782111          0.582578  1.000000   \n",
       "rain                 -0.026901 -0.219374         -0.183682 -0.294689   \n",
       "rooms                 0.250758  0.439946          0.226850  0.562313   \n",
       "\n",
       "                      rain     rooms  \n",
       "Bioindicator     -0.026901  0.250758  \n",
       "Coliform         -0.219374  0.439946  \n",
       "Bioindicator UV+ -0.183682  0.226850  \n",
       "total            -0.294689  0.562313  \n",
       "rain              1.000000 -0.013363  \n",
       "rooms            -0.013363  1.000000  "
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.query('year==17')[main_vars].corr(method ='spearman')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "55a1cfdf",
   "metadata": {},
   "outputs": [
    {
     "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>Bioindicator</th>\n",
       "      <th>Coliform</th>\n",
       "      <th>Bioindicator UV+</th>\n",
       "      <th>total</th>\n",
       "      <th>rain</th>\n",
       "      <th>rooms</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bioindicator</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.282178</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.313524</td>\n",
       "      <td>0.121148</td>\n",
       "      <td>0.261842</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Coliform</th>\n",
       "      <td>0.282178</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.808145</td>\n",
       "      <td>0.159472</td>\n",
       "      <td>0.508954</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bioindicator UV+</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>0.313524</td>\n",
       "      <td>0.808145</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.292094</td>\n",
       "      <td>0.668511</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rain</th>\n",
       "      <td>0.121148</td>\n",
       "      <td>0.159472</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.292094</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.579178</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rooms</th>\n",
       "      <td>0.261842</td>\n",
       "      <td>0.508954</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.668511</td>\n",
       "      <td>0.579178</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  Bioindicator  Coliform  Bioindicator UV+     total  \\\n",
       "Bioindicator          1.000000  0.282178               NaN  0.313524   \n",
       "Coliform              0.282178  1.000000               NaN  0.808145   \n",
       "Bioindicator UV+           NaN       NaN               NaN       NaN   \n",
       "total                 0.313524  0.808145               NaN  1.000000   \n",
       "rain                  0.121148  0.159472               NaN  0.292094   \n",
       "rooms                 0.261842  0.508954               NaN  0.668511   \n",
       "\n",
       "                      rain     rooms  \n",
       "Bioindicator      0.121148  0.261842  \n",
       "Coliform          0.159472  0.508954  \n",
       "Bioindicator UV+       NaN       NaN  \n",
       "total             0.292094  0.668511  \n",
       "rain              1.000000  0.579178  \n",
       "rooms             0.579178  1.000000  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.query('year==20')[main_vars].corr(method ='spearman')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ca2a01b0",
   "metadata": {},
   "outputs": [
    {
     "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>Bioindicator</th>\n",
       "      <th>Coliform</th>\n",
       "      <th>Bioindicator UV+</th>\n",
       "      <th>total</th>\n",
       "      <th>rain</th>\n",
       "      <th>rooms</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bioindicator</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.413961</td>\n",
       "      <td>0.798730</td>\n",
       "      <td>0.479667</td>\n",
       "      <td>0.108622</td>\n",
       "      <td>0.321153</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Coliform</th>\n",
       "      <td>0.413961</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.332692</td>\n",
       "      <td>0.828321</td>\n",
       "      <td>-0.048135</td>\n",
       "      <td>0.240064</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bioindicator UV+</th>\n",
       "      <td>0.798730</td>\n",
       "      <td>0.332692</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.579050</td>\n",
       "      <td>-0.160397</td>\n",
       "      <td>0.228377</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>0.479667</td>\n",
       "      <td>0.828321</td>\n",
       "      <td>0.579050</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.074906</td>\n",
       "      <td>0.223255</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rain</th>\n",
       "      <td>0.108622</td>\n",
       "      <td>-0.048135</td>\n",
       "      <td>-0.160397</td>\n",
       "      <td>-0.074906</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.399218</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rooms</th>\n",
       "      <td>0.321153</td>\n",
       "      <td>0.240064</td>\n",
       "      <td>0.228377</td>\n",
       "      <td>0.223255</td>\n",
       "      <td>0.399218</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  Bioindicator  Coliform  Bioindicator UV+     total  \\\n",
       "Bioindicator          1.000000  0.413961          0.798730  0.479667   \n",
       "Coliform              0.413961  1.000000          0.332692  0.828321   \n",
       "Bioindicator UV+      0.798730  0.332692          1.000000  0.579050   \n",
       "total                 0.479667  0.828321          0.579050  1.000000   \n",
       "rain                  0.108622 -0.048135         -0.160397 -0.074906   \n",
       "rooms                 0.321153  0.240064          0.228377  0.223255   \n",
       "\n",
       "                      rain     rooms  \n",
       "Bioindicator      0.108622  0.321153  \n",
       "Coliform         -0.048135  0.240064  \n",
       "Bioindicator UV+ -0.160397  0.228377  \n",
       "total            -0.074906  0.223255  \n",
       "rain              1.000000  0.399218  \n",
       "rooms             0.399218  1.000000  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[main_vars].corr(method ='spearman') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f32282f0",
   "metadata": {},
   "outputs": [
    {
     "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>cfu</th>\n",
       "      <th>rooms</th>\n",
       "      <th>rain</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>cfu</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.221528</td>\n",
       "      <td>-0.035081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rooms</th>\n",
       "      <td>0.221528</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.345761</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rain</th>\n",
       "      <td>-0.035081</td>\n",
       "      <td>0.345761</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            cfu     rooms      rain\n",
       "cfu    1.000000  0.221528 -0.035081\n",
       "rooms  0.221528  1.000000  0.345761\n",
       "rain  -0.035081  0.345761  1.000000"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A2[['cfu', 'rooms', 'rain']].corr(method ='spearman')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "863f7146",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "9452e227",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bootstrap_ci(df, variable, classes, repetitions = 1000, alpha = 0.05, random_state=None): \n",
    "    \n",
    "    df = df[[variable, classes]]\n",
    "    bootstrap_sample_size = len(df) \n",
    "    \n",
    "    mean_diffs = []\n",
    "    for i in range(repetitions):\n",
    "            bootstrap_sample = df.sample(n = bootstrap_sample_size, replace = True, random_state = random_state)\n",
    "            mean_diff = bootstrap_sample.groupby(classes).mean().iloc[1,0] - bootstrap_sample.groupby(classes).mean().iloc[0,0]\n",
    "            mean_diffs.append(mean_diff)\n",
    "    # confidence interval\n",
    "    left = np.percentile(mean_diffs, alpha/2*100)\n",
    "    right = np.percentile(mean_diffs, 100-alpha/2*100)\n",
    "    # point estimate\n",
    "    point_est = df.groupby(classes).mean().iloc[1,0] - df.groupby(classes).mean().iloc[0,0]\n",
    "    print('Point estimate of difference between means:', round(point_est,2))\n",
    "    print((1-alpha)*100,'%','confidence interval for the difference between means:', (round(left,2), round(right,2)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "hazardous-fashion",
   "metadata": {},
   "outputs": [
    {
     "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>Bioindicator</th>\n",
       "      <th>rain</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>200.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67</th>\n",
       "      <td>100.0</td>\n",
       "      <td>1.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>68</th>\n",
       "      <td>300.0</td>\n",
       "      <td>1.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>69</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>71</th>\n",
       "      <td>200.0</td>\n",
       "      <td>1.8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>216 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    Bioindicator  rain\n",
       "0          200.0   0.0\n",
       "1            0.0   0.0\n",
       "2            0.0   0.0\n",
       "3            0.0   0.0\n",
       "4            0.0   0.0\n",
       "..           ...   ...\n",
       "67         100.0   1.8\n",
       "68         300.0   1.8\n",
       "69           0.0   1.8\n",
       "70           0.0   1.8\n",
       "71         200.0   1.8\n",
       "\n",
       "[216 rows x 2 columns]"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[['Bioindicator','rain']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "8e038d1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import arch.data.frenchdata\n",
    "from arch.bootstrap import StationaryBootstrap,IndependentSamplesBootstrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "sought-integrity",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXDUlEQVR4nO3df5TddX3n8ee7WF2WQQLFvU0D7ZRzoj1C3NRMqefUtTPFbhG7oK2HyiJNlHa0a3f/aM5ZU/VUtx7Poitat1ptFA7QugwuSEkjusUss+iusSYeJCCKYOOWNE0QQsJgyprw3j/uN3oZ7sz9NffON588H+fcM9/7/X7u/b5mcueV73zv9/u9kZlIksryY8sdQJK09Cx3SSqQ5S5JBbLcJalAlrskFchyl6QCWe6SVCDLXapExMsj4v9ExMGIeCwi/ndE/EJEbIiIoxExFxGHIuLrEfHry51XWozlLgER8XxgK/CnwBnAKuA/AU9VQ76cmWPACuDPgJmIWDH6pFJ3LHep6YUAmXljZh7NzMOZ+TeZeU/roMx8GvgL4BRg9TLklLpiuUtNDwBHI+L6iHhVRJzeblBEnAS8EfgB8N1RBpR6YblLQGYeAl4OJPAJ4JGI2BIRjWrIyyLiceCfgA8Ab8jM/csSVuqC5S5VMvP+zNyQmWcB5wE/BfxJtXh7Zq4ATge2AP9qWUJKXbLcpTYy85vAdTRLvnX+HPB7wBUR8fPLEE3qiuUuARHxcxGxMSLOqu6fDVwGbJ8/NjMfAz4J/NFoU0rds9ylpieAXwS+EhFP0iz1e4GNC4z/E+CiiHjJaOJJvQk/rEOSyuOWuyQVyHKXpAJZ7pJUIMtdkgr0nE4DqkPCbgAaNM/e25yZH46IM4CbgHFgN3BpZh6IiAA+DFwEfB/YkJlfW2wdZ555Zo6Pjw/wbfTuySef5JRTThnpOntR53xm60+ds0G985mtvZ07d34vM1/QdmFmLnoDVgIvraZPpXkNjhcD7wc2VfM3Ae+rpi8CPgcE8DLgK53WsW7duhy1O++8c+Tr7EWd85mtP3XOllnvfGZrD9iRC/Rqx90ymbk3qy3vzHwCuJ/m5VAvAa6vhl0PvKaavgS4oVr3dmBFRKzs/v8iSdKgejrOPSLGgbtonpL9f7N5rQ2qXTEHMnNFRGwFrsrML1XLtgFvy8wd855rGpgGaDQa62ZmZgb/bnowNzfH2NjYSNfZizrnM1t/6pwN6p3PbO1NTU3tzMyJtgsX2qSffwPGgJ3Ab1T3H5+3/ED1dSvw8pb524CJxZ7b3TLPVud8ZutPnbNl1juf2dpjkN0yABHx48AtwKcy8zPV7H3HdrdUX49d/nQPcHbLw8+q5kmSRqRjuVe7XK4B7s/MD7Ys2gKsr6bXA7e1zP/taHoZcDAz9y5hZklSBx0PhQR+CbgC2BURd1fz3g5cBXw6Iq6k+Yk0l1bLbqd5xMyDNA+FfONSBpYkddax3LP5xmgssPiCNuMTeOuAuSRJA/AMVUkqkOUuSQXqZp+7dEIb3/TZBZftvurVI0widc8td0kqkOUuSQWy3CWpQJa7JBXIcpekAlnuklQgD4XUCcHDGXWicctdkgpkuUtSgSx3SSqQ5S5JBbLcJalAlrskFchyl6QCWe6SVKBuPiD72ojYHxH3tsy7KSLurm67j322akSMR8ThlmUfH2J2SdICujlD9TrgI8ANx2Zk5m8dm46Iq4GDLeMfysy1S5RPktSHbj4g+66IGG+3LCICuBT4lSXOJUkaQGRm50HNct+amefNm/8K4IOZOdEy7j7gAeAQ8M7M/OICzzkNTAM0Go11MzMz/X8XfZibm2NsbGyk6+xFnfMdj9l27TnYZnTTmlWnLfqcgzy2m2x1Ued8Zmtvampq57H+nW/QC4ddBtzYcn8v8NOZ+WhErAP+KiLOzcxD8x+YmZuBzQATExM5OTk5YJTezM7OMup19qLO+Y7HbBsWu3DY5c8ev1SPbVXnnxvUO5/Zetf30TIR8RzgN4Cbjs3LzKcy89FqeifwEPDCQUNKknozyJb7K4FvZubDx2ZExAuAxzLzaEScA6wGvjNgRmnRS/aCl+2V5uvmUMgbgS8DL4qIhyPiymrR63nmLhmAVwD3VIdG3gy8JTMfW8K8kqQudHO0zGULzN/QZt4twC2Dx5IkDcIzVCWpQJa7JBXIcpekAlnuklQgy12SCmS5S1KBLHdJKpDlLkkFstwlqUCDXhVSqoVj157ZuObIoldxXOyxUkksd9WGJSstHXfLSFKBLHdJKpDlLkkFstwlqUCWuyQVyHKXpAJZ7pJUoG4+Q/XaiNgfEfe2zHt3ROyJiLur20Uty/4wIh6MiG9FxK8NK7gkaWHdbLlfB1zYZv6HMnNtdbsdICJeTPODs8+tHvNnEXHSUoWVJHWnY7ln5l3AY10+3yXATGY+lZl/BzwInD9APklSHyIzOw+KGAe2ZuZ51f13AxuAQ8AOYGNmHoiIjwDbM/Mvq3HXAJ/LzJvbPOc0MA3QaDTWzczMLMX307W5uTnGxsZGus5e1DnfsLLt2nNw4OdonAz7Di9BmC6tWXVa12Pr/G8K9c5ntvampqZ2ZuZEu2X9XlvmY8B7gKy+Xg28qZcnyMzNwGaAiYmJnJyc7DNKf2ZnZxn1OntR53zDytbrBb/a2bjmCFfvGt0lk3ZfPtn12Dr/m0K985mtd30dLZOZ+zLzaGY+DXyCH+162QOc3TL0rGqeJGmE+ir3iFjZcve1wLEjabYAr4+I50XEzwKrgb8dLKIkqVcd/36NiBuBSeDMiHgYeBcwGRFrae6W2Q28GSAz74uITwPfAI4Ab83Mo0NJLklaUMdyz8zL2sy+ZpHx7wXeO0go6XjR6Rr0u6969YiSSM/kGaqSVCDLXZIKZLlLUoEsd0kqkOUuSQWy3CWpQJa7JBXIcpekAlnuklQgy12SCjS6a6PqhNfpVH1JS8ctd0kqkOUuSQVyt4w0RK27ojauOfKsT5vyqpEaFrfcJalAlrskFchyl6QCWe6SVKCO5R4R10bE/oi4t2Xef4mIb0bEPRFxa0SsqOaPR8ThiLi7un18iNklSQvoZsv9OuDCefPuAM7LzJcADwB/2LLsocxcW93esjQxJUm96FjumXkX8Ni8eX+TmUequ9uBs4aQTZLUp8jMzoMixoGtmXlem2V/DdyUmX9ZjbuP5tb8IeCdmfnFBZ5zGpgGaDQa62ZmZvr9HvoyNzfH2NjYSNfZizrn6zfbrj0Hh5DmmRonw77DQ19NX9plW7PqtOUJ00aJr7lRWM5sU1NTOzNzot2ygU5iioh3AEeAT1Wz9gI/nZmPRsQ64K8i4tzMPDT/sZm5GdgMMDExkZOTk4NE6dns7CyjXmcv6pyv32zzT+AZho1rjnD1rnqem9cu2+7LJ5cnTBslvuZGoa7Z+j5aJiI2AL8OXJ7V5n9mPpWZj1bTO4GHgBcuQU5JUg/6KveIuBD4j8DFmfn9lvkviIiTqulzgNXAd5YiqCSpex3/fo2IG4FJ4MyIeBh4F82jY54H3BERANurI2NeAfxxRPwAeBp4S2Y+1vaJJS16GWSvO6NBdCz3zLyszexrFhh7C3DLoKEkSYPxDFVJKpDlLkkFstwlqUCWuyQVyHKXpAJZ7pJUIMtdkgpkuUtSgSx3SSpQPS+fp+PSYqfSSxott9wlqUCWuyQVyHKXpAJZ7pJUIN9QVU/GN32WjWuOjOQj8yT1zy13SSqQ5S5JBeqq3CPi2ojYHxH3tsw7IyLuiIhvV19Pr+ZHRPzXiHgwIu6JiJcOK7wkqb1ut9yvAy6cN28TsC0zVwPbqvsAr6L5wdirgWngY4PHlCT1oqtyz8y7gPkfdH0JcH01fT3wmpb5N2TTdmBFRKxcgqySpC4Nss+9kZl7q+l/BBrV9Crg71vGPVzNkySNSGRmdwMjxoGtmXledf/xzFzRsvxAZp4eEVuBqzLzS9X8bcDbMnPHvOebprnbhkajsW5mZmYJvp3uzc3NMTY2NtJ19qKu+XbtOUjjZNh3eLmTtFdStjWrThtemDbq+poDsy1kampqZ2ZOtFs2yHHu+yJiZWburXa77K/m7wHObhl3VjXvGTJzM7AZYGJiIicnJweI0rvZ2VlGvc5e1DXfhuo496t31fMUiZKy7b58cnhh2qjraw7M1o9BdstsAdZX0+uB21rm/3Z11MzLgIMtu28kSSPQ1WZERNwITAJnRsTDwLuAq4BPR8SVwHeBS6vhtwMXAQ8C3wfeuMSZJUkddFXumXnZAosuaDM2gbcOEkqSNBjPUJWkAlnuklQgy12SCmS5S1KBLHdJKlA9z/aQxHiHD0TZfdWrR5RExyO33CWpQJa7JBXIcpekAlnuklQgy12SCmS5S1KBLHdJKpDHuUvHqcWOg/cYeLnlLkkFstwlqUCWuyQVyHKXpAL1/YZqRLwIuKll1jnAHwErgN8FHqnmvz0zb+93PZKk3vVd7pn5LWAtQEScBOwBbqX5gdgfyswPLEVASVLvlupQyAuAhzLzuxGxRE+p5dDpMrOSjg+RmYM/ScS1wNcy8yMR8W5gA3AI2AFszMwDbR4zDUwDNBqNdTMzMwPn6MXc3BxjY2MjXWcvlivfrj0HO45pnAz7Do8gTB/M1rRm1Wk9P6bOvxNma29qampnZk60WzZwuUfEc4F/AM7NzH0R0QC+ByTwHmBlZr5pseeYmJjIHTt2DJSjV7Ozs0xOTo50nb1YrnzdbLlvXHOEq3fV8/w3szX1cxJTnX8nzNZeRCxY7ktxtMyraG617wPIzH2ZeTQznwY+AZy/BOuQJPVgKcr9MuDGY3ciYmXLstcC9y7BOiRJPRjob8SIOAX4VeDNLbPfHxFrae6W2T1vmSRpBAYq98x8EviJefOuGCiRhs4jYqTyeYaqJBXIcpekAlnuklQgy12SCmS5S1KBLHdJKpDlLkkFstwlqUCWuyQVqJ6Xz5M0kE5nIfdz1UgdX9xyl6QCWe6SVCDLXZIKZLlLUoEsd0kqkOUuSQWy3CWpQJa7JBVo4JOYImI38ARwFDiSmRMRcQZwEzBO83NUL83MA4OuS5LUnaXacp/KzLWZOVHd3wRsy8zVwLbqviRpRIZ1+YFLgMlq+npgFnjbkNZ1wvHUckmdRGYO9gQRfwccABL488zcHBGPZ+aKankAB47db3ncNDAN0Gg01s3MzAyUo1dzc3OMjY2NdJ29WCzfrj0HF33smlWnLbq80+M7aZwM+w4P9BRDY7butHuN1Pl3wmztTU1N7WzZY/IMS1HuqzJzT0T8C+AO4N8DW1rLPCIOZObpCz3HxMRE7tixY6AcvZqdnWVycnKk6+zFYvkG3XLv9PhONq45wtW76nnNObN1p91rpM6/E2ZrLyIWLPeBX2mZuaf6uj8ibgXOB/ZFxMrM3BsRK4H9g65H3Ru0vCUd/wZ6QzUiTomIU49NA/8auBfYAqyvhq0HbhtkPZKk3gy65d4Abm3uVuc5wH/LzM9HxFeBT0fElcB3gUsHXI8kqQcDlXtmfgf4l23mPwpcMMhzS5L65xmqklQgy12SCmS5S1KBLHdJKpDlLkkFstwlqUCWuyQVyHKXpALV4ypGepZdew6ywWvESOqT5S6dgNpdXG7jmiM/3KDwMwGOf+6WkaQCWe6SVCDLXZIKZLlLUoEsd0kqkOUuSQWy3CWpQJa7JBWo75OYIuJs4Aaan6OawObM/HBEvBv4XeCRaujbM/P2QYNKGp12Jzkd4wlOx4dBzlA9AmzMzK9FxKnAzoi4o1r2ocz8wODxJEn96LvcM3MvsLeafiIi7gdWLVUwSVL/IjMHf5KIceAu4DzgD4ANwCFgB82t+wNtHjMNTAM0Go11MzMzA+foxdzcHGNjYyNdZy/2P3aQfYeXO0V7jZMxWx/qnA2WLt+aVacN/iTz1Pn3dTmzTU1N7czMiXbLBi73iBgD/hfw3sz8TEQ0gO/R3A//HmBlZr5pseeYmJjIHTt2DJSjV7Ozs0xOTo50nb3400/dxtW76nldt41rjpitD3XOBqPJ1+/++jr/vi5ntohYsNwHOlomIn4cuAX4VGZ+BiAz92Xm0cx8GvgEcP4g65Ak9a7vco+IAK4B7s/MD7bMX9ky7LXAvf3HkyT1Y5C/wX4JuALYFRF3V/PeDlwWEWtp7pbZDbx5gHVIkvowyNEyXwKizSKPaZekZeYZqpJUIMtdkgpkuUtSgSx3SSqQ5S5JBarv6XIF8Mp6kpaLW+6SVCC33JfJYlv1ABvXjCiIpCK55S5JBXLLXdLIdPqL1feilo7lLqk2Fir/jWuOsGHTZy3/HrhbRpIKZLlLUoHcLSPpuOG5I92z3BfR6c0fSaqrE77cLXBJJXKfuyQV6ITfcpdUBo+hf6ahbblHxIUR8a2IeDAiNg1rPZKkZxvKlntEnAR8FPhV4GHgqxGxJTO/MYz1+Q66pE4GeX9tkB5Zrr8ohrVb5nzgwcz8DkBEzACXAEMpd0kapsUK+tjZs3UTmbn0TxrxOuDCzPyd6v4VwC9m5u+3jJkGpqu7LwK+teRBFncm8L0Rr7MXdc5ntv7UORvUO5/Z2vuZzHxBuwXL9oZqZm4GNi/X+iNiR2ZOLNf6O6lzPrP1p87ZoN75zNa7Yb2hugc4u+X+WdU8SdIIDKvcvwqsjoifjYjnAq8HtgxpXZKkeYayWyYzj0TE7wP/AzgJuDYz7xvGugawbLuEulTnfGbrT52zQb3zma1HQ3lDVZK0vLz8gCQVyHKXpAKdMOUeEWdExB0R8e3q6+ltxqyNiC9HxH0RcU9E/Fad8lXjPh8Rj0fE1hFkWvQSEhHxvIi4qVr+lYgYH3amHrK9IiK+FhFHqvMuRqaLbH8QEd+oXmPbIuJnapTtLRGxKyLujogvRcSLR5Wtm3wt434zIjIiRnYIYhc/uw0R8Uj1s7s7In5nVNnayswT4ga8H9hUTW8C3tdmzAuB1dX0TwF7gRV1yVctuwD4N8DWIec5CXgIOAd4LvB14MXzxvw74OPV9OuBm0b0s+om2zjwEuAG4HUjfJ11k20K+OfV9O/V7Of2/Jbpi4HP1+lnV407FbgL2A5M1CUbsAH4yKh+Xp1uJ8yWO83LH1xfTV8PvGb+gMx8IDO/XU3/A7AfaHv213LkA8jMbcATI8jzw0tIZOb/A45dQqJVa+abgQsiIuqQLTN3Z+Y9wNMjyNNrtjsz8/vV3e00zwOpS7ZDLXdPAUZ5xEU3rzmA9wDvA/6phtlq40Qq90Zm7q2m/xFoLDY4Is6n+T/0Q8MOVukp3wisAv6+5f7D1by2YzLzCHAQ+ImaZFsuvWa7EvjcUBP9SFfZIuKtEfEQzb8m/8OIskEX+SLipcDZmTnqi7l0++/6m9Xutpsj4uw2y0emqOu5R8QXgJ9ss+gdrXcyMyNiwS2SiFgJ/AWwPjOXbMtvqfKpDBHxBmAC+OXlztIqMz8KfDQi/i3wTmD9MkcCICJ+DPggzd0fdfTXwI2Z+VREvJnmX7W/slxhiir3zHzlQssiYl9ErMzMvVV5719g3POBzwLvyMztdcs3Qt1cQuLYmIcj4jnAacCjNcm2XLrKFhGvpPmf+i9n5lN1ytZiBvjYUBM9U6d8pwLnAbPV3r+fBLZExMWZuWOZs5GZra/9T9L8y2fZnEi7Zbbwoy2Q9cBt8wdUl0q4FbghM28eYTboIt+IdXMJidbMrwP+Z1bvLNUg23LpmC0ifh74c+DizBzlf+LdZFvdcvfVwLfrki8zD2bmmZk5npnjNN+vGEWxd8wGP/yL/5iLgftHkGthy/2O7qhuNPcFb6P5Yv0CcEY1fwL4ZDX9BuAHwN0tt7V1yVfd/yLwCHCY5n6/XxtipouAB2i+7/COat4f0/yFAvhnwH8HHgT+FjhnhP+enbL9QvXzeZLmXxP31SjbF4B9La+xLTXK9mHgvirXncC5o8rWTb55Y2cZ0dEyXf7s/nP1s/t69bP7uVH+7ObfvPyAJBXoRNotI0knDMtdkgpkuUtSgSx3SSqQ5S5JBbLcJalAlrskFej/A/1j7le7fdlNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def foi(X):\n",
    "    return X.corr(method ='spearman').iloc[0,1] #['cfu']['rain']\n",
    "    \n",
    "bs = StationaryBootstrap(12, A[['Bioindicator','rain']])\n",
    "results = bs.apply(foi, 2500)\n",
    "SR = pd.DataFrame(results[:, -1:], columns=[\"SR\"])\n",
    "fig = SR.hist(bins=40)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "christian-static",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.12390000604788295"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[['Bioindicator','rain']].corr(method='spearman').iloc[0,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "prerequisite-investor",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4904"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l = SR.SR\n",
    "\n",
    "total_size = len(SR.SR)\n",
    "true = np.sum(l <0.12)\n",
    "\n",
    "true/total_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "4c090135",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAU9ElEQVR4nO3df5DkdX3n8ef7NupxjGEhYGdd8Ear0JSwZuNOiHUxVs9pLgjWoVeJB8VxbORupS7WJZVNVTCm1MSiisuJWommklUoMBoG65C4x3GXECtz5E5Qdz1kF0ECuom7bO0GWHYds0ey6zt/zHeTZpiZ/vH99nT3h+ejqmu6P99vf/r9ZnpffOfb3/5+IzORJJXln4y6AElS8wx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXapExBsj4ksRcTQino6I/xsRPx4RWyPiZEQsRMSxiPh6RLxt1PVKqzHcJSAifhC4C/gd4CxgI/AbwLPVKvdl5hSwHvhdYC4i1q99pVJvDHdp0asBMvO2zDyZmccz808y88HOlTLz+8AfAKcD54+gTqknhru06FHgZETcGhFvjYgzl1spItYBPw/8HfCXa1mg1A/DXQIy8xjwRiCBTwJ/HRE7I6JVrfKGiHgG+P/Ah4F/l5mHR1Ks1APDXapk5sOZuTUzzwUuBF4OfKxafH9mrgfOBHYCPzWSIqUeGe7SMjLzEeAWFkO+c3wB+E/AVRHxYyMoTeqJ4S4BEfEjEbE9Is6tHp8HXAHcv3TdzHwK+BTw/rWtUuqd4S4t+i7wE8CXI+J7LIb6XmD7Cut/DLgkIl63NuVJ/Qkv1iFJ5XHLXZIKZLhLUoEMd0kqkOEuSQX6gVEXAHD22Wfn9PT00Ob/3ve+x+mnnz60+UfFviaLfU2ece9t9+7dT2bmOcstG4twn56eZteuXUObf35+nna7PbT5R8W+Jot9TZ5x7y0iVjy/kbtlJKlAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQGPxDVVplKav+x+rLt93w6VrVInUHLfcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBWoa7hHxM0RcTgi9naM3R4RD1S3fRHxQDU+HRHHO5b93hBrlyStoJdzy9wCfBz49KmBzPy3p+5HxI3A0Y71H8/MzQ3VJ0kaQNdwz8x7I2J6uWUREcA7gX/ZcF1So7qdHEwqTWRm95UWw/2uzLxwyfibgI9k5kzHeg8BjwLHgF/PzD9fYc5twDaAVqu1ZW5ubvAuulhYWGBqampo84+KffVuz4Gj3VdawaaNZzRSg7+vyTPuvc3Ozu4+lb9L1T3l7xXAbR2PDwKvyMynImIL8EcRcUFmHlv6xMzcAewAmJmZyXa7XbOUlc3PzzPM+UfFvnq3tcaW+74r243U4O9r8kxybwMfLRMRPwD8G+D2U2OZ+WxmPlXd3w08Dry6bpGSpP7UORTyLcAjmbn/1EBEnBMR66r7rwLOB75Vr0RJUr96ORTyNuA+4DURsT8irqkWXc5zd8kAvAl4MCK+Dvw34NrMfLrJgiVJ3fVytMwVK4xvXWbsDuCO+mVJkurwG6qSVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKlDdi3VIWkXn5f22bzrxvIuG7Lvh0rUuSS8QbrlLUoEMd0kqkOEuSQVyn7tUw3SNC29Lw+SWuyQVqJdrqN4cEYcjYm/H2Acj4kBEPFDdLulY9t6IeCwivhkRPzOswiVJK+tly/0W4OJlxj+amZur290AEfFaFi+cfUH1nN+NiHVNFStJ6k3XcM/Me4Gne5zvMmAuM5/NzG8DjwEX1ahPkjSAOvvc3xMRD1a7bc6sxjYC3+lYZ381JklaQ5GZ3VeKmAbuyswLq8ct4EkggQ8BGzLzXRHxCeC+zPxMtd5NwN2Zeccyc24DtgG0Wq0tc3NzzXS0jIWFBaampoY2/6jYV+/2HDg68HM3bTyjkXlbp8Gh473PPSlKfR/C+Pc2Ozu7OzNnlls20KGQmXno1P2I+CRwV/VwP3Bex6rnAk+sMMcOYAfAzMxMttvtQUrpyfz8PMOcf1Tsq3dLv/bfj31XthuZd/umE9y457n/5Fabe1KU+j6Eye5toN0yEbGh4+E7gFNH0uwELo+Il0TEK4Hzga/UK1GS1K+uW+4RcRvQBs6OiP3AB4B2RGxmcbfMPuDdAJn5UER8DvgGcAL4hcw8OZTKJUkr6hrumXnFMsM3rbL+9cD1dYqSJNXjN1QlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgbxYh9SFF+TQJDLcVQQDWHoud8tIUoEMd0kqkLtlNDZO7VrZvunE8862uO+GS0dRkjSx3HKXpAK55a6J4AemUn/ccpekAhnuklQgw12SCmS4S1KB/EBVGqHVPij28E/V4Za7JBWoa7hHxM0RcTgi9naM/deIeCQiHoyIOyNifTU+HRHHI+KB6vZ7Q6xdkrSCXrbcbwEuXjJ2D3BhZr4OeBR4b8eyxzNzc3W7tpkyJUn96LrPPTPvjYjpJWN/0vHwfuBnG65LBfKLSNLaiczsvtJiuN+VmRcus+y/A7dn5meq9R5icWv+GPDrmfnnK8y5DdgG0Gq1tszNzQ3aQ1cLCwtMTU0Nbf5RmbS+9hw42tN6rdPg0PEhFzMC/fa1aeMZwyumQZP2PuzHuPc2Ozu7OzNnlltW62iZiHgfcAL4bDV0EHhFZj4VEVuAP4qICzLz2NLnZuYOYAfAzMxMttvtOqWsan5+nmHOPyqT1tfSk4GtZPumE9y4p7wDufrta9+V7eEV06BJex/2Y5J7G/homYi4GngbcGVWm/+Z+WxmPlXd3w08Dry6iUIlSb0bKNwj4mLgV4F/nZl/0zF+TkSsq+6/Cjgf+FYThUqSetf1b8SIuA1oA2dHxH7gAyweHfMS4J6IALi/OjLmTcBvRsQJ4CRwbWY+PaTaJUkr6OVomSuWGb5phXXvAO6oW5Sk7kcX+Q1WrcZvqEpSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKlAv11C9GXgbcDgzL6zGzgJuB6aBfcA7M/NItey9wDUsXkP1P2fmHw+lckkr8hJ96hruwC3Ax4FPd4xdB3wxM2+IiOuqx78aEa8FLgcuAF4O/GlEvDozTzZbtsZRt0CRtHa67pbJzHuBp5cMXwbcWt2/FXh7x/hcZj6bmd8GHgMuaqZUSVKvBt3n3srMgwDVz5dV4xuB73Sst78akyStocjM7itFTAN3dexzfyYz13csP5KZZ0bEJ4D7MvMz1fhNwN2Zeccyc24DtgG0Wq0tc3NzDbSzvIWFBaampoY2/6iMW197DhxtZJ7WaXDoeCNTjZWm+9q08YwVl3X7Xaz23H6N2/uwSePe2+zs7O7MnFluWS/73JdzKCI2ZObBiNgAHK7G9wPndax3LvDEchNk5g5gB8DMzEy22+0BS+lufn6eYc4/KuPW19aG9rlv33SCG/cM+tYcX033te/K9orLuv0uVntuv8btfdikSe5t0N0yO4Grq/tXA1/oGL88Il4SEa8Ezge+Uq9ESVK/ejkU8jagDZwdEfuBDwA3AJ+LiGuAvwJ+DiAzH4qIzwHfAE4Av+CRMpK09rqGe2ZescKiN6+w/vXA9XWKkiTV4zdUJalAhrskFai8QxKkFwi/EazVuOUuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHPLSC9A3c5Ls++GS9eoEg2LW+6SVCDDXZIKZLhLUoEMd0kq0MAfqEbEa4DbO4ZeBbwfWA/8R+Cvq/Ffy8y7B30dSVL/Bg73zPwmsBkgItYBB4A7gZ8HPpqZH26iQElS/5raLfNm4PHM/MuG5pMk1RCZWX+SiJuBr2XmxyPig8BW4BiwC9iemUeWec42YBtAq9XaMjc3V7uOlSwsLDA1NTW0+Udl3Prac+BoI/O0ToNDxxuZaqxMUl+bNp7R87rj9j5s0rj3Njs7uzszZ5ZbVjvcI+LFwBPABZl5KCJawJNAAh8CNmTmu1abY2ZmJnft2lWrjtXMz8/TbreHNv+ojFtfTV2wefumE9y4p7zv101SX/18iWnc3odNGvfeImLFcG9it8xbWdxqPwSQmYcy82Rmfh/4JHBRA68hSepDE+F+BXDbqQcRsaFj2TuAvQ28hiSpD7X+RoyIfwb8NPDujuHfiojNLO6W2bdkmSRpDdQK98z8G+CHloxdVasiSVJtfkNVkgo0GR/da2w0dUSMpOFyy12SCmS4S1KBDHdJKpDhLkkFMtwlqUAeLSPpeVY7KsqLZ08Gt9wlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFcjTD+g5vBiHVIa6F8jeB3wXOAmcyMyZiDgLuB2YZvEC2e/MzCP1ypQk9aOJ3TKzmbk5M2eqx9cBX8zM84EvVo8lSWtoGPvcLwNure7fCrx9CK8hSVpFZObgT474NnAESOD3M3NHRDyTmes71jmSmWcu89xtwDaAVqu1ZW5ubuA6ullYWGBqampo84/KMPrac+Boo/MNonUaHDo+6iqa90Lpa9PGM0ZXTMPGPTtmZ2d3d+w1eY66H6j+ZGY+EREvA+6JiEd6fWJm7gB2AMzMzGS73a5Zysrm5+cZ5vyjMoy+to7BB6rbN53gxj3lfdb/Qulr35Xt0RXTsEnOjlq7ZTLziernYeBO4CLgUERsAKh+Hq5bpCSpPwOHe0ScHhEvPXUf+FfAXmAncHW12tXAF+oWKUnqT52/EVvAnRFxap4/zMz/FRFfBT4XEdcAfwX8XP0y1SSPZZfKN3C4Z+a3gB9dZvwp4M11ipIk1ePpBySpQIa7JBXIcJekApV30K38wFSSW+6SVCLDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBPLeMpEatdm6jfTdcuoaVvLC55S5JBTLcJalAdS6QfV5E/FlEPBwRD0XEL1bjH4yIAxHxQHW7pLlyJUm9qLPP/QSwPTO/FhEvBXZHxD3Vso9m5ofrlydJGkSdC2QfBA5W978bEQ8DG5sqTCs79YHV9k0n2OqFOSQtIzKz/iQR08C9wIXALwNbgWPALha37o8s85xtwDaAVqu1ZW5urnYdK1lYWGBqampo86+1PQeOAtA6DQ4dH3ExQ2Bfk6WfvjZtPGPV5afe24M+v2njnh2zs7O7M3NmuWW1wz0ipoD/DVyfmZ+PiBbwJJDAh4ANmfmu1eaYmZnJXbt21apjNfPz87Tb7aHNv9Y6t9xv3FPe0az2NVn66avboZDdLhG51odSjnt2RMSK4V7raJmIeBFwB/DZzPw8QGYeysyTmfl94JPARXVeQ5LUvzpHywRwE/BwZn6kY3xDx2rvAPYOXp4kaRB1/kb8SeAqYE9EPFCN/RpwRURsZnG3zD7g3TVeQ5I0gDpHy/wfIJZZdPfg5UgqWbd96mqO31CVpAIZ7pJUIMNdkgpkuEtSgQx3SSpQeV+XK4RHFUj9Gbdvt46aW+6SVCDDXZIKZLhLUoEMd0kqkB+oSlIXk/hhreE+Ih4NI2mY3C0jSQVyy32I3DqXmlXn39Rqzx3H3Sp1ueUuSQUy3CWpQIa7JBXIfe6SXvBW2h+/fdMJtk7oZ2duuUtSgYYW7hFxcUR8MyIei4jrhvU6kqTnG8pumYhYB3wC+GlgP/DViNiZmd8Yxut1OzxqWH9alXj4lKQyDGuf+0XAY5n5LYCImAMuA4YS7qPiceySoF4WDGsjMTKz+Ukjfha4ODP/Q/X4KuAnMvM9HetsA7ZVD18DfLPxQv7R2cCTQ5x/VOxrstjX5Bn33v55Zp6z3IJhbbnHMmPP+b9IZu4Adgzp9Z9bTMSuzJxZi9daS/Y1Wexr8kxyb8P6QHU/cF7H43OBJ4b0WpKkJYYV7l8Fzo+IV0bEi4HLgZ1Dei1J0hJD2S2TmSci4j3AHwPrgJsz86FhvFaP1mT3zwjY12Sxr8kzsb0N5QNVSdJo+Q1VSSqQ4S5JBSoy3CPirIi4JyL+ovp55jLr/NOI+EpEfD0iHoqI3xhFrf3osa/zIuLPIuLhqq9fHEWt/eilr2q9myPicETsXesa+9Ht1Bux6Ler5Q9GxOtHUWe/eujrRyLivoh4NiJ+ZRQ1DqKHvq6sfk8PRsSXIuJHR1Fn3zKzuBvwW8B11f3rgP+yzDoBTFX3XwR8GXjDqGtvoK8NwOur+y8FHgVeO+ra6/ZVLXsT8Hpg76hrXqWXdcDjwKuAFwNfX/rfH7gE+J/Ve/ANwJdHXXdDfb0M+HHgeuBXRl1zg339C+DM6v5bJ+H3lZllbrmzeKqDW6v7twJvX7pCLlqoHr6ouo37p8u99HUwM79W3f8u8DCwca0KHFDXvgAy817g6TWqaVD/cOqNzPxb4NSpNzpdBny6eg/eD6yPiA1rXWifuvaVmYcz86vA342iwAH10teXMvNI9fB+Fr+3M/ZKDfdWZh6ExbBjcYvieSJiXUQ8ABwG7snML69diQPpqa9TImIa+DEW/yoZZ331NeY2At/peLyf5//PtZd1xs0k1tyLfvu6hsW/usbexF6sIyL+FPjhZRa9r9c5MvMksDki1gN3RsSFmTnS/blN9FXNMwXcAfxSZh5rorY6muprAnQ99UaP64ybSay5Fz33FRGzLIb7G4daUUMmNtwz8y0rLYuIQxGxITMPVn/uHu4y1zMRMQ9cDIw03JvoKyJexGKwfzYzPz+kUvvS5O9rzPVy6o1JPD3HJNbci576iojXAZ8C3pqZT61RbbWUultmJ3B1df9q4AtLV4iIc6otdiLiNOAtwCNrVeCAeukrgJuAhzPzI2tYWx1d+5ogvZx6Yyfw76ujZt4AHD21W2qMlXpKka59RcQrgM8DV2XmoyOocTCj/kR3GDfgh4AvAn9R/TyrGn85cHd1/3XA/wMeZHFr/f2jrruhvt7I4p+VDwIPVLdLRl173b6qx7cBB1n8wG4/cM2oa1+hn0tYPErpceB91di1wLXV/WDxYjaPA3uAmVHX3FBfP1z9Xo4Bz1T3f3DUdTfQ16eAIx3/nnaNuuZebp5+QJIKVOpuGUl6QTPcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoH+HjDmCN8zq+yZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def foi(X):\n",
    "    return X[['cfu', 'rooms', 'rain']].corr(method ='spearman')['cfu']['rain']\n",
    "    \n",
    "bs = StationaryBootstrap(12, A2)\n",
    "results = bs.apply(foi, 2500)\n",
    "SR = pd.DataFrame(results[:, -1:], columns=[\"SR\"])\n",
    "fig = SR.hist(bins=40)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8626f555",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/kato/anaconda3/lib/python3.7/site-packages/arch/bootstrap/base.py:685: RuntimeWarning: The conditions to reuse the previous bootstrap has not been satisfied. A new bootstrap will be used.\n",
      "  warnings.warn(warn, RuntimeWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[-0.20008459],\n",
       "       [ 0.1353287 ]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ci = bs.conf_int(foi, 1000, method=\"percentile\", reuse=True)\n",
    "ci\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "consistent-shore",
   "metadata": {},
   "source": [
    "## Bootstrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5439dd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def boot_matrix(z, B):\n",
    "    \"\"\"Bootstrap sample\n",
    "    \n",
    "    Returns all bootstrap samples in a matrix\"\"\"\n",
    "    \n",
    "    n = len(z)  # sample size\n",
    "    idz = np.random.randint(0, n, size=(B, n))  # indices to pick for all boostrap samples\n",
    "    return z[idz]\n",
    "\n",
    "def bootstrap_mean(x, B=10000, alpha=0.05, plot=False):\n",
    "    \"\"\"Bootstrap standard error and (1-alpha)*100% c.i. for the population mean\n",
    "    \n",
    "    Returns bootstrapped standard error and different types of confidence intervals\"\"\"\n",
    "   \n",
    "    # Deterministic things\n",
    "    n = len(x)  # sample size\n",
    "    orig = x.mean()  # sample mean\n",
    "    se_mean = x.std()/np.sqrt(n) # standard error of the mean\n",
    "    qt = stats.t.ppf(q=1 - alpha/2, df=n - 1) # Student quantile\n",
    "    \n",
    "    # Generate boostrap distribution of sample mean\n",
    "    xboot = boot_matrix(x, B=B)\n",
    "    sampling_distribution = xboot.mean(axis=1)\n",
    "   \n",
    "   # Standard error and sample quantiles\n",
    "    se_mean_boot = sampling_distribution.std()\n",
    "    quantile_boot = np.percentile(sampling_distribution, q=(100*alpha/2, 100*(1-alpha/2)))\n",
    " \n",
    "    # RESULTS\n",
    "    print(\"Estimated mean:\", orig)\n",
    "    print(\"Classic standard error:\", se_mean)\n",
    "    print(\"Classic student c.i.:\", orig + np.array([-qt, qt])*se_mean)\n",
    "    print(\"\\nBootstrap results:\")\n",
    "    print(\"Standard error:\", se_mean_boot)\n",
    "    print(\"t-type c.i.:\", orig + np.array([-qt, qt])*se_mean_boot)\n",
    "    print(\"Percentile c.i.:\", quantile_boot)\n",
    "    print(\"Basic c.i.:\", 2*orig - quantile_boot[::-1])\n",
    "\n",
    "    if plot:\n",
    "        plt.hist(sampling_distribution, bins=\"fd\")\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b169f30e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "10e19f63",
   "metadata": {},
   "outputs": [
    {
     "ename": "StudentizationError",
     "evalue": "\nThe estimated covariance computed in the studentication is numerically 0.\nThis might occur if your statistic has no variation. It is not possible to\napply the studentized bootstrap if any of the variances the values returned\nby func have not variability when resampled. The estimated covariance\nis:\n\n [[6.09003938e+05 0.00000000e+00 4.05197396e+02 1.44864454e+04\n  3.80442704e+04 2.58022158e+05 4.02730564e+02            nan]\n [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n  0.00000000e+00 0.00000000e+00 0.00000000e+00            nan]\n [4.05197396e+02 0.00000000e+00 1.20137112e+00 1.89703409e+01\n  3.85232060e-01 1.19202958e+02 1.74850725e-01            nan]\n [1.44864454e+04 0.00000000e+00 1.89703409e+01 2.85057591e+03\n  2.97990400e+03 1.27440394e+04 2.54599447e+00            nan]\n [3.80442704e+04 0.00000000e+00 3.85232060e-01 2.97990400e+03\n  2.19226485e+04 7.36443613e+04 1.36616522e+01            nan]\n [2.58022158e+05 0.00000000e+00 1.19202958e+02 1.27440394e+04\n  7.36443613e+04 3.70711246e+05 9.70945622e+01            nan]\n [4.02730564e+02 0.00000000e+00 1.74850725e-01 2.54599447e+00\n  1.36616522e+01 9.70945622e+01 3.78825688e-01            nan]\n [           nan            nan            nan            nan\n             nan            nan            nan            nan]]\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mStudentizationError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-23-adc3cdb23f1e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mbs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mIndependentSamplesBootstrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'year'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m16\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mA\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'year'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m20\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconf_int\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmean_diff\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"studentized\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/arch/bootstrap/base.py\u001b[0m in \u001b[0;36mconf_int\u001b[0;34m(self, func, reps, method, size, tail, extra_kwargs, reuse, sampling, std_err_func, studentize_reps)\u001b[0m\n\u001b[1;32m    690\u001b[0m                 \u001b[0mstd_err_func\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mstd_err_func\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    691\u001b[0m                 \u001b[0mstudentize_reps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mstudentize_reps\u001b[0m\u001b[0;34m,\u001b[0m  \u001b[0;31m# noqa\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 692\u001b[0;31m                 \u001b[0msampling\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msampling\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    693\u001b[0m             )\n\u001b[1;32m    694\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/arch/bootstrap/base.py\u001b[0m in \u001b[0;36m_construct_bootstrap_estimates\u001b[0;34m(self, func, reps, extra_kwargs, std_err_func, studentize_reps, sampling)\u001b[0m\n\u001b[1;32m    943\u001b[0m                 \u001b[0merr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mresults\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mcount\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mbase\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    944\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0many\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstd_err\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0meps\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 945\u001b[0;31m                     \u001b[0;32mraise\u001b[0m \u001b[0mStudentizationError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstudentization_error\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcov\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcov\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    946\u001b[0m                 \u001b[0mstudentized_results\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mcount\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0merr\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mstd_err\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    947\u001b[0m             \u001b[0mcount\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mStudentizationError\u001b[0m: \nThe estimated covariance computed in the studentication is numerically 0.\nThis might occur if your statistic has no variation. It is not possible to\napply the studentized bootstrap if any of the variances the values returned\nby func have not variability when resampled. The estimated covariance\nis:\n\n [[6.09003938e+05 0.00000000e+00 4.05197396e+02 1.44864454e+04\n  3.80442704e+04 2.58022158e+05 4.02730564e+02            nan]\n [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n  0.00000000e+00 0.00000000e+00 0.00000000e+00            nan]\n [4.05197396e+02 0.00000000e+00 1.20137112e+00 1.89703409e+01\n  3.85232060e-01 1.19202958e+02 1.74850725e-01            nan]\n [1.44864454e+04 0.00000000e+00 1.89703409e+01 2.85057591e+03\n  2.97990400e+03 1.27440394e+04 2.54599447e+00            nan]\n [3.80442704e+04 0.00000000e+00 3.85232060e-01 2.97990400e+03\n  2.19226485e+04 7.36443613e+04 1.36616522e+01            nan]\n [2.58022158e+05 0.00000000e+00 1.19202958e+02 1.27440394e+04\n  7.36443613e+04 3.70711246e+05 9.70945622e+01            nan]\n [4.02730564e+02 0.00000000e+00 1.74850725e-01 2.54599447e+00\n  1.36616522e+01 9.70945622e+01 3.78825688e-01            nan]\n [           nan            nan            nan            nan\n             nan            nan            nan            nan]]\n"
     ]
    }
   ],
   "source": [
    "def mean_diff(x, y):\n",
    "    return x.mean() - y.mean()\n",
    "\n",
    "bs = IndependentSamplesBootstrap(A.loc[A['year'] == 16], A.loc[A['year'] == 20])\n",
    "print(bs.conf_int(mean_diff, method=\"studentized\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "346ea7e2",
   "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}