Research & experiments with data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "50510fd9-1698-4b72-abc5-585d1a96a1a2",
   "metadata": {},
   "source": [
    "I decided to review article Kjærgaard J, Stensballe LG, Birk NM, et al. Lack of a Negative Effect of BCG-Vaccination on Child Psychomotor Development: Results from the Danish Calmette Study - A Randomised Clinical Trial. PLoS One. 2016;11(4):e0154541. Published 2016 Apr 28. doi:10.1371/journal.pone.0154541. In this paper researchears tried to assess the non-specific effect of Bacillus Calmette-Guérin (BCG) vaccination at birth on psychomotor development of children. The samples were collected in maternity units and pediatric wards at three university hospitals in Denmark. Randomization was stratified by pre-maturity (gestational age [GA] < 37 weeks or ≥ 37 weeks) and children were allocated 1:1 in permuting blocks of 2-4-6 using an online system."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6696743b-4c9c-4631-a9c7-e651816e676f",
   "metadata": {},
   "source": [
    "Firstly i tried to make synthetic data according to table 1. I focused only for participants who filled out the entire questionnaire (not partialy or non-responders). Original distribution between control and BCG group is 1647 for control group and 1779 responders for BCG group. My set consist from 10000 responders (5000 control and 5000 BCG group)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1c0d0c12-8f1a-45ad-ab49-f63a8fed439c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import libraries\n",
    "from matplotlib import pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import random, uuid\n",
    "from scipy.stats import mannwhitneyu\n",
    "from matplotlib import pyplot\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeff4ba7-0876-45cf-94db-f0d46f9c4c4c",
   "metadata": {},
   "source": [
    "Distribution in table 1 is represtened mostly in percentage. So i made function for generating data based on percentages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ad349381-1dbd-45b9-8414-ec775137a28c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate data based on percentage distribution\n",
    "def feat_gen(percentage, data: pd.DataFrame, col_name: str):\n",
    "    dis = [True] * int((percentage * len(data)) / 100) + [False] * int((len(data) - int((percentage * len(data)) / 100)))\n",
    "    random.shuffle(dis)\n",
    "    result = pd.DataFrame(dis, columns = [col_name])\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0200b963-56a8-4823-b379-8827aa088e7b",
   "metadata": {},
   "source": [
    "Now i can generate patient population with selected features: population size, age, smoking, premature birth, sex, older biological siblings and multiple births"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a52d547b-410a-4831-8439-b8ade7c00f3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function for generating mothers population\n",
    "def create_patient_population(sample_size: float,age_loc: float, age_scale: float, BCG_loc: float, BCG_scale: float, c_smoke: float, b_smoke: float,\n",
    "                             c_premature: float, b_premature: float, c_sex: float, b_sex: float, c_siblings: float, b_siblings: float, c_multiple: float,\n",
    "                             b_multiple: float):\n",
    "    # population size\n",
    "    population_control = pd.DataFrame(index=np.arange(int(sample_size/2)))\n",
    "    population_control.set_index([['control']*len(population_control)], inplace=True)\n",
    "    population_BCG = pd.DataFrame(index=np.arange(int(sample_size/2)))\n",
    "    population_BCG.set_index([['BCG']*len(population_BCG)], inplace=True)\n",
    "    \n",
    "    control = population_control.reset_index()\n",
    "    BCG = population_BCG.reset_index()\n",
    "    \n",
    "    # generate age from normal distribution\n",
    "    age_control = np.random.normal(loc = age_loc, scale = age_scale, size = int(sample_size/2)).astype(float)\n",
    "    age_control_df = pd.DataFrame(age_control, columns = ['age'])\n",
    "    age_BCG = np.random.normal(loc = BCG_loc, scale = BCG_scale, size = int(sample_size/2)).astype(float)\n",
    "    age_BCG_df = pd.DataFrame(age_BCG, columns = ['age'])\n",
    "    \n",
    "    # generate fetures based on table 1 percentage with function feat_gen\n",
    "    con_smoke = feat_gen(c_smoke,control,'smoke')\n",
    "    bcg_smoke = feat_gen(b_smoke,BCG,'smoke')\n",
    "    \n",
    "    con_premature = feat_gen(c_premature,control,'premature')\n",
    "    bcg_premature = feat_gen(b_premature,BCG,'premature')\n",
    "    \n",
    "    con_sex = feat_gen(c_sex,control,'sex')\n",
    "    bcg_sex = feat_gen(b_sex,BCG,'sex')\n",
    "    \n",
    "    con_siblings = feat_gen(c_siblings,control,'siblings')\n",
    "    bcg_siblings = feat_gen(b_siblings,BCG,'siblings')\n",
    "    \n",
    "    con_multiple = feat_gen(c_multiple,control,'multiple')\n",
    "    bcg_multiple = feat_gen(b_multiple,BCG,'multiple')\n",
    "    \n",
    "    dfs_control = [control, age_control_df, con_smoke, con_premature, con_sex, con_siblings, con_multiple]\n",
    "    control_fin = dfs_control[0].join(dfs_control[1:])\n",
    "    \n",
    "    # features for groups together\n",
    "    dfs_BCG = [BCG, age_BCG_df, bcg_smoke, bcg_premature, bcg_sex, bcg_siblings, bcg_multiple]\n",
    "    BCG_fin = dfs_BCG[0].join(dfs_BCG[1:])\n",
    "    \n",
    "    # concat groups\n",
    "    final_df = pd.concat([control_fin, BCG_fin], axis=0).reset_index(drop = True)\n",
    "    \n",
    "    return final_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c6705c4f-8203-4e15-a1aa-f0aa3444fb04",
   "metadata": {},
   "outputs": [],
   "source": [
    "# data based on table 1\n",
    "sample_size = 10000\n",
    "c_age_loc = 31.6\n",
    "c_age_scale = 4.3\n",
    "b_age_loc = 31.6\n",
    "b_age_scale = 4.5\n",
    "c_smoke = 9.3\n",
    "b_smoke = 9.4\n",
    "c_premature = 3.5\n",
    "b_premature = 3.2\n",
    "c_sex = 53.1\n",
    "b_sex = 52.2\n",
    "c_siblings = 39.7 \n",
    "b_siblings = 40.9\n",
    "c_multiple = 3.7\n",
    "b_multiple = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "71344ec5-2260-4806-bf77-0b8267525b80",
   "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>index</th>\n",
       "      <th>age</th>\n",
       "      <th>smoke</th>\n",
       "      <th>premature</th>\n",
       "      <th>sex</th>\n",
       "      <th>siblings</th>\n",
       "      <th>multiple</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>control</td>\n",
       "      <td>31.557752</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>control</td>\n",
       "      <td>25.057687</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>control</td>\n",
       "      <td>33.321448</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>control</td>\n",
       "      <td>35.037117</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>control</td>\n",
       "      <td>36.864770</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9995</th>\n",
       "      <td>BCG</td>\n",
       "      <td>42.995176</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9996</th>\n",
       "      <td>BCG</td>\n",
       "      <td>38.251725</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9997</th>\n",
       "      <td>BCG</td>\n",
       "      <td>36.881350</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9998</th>\n",
       "      <td>BCG</td>\n",
       "      <td>40.814282</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999</th>\n",
       "      <td>BCG</td>\n",
       "      <td>25.290442</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10000 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        index        age  smoke  premature    sex  siblings  multiple\n",
       "0     control  31.557752  False      False   True      True     False\n",
       "1     control  25.057687  False      False  False     False      True\n",
       "2     control  33.321448  False      False   True      True     False\n",
       "3     control  35.037117  False       True   True      True     False\n",
       "4     control  36.864770  False      False   True     False     False\n",
       "...       ...        ...    ...        ...    ...       ...       ...\n",
       "9995      BCG  42.995176  False      False   True      True     False\n",
       "9996      BCG  38.251725  False      False  False      True     False\n",
       "9997      BCG  36.881350  False      False   True     False     False\n",
       "9998      BCG  40.814282  False      False  False     False     False\n",
       "9999      BCG  25.290442  False      False  False     False     False\n",
       "\n",
       "[10000 rows x 7 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# generate population\n",
    "population = create_patient_population(sample_size, c_age_loc, c_age_scale, b_age_loc, b_age_scale, c_smoke, b_smoke, \n",
    "                              c_premature, b_premature, c_sex, b_sex, c_siblings, b_siblings, c_multiple, b_multiple)\n",
    "population"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c99afd43-5c5f-41d5-a4e4-248a5f09ec67",
   "metadata": {},
   "source": [
    "In table 1 is also information about education level so I add this. Function is based on percentage fractions, like in table 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5fd1cccf-e6ae-40ed-b198-d232bb580319",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function\n",
    "def generate_school(group_size: int, fracs:list, random_state:int=42):\n",
    "    \n",
    "    # empty df with defined nuber of rows\n",
    "    df = pd.DataFrame(index=range(group_size))\n",
    "    \n",
    "    # generate education based on fractions\n",
    "    remain = df.index.copy().to_frame()\n",
    "    res = []\n",
    "    for i in range(len(fracs)):\n",
    "        fractions_sum=sum(fracs[i:])\n",
    "        frac = fracs[i]/fractions_sum\n",
    "        idxs = remain.sample(frac=frac, random_state=random_state).index\n",
    "        remain=remain.drop(idxs)\n",
    "        res.append(idxs)\n",
    "        \n",
    "    samples = [df.loc[idxs] for idxs in res]\n",
    "    basic = samples[0].index = ['Basic schooling and non-theoretical education'] * len(samples[0])\n",
    "    basic = pd.DataFrame(basic, columns = ['school'])\n",
    "    BA = samples[1].index = ['Theoretical education incl. BA level'] * len(samples[1])\n",
    "    BA = pd.DataFrame(BA, columns = ['school'])\n",
    "    master = samples[2].index = ['Master level or more'] * len(samples[2])\n",
    "    master = pd.DataFrame(master, columns = ['school'])\n",
    "    \n",
    "    final_df = pd.concat([basic, BA, master], axis=0).reset_index(drop = True)\n",
    "    \n",
    "    return final_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b17f853a-568b-4878-b341-4e9b56d5bbac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate education\n",
    "control_school = generate_school(5000, [0.196,0.444,0.36])\n",
    "BCG_school = generate_school(5000, [0.206,0.445,0.349])\n",
    "\n",
    "# concat tables\n",
    "school = pd.concat([control_school, BCG_school], axis=0).reset_index(drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "dbceb75e-2407-4926-b83d-38a1d0654ebd",
   "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>index</th>\n",
       "      <th>age</th>\n",
       "      <th>smoke</th>\n",
       "      <th>premature</th>\n",
       "      <th>sex</th>\n",
       "      <th>siblings</th>\n",
       "      <th>multiple</th>\n",
       "      <th>school</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>control</td>\n",
       "      <td>31.557752</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>Basic schooling and non-theoretical education</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>control</td>\n",
       "      <td>25.057687</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>Basic schooling and non-theoretical education</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>control</td>\n",
       "      <td>33.321448</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>Basic schooling and non-theoretical education</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>control</td>\n",
       "      <td>35.037117</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>Basic schooling and non-theoretical education</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>control</td>\n",
       "      <td>36.864770</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>Basic schooling and non-theoretical education</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9995</th>\n",
       "      <td>BCG</td>\n",
       "      <td>42.995176</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>Master level or more</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9996</th>\n",
       "      <td>BCG</td>\n",
       "      <td>38.251725</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>Master level or more</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9997</th>\n",
       "      <td>BCG</td>\n",
       "      <td>36.881350</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>Master level or more</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9998</th>\n",
       "      <td>BCG</td>\n",
       "      <td>40.814282</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>Master level or more</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999</th>\n",
       "      <td>BCG</td>\n",
       "      <td>25.290442</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>Master level or more</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10000 rows × 8 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        index        age  smoke  premature    sex  siblings  multiple  \\\n",
       "0     control  31.557752  False      False   True      True     False   \n",
       "1     control  25.057687  False      False  False     False      True   \n",
       "2     control  33.321448  False      False   True      True     False   \n",
       "3     control  35.037117  False       True   True      True     False   \n",
       "4     control  36.864770  False      False   True     False     False   \n",
       "...       ...        ...    ...        ...    ...       ...       ...   \n",
       "9995      BCG  42.995176  False      False   True      True     False   \n",
       "9996      BCG  38.251725  False      False  False      True     False   \n",
       "9997      BCG  36.881350  False      False   True     False     False   \n",
       "9998      BCG  40.814282  False      False  False     False     False   \n",
       "9999      BCG  25.290442  False      False  False     False     False   \n",
       "\n",
       "                                             school  \n",
       "0     Basic schooling and non-theoretical education  \n",
       "1     Basic schooling and non-theoretical education  \n",
       "2     Basic schooling and non-theoretical education  \n",
       "3     Basic schooling and non-theoretical education  \n",
       "4     Basic schooling and non-theoretical education  \n",
       "...                                             ...  \n",
       "9995                           Master level or more  \n",
       "9996                           Master level or more  \n",
       "9997                           Master level or more  \n",
       "9998                           Master level or more  \n",
       "9999                           Master level or more  \n",
       "\n",
       "[10000 rows x 8 columns]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# add education level to population table\n",
    "population_2 = pd.concat([population, school], axis=1)\n",
    "population_2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b96f53bf-e2ed-4318-a794-551933d77ca3",
   "metadata": {},
   "source": [
    "I have artificial population so I can generate their ASQ score. I work with Table 2 - Mean ASQ score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3fd19bf6-dedf-493e-ba6c-fcee17b19580",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function for generating ASQ score based on Table2 (mean and sd)\n",
    "def gen_score(pop_size: int, com_loc: float, com_scale: float, fin_loc: float, fin_scale: float, gross_loc: float, gross_scale: float, problem_loc:float,\n",
    "              problem_scale: float, person_loc: float, person_scale):\n",
    "    \n",
    "    # generate ASQ score from normal distribution\n",
    "    Communication = np.random.normal(loc = com_loc, scale = com_scale, size = pop_size).astype(int)\n",
    "    Fine_motor = np.random.normal(loc = fin_loc, scale = fin_scale, size = pop_size).astype(int)\n",
    "    Gross_motor = np.random.normal(loc = gross_loc, scale = gross_scale, size = pop_size).astype(int)\n",
    "    Problem_solving = np.random.normal(loc = problem_loc, scale = problem_scale, size = pop_size).astype(int)\n",
    "    Personal_social = np.random.normal(loc = person_loc, scale = person_scale, size = pop_size).astype(int)\n",
    "    \n",
    "    # convert numpy to df\n",
    "    df1 = pd.DataFrame(Communication, columns = ['Communication'])\n",
    "    df2 = pd.DataFrame(Fine_motor, columns = ['Fine_motor'])\n",
    "    df3 = pd.DataFrame(Gross_motor, columns = ['Gross_motor'])\n",
    "    df4 = pd.DataFrame(Problem_solving, columns = ['Problem_solving'])\n",
    "    df5 = pd.DataFrame(Personal_social, columns = ['Personal_social'])\n",
    "    \n",
    "    # concat dfs\n",
    "    df_all = pd.concat([df1, df2, df3, df4, df5], axis=1, join='inner')\n",
    "\n",
    "    return df_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "91ef42ca-041f-45f2-b226-6e3a24051bf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CONTROL GROUP ASQ score (table 2)\n",
    "pop_size = int(sample_size/2)\n",
    "\n",
    "com_loc = 29.7\n",
    "com_scale = 13.3\n",
    "\n",
    "fin_loc = 38.0\n",
    "fin_scale = 19.2\n",
    "\n",
    "gross_loc = 38.7\n",
    "gross_scale = 12.7\n",
    "\n",
    "problem_loc = 36.7\n",
    "problem_scale = 14.1\n",
    "\n",
    "person_loc = 36.7\n",
    "person_scale = 13.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4bc1fcd8-9e60-409e-9e6a-a166ffe9deea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# use function - control group\n",
    "control_score = gen_score(pop_size, com_loc, com_scale, fin_loc, fin_scale, gross_loc, gross_scale, problem_loc,\n",
    "                    problem_scale, person_loc, person_scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "72a135dc-c708-4f28-8ea7-70674542162d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# BCG GROUP ASQ score (table 2)\n",
    "pop_size = int(sample_size/2)\n",
    "\n",
    "com_loc = 29.2\n",
    "com_scale = 13.1\n",
    "\n",
    "fin_loc = 37.8\n",
    "fin_scale = 19.2\n",
    "\n",
    "gross_loc = 38.2\n",
    "gross_scale = 12.6\n",
    "\n",
    "problem_loc = 36.6\n",
    "problem_scale = 13.6\n",
    "\n",
    "person_loc = 36.2\n",
    "person_scale = 13.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ce94378a-9d9f-4aac-8b58-3095beea9ab7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# use function - BCG group\n",
    "BCG_score = gen_score(pop_size, com_loc, com_scale, fin_loc, fin_scale, gross_loc, gross_scale, problem_loc,\n",
    "                    problem_scale, person_loc, person_scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c7e64c90-0456-4da2-b159-e04af6fafe72",
   "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>Communication</th>\n",
       "      <th>Fine_motor</th>\n",
       "      <th>Gross_motor</th>\n",
       "      <th>Problem_solving</th>\n",
       "      <th>Personal_social</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>32</td>\n",
       "      <td>52</td>\n",
       "      <td>56</td>\n",
       "      <td>32</td>\n",
       "      <td>47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>15</td>\n",
       "      <td>61</td>\n",
       "      <td>47</td>\n",
       "      <td>28</td>\n",
       "      <td>42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>28</td>\n",
       "      <td>31</td>\n",
       "      <td>45</td>\n",
       "      <td>39</td>\n",
       "      <td>47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>36</td>\n",
       "      <td>31</td>\n",
       "      <td>53</td>\n",
       "      <td>21</td>\n",
       "      <td>56</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>11</td>\n",
       "      <td>31</td>\n",
       "      <td>33</td>\n",
       "      <td>39</td>\n",
       "      <td>42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9995</th>\n",
       "      <td>51</td>\n",
       "      <td>30</td>\n",
       "      <td>23</td>\n",
       "      <td>48</td>\n",
       "      <td>26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9996</th>\n",
       "      <td>46</td>\n",
       "      <td>70</td>\n",
       "      <td>57</td>\n",
       "      <td>35</td>\n",
       "      <td>55</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9997</th>\n",
       "      <td>32</td>\n",
       "      <td>18</td>\n",
       "      <td>52</td>\n",
       "      <td>23</td>\n",
       "      <td>59</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9998</th>\n",
       "      <td>45</td>\n",
       "      <td>13</td>\n",
       "      <td>29</td>\n",
       "      <td>51</td>\n",
       "      <td>32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999</th>\n",
       "      <td>31</td>\n",
       "      <td>57</td>\n",
       "      <td>51</td>\n",
       "      <td>48</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10000 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      Communication  Fine_motor  Gross_motor  Problem_solving  Personal_social\n",
       "0                32          52           56               32               47\n",
       "1                15          61           47               28               42\n",
       "2                28          31           45               39               47\n",
       "3                36          31           53               21               56\n",
       "4                11          31           33               39               42\n",
       "...             ...         ...          ...              ...              ...\n",
       "9995             51          30           23               48               26\n",
       "9996             46          70           57               35               55\n",
       "9997             32          18           52               23               59\n",
       "9998             45          13           29               51               32\n",
       "9999             31          57           51               48               30\n",
       "\n",
       "[10000 rows x 5 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# concat score for both groups\n",
    "score = pd.concat([control_score, BCG_score], axis=0).reset_index(drop = True)\n",
    "score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9ef60bfe-aa2e-44f9-a347-378e052383af",
   "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>index</th>\n",
       "      <th>age</th>\n",
       "      <th>smoke</th>\n",
       "      <th>premature</th>\n",
       "      <th>sex</th>\n",
       "      <th>siblings</th>\n",
       "      <th>multiple</th>\n",
       "      <th>school</th>\n",
       "      <th>Communication</th>\n",
       "      <th>Fine_motor</th>\n",
       "      <th>Gross_motor</th>\n",
       "      <th>Problem_solving</th>\n",
       "      <th>Personal_social</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>control</td>\n",
       "      <td>31.557752</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>Basic schooling and non-theoretical education</td>\n",
       "      <td>32</td>\n",
       "      <td>52</td>\n",
       "      <td>56</td>\n",
       "      <td>32</td>\n",
       "      <td>47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>control</td>\n",
       "      <td>25.057687</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>Basic schooling and non-theoretical education</td>\n",
       "      <td>15</td>\n",
       "      <td>61</td>\n",
       "      <td>47</td>\n",
       "      <td>28</td>\n",
       "      <td>42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>control</td>\n",
       "      <td>33.321448</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>Basic schooling and non-theoretical education</td>\n",
       "      <td>28</td>\n",
       "      <td>31</td>\n",
       "      <td>45</td>\n",
       "      <td>39</td>\n",
       "      <td>47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>control</td>\n",
       "      <td>35.037117</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>Basic schooling and non-theoretical education</td>\n",
       "      <td>36</td>\n",
       "      <td>31</td>\n",
       "      <td>53</td>\n",
       "      <td>21</td>\n",
       "      <td>56</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>control</td>\n",
       "      <td>36.864770</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>Basic schooling and non-theoretical education</td>\n",
       "      <td>11</td>\n",
       "      <td>31</td>\n",
       "      <td>33</td>\n",
       "      <td>39</td>\n",
       "      <td>42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9995</th>\n",
       "      <td>BCG</td>\n",
       "      <td>42.995176</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>Master level or more</td>\n",
       "      <td>51</td>\n",
       "      <td>30</td>\n",
       "      <td>23</td>\n",
       "      <td>48</td>\n",
       "      <td>26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9996</th>\n",
       "      <td>BCG</td>\n",
       "      <td>38.251725</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>Master level or more</td>\n",
       "      <td>46</td>\n",
       "      <td>70</td>\n",
       "      <td>57</td>\n",
       "      <td>35</td>\n",
       "      <td>55</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9997</th>\n",
       "      <td>BCG</td>\n",
       "      <td>36.881350</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>Master level or more</td>\n",
       "      <td>32</td>\n",
       "      <td>18</td>\n",
       "      <td>52</td>\n",
       "      <td>23</td>\n",
       "      <td>59</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9998</th>\n",
       "      <td>BCG</td>\n",
       "      <td>40.814282</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>Master level or more</td>\n",
       "      <td>45</td>\n",
       "      <td>13</td>\n",
       "      <td>29</td>\n",
       "      <td>51</td>\n",
       "      <td>32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999</th>\n",
       "      <td>BCG</td>\n",
       "      <td>25.290442</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>Master level or more</td>\n",
       "      <td>31</td>\n",
       "      <td>57</td>\n",
       "      <td>51</td>\n",
       "      <td>48</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10000 rows × 13 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        index        age  smoke  premature    sex  siblings  multiple  \\\n",
       "0     control  31.557752  False      False   True      True     False   \n",
       "1     control  25.057687  False      False  False     False      True   \n",
       "2     control  33.321448  False      False   True      True     False   \n",
       "3     control  35.037117  False       True   True      True     False   \n",
       "4     control  36.864770  False      False   True     False     False   \n",
       "...       ...        ...    ...        ...    ...       ...       ...   \n",
       "9995      BCG  42.995176  False      False   True      True     False   \n",
       "9996      BCG  38.251725  False      False  False      True     False   \n",
       "9997      BCG  36.881350  False      False   True     False     False   \n",
       "9998      BCG  40.814282  False      False  False     False     False   \n",
       "9999      BCG  25.290442  False      False  False     False     False   \n",
       "\n",
       "                                             school  Communication  \\\n",
       "0     Basic schooling and non-theoretical education             32   \n",
       "1     Basic schooling and non-theoretical education             15   \n",
       "2     Basic schooling and non-theoretical education             28   \n",
       "3     Basic schooling and non-theoretical education             36   \n",
       "4     Basic schooling and non-theoretical education             11   \n",
       "...                                             ...            ...   \n",
       "9995                           Master level or more             51   \n",
       "9996                           Master level or more             46   \n",
       "9997                           Master level or more             32   \n",
       "9998                           Master level or more             45   \n",
       "9999                           Master level or more             31   \n",
       "\n",
       "      Fine_motor  Gross_motor  Problem_solving  Personal_social  \n",
       "0             52           56               32               47  \n",
       "1             61           47               28               42  \n",
       "2             31           45               39               47  \n",
       "3             31           53               21               56  \n",
       "4             31           33               39               42  \n",
       "...          ...          ...              ...              ...  \n",
       "9995          30           23               48               26  \n",
       "9996          70           57               35               55  \n",
       "9997          18           52               23               59  \n",
       "9998          13           29               51               32  \n",
       "9999          57           51               48               30  \n",
       "\n",
       "[10000 rows x 13 columns]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# add score to table with population data\n",
    "population_3 = pd.concat([population_2, score], axis=1)\n",
    "population_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "3b113ecf-8ad2-42af-aeff-4f7941baa0cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now i can calculate mean score for each subgroup of questions and mean SUM ASQ score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "6b15f84f-16f0-47b8-b973-26323784adfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function for calculating mean, sum and sd\n",
    "def ASQ_score_means(ASQ_score: pd.DataFrame):\n",
    "    \n",
    "    # sum ASQ score\n",
    "    summary = ASQ_score.sum(1)\n",
    "    \n",
    "    # mean score\n",
    "    mean = summary.mean()\n",
    "    mean_frame = pd.DataFrame([mean], columns = [\"ASQ_score_mean\"])\n",
    "    \n",
    "    # standard deviation\n",
    "    sd = summary.std()\n",
    "    sd_frame = pd.DataFrame([sd], columns = [\"ASQ_score_std\"])\n",
    "    \n",
    "    # to df\n",
    "    summary2 = ASQ_score.mean().to_frame('ASQ_score_mean')\n",
    "    \n",
    "    summary2_sd = ASQ_score.std().to_frame('ASQ_score_std')\n",
    "    \n",
    "    # concat \n",
    "    df_all = pd.concat([mean_frame, sd_frame, summary2, summary2_sd], axis=1)\n",
    "    df_all = df_all.fillna(0)\n",
    "    df_all2 = df_all.groupby(lambda x:x, axis=1).sum().rename(index={0: 'sum'})\n",
    "    \n",
    "    return df_all2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "30a7c213-a0ea-4807-b059-55d21fc5625b",
   "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>ASQ_score_mean</th>\n",
       "      <th>ASQ_score_std</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sum</th>\n",
       "      <td>177.7228</td>\n",
       "      <td>33.046367</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Communication</th>\n",
       "      <td>28.8926</td>\n",
       "      <td>13.472452</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Fine_motor</th>\n",
       "      <td>37.7380</td>\n",
       "      <td>19.045364</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gross_motor</th>\n",
       "      <td>38.2864</td>\n",
       "      <td>12.801592</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Problem_solving</th>\n",
       "      <td>36.2954</td>\n",
       "      <td>13.985352</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Personal_social</th>\n",
       "      <td>36.5104</td>\n",
       "      <td>13.403157</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 ASQ_score_mean  ASQ_score_std\n",
       "sum                    177.7228      33.046367\n",
       "Communication           28.8926      13.472452\n",
       "Fine_motor              37.7380      19.045364\n",
       "Gross_motor             38.2864      12.801592\n",
       "Problem_solving         36.2954      13.985352\n",
       "Personal_social         36.5104      13.403157"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# score - control group\n",
    "fin_score_control = ASQ_score_means(control_score)\n",
    "fin_score_control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "ae7a19a4-e586-4463-ad7d-dea88f0ae1ee",
   "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>ASQ_score_mean</th>\n",
       "      <th>ASQ_score_std</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sum</th>\n",
       "      <td>175.9362</td>\n",
       "      <td>32.644118</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Communication</th>\n",
       "      <td>29.1678</td>\n",
       "      <td>13.000009</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Fine_motor</th>\n",
       "      <td>37.1996</td>\n",
       "      <td>19.147112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gross_motor</th>\n",
       "      <td>38.1682</td>\n",
       "      <td>12.512331</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Problem_solving</th>\n",
       "      <td>35.8478</td>\n",
       "      <td>13.336544</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Personal_social</th>\n",
       "      <td>35.5528</td>\n",
       "      <td>13.298277</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 ASQ_score_mean  ASQ_score_std\n",
       "sum                    175.9362      32.644118\n",
       "Communication           29.1678      13.000009\n",
       "Fine_motor              37.1996      19.147112\n",
       "Gross_motor             38.1682      12.512331\n",
       "Problem_solving         35.8478      13.336544\n",
       "Personal_social         35.5528      13.298277"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# score - BCG group\n",
    "fin_score_BCG = ASQ_score_means(BCG_score)\n",
    "fin_score_BCG"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0f6f293-ffdb-491a-9cba-5156e2b3e171",
   "metadata": {},
   "source": [
    "As you can see means are relativelly simillar between control and BCG groups. They are also close to score values in paper - Table 2. I want to see distributions and calcuate p value. I used Mann-Whitney U test, where H0 tells the both samples comes from same distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "26ad28c9-038f-4d6b-a454-638b0c2fa6ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# I decided to generate new data based on ASQ sum mean and sd. First i select loc and scale values from generated tables.\n",
    "\n",
    "# mean and sd \n",
    "loc_ct = fin_score_control.iloc[0:1,0:1].to_numpy().reshape(1)\n",
    "scale_ct = fin_score_control.iloc[0:1,1:2].to_numpy().reshape(1)\n",
    "\n",
    "loc_BCG = fin_score_BCG.iloc[0:1,0:1].to_numpy().reshape(1)\n",
    "scale_BCG = fin_score_BCG.iloc[0:1,1:2].to_numpy().reshape(1)\n",
    "\n",
    "# generate sum score from normal distribution\n",
    "control_sum = np.random.normal(loc= loc_ct, scale = scale_ct, size=5000).astype(int)\n",
    "BCG_sum = np.random.normal(loc= loc_BCG, scale= scale_BCG, size=5000).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "2eeec36c-b13f-4591-aa2d-f92efb9826c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# optimalization of bins size\n",
    "bins_control = np.histogram_bin_edges(control_sum, bins='auto', range=(control_sum.min(), control_sum.max()))\n",
    "bins_BCG = np.histogram_bin_edges(BCG_sum, bins='auto', range=(BCG_sum.min(), BCG_sum.max()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "87d4531e-1351-4c43-8c30-2e810db1d58c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6NElEQVR4nO3de3gU1cHH8d8mbDZZEnIl2aSEGCWBQgAVLIIWkZuk5SZWaKEKfSmVCrQRqBatGrUSxFfAFwptLQUUEB/fEusVCXIrIghRFNACvo0SamIQciVhE8i8f0S2Lrntktsk+X6eZx6zM2dnzowD+XHmzDkWwzAMAQAAmIhPS1cAAADgcgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOh1augJXorKyUl9++aWCgoJksVhaujoAAMADhmGouLhYMTEx8vGpu42kVQaUL7/8UrGxsS1dDQAAcAWys7PVpUuXOsu0yoASFBQkqeoEO3Xq1MK1AQAAnigqKlJsbKzr93hdWmVAufRYp1OnTgQUAABaGU+6Z9BJFgAAmA4BBQAAmA4BBQAAmE6r7IMCADAPwzB04cIFXbx4saWrAhOwWq3y9fVt8H4IKACAK1ZeXq6cnByVlpa2dFVgEhaLRV26dFFgYGCD9kNAAQBckcrKSmVlZcnX11cxMTHy8/Nj8Mx2zjAMnT59WqdOnVJCQkKDWlIIKACAK1JeXq7KykrFxsbKbre3dHVgEp07d9bnn3+uioqKBgUUOskCABqkviHL0b40VisadxUAADAdrwLKqlWr1KdPH9cIrgMHDtRbb73l2j5t2jRZLBa35cYbb3Tbh9Pp1Jw5cxQREaGOHTtq7NixOnXqVOOcDQAAaBO86oPSpUsXLVq0SN26dZMkrVu3TuPGjdOHH36oXr16SZJGjRqlNWvWuL7j5+fnto+UlBS99tpr2rRpk8LDwzVv3jyNHj1amZmZjfJaEgCgZS3NON6sx7tvRGKzHq+5rV27VikpKSooKGjpqjQrr1pQxowZox/84AdKTExUYmKinnzySQUGBmrfvn2uMjabTQ6Hw7WEhYW5thUWFmr16tV65plnNHz4cF133XVav369Dh8+rG3btjXeWQEA0IKuuuoqLVu2rKWr0apdcR+UixcvatOmTTp37pwGDhzoWr9z505FRkYqMTFRM2bMUF5enmtbZmamKioqNHLkSNe6mJgYJSUlae/evbUey+l0qqioyG0BAKA1u3jxoiorK1u6GqbldUA5fPiwAgMDZbPZNHPmTKWnp6tnz56SpOTkZG3YsEHbt2/XM888owMHDmjo0KFyOp2SpNzcXPn5+Sk0NNRtn1FRUcrNza31mGlpaQoODnYtsbGx3lYbAACXyspKPfXUU+rWrZtsNpu6du2qJ598UlLV77mhQ4cqICBA4eHh+sUvfqGSkhLXd6dNm6bx48frv//7vxUdHa3w8HDNmjVLFRUVkqQhQ4boiy++0H333efqjylVPaoJCQnR66+/rp49e8pms+mLL75Qfn6+7r77boWGhsputys5OVknTpy4ovP66KOPdOuttyooKEidOnVSv379dPDgQUlSamqqrr32Wrfyy5Yt01VXXVXt3BYuXKioqCiFhIToscce04ULF/Sb3/xGYWFh6tKli/76179eUf284fU4KN27d9ehQ4dUUFCgv/3tb5o6dap27dqlnj17atKkSa5ySUlJ6t+/v+Li4vTGG29owoQJte7TMIw6X0tasGCB5s6d6/pcVFRESAHaEE/6LLT1fgZoXgsWLNBzzz2npUuX6uabb1ZOTo7++c9/qrS0VKNGjdKNN96oAwcOKC8vTz//+c81e/ZsrV271vX9HTt2KDo6Wjt27NBnn32mSZMm6dprr9WMGTO0efNm9e3bV7/4xS80Y8YMt+OWlpYqLS1Nf/nLXxQeHq7IyEhNnjxZJ06c0KuvvqpOnTrpgQce0A9+8AN98sknslqtXp3XlClTdN1112nVqlXy9fXVoUOHvN7H9u3b1aVLF+3evVvvvvuupk+frvfee0+DBw/W/v379dJLL2nmzJkaMWJEk/4u9jqg+Pn5uTrJ9u/fXwcOHNCzzz6rP/3pT9XKRkdHKy4uzpUEHQ6HysvLlZ+f79aKkpeXp0GDBtV6TJvNJpvN5m1VAQCopri4WM8++6xWrFihqVOnSpKuueYa3XzzzXruuedUVlam559/Xh07dpQkrVixQmPGjNFTTz2lqKgoSVJoaKhWrFghX19f9ejRQz/84Q/1zjvvaMaMGQoLC5Ovr6+CgoLkcDjcjl1RUaGVK1eqb9++kuQKJu+++67r9+CGDRsUGxurV155RXfeeadX53by5En95je/UY8ePSRJCQkJXl+fsLAw/c///I98fHzUvXt3LV68WKWlpXrwwQclVYW7RYsW6d1339WPf/xjr/fvqQaPg2IYhusRzuXOnDmj7OxsRUdHS5L69esnq9WqjIwMV5mcnBwdOXKkzoACAEBj+fTTT+V0OjVs2LAat/Xt29cVTiTppptuUmVlpY4dO+Za16tXL7c3T6Ojo936XNbGz89Pffr0cTtehw4dNGDAANe68PBwde/eXZ9++qnX5zZ37lz9/Oc/1/Dhw7Vo0SL93//9n9f76NWrl9vge1FRUerdu7frs6+vr8LDwz0634bwKqA8+OCD+sc//qHPP/9chw8f1kMPPaSdO3dqypQpKikp0fz58/Xee+/p888/186dOzVmzBhFRETo9ttvlyQFBwdr+vTpmjdvnt555x19+OGH+ulPf6revXtr+PDhTXKCAAB8W0BAQK3b6upy8O31lz82sVgsHnV4DQgIcNuPYRhe16MuqampOnr0qH74wx9q+/bt6tmzp9LT0yVVjfh7+fEu9Zv5tprO7UrPtyG8CihfffWV7rrrLnXv3l3Dhg3T/v37tWXLFo0YMUK+vr46fPiwxo0bp8TERE2dOlWJiYl67733FBQU5NrH0qVLNX78eE2cOFE33XST7Ha7XnvtNcZAAQA0i4SEBAUEBOidd96ptq1nz546dOiQzp0751r37rvvysfHR4mJnveD8vPz08WLF+st17NnT124cEH79+93rTtz5oyOHz+u7373ux4f79sSExN13333aevWrZowYYJrbLLOnTsrNzfXLaQcOnToio7RHLzqg7J69epatwUEBOjtt9+udx/+/v5avny5li9f7s2hAQBoFP7+/nrggQd0//33y8/PTzfddJNOnz6to0ePasqUKXr00Uc1depUpaam6vTp05ozZ47uuusuV/8TT1x11VXavXu3fvzjH8tmsykiIqLGcgkJCRo3bpxmzJihP/3pTwoKCtJvf/tbfec739G4ceO8Oq+ysjL95je/0Y9+9CPFx8fr1KlTOnDggO644w5JVW8XnT59WosXL9aPfvQjbdmyRW+99ZY6derk1XGaC7MZAwAaVWt44+rhhx9Whw4d9Mgjj+jLL79UdHS0Zs6cKbvdrrffflu//vWvdcMNN8hut+uOO+7QkiVLvNr/448/rnvuuUfXXHONnE5nrY9yJGnNmjX69a9/rdGjR6u8vFyDBw/Wm2++6fXbN76+vjpz5ozuvvtuffXVV4qIiNCECRP02GOPSZK++93vauXKlVq4cKGeeOIJ3XHHHZo/f77+/Oc/e3Wc5mIx6rpqJlVUVKTg4GAVFhaaNvkB8ByvGbdO58+fV1ZWluLj4+Xv79/S1YFJ1HVfePP7m9mMAQCA6RBQAABoJXr16qXAwMAalw0bNrR09RoVfVAAAGgl3nzzzRpfDZbkVSfe1oCAAgBAKxEXF9fSVWg2POIBAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmw1s8AIDGVZAtlZ5pvuPZw6WQ2OY7HpoFAQUA0HgKsqU/fE+qKG2+Y1rt0qz3vQop06ZN07p161yfw8LCdMMNN2jx4sXq06ePJMkwDD333HNavXq1jh49qg4dOqhbt2766U9/ql/84hey2+2SqoZvf/rpp7V582b961//kt1u19VXX60777xTM2bMUGhoaOOebztBQAEANJ7SM1XhZMJzUkQzzJ/09XFp84yq43rZijJq1CitWbNGkpSbm6vf/e53Gj16tE6ePClJuuuuu7R582b97ne/04oVK9S5c2d99NFHWrZsma666iqNHz9eZ8+e1c0336yioiI98cQT6tevn/z8/PTZZ59p48aN2rhxo2bNmtXop90eEFAAAI0vIlGKubala1Enm80mh8MhSXI4HHrggQc0ePBgnT59Wjt27NCGDRv0yiuvaNy4ca7vXHXVVRo7dqyKiookSQ8++KBOnjypY8eO6Tvf+Y6rXI8ePTR69Og6ZzH+tpUrV2rp0qXKzs5WcHCwvv/97+t///d/XcdMSUlRSkqKq/y1116r8ePHKzU1VZJksVj0xz/+Ua+99pq2b9+uuLg4/fWvf1Xnzp3185//XAcOHFCfPn20fv16XXPNNQ25bM2GTrIAgHavpKREGzZsULdu3RQeHq4NGzaoe/fubuHkEovFouDgYFVWVuqll17ST3/6U7dwcnnZ+hw8eFC/+tWv9Pjjj+vYsWPasmWLBg8e7PU5PPHEE7r77rt16NAh9ejRQ5MnT9Y999yjBQsW6ODBg5Kk2bNne73flkILCoCW8a2OlJElJ2ssUmYNUbHN0Zy1Qjvy+uuvKzAwUJJ07tw5RUdH6/XXX5ePj49OnDih7t271/n906dPq6CgoFq5fv366dixY5KkMWPG6MUXX6xzPydPnlTHjh01evRoBQUFKS4uTtddd53X5/Ozn/1MEydOlCQ98MADGjhwoB5++GHddtttkqRf//rX+tnPfub1flsKAQVA87usI+WUWopV+Phr3fUvE1LQJG699VatWrVKknT27FmtXLlSycnJev/992UYhketH1L1VpL09HSVl5frgQceUFlZWb3fHzFihOLi4nT11Vdr1KhRGjVqlG6//XZXJ1xPXercK/1n4sDevXu7rTt//ryKiorUqVMnr/bdEggoAJrfZR0pN+yv3oISVpql5BOPKKCigICCJtGxY0d169bN9blfv34KDg7Wc889p8TERH366ad1fr9z584KCQnRP//5T7f1Xbt2lSQFBQWpoKCg3noEBQXpgw8+0M6dO7V161Y98sgjSk1N1YEDBxQSEiIfH59qfVlqmtHYarW6fr4UmmpaV1lZWW+dzIA+KABazjcdKfMCe1RbztrjW7p2aGcsFot8fHxUVlamyZMn6/jx4/r73/9erZxhGCosLJSPj48mTpyo9evX69///neDjt2hQwcNHz5cixcv1scff6zPP/9c27dvl1QVhHJyclxli4qKlJWV1aDjtQa0oAAAGt/Xx01/HKfTqdzcXElSfn6+VqxYoZKSEo0ZM0a33HKL0tPT9ZOf/EQPP/ywRowYoc6dO+vw4cNaunSp5syZo/Hjx2vhwoXauXOnBgwYoMcff1z9+/dXx44d9fHHH+u9995TUlJSvfV4/fXX9a9//UuDBw9WaGio3nzzTVVWVrr6tgwdOlRr167VmDFjFBoaqocffli+vr5XfN6tBQEFQJNZmlHzL4/IkpOaImnD/pPKC/TuOTtMzh5eNXDa5hnNd0yrveq4XtqyZYuio6MlVT1m6dGjh15++WUNGTJEkrRx40b9+c9/1l//+lf9/ve/V4cOHZSQkKC7777b1fE0PDxc77//vp566ik9/fTTysrKko+PjxISEjRp0iS3V4NrExISos2bNys1NVXnz59XQkKCXnzxRfXq1UuStGDBAv3rX//S6NGjFRwcrCeeeKJdtKBYDE9f0jaRoqIiBQcHq7CwsFV09AHaq9oDyj815aO7tKHvC8oL7OFRmftGNMOgX/DK+fPnlZWVpfj4ePn7+/9nA0Pdt2u13hfy7vc3LSgAgMYVEktgQIPRSRYAgCb0j3/8Q4GBgbUuqBktKAAANKH+/fvr0KFDLV2NVoeAAgBAEwoICHAbbwWe4REPAKBBWuG7FmhCjXU/EFAAAFfk0iilpaWlLVwTmEl5ebkkNXisFh7xAKjiyauhnr7O+c2+apsEMKy07Y/h0B74+voqJCREeXl5kiS73e7x/DVomyorK3X69GnZ7XZ16NCwiEFAAVBt8r5aWe3SrPfrDinf2ldtkwBKVRMBlllDrqS2MBGHo2qepEshBfDx8VHXrl0bHFYJKACqTd5Xo6+PV40OWnqm7oDyrX1t+D//WouVWUOYBLANsFgsio6OVmRkZI0T2KH98fPzk49Pw3uQEFAA/Mc3k/c11r7yvmIY+/bC19e3XcwPg+ZDJ1kAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6vGYMoFVYmnHco3L3jahlHBcArYpXLSirVq1Snz591KlTJ3Xq1EkDBw7UW2+95dpuGIZSU1MVExOjgIAADRkyREePHnXbh9Pp1Jw5cxQREaGOHTtq7NixOnXqVOOcDQAAaBO8CihdunTRokWLdPDgQR08eFBDhw7VuHHjXCFk8eLFWrJkiVasWKEDBw7I4XBoxIgRKi4udu0jJSVF6enp2rRpk/bs2aOSkhKNHj1aFy9ebNwzAwAArZZXAWXMmDH6wQ9+oMTERCUmJurJJ59UYGCg9u3bJ8MwtGzZMj300EOaMGGCkpKStG7dOpWWlmrjxo2SpMLCQq1evVrPPPOMhg8fruuuu07r16/X4cOHtW3btiY5QQAA0PpccSfZixcvatOmTTp37pwGDhyorKws5ebmauTIka4yNptNt9xyi/bu3StJyszMVEVFhVuZmJgYJSUlucrUxOl0qqioyG0BAABtl9cB5fDhwwoMDJTNZtPMmTOVnp6unj17Kjc3V5IUFRXlVj4qKsq1LTc3V35+fgoNDa21TE3S0tIUHBzsWmJjPZjuHQAAtFpeB5Tu3bvr0KFD2rdvn375y19q6tSp+uSTT1zbL59e2TCMeqdcrq/MggULVFhY6Fqys7O9rTYAAGhFvA4ofn5+6tatm/r376+0tDT17dtXzz77rByOqmnTL28JycvLc7WqOBwOlZeXKz8/v9YyNbHZbK43hy4tAACg7WrwQG2GYcjpdCo+Pl4Oh0MZGRmubeXl5dq1a5cGDRokSerXr5+sVqtbmZycHB05csRVBgAAwKuB2h588EElJycrNjZWxcXF2rRpk3bu3KktW7bIYrEoJSVFCxcuVEJCghISErRw4ULZ7XZNnjxZkhQcHKzp06dr3rx5Cg8PV1hYmObPn6/evXtr+PDhTXKCAACg9fEqoHz11Ve66667lJOTo+DgYPXp00dbtmzRiBEjJEn333+/ysrKdO+99yo/P18DBgzQ1q1bFRQU5NrH0qVL1aFDB02cOFFlZWUaNmyY1q5dK19f38Y9MwAA0Gp5FVBWr15d53aLxaLU1FSlpqbWWsbf31/Lly/X8uXLvTk0AABoR5gsEAAAmA4BBQAAmA6zGQPwWl0zC0eWnNQUSRv2n5QCezRfpQC0KQQUAN75+rgiS/xr3RxWmtWMlQHQVhFQAHjGHi5Z7dLmGZpST9EKH3+VWUOao1YA2igCCgA3dT2+Cer7kgIqCurdR5k1RMU2RyPWCkB7Q0AB4LFim4PgAaBZ8BYPAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHV4zBtCqBTlz3cdm+bK05oL2cCkktlnqBKDhCChAO9EW588JcuZq6gd3ylp5/j8rP6qlsNUuzXqfkAK0EgQUAK1WQEWBrJXn9VbC4zprj5ckTRnQtXrBr49Lm2dIpWcIKEArQUAB0Oqdtccr71LLT0xiy1YGQKOgkywAADAdWlAAtCk19bX5dh+bvEC77htBKwtgdrSgAAAA0yGgAAAA0yGgAAAA0yGgAAAA06GTLNAG1DUIm1Q1oFnkt0dbvUxYaVYj1wgAGoaAArRxNY62WoMKH3+VWUOap1IAUA8CCtDG1TTaak3KrCEqtjmasWYAUDsCCtBOuI22CgAmR0ABWruCbEWW/LPWza29f0ld9W/t5wagdgQUoDUryJb+8D1NqSits1hr7F9SZg1RhY+/kk88Ume51nhuAOpHQAFas9IzUkVpm+xfUmxzaN31LyugjrePpNZ5bgDqR0AB2oC22r+k2OYgfADtFAO1AQAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0/EqoKSlpemGG25QUFCQIiMjNX78eB07dsytzLRp02SxWNyWG2+80a2M0+nUnDlzFBERoY4dO2rs2LE6depUw88GAAC0CV4FlF27dmnWrFnat2+fMjIydOHCBY0cOVLnzp1zKzdq1Cjl5OS4ljfffNNte0pKitLT07Vp0ybt2bNHJSUlGj16tC5evNjwMwIAAK2eV0Pdb9myxe3zmjVrFBkZqczMTA0ePNi13mazyeGoeXjqwsJCrV69Wi+88IKGDx8uSVq/fr1iY2O1bds23Xbbbd6eAwAAaGMa1AelsLBQkhQWFua2fufOnYqMjFRiYqJmzJihvLw817bMzExVVFRo5MiRrnUxMTFKSkrS3r17azyO0+lUUVGR2wIAANquKw4ohmFo7ty5uvnmm5WUlORan5ycrA0bNmj79u165plndODAAQ0dOlROp1OSlJubKz8/P4WGhrrtLyoqSrm5uTUeKy0tTcHBwa4lNjb2SqsNAABagSuezXj27Nn6+OOPtWfPHrf1kyZNcv2clJSk/v37Ky4uTm+88YYmTJhQ6/4Mw5DFYqlx24IFCzR37lzX56KiIkIKAABt2BUFlDlz5ujVV1/V7t271aVLlzrLRkdHKy4uTidOnJAkORwOlZeXKz8/360VJS8vT4MGDapxHzabTTab7UqqCrRuBdlS6Znat399vPnqAgDNyKuAYhiG5syZo/T0dO3cuVPx8fH1fufMmTPKzs5WdHS0JKlfv36yWq3KyMjQxIkTJUk5OTk6cuSIFi9efAWnALRRBdnSH74nVZTWXc5qV5k1pFmq1NqFlWZV/fBlLdfUHi6F0DoLmIFXAWXWrFnauHGj/v73vysoKMjVZyQ4OFgBAQEqKSlRamqq7rjjDkVHR+vzzz/Xgw8+qIiICN1+++2ustOnT9e8efMUHh6usLAwzZ8/X71793a91QNAVS0nFaXShOekiMTay9nDVXygrPnq1QqVWUNU4eOv5BOPVK34qJaCVrs0631CCmACXgWUVatWSZKGDBnitn7NmjWaNm2afH19dfjwYT3//PMqKChQdHS0br31Vr300ksKCgpylV+6dKk6dOigiRMnqqysTMOGDdPatWvl6+vb8DMC2pqIRCnm2noK8ainLsU2h9Zd/7ICKgokSVMGdK1e6Ovj0uYZVcGQgAK0OK8f8dQlICBAb7/9dr378ff31/Lly7V8+XJvDg8AV6zY5lCx7ZvxmWLqaJECYArMxQMAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEynQ0tXAACa29KM49XWRZac1BRJG/afVF6gXZJ034jEZq4ZgEtoQQEAAKZDQAEAAKbDIx7A5L79yAEA2gtaUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOnwFg/QAmoaKOxylwYOQ8vx5P+TxIBuQFOgBQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJiOVwElLS1NN9xwg4KCghQZGanx48fr2LFjbmUMw1BqaqpiYmIUEBCgIUOG6OjRo25lnE6n5syZo4iICHXs2FFjx47VqVOnGn42AACgTfAqoOzatUuzZs3Svn37lJGRoQsXLmjkyJE6d+6cq8zixYu1ZMkSrVixQgcOHJDD4dCIESNUXFzsKpOSkqL09HRt2rRJe/bsUUlJiUaPHq2LFy823pkBAIBWy6uRZLds2eL2ec2aNYqMjFRmZqYGDx4swzC0bNkyPfTQQ5owYYIkad26dYqKitLGjRt1zz33qLCwUKtXr9YLL7yg4cOHS5LWr1+v2NhYbdu2TbfddlsjnRoAAGitGtQHpbCwUJIUFhYmScrKylJubq5GjhzpKmOz2XTLLbdo7969kqTMzExVVFS4lYmJiVFSUpKrDAAAaN+ueC4ewzA0d+5c3XzzzUpKSpIk5ebmSpKioqLcykZFRemLL75wlfHz81NoaGi1Mpe+fzmn0ymn0+n6XFRUdKXVBpqcp/O3AABqd8UtKLNnz9bHH3+sF198sdo2i8Xi9tkwjGrrLldXmbS0NAUHB7uW2NjYK602AABoBa4ooMyZM0evvvqqduzYoS5durjWOxwOSarWEpKXl+dqVXE4HCovL1d+fn6tZS63YMECFRYWupbs7OwrqTYAAGglvAoohmFo9uzZ2rx5s7Zv3674+Hi37fHx8XI4HMrIyHCtKy8v165duzRo0CBJUr9+/WS1Wt3K5OTk6MiRI64yl7PZbOrUqZPbAgAA2i6v+qDMmjVLGzdu1N///ncFBQW5WkqCg4MVEBAgi8WilJQULVy4UAkJCUpISNDChQtlt9s1efJkV9np06dr3rx5Cg8PV1hYmObPn6/evXu73uoBAADtm1cBZdWqVZKkIUOGuK1fs2aNpk2bJkm6//77VVZWpnvvvVf5+fkaMGCAtm7dqqCgIFf5pUuXqkOHDpo4caLKyso0bNgwrV27Vr6+vg07GwAA0CZ4FVAMw6i3jMViUWpqqlJTU2st4+/vr+XLl2v58uXeHB4AALQTzMUDAABMh4ACAABMh4ACAABM54pHkgWAtiisNKvO7WXWEBXbHM1UG6D9IqAALSDImauAioI6y9T3ixKNq8waogoffyWfeKTOchU+/lp3/cuEFKCJEVCAZhbkzNXUD+6UtfJ8vWUrfPxVZg1p+kpBxTaH1l3/cp3BMaw0S8knHlFARQEBBWhiBBSgmQVUFMhaeV5vJTyus/b4OsvyOKF5FdscXG/AJAgoQAs5a49XXmCPlq4GAJgSb/EAAADToQUFABpbQbZUeqb+cvZwKSS26esDtEIEFABoTAXZ0h++J1WU1l/WapdmvU9IAWpAQAGAxlR6piqcTHhOikisvdzXx6XNM6rKE1CAaggoANAUIhKlmGtbuhZAq0VAAQAvVRtE78tvPc75+njzVgZoowgoAOChWkeb/eiyglZ7VQdYAFeMgAIAHqpttNkpA7q6F+TtHKDBCCgA4IUaR5uNqaMzLIArwkBtAADAdAgoAADAdAgoAADAdOiDAgANtDSj/leL7xtBPxXAG7SgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0+EtHgBoBpe/6RNZclJTJG3Yf1J5gXbXet72AarQggIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHuXiARhbkzFVARUGt28NKs5qvMgDQShFQgEsKsqXSM3WXsYdLIbF17mPqB3fKWnm+zt1U+PirzBrifR0BoJ3wOqDs3r1bTz/9tDIzM5WTk6P09HSNHz/etX3atGlat26d23cGDBigffv2uT47nU7Nnz9fL774osrKyjRs2DCtXLlSXbp0ufIzARqiIFv6w/ekitK6y1nt0qz3aw8ppWdkrTyvtxIe11l7fK27KbOGqNjmaECFAaBt8zqgnDt3Tn379tXPfvYz3XHHHTWWGTVqlNasWeP67Ofn57Y9JSVFr732mjZt2qTw8HDNmzdPo0ePVmZmpnx9fb2tEtBwpWeqwsmE56SIWqa7//q4tHlGVdm6WlEknbXHKy+wRxNUFADaB68DSnJyspKTk+ssY7PZ5HDU/K/DwsJCrV69Wi+88IKGDx8uSVq/fr1iY2O1bds23Xbbbd5WCWg8EYlSzLUtXQsAaPea5C2enTt3KjIyUomJiZoxY4by8vJc2zIzM1VRUaGRI0e61sXExCgpKUl79+6tcX9Op1NFRUVuCwAAaLsavZNscnKy7rzzTsXFxSkrK0sPP/ywhg4dqszMTNlsNuXm5srPz0+hoaFu34uKilJubm6N+0xLS9Njjz3W2FUFgBZX7a2uLy/rB1Vfx2ygjWr0gDJp0iTXz0lJSerfv7/i4uL0xhtvaMKECbV+zzAMWSyWGrctWLBAc+fOdX0uKipSbCx/YAG0XmXWEFX4+Cv5xCPuGz66rGB9HbOBNqrJXzOOjo5WXFycTpw4IUlyOBwqLy9Xfn6+WytKXl6eBg0aVOM+bDabbDZbU1cVAJpNsc2hdde/XG3MnCkDuv7ngxcds4G2pskDypkzZ5Sdna3o6GhJUr9+/WS1WpWRkaGJEydKknJycnTkyBEtXry4qasDAKZRbHNUf908ppa3yIB2xuuAUlJSos8++8z1OSsrS4cOHVJYWJjCwsKUmpqqO+64Q9HR0fr888/14IMPKiIiQrfffrskKTg4WNOnT9e8efMUHh6usLAwzZ8/X71793a91QMAANo3rwPKwYMHdeutt7o+X+obMnXqVK1atUqHDx/W888/r4KCAkVHR+vWW2/VSy+9pKCgINd3li5dqg4dOmjixImugdrWrl3LGCgAAEDSFQSUIUOGyDCMWre//fbb9e7D399fy5cv1/Lly709PAAAaAeYiwfw0Ib9JzXlm//mBdprLBNZUlUGANAwTTJQGwAAQEPQggJ4qdrAWh5uAwB4joACeKjWgbUuU+HjrzJrSPNUCgDaKAIK4KHaBta6XJk1pPrYFgAArxBQAC/UOLAW0MyWZhz3qNx9Ixj0Da0XnWQBAIDpEFAAAIDpEFAAAIDpEFAAAIDpEFAAAIDpEFAAAIDpEFAAAIDpEFAAAIDpEFAAAIDpEFAAAIDpEFAAAIDpEFAAAIDpEFAAAIDpEFAAAIDpdGjpCgBmsGH/SU355r95gfaWrg4AtHu0oAAAANMhoAAAANMhoAAAANMhoAAAANOhkywAmMjSjOOunyNL6LyN9osWFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDpeB5Tdu3drzJgxiomJkcVi0SuvvOK23TAMpaamKiYmRgEBARoyZIiOHj3qVsbpdGrOnDmKiIhQx44dNXbsWJ06dapBJwIAANoOrwPKuXPn1LdvX61YsaLG7YsXL9aSJUu0YsUKHThwQA6HQyNGjFBxcbGrTEpKitLT07Vp0ybt2bNHJSUlGj16tC5evHjlZwIAANqMDt5+ITk5WcnJyTVuMwxDy5Yt00MPPaQJEyZIktatW6eoqCht3LhR99xzjwoLC7V69Wq98MILGj58uCRp/fr1io2N1bZt23Tbbbc14HQAAEBb0Kh9ULKyspSbm6uRI0e61tlsNt1yyy3au3evJCkzM1MVFRVuZWJiYpSUlOQqczmn06mioiK3BQAAtF2NGlByc3MlSVFRUW7ro6KiXNtyc3Pl5+en0NDQWstcLi0tTcHBwa4lNja2MasNAABMpkne4rFYLG6fDcOotu5ydZVZsGCBCgsLXUt2dnaj1RUAAJiP131Q6uJwOCRVtZJER0e71ufl5blaVRwOh8rLy5Wfn+/WipKXl6dBgwbVuF+bzSabzdaYVQWANiPImauAioLqG74s/c/P9nAphNZntB6NGlDi4+PlcDiUkZGh6667TpJUXl6uXbt26amnnpIk9evXT1arVRkZGZo4caIkKScnR0eOHNHixYsbszoA0OYFOXM19YM7Za08X33jR9/62WqXZr1PSEGr4XVAKSkp0Weffeb6nJWVpUOHDiksLExdu3ZVSkqKFi5cqISEBCUkJGjhwoWy2+2aPHmyJCk4OFjTp0/XvHnzFB4errCwMM2fP1+9e/d2vdUDeKwgWyo90+DdhJVmNUJlgOYXUFEga+V5vZXwuM7a4922TRnQteqHr49Lm2dU/VkhoKCV8DqgHDx4ULfeeqvr89y5cyVJU6dO1dq1a3X//ferrKxM9957r/Lz8zVgwABt3bpVQUFBru8sXbpUHTp00MSJE1VWVqZhw4Zp7dq18vX1bYRTQrtRkC394XtSRWn9ZeuRLKnCx19l1pAG7wtoCWft8coL7OG+MiaxZSoDNAKvA8qQIUNkGEat2y0Wi1JTU5WamlprGX9/fy1fvlzLly/39vDAf5SeqQonE56TIhr2F/GG/SdVZg1Rsc3RSJUDADREo/ZBAVpERKIUc22DdpEXaG+cugAAGgUBBW3a0ozjLV0FAMAVYDZjAABgOgQUAABgOgQUAABgOvRBAQCTq2ucHsbwQVtFQAEAkyqzhqjCx1/JJx6ps1xtY/hc6iQeWXJSU1T1Ov3lb6zdN4KxUmBOBBQAMKlim0Prrn+55nl2voUxfNAWEVAAwMSKbQ7CB9olOskCAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTYbJAtEpLM47XOYU8AKB1owUFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDgEFAACYDiPJAkA7tjTjuEfl7huR2MQ1AdzRggIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHTrIA0E6ElWbVW6bMGqJim6MZagPUjYAC8yrIlkrP1LgpsuSkR3/ZAqgKHRU+/ko+8Ui9ZSt8/LXu+pcJKWhxBBSYU0G29IfvSRWlNW6e8s1/K3z8VWYNabZqAa1Rsc2hdde/rICKgjrLhZVmKfnEIwqoKCCgoMURUGBOpWeqwsmE56SI6uMvbNh/UhLN0YCnim2OBv1ZYbwUNLdG7ySbmpoqi8Xitjgc//lDYRiGUlNTFRMTo4CAAA0ZMkRHjx5t7GqgrYhIlGKurbbkBfZQXmAPwgkAtFFN8hZPr169lJOT41oOHz7s2rZ48WItWbJEK1as0IEDB+RwODRixAgVFxc3RVUAAEAr1CSPeDp06ODWanKJYRhatmyZHnroIU2YMEGStG7dOkVFRWnjxo265557mqI6MKM6OsBKkr72rDkZANA2NUlAOXHihGJiYmSz2TRgwAAtXLhQV199tbKyspSbm6uRI0e6ytpsNt1yyy3au3dvrQHF6XTK6XS6PhcVFTVFtdFc6ukA62K1S/bw5qkTAMBUGj2gDBgwQM8//7wSExP11Vdf6fe//70GDRqko0ePKjc3V5IUFRXl9p2oqCh98cUXte4zLS1Njz32WGNXFS2lng6wLvZwKSS2+eoFADCNRg8oycnJrp979+6tgQMH6pprrtG6det04403SpIsFovbdwzDqLbu2xYsWKC5c+e6PhcVFSk2ll9crd6lDrCX+c/bAmWSeNQDAO1Rk79m3LFjR/Xu3VsnTpzQ+PHjJUm5ubmKjo52lcnLy6vWqvJtNptNNputqasKAFD9I87yej+aQ5MHFKfTqU8//VTf//73FR8fL4fDoYyMDF133XWSpPLycu3atUtPPfVUU1cFAFAHT0ecrfDx12s9FqvMGlp945ff9C3jES0aqNEDyvz58zVmzBh17dpVeXl5+v3vf6+ioiJNnTpVFotFKSkpWrhwoRISEpSQkKCFCxfKbrdr8uTJjV0VAIAXPBlxNqAiX2P+eb8mfPKrmgt89M1/rXZp1vuEFFyxRg8op06d0k9+8hN9/fXX6ty5s2688Ubt27dPcXFxkqT7779fZWVluvfee5Wfn68BAwZo69atCgoKauyqAAC85MmIs3WFmCkDulYNE7B5RlWHeAIKrlCjB5RNmzbVud1isSg1NVWpqamNfWgAQDOoM8TEMNQ9Ggdz8aBReDJPB3N0AAA8RUBBs7kUYiJLTmqKqib8ywu0t2ylAACmREABADSapRnH6/1HCK2p8AQBBZ6rY/6cyJKTkhgfAQDQOAgo8Ew98+dM+ea/FT7+Wnf9y4QUAECDEFDgmXrmz9mw/6TCSrOUfOIRBVQUEFAA1D4iLYO5wQMEFHinlvlz6OwK4JJ6R6T9ZjC3+lpc6avSvhFQAACNypMRaWlxRX0IKACARufJiLRAXQgoaHT1zYRa33YAAAgoaDSezoQqVT17LrOGNH2lAACtEgEFjcaT586XMF4KAKAuBBQ0Kp47AwAag09LVwAAAOByBBQAAGA6POIBALSYOt/q+7KU0WbbMQIKAKDZefTW30eSrHZp1vuElHaIgAIAaHaevPU35Zrz0uYZVXOBEVDaHQIKAKBF1PvWX0TNs6ejfaCTLAAAMB0CCgAAMB0CCgAAMB36oKBOSzOOS5IiS05qiqQN+08qL9DespUCALR5tKAAAADTIaAAAADTIaAAAADTIaAAAADToZMsAMCUNuyvv3P+fSMSm7dSaDYEFFQpyK4aTvoykSUnJdUzoRcAAI2MgIKqcPKH70kV1YeVnvKtnyt8/FVmDWm2agGA5MGMxxKzHrdBBJR26tL4JpIUWfJPTako1VsJj+usPb7W75RZQ+qeNwMAGpHHMx5LzHrcBhFQ2qBvhw9JCnLmVpsxNPJbP1/618lZe7zyAns0ce0AwDOezHgsVf0dlnziEW3Y8aHyAstqLENfldaHgNLGBTlzNfWDO2WtPF9nOR7fADCjemc8RptFQGnjAioKZK08z+MbAO3a5S3LtaGlxTwIKO0Ej28AAK0JA7UBAADToQXFrGoZl8QTl8YukRi/BED7Ud/fdzzKbl0IKGZUx7gknphy2Wc6wAJoyzx6HVlVfxeuu/7lZgspnvR7oc9L7QgoZlR6piqcTHhOivD+5t2w/6TbZ/7VAKAt8+R15EuvIgdUFNT596FHoeKGAI9auIOcRfzd2wAtGlBWrlypp59+Wjk5OerVq5eWLVum73//+y1ZpeZR3+Obr7/5AxKRKMVc61rtaS900RkWQDvj6evIDX3sHVCRL/3htx61cE/18ddrPRarzBpae6GCAAaXq0WLBZSXXnpJKSkpWrlypW666Sb96U9/UnJysj755BN17dq1paolqfYgUNOAZ5e7vLWiWvOdh49vKnz8te6DIhUf9TCUAABq5eljII9Y7dJP/ybZI2ovU/q1tHGKJnzyq7r39YldmvRC3ftqp8P4WwzDMFriwAMGDND111+vVatWudZ997vf1fjx45WWllbnd4uKihQcHKzCwkJ16tSp0etWU0DxZsCzbz/jrBZQvjwk/fmWao9veCwDAE3Lk39kemLKrdd5FBj+8vruOo8XUJGvCSc8aI3xZBh/T1+saOGw483v7xZpQSkvL1dmZqZ++9vfuq0fOXKk9u7dW6280+mU0+l0fS4sLJRUdaJN4fy5kmrrAs99qbKyMr169UM6a4+r8XthpV9oxL+eVGjOezK+KVN0LM+90JnPJKch2WKkwKtdq09aKt3LXZB0oXo9AABX5rwCJUtgg/dT5BMsefD75/SFeo7n10VFd70jlZ2tvcyZz6TXfiV98o4U3q3mMqVnpM2/kC7UPMy/mw4BWtd7rUr8ouotOmtoLcdrgEu/tz1qGzFawL///W9DkvHuu++6rX/yySeNxMTEauUfffRRQxILCwsLCwtLG1iys7PrzQot2knWYrG4fTYMo9o6SVqwYIHmzp3r+lxZWamzZ88qPDy8xvJoXEVFRYqNjVV2dnaTPFJDzbjuzY9r3jK47s2vpa65YRgqLi5WTExMvWVbJKBERETI19dXubm5buvz8vIUFVW92clms8lms7mtCwkJacoqogadOnXiL48WwHVvflzzlsF1b34tcc2Dg4M9KtciQ937+fmpX79+ysjIcFufkZGhQYMGtUSVAACAibTYI565c+fqrrvuUv/+/TVw4ED9+c9/1smTJzVz5syWqhIAADCJFgsokyZN0pkzZ/T4448rJydHSUlJevPNNxUXF9dSVUItbDabHn300WqP2dC0uO7Nj2veMrjuza81XPMWGwcFAACgNi3SBwUAAKAuBBQAAGA6BBQAAGA6BBQAAGA6BJR2avfu3RozZoxiYmJksVj0yiuvuG03DEOpqamKiYlRQECAhgwZoqNHj7qVcTqdmjNnjiIiItSxY0eNHTtWp06dasazaH3qu+7Tpk2TxWJxW2688Ua3Mlx376SlpemGG25QUFCQIiMjNX78eB07dsytDPd74/LkmnOvN75Vq1apT58+rsHXBg4cqLfeesu1vbXd5wSUdurcuXPq27evVqxYUeP2xYsXa8mSJVqxYoUOHDggh8OhESNGqLi42FUmJSVF6enp2rRpk/bs2aOSkhKNHj1aFy9ebK7TaHXqu+6SNGrUKOXk5LiWN99802071907u3bt0qxZs7Rv3z5lZGTowoULGjlypM6dO+cqw/3euDy55hL3emPr0qWLFi1apIMHD+rgwYMaOnSoxo0b5wohre4+b/jUf2jtJBnp6emuz5WVlYbD4TAWLVrkWnf+/HkjODjY+OMf/2gYhmEUFBQYVqvV2LRpk6vMv//9b8PHx8fYsmVLs9W9Nbv8uhuGYUydOtUYN25crd/hujdcXl6eIcnYtWuXYRjc783h8mtuGNzrzSU0NNT4y1/+0irvc1pQUE1WVpZyc3M1cuRI1zqbzaZbbrlFe/fulSRlZmaqoqLCrUxMTIySkpJcZXBldu7cqcjISCUmJmrGjBnKy8tzbeO6N1xhYaEkKSwsTBL3e3O4/Jpfwr3edC5evKhNmzbp3LlzGjhwYKu8zwkoqObSJI6XT9wYFRXl2pabmys/Pz+FhobWWgbeS05O1oYNG7R9+3Y988wzOnDggIYOHSqn0ymJ695QhmFo7ty5uvnmm5WUlCSJ+72p1XTNJe71pnL48GEFBgbKZrNp5syZSk9PV8+ePVvlfd5iQ93D/CwWi9tnwzCqrbucJ2VQu0mTJrl+TkpKUv/+/RUXF6c33nhDEyZMqPV7XHfPzJ49Wx9//LH27NlTbRv3e9Oo7ZpzrzeN7t2769ChQyooKNDf/vY3TZ06Vbt27XJtb033OS0oqMbhcEhStcScl5fnSt8Oh0Pl5eXKz8+vtQwaLjo6WnFxcTpx4oQkrntDzJkzR6+++qp27NihLl26uNZzvzed2q55TbjXG4efn5+6deum/v37Ky0tTX379tWzzz7bKu9zAgqqiY+Pl8PhUEZGhmtdeXm5du3apUGDBkmS+vXrJ6vV6lYmJydHR44ccZVBw505c0bZ2dmKjo6WxHW/EoZhaPbs2dq8ebO2b9+u+Ph4t+3c742vvmteE+71pmEYhpxOZ+u8z5u9Wy5Mobi42Pjwww+NDz/80JBkLFmyxPjwww+NL774wjAMw1i0aJERHBxsbN682Th8+LDxk5/8xIiOjjaKiopc+5g5c6bRpUsXY9u2bcYHH3xgDB061Ojbt69x4cKFljot06vruhcXFxvz5s0z9u7da2RlZRk7duwwBg4caHznO9/hujfAL3/5SyM4ONjYuXOnkZOT41pKS0tdZbjfG1d915x7vWksWLDA2L17t5GVlWV8/PHHxoMPPmj4+PgYW7duNQyj9d3nBJR2aseOHYakasvUqVMNw6h69fLRRx81HA6HYbPZjMGDBxuHDx9220dZWZkxe/ZsIywszAgICDBGjx5tnDx5sgXOpvWo67qXlpYaI0eONDp37mxYrVaja9euxtSpU6tdU667d2q63pKMNWvWuMpwvzeu+q4593rT+K//+i8jLi7O8PPzMzp37mwMGzbMFU4Mo/Xd5xbDMIzma68BAACoH31QAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6fw/7j2lYopTUVsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# histogram of distributions\n",
    "pyplot.hist(control_sum, bins = bins_control, alpha=0.5, label='control_sum')\n",
    "pyplot.hist(BCG_sum, bins = bins_BCG, alpha=1, label='BCG_sum', histtype=u'step')\n",
    "pyplot.legend(loc='upper right')\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "cbfb9b92-9d23-4f1a-8375-51f40aaf5f2f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.006161525142376038"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# p value\n",
    "U1, p = mannwhitneyu(control_sum, BCG_sum)\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd2b65c9-272a-4e20-ba32-bdf224db8fa4",
   "metadata": {},
   "source": [
    "p value is <0.01, so we can reject H0 hypothesis which means, samples comes from two different distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dd27c52-f702-47f4-8d25-3cc5c1617c67",
   "metadata": {},
   "source": [
    "#### Data from paper are collected at three university hospitals in Denmark. Now I split my artificial data set randomly to 3 subsets. And use only fraction of data from all 3 sets. Let's say 20%. In reality it can be economic nonsense and waste of precious data, but I want to see how big can be effect of this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "505870a9-a568-49f7-a2e6-e8e7800c05ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# split data randomly to 3 subsets\n",
    "ppp = population_3.set_index('index')\n",
    "\n",
    "# split to groups\n",
    "p_con = ppp.loc['control']\n",
    "p_BCG = ppp.loc['BCG']\n",
    "\n",
    "# shuffle\n",
    "shuffled_con = p_con.sample(frac=1)\n",
    "shuffled_BCG = p_BCG.sample(frac=1)\n",
    "\n",
    "# split\n",
    "result_con = np.array_split(shuffled_con, 3)\n",
    "result_BCG = np.array_split(shuffled_BCG, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "0dc544b4-e43a-4379-b538-5c210023a23a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define fraction of subsets which we want to use \n",
    "def frac(df_list: list, fraction: float):\n",
    "    lst = []\n",
    "    for DataFrame in df_list:\n",
    "        df_percent = DataFrame.sample(frac = fraction)\n",
    "        lst.append(df_percent)\n",
    "\n",
    "    select = pd.concat(lst)\n",
    "    \n",
    "    return select"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "b3c4276c-0c94-4cb2-b13f-6569925ddcb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select let's say 20% from each subset\n",
    "select_control = frac(result_con, 0.2)\n",
    "select_BCG = frac(result_BCG, 0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9701265-4b82-4869-9b08-04eed01fcd1b",
   "metadata": {},
   "source": [
    "Now, I will test distributions of several features which represent mothers "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9171a398-6bdd-4834-ad9d-24665025eba1",
   "metadata": {},
   "source": [
    "#### age"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "f3db17fc-4856-4c45-adff-4681f6819a87",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select age values\n",
    "age_control = select_control.iloc[:,0:1].to_numpy().reshape(999)\n",
    "age_BCG = select_BCG.iloc[:,0:1].to_numpy().reshape(999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "93adcf7b-326f-4930-9a28-d0c491855e53",
   "metadata": {},
   "outputs": [],
   "source": [
    "# bins optimization\n",
    "bins_control_age = np.histogram_bin_edges(age_control, bins='auto', range=(age_control.min(), age_control.max()))\n",
    "bins_BCG_age = np.histogram_bin_edges(age_BCG, bins='auto', range=(age_BCG.min(), age_BCG.max()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "236372ca-a022-4cd2-be51-2637e388f970",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzoElEQVR4nO3de1hVZf7//9dGcMNWBEFhQ6FSomaW5SHTGrVMzNJ0rKy00qlBS3Mim9HMDmQKqV1KSVlan7RRp67mo46fvjZJjeKUM3k2szxUJDTJRR44uRFI1u8Pfu5pyxn2Zi/Yz8d1reti3+tea7/37brg5b3XwWIYhiEAAAAT8fN2AQAAABcjoAAAANMhoAAAANMhoAAAANMhoAAAANMhoAAAANMhoAAAANMhoAAAANPx93YBDVFeXq6ffvpJwcHBslgs3i4HAADUgWEYKiwsVHR0tPz8ap4jaZYB5aefflJMTIy3ywAAAA2QnZ2tSy+9tMY+zTKgBAcHS6r4gO3atfNyNQAAoC4KCgoUExPj/Dtek2YZUC58rdOuXTsCCgAAzUxdTs/gJFkAAGA6BBQAAGA6BBQAAGA6zfIcFACAORmGoV9++UXnz5/3dinwkoCAALVq1arR+yGgAADcorS0VCdOnJDD4fB2KfAii8WiSy+9VG3btm3UfggoAIBGKy8vV2Zmplq1aqXo6Gi1bt2aG2n6IMMw9PPPP+vHH39UXFxco2ZSCCgAgEYrLS1VeXm5YmJiZLPZvF0OvKhjx4764YcfVFZW1qiAwkmyAAC3qe325Wj53DVzxpEEAABMh4ACAABMh3NQAAAetTT9aJO91xPDuzXZe3nDqlWrlJiYqLy8PG+X4nHMoAAA4EFdunRRamqqt8todggoAAB42fnz51VeXu7tMkyl3gFl+/btGj16tKKjo2WxWLRx40bnurKyMs2ePVtXXXWV2rRpo+joaD344IP66aefXPZRUlKiGTNmqEOHDmrTpo3uuOMO/fjjj43+MAAA1Fd5ebkWLlyorl27ymq1qlOnTlqwYIEk6eDBg7r55psVFBSk8PBwTZkyRUVFRc5tJ0+erLFjx+rll19WVFSUwsPDNX36dJWVlUmShg4dquPHj+uJJ56QxWJxXuGyatUqhYaG6sMPP1TPnj1ltVp1/PhxnTlzRg8++KDat28vm82mkSNH6tixYw36XN99953GjBmjyMhItW3bVv3799cnn3zi0ufEiRO6/fbbFRQUpNjYWK1bt67SjE9+fr6mTJmiiIgItWvXTjfffLMOHDjQoJrqo97noJw9e1a9e/fW7373O915550u6xwOh/bu3atnn31WvXv31pkzZ5SYmKg77rhDu3fvdvZLTEzU//3f/+m9995TeHi4nnzySY0aNUp79uxxy+1xATQdd5xf0NLPG4C5zZkzRytXrtTSpUt144036sSJEzp8+LAcDoduvfVWXX/99dq1a5dyc3P1+9//Xo899phWrVrl3H7r1q2KiorS1q1b9e233+qee+7RNddco4SEBK1fv169e/fWlClTlJCQ4PK+DodDKSkpeuuttxQeHq6IiAhNmDBBx44d06ZNm9SuXTvNnj1bt912m77++msFBATU63MVFRXptttu0/z58xUYGKjVq1dr9OjROnLkiDp16iRJevDBB3Xy5Elt27ZNAQEBmjlzpnJzc537MAxDt99+u8LCwrR582aFhITozTff1LBhw3T06FGFhYU1fOBrUe+AMnLkSI0cObLKdSEhIUpPT3dpW7Zsma677jplZWWpU6dOys/P19tvv60///nPuuWWWyRJa9asUUxMjD755BONGDGiAR8DAID6Kyws1CuvvKK0tDRNmjRJknT55Zfrxhtv1MqVK1VcXKx3331Xbdq0kSSlpaVp9OjRWrhwoSIjIyVJ7du3V1pamlq1aqUePXro9ttv16effqqEhASFhYWpVatWCg4Olt1ud3nvsrIyvf766+rdu7ckOYPJ559/rkGDBkmS1q5dq5iYGG3cuFF33313vT5b7969nfuWpPnz52vDhg3atGmTHnvsMR0+fFiffPKJdu3apX79+kmS3nrrLcXFxTm32bp1qw4ePKjc3FxZrVZJ0ssvv6yNGzfqr3/9q6ZMmVKvmurD4+eg5Ofny2KxKDQ0VJK0Z88elZWVKT4+3tknOjpavXr10o4dO6rcR0lJiQoKClwWAAAa65tvvlFJSYmGDRtW5brevXs7w4kk3XDDDSovL9eRI0ecbVdeeaXL7H9UVJTLLER1Wrdurauvvtrl/fz9/TVgwABnW3h4uLp3765vvvmm3p/t7NmzmjVrlnr27KnQ0FC1bdtWhw8fVlZWliTpyJEj8vf3V58+fZzbdO3aVe3bt3e+3rNnj4qKihQeHq62bds6l8zMTH333Xf1rqk+PHqZ8blz5/TUU09pwoQJateunSQpJydHrVu3dhkASYqMjFROTk6V+0lJSdELL7zgyVIBAD4oKCio2nWGYVR7V9Rft1/81YvFYqnTCa9BQUEu+zEMo9511ORPf/qTPv74Y7388svq2rWrgoKCdNddd6m0tLTW97ugvLxcUVFR2rZtW6V+FyYePMVjMyhlZWW69957VV5ertdff73W/jX9A8yZM0f5+fnOJTs7293lAgB8UFxcnIKCgvTpp59WWtezZ0/t379fZ8+edbZ9/vnn8vPzU7dudT9vqnXr1jp//nyt/Xr27KlffvlFX3zxhbPt1KlTOnr0qK644oo6v98F//znPzV58mT99re/1VVXXSW73a4ffvjBub5Hjx765ZdftG/fPmfbt99+63KPlT59+ignJ0f+/v7q2rWry9KhQ4d611QfHgkoZWVlGj9+vDIzM5Wenu6cPZEku92u0tJSnTlzxmWb3Nxc5/d5F7NarWrXrp3LAgBAYwUGBmr27NmaNWuW3n33XX333Xf697//rbffflsTJ05UYGCgJk2apK+++kpbt27VjBkz9MADD1T796oqXbp00fbt2/Wf//xHJ0+erLZfXFycxowZo4SEBH322Wc6cOCA7r//fl1yySUaM2ZMvT9b165dtX79eu3fv18HDhzQhAkTXGZ2evTooVtuuUVTpkzRzp07tW/fPk2ZMsVlZueWW27RwIEDNXbsWH388cf64YcftGPHDj3zzDMuF794gtu/4rkQTo4dO6atW7cqPDzcZX3fvn0VEBCg9PR0jR8/XlLFZU5fffWVFi1a5O5yAABeZvartJ599ln5+/vrueee008//aSoqCg98sgjstls+vjjj/X444+rf//+stlsuvPOO7VkyZJ67X/evHmaOnWqLr/8cpWUlFT71YokvfPOO3r88cc1atQolZaWavDgwdq8eXO9r+CRpKVLl+qhhx7SoEGD1KFDB82ePbvSOZzvvvuuHn74YQ0ePFh2u10pKSk6dOiQAgMDJVV8XbV582bNnTtXDz30kH7++WfZ7XYNHjy4XiGtISxGTSNVhaKiIn377beSpGuvvVZLlizRTTfdpLCwMEVHR+vOO+/U3r179eGHH7oUHxYWptatW0uSHn30UX344YdatWqVwsLC9Mc//lGnTp2q82XGBQUFCgkJUX5+PrMpgJdxmTGkinMOMzMzFRsb6/zjhubnxx9/dF5VW9WJw3VR07FQn7/f9Z5B2b17t2666Sbn65kzZ0qSJk2apKSkJG3atEmSdM0117hst3XrVg0dOlRSRarz9/fX+PHjVVxcrGHDhmnVqlXcAwUAgCb0j3/8Q0VFRbrqqqt04sQJzZo1S126dNHgwYO9XVr9A8rQoUNrnJ6qy4RMYGCgli1bpmXLltX37QEAgCoubz5+/HiV6958801NnDix1n2UlZXp6aef1vfff6/g4GANGjRIa9eubdBXSu7G04wBAGiGNm/e7Lyl/sXqen7IiBEjTHuDVAIKAADNUOfOnb1dgkfxNGMAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6XMUDAPCsvGzJcapp3ssWLoXGNM17waMIKAAAz8nLll67TipzNM37Bdik6TvrFVImT56s1atXO1+HhYWpf//+WrRoka6++mpJFTchXblypd5++20dOnTI+XTf+++/X1OmTJHNZpNUcSv3xYsXa/369fr+++9ls9l02WWX6e6771ZCQoLat2/v3s/bghFQAACe4zhVEU7GrZQ6ePiZSyePSusTKt6znrMot956q9555x1JUk5Ojp555hmNGjVKWVlZkqQHHnhA69ev1zPPPKO0tDR17NhRBw4cUGpqqrp06aKxY8fq9OnTuvHGG1VQUKAXX3xRffv2VevWrfXtt99q3bp1WrdunaZPn+72j91SEVAAAJ7XoZsUfY23q6iW1WqV3W6XJNntds2ePVuDBw/Wzz//rK1bt2rt2rXauHGjxowZ49ymS5cuuuOOO5xPCH766aeVlZWlI0eO6JJLLnH269Gjh0aNGlWnR8FI0po1a5SamqojR46oTZs2uvnmm5WamqqIiAhnn02bNunJJ5/Ujz/+qOuvv16TJ0/W5MmTdebMGYWGhkqSduzYoaeeekq7du1Shw4d9Nvf/lYpKSlq06ZNY4erSXCSLAAAv1JUVKS1a9eqa9euCg8P19q1a9W9e3eXcHKBxWJRSEiIysvL9f777+v+++93CScX962L0tJSvfjiizpw4IA2btyozMxMTZ482bn+hx9+0F133aWxY8dq//79mjp1qubOneuyj4MHD2rEiBEaN26cvvzyS73//vv67LPP9Nhjj9V9ILyMGRQAgM/78MMP1bZtW0nS2bNnFRUVpQ8//FB+fn46duyYunfvXuP2P//8s/Ly8ir169u3r44cOSJJGj16tP7yl7/UWstDDz3k/Pmyyy7Tq6++quuuu05FRUVq27at3njjDXXv3l2LFy+WJHXv3l1fffWVFixY4Nxu8eLFmjBhghITEyVJcXFxevXVVzVkyBAtX75cgYGBtQ+KlxFQAJhTY6/84GoO1MNNN92k5cuXS5JOnz6t119/XSNHjtTOnTtlGEadZz8u7rdhwwaVlpZq9uzZKi4urtM+9u3bp6SkJO3fv1+nT59WeXm5JCkrK0s9e/bUkSNH1L9/f5dtrrvuOpfXe/bs0bfffqu1a9c62wzDUHl5uTIzM3XFFVfUqRZvIqAAMB93XPnRgKs54LvatGmjrl27Ol/37dtXISEhWrlypbp166Zvvvmmxu07duyo0NBQHT582KW9U6dOkqTg4GDl5eXVWsfZs2cVHx+v+Ph4rVmzRh07dlRWVpZGjBih0tJSSaoyMF18fkt5ebmmTp2qP/zhD5Xe40JNZkdAAWA+jb3yoxFXcwBSxUyIn5+fiouLNWHCBN17773629/+Vuk8FMMwVFBQoJCQEI0fP15r1qzRs88+W+15KLU5fPiwTp48qZdeekkxMRXH7u7du1369OjRQ5s3b3Zpu7hPnz59dOjQIZfQ1dwQUACYl8mv/EA9nDxq6vcoKSlRTk6OJOnMmTNKS0tTUVGRRo8erSFDhmjDhg2677779Oyzz2r48OHq2LGjDh48qKVLl2rGjBkaO3askpOTtW3bNg0YMEDz5s1Tv3791KZNG3355Zf617/+pV69etVaR6dOndS6dWstW7ZMjzzyiL766iu9+OKLLn2mTp2qJUuWaPbs2Xr44Ye1f/9+rVq1StJ/v2KaPXu2rr/+ek2fPl0JCQlq06aNvvnmG6Wnp2vZsmUNHqemREABAHiOLbzi67b1CU3zfgG2ivesp7///e+KioqSVPF1TI8ePfTBBx9o6NChkqR169ZpxYoV+p//+R/Nnz9f/v7+iouL04MPPqgRI0ZIksLDw7Vz504tXLhQixcvVmZmpvz8/BQXF6d77rnHecJqTTp27KhVq1bp6aef1quvvqo+ffro5Zdf1h133OHsExsbq7/+9a968skn9corr2jgwIGaO3euHn30UVmtVknS1VdfrYyMDM2dO1e/+c1vZBiGLr/8ct1zzz31HhtvsRh1vTDbRC5Mp+Xn56tdu3beLgfwaUvTG/8/4yeGX/Q1zk/7pRVDpCkZDZtBaez2qLdz584pMzNTsbGxla8Q4Vb3HrdgwQK98cYbys7O9nYpNR4L9fn7zQwKAMCzQmN8MjR40uuvv67+/fsrPDxcn3/+uRYvXtys7nFSFwQUAACayD//+U+NHDmy2vVFRUV12s+xY8c0f/58nT59Wp06ddKTTz6pOXPmuKtMUyCgAADQRPr166f9+/c3ej9Lly7V0qVLG1+QiRFQAABoIkFBQc360t+mxLN4AABu0wyvu4CbuesYIKAAABotICBAkuRwNOLuv2gRLtzxtlWrVo3aD1/xAAAarVWrVgoNDVVubq4kyWaz1fn5NWg5ysvL9fPPP8tms8nfv3ERg4ACAHALu90uSc6QAt/k5+enTp06NTqgElAAeERwSY6CyvLq1vmni74WaIrbosPtLBaLoqKiFBERobKyMm+XAy9p3bq1/PwafwYJAQWA2wWX5GjS3rsVUH6ubhscqKKtgbcsh/e1atWq0ecfAAQUAG4XVJangPJz+ihunk7bYmvtP3FAFY9/99FblgOoQEAB4DGnbbHKbduj9o7R3WrvA8CnEFAAH+eOh/0BgLtxHxQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6/t4uAACWph9t1PZPDO/mpkoAmEW9Z1C2b9+u0aNHKzo6WhaLRRs3bnRZbxiGkpKSFB0draCgIA0dOlSHDh1y6VNSUqIZM2aoQ4cOatOmje644w79+OOPjfogAACg5ah3QDl79qx69+6ttLS0KtcvWrRIS5YsUVpamnbt2iW73a7hw4ersLDQ2ScxMVEbNmzQe++9p88++0xFRUUaNWqUzp8/3/BPAgAAWox6f8UzcuRIjRw5ssp1hmEoNTVVc+fO1bhx4yRJq1evVmRkpNatW6epU6cqPz9fb7/9tv785z/rlltukSStWbNGMTEx+uSTTzRixIhGfBwAANASuPUclMzMTOXk5Cg+Pt7ZZrVaNWTIEO3YsUNTp07Vnj17VFZW5tInOjpavXr10o4dO6oMKCUlJSopKXG+LigocGfZQMuTly05TtWpa0RRVqW24oBQFVrt7q4KAOrMrQElJydHkhQZGenSHhkZqePHjzv7tG7dWu3bt6/U58L2F0tJSdELL7zgzlKBlisvW3rtOqnMUafuE6toK/ML1Oo+HxBSAHiNR67isVgsLq8Nw6jUdrGa+syZM0czZ850vi4oKFBMTEzjCwVaIsepinAybqXUofarW9Z+4TqDEubI1MhjzymoLI+AAsBr3BpQ7PaKX2Y5OTmKiopytufm5jpnVex2u0pLS3XmzBmXWZTc3FwNGjSoyv1arVZZrVZ3lgq0fB26SdHX1Nott63N87UAQD259UZtsbGxstvtSk9Pd7aVlpYqIyPDGT769u2rgIAAlz4nTpzQV199VW1AAQAAvqXeMyhFRUX69ttvna8zMzO1f/9+hYWFqVOnTkpMTFRycrLi4uIUFxen5ORk2Ww2TZgwQZIUEhKihx9+WE8++aTCw8MVFhamP/7xj7rqqqucV/UAAADfVu+Asnv3bt10003O1xfODZk0aZJWrVqlWbNmqbi4WNOmTdOZM2c0YMAAbdmyRcHBwc5tli5dKn9/f40fP17FxcUaNmyYVq1apVatWrnhIwEAgOau3gFl6NChMgyj2vUWi0VJSUlKSkqqtk9gYKCWLVumZcuW1fftAQCAD+BhgQAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHTq/SweAGhplqYfbfQ+nhjezQ2VALiAGRQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6XMUDNGNVXX0SUZSliZLWfpGl3La2pi8KANyAGRQAAGA6zKAAaLlO1u3+JhFFWZXaigNCVWi1u7siAHVEQAHQ8tjCpQCbtD6hTt0nVtFW5heo1X0+IKQAXkJAAdDyhMZI03dKjlN16r72C9cZlDBHpkYee05BZXkEFMBLCCgAWqbQmIqlDjiZGDAfTpIFAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACmQ0ABAACm4/aA8ssvv+iZZ55RbGysgoKCdNlll2nevHkqLy939jEMQ0lJSYqOjlZQUJCGDh2qQ4cOubsUAADQTLk9oCxcuFBvvPGG0tLS9M0332jRokVavHixli1b5uyzaNEiLVmyRGlpadq1a5fsdruGDx+uwsJCd5cDAACaIbcHlH/9618aM2aMbr/9dnXp0kV33XWX4uPjtXv3bkkVsyepqamaO3euxo0bp169emn16tVyOBxat26du8sBAADNkNsDyo033qhPP/1UR48elSQdOHBAn332mW677TZJUmZmpnJychQfH+/cxmq1asiQIdqxY4e7ywEAAM2Qv7t3OHv2bOXn56tHjx5q1aqVzp8/rwULFui+++6TJOXk5EiSIiMjXbaLjIzU8ePHq9xnSUmJSkpKnK8LCgrcXTYAADARt8+gvP/++1qzZo3WrVunvXv3avXq1Xr55Ze1evVql34Wi8XltWEYldouSElJUUhIiHOJiYlxd9kAAMBE3B5Q/vSnP+mpp57Svffeq6uuukoPPPCAnnjiCaWkpEiS7Ha7pP/OpFyQm5tbaVblgjlz5ig/P9+5ZGdnu7tsAABgIm4PKA6HQ35+rrtt1aqV8zLj2NhY2e12paenO9eXlpYqIyNDgwYNqnKfVqtV7dq1c1kAAEDL5fZzUEaPHq0FCxaoU6dOuvLKK7Vv3z4tWbJEDz30kKSKr3YSExOVnJysuLg4xcXFKTk5WTabTRMmTHB3OQAAoBlye0BZtmyZnn32WU2bNk25ubmKjo7W1KlT9dxzzzn7zJo1S8XFxZo2bZrOnDmjAQMGaMuWLQoODnZ3OQAAoBlye0AJDg5WamqqUlNTq+1jsViUlJSkpKQkd789AB+0NP2ot0sA4GY8iwcAAJiO22dQALQMYY5Mr2wLABIBBcBFigNCVeYXqJHHnqu9cw3K/AJVHBDqnqIA+BwCCgAXhVa7Vvf5QEFleY3aT3FAqAqtdvcUBcDnEFAAVFJotRMuAHgVJ8kCAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADTIaAAAADT4UZtAOAJedmS41TDt7eFS6Ex7qsHaGYIKADgbnnZ0mvXSWWOhu8jwCZN30lIgc8ioACAuzlOVYSTcSulDt3qv/3Jo9L6hIr9EFDgowgoAOApHbpJ0dd4uwqgWeIkWQAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDoEFAAAYDr+3i4AQBXysiXHqVq7RRRlVWoLc2R6oiIAaFIEFMBs8rKl166Tyhy1dp1YTXuZX6CKA0LdWhZqtjT9qPPniKIsTZS09oss5ba11Wn7J4Z381BlQPNEQAHMxnGqIpyMWyl1qPmP1tovKs+gSFJxQKgKrXZPVAcATYKAAphVh25S9DU1dqnr/84BoLnhJFkAAGA6BBQAAGA6Hgko//nPf3T//fcrPDxcNptN11xzjfbs2eNcbxiGkpKSFB0draCgIA0dOlSHDh3yRCkA0GBhjkxFFB2u98KVVEDjuf0clDNnzuiGG27QTTfdpI8++kgRERH67rvvFBoa6uyzaNEiLVmyRKtWrVK3bt00f/58DR8+XEeOHFFwcLC7SwKAeikOCFWZX6BGHnuuwfvgSiqgcdweUBYuXKiYmBi98847zrYuXbo4fzYMQ6mpqZo7d67GjRsnSVq9erUiIyO1bt06TZ061d0lAUC9FFrtWt3nAwWV5TV4H1xJBTSO2wPKpk2bNGLECN19993KyMjQJZdcomnTpikhIUGSlJmZqZycHMXHxzu3sVqtGjJkiHbs2FFlQCkpKVFJSYnzdUFBgbvLBgAXhVY7AQPwIrefg/L9999r+fLliouL08cff6xHHnlEf/jDH/Tuu+9KknJyciRJkZGRLttFRkY6110sJSVFISEhziUmJsbdZQMAABNxe0ApLy9Xnz59lJycrGuvvVZTp05VQkKCli9f7tLPYrG4vDYMo1LbBXPmzFF+fr5zyc7OdnfZAADARNweUKKiotSzZ0+XtiuuuEJZWRV3vLTbK6ZML54tyc3NrTSrcoHValW7du1cFgAA0HK5PaDccMMNOnLkiEvb0aNH1blzZ0lSbGys7Ha70tPTnetLS0uVkZGhQYMGubscAADQDLn9JNknnnhCgwYNUnJyssaPH6+dO3dqxYoVWrFihaSKr3YSExOVnJysuLg4xcXFKTk5WTabTRMmTHB3OQAAoBlye0Dp37+/NmzYoDlz5mjevHmKjY1VamqqJk7873NXZ82apeLiYk2bNk1nzpzRgAEDtGXLFu6BAgAAJHnoYYGjRo3SqFGjql1vsViUlJSkpKQkT7w9AABo5ngWDwAAMB2PzKAAqJul6UcrtUUUZWmipLVfZCm3ra3pi4JX/PpYaOgx8MTwbh6oDPAOZlAAAIDpMIMCeEJetuQ4VWu3iKKsSm08CRcACCiA++VlS69dJ5U5au06sZp2noQLwNcRUAB3c5yqCCfjVkodaj4nYO0XlWdQJJ6ECwAEFMBTOnSToq+psQsnwQJA1ThJFgAAmA4BBQAAmA4BBQAAmA4BBQAAmA4BBQAAmA4BBQAAmA4BBQAAmA4BBQAAmA4BBQAAmA4BBQAAmA4BBQAAmA7P4gEAkwpzZNZvg59+9QRtW7gUGuPegoAmREABAJMpDghVmV+gRh57rn4bHvjVzwE2afpOQgqaLQIKAJhModWu1X0+UFBZXr22mzigU8UPJ49K6xMkxykCCpotAgoAmFCh1a5Cq71+G0V380wxgBdwkiwAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdnmYMVCUvu+JR9Q1x8qh7awEAH0RAAS6Wly29dp1U5mj4PgJski3cfTUBgI8hoAAXc5yqCCfjVkodujVsH7ZwKTTGvXUBgA8hoADV6dBNir7G21UAgE8ioAANtDSdc00AwFO4igcAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOV/EAQAtx4cqyiKIsTZS09oss5ba11Xn7J4Y38L4/gAcwgwIAAEzH4wElJSVFFotFiYmJzjbDMJSUlKTo6GgFBQVp6NChOnTokKdLAQAAzYRHA8quXbu0YsUKXX311S7tixYt0pIlS5SWlqZdu3bJbrdr+PDhKiws9GQ5AACgmfBYQCkqKtLEiRO1cuVKtW/f3tluGIZSU1M1d+5cjRs3Tr169dLq1avlcDi0bt06T5UDAACaEY8FlOnTp+v222/XLbfc4tKemZmpnJwcxcfHO9usVquGDBmiHTt2VLmvkpISFRQUuCwAAKDl8shVPO+995727t2rXbt2VVqXk5MjSYqMjHRpj4yM1PHjx6vcX0pKil544QX3FwoAAEzJ7TMo2dnZevzxx7VmzRoFBgZW289isbi8NgyjUtsFc+bMUX5+vnPJzs52a80AAMBc3D6DsmfPHuXm5qpv377OtvPnz2v79u1KS0vTkSNHJFXMpERFRTn75ObmVppVucBqtcpqtbq7VABo0cIcmfXb4CfHf3+2hUuhMe4tCKgHtweUYcOG6eDBgy5tv/vd79SjRw/Nnj1bl112mex2u9LT03XttddKkkpLS5WRkaGFCxe6uxwA8DnFAaEq8wvUyGPP1W/DA7/6OcAmTd9JSIHXuD2gBAcHq1evXi5tbdq0UXh4uLM9MTFRycnJiouLU1xcnJKTk2Wz2TRhwgR3lwMAPqfQatfqPh8oqCyvXttNHNCp4oeTR6X1CZLjFAEFXuOVW93PmjVLxcXFmjZtms6cOaMBAwZoy5YtCg4O9kY5ANDiFFrtKrTa67dRNLe6h3k0SUDZtm2by2uLxaKkpCQlJSU1xdsDAIBmhmfxAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA02mSpxkDZrQ0/WiV7RFFWZooae0XWcpta2vaogAAkphBAQAAJkRAAQAApkNAAQAApkNAAQAApkNAAQAApkNAAQAApkNAAQAApsN9UNAy5WVLjlM1dokoyqqyPcyR6YmKAAD1QEBBs1XdjdaCS3I0ae/dCig/V+P2E2tYV+YXqOKA0IYXBwBoFAIKWpygsjwFlJ/TR3HzdNoW26B9FAeEqtBqd3NlAIC6IqCgxTpti1Vu2x7eLgMA0ACcJAsAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEzH3907TElJ0fr163X48GEFBQVp0KBBWrhwobp37+7sYxiGXnjhBa1YsUJnzpzRgAED9Nprr+nKK690dzkAgIY6ebRx29vCpdAY99QCn+P2gJKRkaHp06erf//++uWXXzR37lzFx8fr66+/Vps2bSRJixYt0pIlS7Rq1Sp169ZN8+fP1/Dhw3XkyBEFBwe7uyQAQH3YwqUAm7Q+oXH7CbBJ03cSUtAgbg8of//7311ev/POO4qIiNCePXs0ePBgGYah1NRUzZ07V+PGjZMkrV69WpGRkVq3bp2mTp3q7pIAAPURGlMRLBynGr6Pk0crAo7jFAEFDeL2gHKx/Px8SVJYWJgkKTMzUzk5OYqPj3f2sVqtGjJkiHbs2FFlQCkpKVFJSYnzdUFBgYerBgAfFxpDsIBXeTSgGIahmTNn6sYbb1SvXr0kSTk5OZKkyMhIl76RkZE6fvx4lftJSUnRCy+84MlSAcDnLU1v5Dknkp4Y3s0NlQAevornscce05dffqm//OUvldZZLBaX14ZhVGq7YM6cOcrPz3cu2dnZHqkXAACYg8dmUGbMmKFNmzZp+/btuvTSS53tdrtdUsVMSlRUlLM9Nze30qzKBVarVVar1VOlwozysmv9/juiKKvK9jBHpicqAgA0IbcHFMMwNGPGDG3YsEHbtm1TbGysy/rY2FjZ7Xalp6fr2muvlSSVlpYqIyNDCxcudHc5aI7ysqXXrpPKHDV2m1jDujK/QBUHhLq1LAC1u/A1UURRliZKWvtFlnLb2uq8PV8R4QK3B5Tp06dr3bp1+tvf/qbg4GDnOSchISEKCgqSxWJRYmKikpOTFRcXp7i4OCUnJ8tms2nChAnuLgfNkeNURTgZt1LqUP0vq7VfVD2DIknFAaEqtNo9UR0AoAm4PaAsX75ckjR06FCX9nfeeUeTJ0+WJM2aNUvFxcWaNm2a80ZtW7Zs4R4ocNWhmxR9TbWr6/O/MgBA8+KRr3hqY7FYlJSUpKSkJHe/PQAAaAF4Fg8AADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdjz8sEKhKTc/8aOgNngAALQczKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQ4SRaekZdd8dC/akQUVf+gvzBHpicqAgA0IwQUuF9etvTadRVPJK7GxFp2UeYXqOKAULeWBQBoPggocD/HqYpwMm5lxROJq7D2i+pnUCSpOCBUhVa7J6oDADQDBBR4ToduUvQ1Va7i/iYAgJpwkiwAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAd7oMCADCnWh6ZUStbuBQa47560KQIKAAA86nDIzNqFWCTpu8kpDRTBBQAgPnU4ZEZNTp5VFqfULEfAkqzREBBZY2dVj151H21AGjW6v108p/+/xmTC79HanhkBlo2AgpcuWNaVaqYWrWFu6cmAM1OcUCoyvwCNfLYc/Xb8MCvfub3iE8joMBVY6dVL+DkNMCnFVrtWt3nAwWV5dVru4kDOv33Bb9HfBoBBVVjWhVAIxVa7Sq02uu3UXQj/mOEFoX7oAAAANMhoAAAANMhoAAAANPhHBQ0yNJ0LiUGAHgOMygAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0uIoHAGAa7rhC8Inh3I22JWAGBQAAmA4zKC1RXnbFQ/8a4iT3NwEAeB8BpaXJy5Zeu67iicQNVOYXqNV7C1R4iLACAPAOAkpL4zhVEU7Grax4InEV1n6RVeMuigNC6/8EUgAA3IiA0lJ16CZFX1Plqty2tqatBQC8pTFfW9vCpdAY99WCeiGgNDHOUAeAJmALlwJs0vqEhu8jwCZN30lI8RICigkFl+QoqCyv+g4/1XB+CSe5AkBFqJi+s3EXDKxPqNiegOIVBBSTCS7J0aS9dyug/Fz1nQ7UvA9OcgUAVQQLwkWzRUAxmaCyPAWUn9NHcfN02hbboH1wkisAoLnzakB5/fXXtXjxYp04cUJXXnmlUlNT9Zvf/MabJVVozH1EJLecWHXaFqvctj0atQ8AAJorrwWU999/X4mJiXr99dd1ww036M0339TIkSP19ddfq1OnTt4qyy33EeHEKgDwHndcjBBRlKWJqrgtQ32vfDTDhQwt4YIMrwWUJUuW6OGHH9bvf/97SVJqaqo+/vhjLV++XCkpKd4qq073EakRJ1YBANBoXgkopaWl2rNnj5566imX9vj4eO3YsaNS/5KSEpWUlDhf5+fnS5IKCgrcX1xhkVRiSNZoqe1lDd/+h/0VP1+kbe6PNW4e6DiughJDDodD5yyVtwcAeJ7D4Wjw72KP/G2qp3NnG//3wxOf48I+DcOovbPhBf/5z38MScbnn3/u0r5gwQKjW7dulfo///zzhiQWFhYWFhaWFrBkZ2fXmhW8epKsxWJxeW0YRqU2SZozZ45mzpzpfF1eXq7Tp08rPDzcpX9BQYFiYmKUnZ2tdu3aea7wZoixqRrjUj3GpnqMTfUYm+oxNhV/5wsLCxUdHV1rX68ElA4dOqhVq1bKyclxac/NzVVkZGSl/larVVar1aUtNDS02v23a9fOZ//xa8PYVI1xqR5jUz3GpnqMTfV8fWxCQkLq1M/Pw3VUqXXr1urbt6/S09Nd2tPT0zVo0CBvlAQAAEzEa1/xzJw5Uw888ID69eungQMHasWKFcrKytIjjzzirZIAAIBJeC2g3HPPPTp16pTmzZunEydOqFevXtq8ebM6d+7c4H1arVY9//zzlb4OAmNTHcaleoxN9Rib6jE21WNs6sdiGHW51gcAAKDpeOUcFAAAgJoQUAAAgOkQUAAAgOkQUAAAgOk0u4Cyfft2jR49WtHR0bJYLNq4caPL+smTJ8tisbgs119/vXeKbWIpKSnq37+/goODFRERobFjx+rIkSMufQzDUFJSkqKjoxUUFKShQ4fq0KFDXqq46dRlbHz12Fm+fLmuvvpq582jBg4cqI8++si53lePGan2sfHVY+ZiKSkpslgsSkxMdLb58nHza1WNDcdN3TS7gHL27Fn17t1baWlp1fa59dZbdeLECeeyefPmJqzQezIyMjR9+nT9+9//Vnp6un755RfFx8fr7Nmzzj6LFi3SkiVLlJaWpl27dslut2v48OEqLCz0YuWeV5exkXzz2Ln00kv10ksvaffu3dq9e7duvvlmjRkzxvnHxFePGan2sZF885j5tV27dmnFihW6+uqrXdp9+bi5oLqxkThu6qTxj/7zHknGhg0bXNomTZpkjBkzxiv1mE1ubq4hycjIyDAMwzDKy8sNu91uvPTSS84+586dM0JCQow33njDW2V6xcVjYxgcO7/Wvn1746233uKYqcKFsTEMjpnCwkIjLi7OSE9PN4YMGWI8/vjjhmHwu8Ywqh8bw+C4qatmN4NSF9u2bVNERIS6deumhIQE5ebmerskr8jPz5ckhYWFSZIyMzOVk5Oj+Ph4Zx+r1aohQ4Zox44dXqnRWy4emwt8/dg5f/683nvvPZ09e1YDBw7kmPmVi8fmAl8+ZqZPn67bb79dt9xyi0s7x031Y3OBLx83deXVpxl7wsiRI3X33Xerc+fOyszM1LPPPqubb75Ze/bs8am79xmGoZkzZ+rGG29Ur169JMn5cMaLH8gYGRmp48ePN3mN3lLV2Ei+fewcPHhQAwcO1Llz59S2bVtt2LBBPXv2dP4x8eVjprqxkXz7mHnvvfe0d+9e7dq1q9I6X/9dU9PYSL593NRHiwso99xzj/PnXr16qV+/furcubP+3//7fxo3bpwXK2tajz32mL788kt99tlnldZZLBaX14ZhVGpryaobG18+drp37679+/crLy9P//u//6tJkyYpIyPDud6Xj5nqxqZnz54+e8xkZ2fr8ccf15YtWxQYGFhtP188buoyNr563NRXi/yK59eioqLUuXNnHTt2zNulNJkZM2Zo06ZN2rp1qy699FJnu91ul/Tf/91ckJubW+l/Oi1VdWNTFV86dlq3bq2uXbuqX79+SklJUe/evfXKK69wzKj6samKrxwze/bsUW5urvr27St/f3/5+/srIyNDr776qvz9/Z3Hhi8eN7WNzfnz5ytt4yvHTX21+IBy6tQpZWdnKyoqytuleJxhGHrssce0fv16/eMf/1BsbKzL+tjYWNntdqWnpzvbSktLlZGRoUGDBjV1uU2qtrGpii8dOxczDEMlJSU+fcxU58LYVMVXjplhw4bp4MGD2r9/v3Pp16+fJk6cqP379+uyyy7z2eOmtrFp1apVpW185bipN6+dnttAhYWFxr59+4x9+/YZkowlS5YY+/btM44fP24UFhYaTz75pLFjxw4jMzPT2Lp1qzFw4EDjkksuMQoKCrxdusc9+uijRkhIiLFt2zbjxIkTzsXhcDj7vPTSS0ZISIixfv164+DBg8Z9991nREVFtfjxqW1sfPnYmTNnjrF9+3YjMzPT+PLLL42nn37a8PPzM7Zs2WIYhu8eM4ZR89j48jFTlYuvVPHl4+Zivx4bjpu6a3YBZevWrYakSsukSZMMh8NhxMfHGx07djQCAgKMTp06GZMmTTKysrK8XXaTqGpcJBnvvPOOs095ebnx/PPPG3a73bBarcbgwYONgwcPeq/oJlLb2PjysfPQQw8ZnTt3Nlq3bm107NjRGDZsmDOcGIbvHjOGUfPY+PIxU5WLA4ovHzcX+/XYcNzUncUwDKOpZ20AAABq0uLPQQEAAM0PAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJgOAQUAAJjO/wc/Z8sZZ4W01AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot data\n",
    "pyplot.hist(age_control, bins = bins_control_age, alpha=0.5, label='control_age')\n",
    "pyplot.hist(age_BCG, bins = bins_BCG_age, alpha=1, label='BCG_age', histtype=u'step')\n",
    "pyplot.legend(loc='upper right')\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "0a9aa299-a6c3-445c-ba7d-12eaec6ee4fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8683017529875808"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# p value\n",
    "U1, p = mannwhitneyu(age_control, age_BCG)\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e20103e-9c59-4fd6-b07a-aa3544019d58",
   "metadata": {},
   "source": [
    "H0 can't be rejected so both groups (control and BCG) are represented by mothers with same age distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7e5c1f8-ec1a-441d-b4ec-4109c8a02542",
   "metadata": {},
   "source": [
    "#### Now i want to see few binary stats from table 1:  smoking, premature birth, sex, older biological siblings and multiple births"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "37bc6688-1761-40a7-91c1-5b11b7551896",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select binary stats\n",
    "control_bin = select_control.iloc[:,1:6]\n",
    "BCG_bin = select_BCG.iloc[:,1:6]\n",
    "length = len(BCG_bin.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "5162c15e-0a7b-4c81-bb00-8150b3221dba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function for calculating percentage distribution\n",
    "def binary_stats(dataframe: pd.DataFrame, length:int, group_name: str):\n",
    "    \n",
    "    # proportion of true values in set\n",
    "    stats = []\n",
    "    for columns in dataframe:\n",
    "        data = (100 / length) * dataframe[columns].sum()\n",
    "        stats.append(data)\n",
    "    \n",
    "    stats_group = pd.DataFrame(stats)\n",
    "\n",
    "    BCG_bin_stats = pd.DataFrame(data=stats_group.T.values, columns=dataframe.columns).T.rename(columns={0: group_name})\n",
    "    \n",
    "    return BCG_bin_stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "6cb90775-4457-496f-9d6f-e520c8b74087",
   "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>control</th>\n",
       "      <th>BCG</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>smoke</th>\n",
       "      <td>10.010010</td>\n",
       "      <td>9.509510</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>premature</th>\n",
       "      <td>2.802803</td>\n",
       "      <td>3.203203</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sex</th>\n",
       "      <td>53.953954</td>\n",
       "      <td>55.655656</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>siblings</th>\n",
       "      <td>38.738739</td>\n",
       "      <td>40.140140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>multiple</th>\n",
       "      <td>3.503504</td>\n",
       "      <td>2.502503</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             control        BCG\n",
       "smoke      10.010010   9.509510\n",
       "premature   2.802803   3.203203\n",
       "sex        53.953954  55.655656\n",
       "siblings   38.738739  40.140140\n",
       "multiple    3.503504   2.502503"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# use function\n",
    "control_bin_stats = binary_stats(control_bin, length, 'control')\n",
    "BCG_bin_stats = binary_stats(BCG_bin, length, 'BCG')\n",
    "\n",
    "stats = pd.concat([control_bin_stats, BCG_bin_stats], axis=1)\n",
    "stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c03f6a80-401f-4d5d-b376-4cbad82ff9eb",
   "metadata": {},
   "source": [
    "Expect smoking, all other features follow original data representation from (paper: table 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90f3b4f4-7523-484f-8d88-e0f61dfe787c",
   "metadata": {},
   "source": [
    "#### ASQ score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "c0a86f0f-b4b8-4fe0-9667-8d9b5bd72a40",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select score data\n",
    "score_control_sel = select_control.iloc[:,7:]\n",
    "score_BCG_sel = select_BCG.iloc[:,7:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "29b75f9a-558a-41d1-8305-9aaa69a22e7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate means, sd and SUM ASQ score\n",
    "fin_score_control_sel = ASQ_score_means(score_control_sel)\n",
    "fin_score_BCG_sel = ASQ_score_means(score_BCG_sel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "46a67225-c5bf-4600-8aa2-ad033967acbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select loc and scale for both, control and BCG group\n",
    "loc_ct_sele = fin_score_control_sel.iloc[0:1,0:1].to_numpy().reshape(1)\n",
    "scale_ct_sele = fin_score_control_sel.iloc[0:1,1:2].to_numpy().reshape(1)\n",
    "\n",
    "loc_BCG_sele = fin_score_BCG_sel.iloc[0:1,0:1].to_numpy().reshape(1)\n",
    "scale_BCG_sele = fin_score_BCG_sel.iloc[0:1,1:2].to_numpy().reshape(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "e43a9f70-5887-4792-979a-9e5b0e773ca3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate new dataset based based on normal distribution and examined loc and scale values\n",
    "control_sele_sum = np.random.normal(loc = loc_ct_sele, scale = scale_ct_sele, size=5000).astype(int)\n",
    "BCG_sele_sum = np.random.normal(loc= loc_BCG_sele, scale = scale_BCG_sele, size=5000).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "f072d65d-38b4-4347-b341-2c8cf54f5692",
   "metadata": {},
   "outputs": [],
   "source": [
    "# bin value optimization\n",
    "bins_control_sele = np.histogram_bin_edges(control_sele_sum, bins='auto', range=(control_sele_sum.min(), control_sele_sum.max()))\n",
    "bins_BCG_sele = np.histogram_bin_edges(BCG_sele_sum, bins='auto', range=(BCG_sele_sum.min(), BCG_sele_sum.max()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "29339181-ced4-4bf1-9700-f37b1442ac5c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAABASElEQVR4nO3de3QU5eH/8c/mtkkgF5KQbFJCREhACKCA5VLlHiBfQTF+hRZUsBTxK1AjUClYbVAkii2gIFQpJSggHluwXhAJBeIPUxSiKFALsUaJNTGKuZKwiWR+f6RsXXK/T5L365w9h515duaZcZSPzzwXi2EYhgAAAEzEpbUrAAAAcCUCCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB231q5AQ5SXl+urr76Sj4+PLBZLa1cHAADUgWEYKiwsVFhYmFxcam4jaZMB5auvvlJ4eHhrVwMAADRAZmamunXrVmOZNhlQfHx8JFVcoK+vbyvXBgAA1EVBQYHCw8Mdf4/XpE0GlMuvdXx9fQkoAAC0MXXpnkEnWQAAYDoEFAAAYDoEFAAAYDqN6oOSmJio5cuX6/7779e6deskVQwhWrFihZ5//nnl5uZq6NChevbZZ9WvXz/H7+x2u5YsWaKXXnpJJSUlGjdunDZu3Fhrj14AQAXDMPT999/r0qVLrV0VwIm7u7tcXV0bfZwGB5Rjx47p+eef14ABA5y2r169WmvWrFFSUpKioqK0cuVKxcTE6MyZM45eu/Hx8Xr99de1a9cuBQYGavHixZo8ebLS0tKa5KIAoD0rLS1VVlaWiouLW7sqQCUWi0XdunVT586dG3ccwzCM+v6oqKhIgwYN0saNG7Vy5Upde+21WrdunQzDUFhYmOLj47V06VJJFa0lISEhevLJJzVv3jzl5+era9euevHFFzV9+nRJ/53XZO/evZo4cWKt5y8oKJCfn5/y8/MZxQOgQykvL1d6erpcXV3VtWtXeXh4MGElTMMwDH3zzTcqLi5WZGRkpUaH+vz93aAWlPnz5+umm27S+PHjtXLlSsf2jIwMZWdna8KECY5tVqtVo0aNUmpqqubNm6e0tDSVlZU5lQkLC1N0dLRSU1OrDCh2u112u93pAgGgIyotLVV5ebnCw8Pl7e3d2tUBKunatas+//xzlZWVNeqtSL0Dyq5du/TBBx/o2LFjlfZlZ2dLkkJCQpy2h4SE6IsvvnCU8fDwUJcuXSqVufz7KyUmJmrFihX1rSoAtFu1TRMOtJamatGr1xOemZmp+++/X9u3b5enp2e15a6snGEYtVa4pjLLli1Tfn6+45OZmVmfagMAgDamXgElLS1NOTk5Gjx4sNzc3OTm5qaUlBQ988wzcnNzc7ScXNkSkpOT49hns9lUWlqq3NzcastcyWq1OmaNZfZYAADav3q94hk3bpxOnjzptO3uu+9Wnz59tHTpUl199dWy2WxKTk7WddddJ6nifWlKSoqefPJJSdLgwYPl7u6u5ORkTZs2TZKUlZWlU6dOafXq1U1xTQDQ4axNPtui53sgJqpFz9fSkpKSFB8fr7y8vGY/11VXXaX4+HjFx8c3+7naknoFFB8fH0VHRztt69SpkwIDAx3b4+PjtWrVKkVGRioyMlKrVq2St7e3ZsyYIUny8/PTnDlztHjxYgUGBiogIEBLlixR//79NX78+Ca6LABAR8Nf9O1Lky8W+OCDD6qkpET33XefY6K2/fv3O61cuHbtWrm5uWnatGmOidqSkpKYAwUA0KwuXboki8VCJ+M2oNH/hA4fPuyYRVaq6CCbkJCgrKwsXbx4USkpKZVaXTw9PbV+/XqdP39excXFev311xUeHt7YqgAATKy8vFxPPvmkevXqJavVqu7du+vxxx+XJJ08eVJjx46Vl5eXAgMDdc8996ioqMjx29mzZ2vq1Kn63e9+p9DQUAUGBmr+/PkqKyuTJI0ePVpffPGFHnjgAVksFsegi6SkJPn7++uNN95Q3759ZbVa9cUXXyg3N1d33XWXunTpIm9vb8XGxio9Pb1B1/XRRx9pzJgx8vHxka+vrwYPHqzjx4879qempmrkyJHy8vJSeHi4fvnLX+rChQvVHi8/P1/33HOPgoOD5evrq7Fjx+qjjz5qdF0SEhJ07bXXOpVft26drrrqKsf3y/d51apVCgkJkb+/v1asWKHvv/9ev/rVrxQQEKBu3brpT3/6U91vUAM1eQsKADRWXftTtPd+EO3NsmXLtHnzZq1du1Y33HCDsrKy9M9//lPFxcWaNGmShg0bpmPHjiknJ0e/+MUvtGDBAiUlJTl+f+jQIYWGhurQoUP69NNPNX36dF177bWaO3eudu/erYEDB+qee+7R3Llznc5bXFysxMRE/fGPf1RgYKCCg4M1Y8YMpaen67XXXpOvr6+WLl2q//mf/9E//vEPubu71+u6Zs6cqeuuu06bNm2Sq6urTpw44TjGyZMnNXHiRD322GPasmWLvvnmGy1YsEALFizQ1q1bKx3LMAzddNNNCggI0N69e+Xn56fnnntO48aN09mzZxUQENDgutTVwYMH1a1bN73zzjt69913NWfOHP3973/XyJEj9d577+nll1/Wvffeq5iYmGZtXCCgAACaXWFhoZ5++mlt2LBBs2bNkiT17NlTN9xwgzZv3qySkhK98MIL6tSpkyRpw4YNmjJlip588knHCM8uXbpow4YNcnV1VZ8+fXTTTTfpb3/7m+bOnauAgAC5urrKx8dHNpvN6dxlZWXauHGjBg4cKEmOYPLuu+9qxIgRkqQdO3YoPDxcr776qm6//fZ6Xdu5c+f0q1/9Sn369JEkRUZGOvY99dRTmjFjhqNfTGRkpJ555hmNGjVKmzZtqjRlx6FDh3Ty5Enl5OTIarVKkn73u9/p1Vdf1Z///Gfdc889Da5LXQUEBOiZZ56Ri4uLevfurdWrV6u4uFjLly+XVBE0n3jiCb377rv66U9/Wu/j1xUv4QAAze6TTz6R3W7XuHHjqtw3cOBARziRpJ/85CcqLy/XmTNnHNv69evn1FcxNDRUOTk5tZ7bw8PDad24Tz75RG5ubho6dKhjW2BgoHr37q1PPvmk3te2aNEi/eIXv9D48eP1xBNP6F//+pdjX1pampKSktS5c2fHZ+LEiSovL1dGRkalY6WlpamoqEiBgYFOv8nIyHA6bkPqUlf9+vVz6qMTEhKi/v37O767uroqMDCwTve+MQgoAIBm5+XlVe2+mibq/OH2K19VWCwWlZeX1+ncPzxOdUvQ1WVS0aokJCTo9OnTuummm3Tw4EH17dtXe/bskVTR72bevHk6ceKE4/PRRx8pPT1dPXv2rHSs8vJyhYaGOpU/ceKEzpw5o1/96leNqouLi0ula7/ch+eHqrrPDb33jUFAAQA0u8jISHl5eelvf/tbpX19+/bViRMnnDqOvvvuu3JxcVFUVN37GXl4eOjSpUu1luvbt6++//57vffee45t58+f19mzZ3XNNdfU+Xw/FBUVpQceeED79+9XXFyco3/JoEGDdPr0afXq1avSx8PDo9JxBg0apOzsbLm5uVUqHxQU1Ki6dO3aVdnZ2U4h5cSJEw263pZAQAEANDtPT08tXbpUDz74oF544QX961//0tGjR7VlyxbNnDlTnp6emjVrlk6dOqVDhw5p4cKFuvPOO6udYbwqV111ld555x39+9//1rffflttucjISN1yyy2aO3eujhw5oo8++kh33HGHfvSjH+mWW26p13WVlJRowYIFOnz4sL744gu9++67OnbsmCPoLF26VH//+981f/58nThxwtH/ZeHChVUeb/z48Ro+fLimTp2qt99+W59//rlSU1P1m9/8xmlkUEPqMnr0aH3zzTdavXq1/vWvf+nZZ5/VW2+9Va/rbUl0kgWAdqAtjGh6+OGH5ebmpkceeURfffWVQkNDde+998rb21tvv/227r//fl1//fXy9vbWbbfdpjVr1tTr+I8++qjmzZunnj17ym63V/sqR5K2bt2q+++/X5MnT1ZpaalGjhypvXv31nvEi6urq86fP6+77rpLX3/9tYKCghQXF+dY4HbAgAFKSUnRQw89pBtvvFGGYahnz56aPn16lcezWCzau3evHnroIf385z/XN998I5vNppEjR9Ya1mqryzXXXKONGzdq1apVeuyxx3TbbbdpyZIlev755+t1zS3FYtT0T9CkCgoK5Ofnp/z8fNblAdohhhlX7+LFi8rIyFCPHj1qXLQVaC01PaP1+fubVzwAAMB0CCgAANSgX79+TkN+f/jZsWNHh61Lc6MPCgAANdi7d2+Vw3El1asTb3urS3MjoAAAUIOIiIjWroKDmerS3HjFAwAATIeAAgAATIeAAgAATIeAAgAATIeAAgAATIdRPADQHuRlSsXnW+583oGSf3jLnc9kDh8+rDFjxig3N1f+/v6tXZ12iYACAG1dXqb07I+lsuKWO6e7tzT//XqFlNmzZ2vbtm2O7wEBAbr++uu1evVqDRgwQJJkGIY2b96sLVu26PTp045Vfe+44w7dc8898vb2llQxZfpTTz2l3bt367PPPpO3t7euvvpq3X777Zo7d666dOnStNeLFkdAAYC2rvh8RTiJ2ywFtcD6RN+elXbPrThvPVtRJk2apK1bt0qSsrOz9Zvf/EaTJ0/WuXPnJEl33nmndu/erd/85jfasGGDunbtqo8++kjr1q3TVVddpalTp+q7777TDTfcoIKCAj322GMaPHiwPDw89Omnn2rnzp3auXOn5s+f3+SXjZZFQAGA9iIoSgq7trVrUSOr1SqbzSZJstlsWrp0qUaOHKlvvvlGhw4d0o4dO/Tqq6/qlltucfzmqquu0s0336yCggJJ0vLly3Xu3DmdOXNGP/rRjxzl+vTpo8mTJ9e4ivEPbdy4UWvXrlVmZqb8/Px044036s9//rOkipacp556Sn/4wx+UlZWlqKgoPfzww/rf//3fao+XmpqqX//61zp27JiCgoJ06623KjExUZ06dWpUXa666irFx8crPj7eUf7aa6/V1KlTlZCQIKliFeQ//OEPev3113Xw4EFFREToT3/6k7p27apf/OIXOnbsmAYMGKDt27erZ8+edbo/rY1OsgCAVlFUVKQdO3aoV69eCgwM1I4dO9S7d2+ncHKZxWKRn5+fysvL9fLLL+uOO+5wCidXlq3N8ePH9ctf/lKPPvqozpw5o3379mnkyJGO/b/5zW+0detWbdq0SadPn9YDDzygO+64QykpKVUe7+TJk5o4caLi4uL08ccf6+WXX9aRI0e0YMGCRtelrh577DHdddddOnHihPr06aMZM2Zo3rx5WrZsmY4fPy5JdaqPWdCCAgBoMW+88YY6d+4sSbpw4YJCQ0P1xhtvyMXFRenp6erdu3eNv//mm2+Ul5dXqdzgwYN15swZSdKUKVP00ksv1Xicc+fOqVOnTpo8ebJ8fHwUERGh6667zlGvNWvW6ODBgxo+fLgk6eqrr9aRI0f03HPPadSoUZWO99RTT2nGjBmOVo7IyEg988wzGjVqlDZt2iRPT88G1aU+7r77bk2bNk2StHTpUg0fPlwPP/ywJk6cKEm6//77dffdd9f7uK2FgAIAaDFjxozRpk2bJEnfffedNm7cqNjYWL3//vsyDKNOrR9S5VaSPXv2qLS0VEuXLlVJSUmtv4+JiVFERISuvvpqTZo0SZMmTdKtt94qb29v/eMf/9DFixcVExPj9JvS0tJqg0NaWpo+/fRTpxWFDcNQeXm5MjIydM011zSoLvVxuaOx9N+FA/v37++07eLFiyooKJCvr2+9jt0aCCgAgBbTqVMn9erVy/F98ODB8vPz0+bNmxUVFaVPPvmkxt937dpV/v7++uc//+m0vXv37pIkHx8f5eXl1VoPHx8fffDBBzp8+LD279+vRx55RAkJCTp27JjKy8slSW+++Wal10hWq7XK45WXl2vevHn65S9/WWnf5bo1pC7+/v5ycXGp1K+mqhWN3d3dHX++HOCq2nb5+syOPigAgFZjsVjk4uKikpISzZgxQ2fPntVf//rXSuUMw1B+fr5cXFw0bdo0bd++Xf/+978bdW43NzeNHz9eq1ev1scff6zPP/9cBw8eVN++fWW1WnXu3Dn16tXL6RMeXvWopUGDBun06dOVyvfq1UseHh4NrotUEcqysrIcZQsKCpSRkdGoa28LaEEBgPbi27OmP4/dbld2drYkKTc3Vxs2bFBRUZGmTJmiUaNGac+ePfrZz36mhx9+WDExMeratatOnjyptWvXauHChZo6dapWrVqlw4cPa+jQoXr00Uc1ZMgQderUSR9//LH+/ve/Kzo6utZ6vPHGG/rss880cuRIdenSRXv37lV5ebl69+4tHx8fLVmyRA888IDKy8sdQ5pTU1PVuXNnzZo1q9Lxli5dqmHDhmn+/PmaO3euOnXqpE8++UTJyclav359g+siSWPHjlVSUpKmTJmiLl266OGHH5arq2sD7n7bQkABgLbOO7Bi4rTdc1vunO7eFeetp3379ik0NFRSxauNPn366JVXXtHo0aMlSTt37tTzzz+vP/3pT1q5cqXc3NwUGRmpu+66y9HZMzAwUO+//76efPJJPfXUU8rIyJCLi4siIyM1ffp0p+G41fH399fu3buVkJCgixcvKjIyUi+99JL69esnqWJETHBwsBITE/XZZ5/J399fgwYN0vLly6s83oABA5SSkqKHHnpIN954owzDUM+ePTV9+vRG12XZsmX67LPPNHnyZPn5+emxxx7rEC0oFqOuA8ZNpKCgQH5+fsrPz28THX0A1M/a5Lr9H/oDMS0wKZnJXLx4URkZGerRo4fzyBCmuodJVPuMqn5/f9OCAgDtgX84gQHtCp1kAQDtzv/7f/9PnTt3rvbTUevSltCCAgBod4YMGaITJ060djUkmasubQkBBQDQ7nh5eTnNt9KazFSXtoRXPADQBrXB8Q3oIJrq2axXQNm0aZMGDBggX19f+fr6avjw4Xrrrbcc+2fPni2LxeL0GTZsmNMx7Ha7Fi5cqKCgIHXq1Ek333yzvvzyyya5GABo7y7PDFpcXNzKNQGqVlpaKkmNnqulXq94unXrpieeeMLRVLVt2zbdcsst+vDDDx3jtSdNmqStW7c6fnPlDHrx8fF6/fXXtWvXLgUGBmrx4sWaPHmy0tLSOsTEMwDQGK6urvL391dOTo4kydvbu87r1wDNrby8XN988428vb3l5ta4XiT1+vWUKVOcvj/++OPatGmTjh496ggoVqtVNputyt/n5+dry5YtevHFFzV+/HhJ0vbt2xUeHq4DBw44JuEBAFTv8n9jL4cUwExcXFzUvXv3RgfnBsebS5cu6ZVXXtGFCxccy1FL0uHDhxUcHCx/f3+NGjVKjz/+uIKDgyVVrPZYVlamCRMmOMqHhYUpOjpaqamp1QYUu90uu93u+F5QUNDQagNAm2exWBQaGqrg4OAqF40DWpOHh4dcXBrfxbXeAeXkyZMaPny4Ll68qM6dO2vPnj3q27evJCk2Nla33367IiIilJGRoYcfflhjx45VWlqarFarsrOz5eHhoS5dujgdMyQkxLE2Q1USExO1YsWK+lYVANo1V1dXXo2j3ap3QOndu7dOnDihvLw8/eUvf9GsWbOUkpKivn37Oq05EB0drSFDhigiIkJvvvmm4uLiqj2mYRg1NgUtW7ZMixYtcnwvKCiodkVJAADQ9tU7oHh4eDg6yQ4ZMkTHjh3T008/reeee65S2dDQUEVERCg9PV1SxXvT0tJS5ebmOrWi5OTkaMSIEdWe02q1ymq11reqAACgjWr0SyLDMJz6h/zQ+fPnlZmZ6Vi5cvDgwXJ3d1dycrKjTFZWlk6dOlVjQAEAAB1LvVpQli9frtjYWIWHh6uwsFC7du3S4cOHtW/fPhUVFSkhIUG33XabQkND9fnnn2v58uUKCgrSrbfeKkny8/PTnDlztHjxYgUGBiogIEBLlixR//79HaN6AAAA6hVQvv76a915553KysqSn5+fBgwYoH379ikmJkYlJSU6efKkXnjhBeXl5Sk0NFRjxozRyy+/LB8fH8cx1q5dKzc3N02bNk0lJSUaN26ckpKS6OgFAAAcLEYbnC+5oKBAfn5+ys/Pl6+vb2tXB0ATW5t8tk7lHoiJauaaAGhK9fn7m7V4AACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6RBQAACA6bi1dgUAdEB5mVLx+Wp3BxedU4m7vwqtthasFAAzIaAAaFl5mdKzP5bKiqstMlNSmYunXu+zWiXuXao/1lfFkneg5B/e9PUE0KoIKABaVvH5inASt1kKiqqyyO4jH2nKPx9U3D9+WfOxPpLk7i3Nf5+QArQzBBQArSMoSgq7tspdX3Tx1rZBr8irLK/GQ8zseVHaPbci9BBQgHaFgAKgSaxNPluncg/0q9vxCq222vugBFX/mghA28YoHgAAYDoEFAAAYDoEFAAAYDr0QQHQ9n1bS/8XhiIDbQ4BBUDb5R1YMcx499yayzEUGWhzCCgA2i7/8IrgUcOstPr2LEORgTaIgAKgbfMPJ3gA7VC9Oslu2rRJAwYMkK+vr3x9fTV8+HC99dZbjv2GYSghIUFhYWHy8vLS6NGjdfr0aadj2O12LVy4UEFBQerUqZNuvvlmffnll01zNQAAoF2oV0Dp1q2bnnjiCR0/flzHjx/X2LFjdcsttzhCyOrVq7VmzRpt2LBBx44dk81mU0xMjAoLCx3HiI+P1549e7Rr1y4dOXJERUVFmjx5si5dutS0VwagVfjYsxVc9M9qP7V2aAUA1fMVz5QpU5y+P/7449q0aZOOHj2qvn37at26dXrooYcUFxcnSdq2bZtCQkK0c+dOzZs3T/n5+dqyZYtefPFFjR8/XpK0fft2hYeH68CBA5o4cWITXRaA1uBjz9asD26Xe/nF6gtdXj/HO7DF6gWg7WlwH5RLly7plVde0YULFzR8+HBlZGQoOztbEyZMcJSxWq0aNWqUUlNTNW/ePKWlpamsrMypTFhYmKKjo5WamlptQLHb7bLb7Y7vBQUFDa02gKrkZdbc0VSq01Bdr7I8uZdf1FuRj+o77x5Vlpk5tDvDfgHUqt4B5eTJkxo+fLguXryozp07a8+ePerbt69SU1MlSSEhIU7lQ0JC9MUXX0iSsrOz5eHhoS5dulQqk52dXe05ExMTtWLFivpWFUBd5GVKz/64YoXhmtRjqO533j2U07lP1TvDql7BGAB+qN4BpXfv3jpx4oTy8vL0l7/8RbNmzVJKSopjv8VicSpvGEalbVeqrcyyZcu0aNEix/eCggKFh/N/X0CTKD5fEU7iNlesMFwVhuoCaGH1DigeHh7q1auXJGnIkCE6duyYnn76aS1dulRSRStJaGioo3xOTo6jVcVms6m0tFS5ublOrSg5OTkaMWJEtee0Wq2yWq31rSqA+giKksKube1atKo6r8gcQysQ0NwaPQ+KYRiy2+3q0aOHbDabkpOTdd1110mSSktLlZKSoieffFKSNHjwYLm7uys5OVnTpk2TJGVlZenUqVNavXp1Y6sCoIOpS6AILjqnmS1QFwBNq14BZfny5YqNjVV4eLgKCwu1a9cuHT58WPv27ZPFYlF8fLxWrVqlyMhIRUZGatWqVfL29taMGTMkSX5+fpozZ44WL16swMBABQQEaMmSJerfv79jVA+A9q2urRQAOrZ6BZSvv/5ad955p7KysuTn56cBAwZo3759iomJkSQ9+OCDKikp0X333afc3FwNHTpU+/fvl4+Pj+MYa9eulZubm6ZNm6aSkhKNGzdOSUlJcnV1bdorAwAAbVa9AsqWLVtq3G+xWJSQkKCEhIRqy3h6emr9+vVav359fU4NAAA6kHrNJAsAANASWCwQQN3VME19cNE5BRRntGBlALRnBBQAtfMOrJiobffcaotcHilT5uKpEnf/FqkWgPaLgAKgdv7hFbPI1jAd/o73zkmSStz9VWi1tVTNALRTBBQAdeMfXuMssjmdvVuwMgDaOzrJAgAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA0yGgAAAA02GYMdBBXV5VOLjonGaqYh6TqoYKPxAT1cI1AwBaUAAAgAkRUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkQUAAAgOkwzBhAx/Dt2Zr3ewe2TD0A1AkBBUC7VuLuL7l7S7vn1lzQ3Vs+A19WodVW6zEvzyFTE+aPARqHgAKgXSu02qT570vF56sv9O1ZafdceZXl1SmgAGh+BBQA7Z9/eMUHQJtBJ1kAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6jOIB0O7VNm9JcNE5zWyhugCoG1pQAACA6RBQAACA6RBQAACA6dAHBQD+I6A4o8b9Je7+TIUPtBACCoAOr8TdX2UunopNf6TGcmUunto26BVCCtACCCgAOrxCq03bBr0ir7K8assEFGcoNv0RFhQEWki9+qAkJibq+uuvl4+Pj4KDgzV16lSdOXPGqczs2bNlsVicPsOGDXMqY7fbtXDhQgUFBalTp066+eab9eWXXzb+agCggQqtNuV07lPt5zvvHq1dRaBDqVdASUlJ0fz583X06FElJyfr+++/14QJE3ThwgWncpMmTVJWVpbjs3fvXqf98fHx2rNnj3bt2qUjR46oqKhIkydP1qVLlxp/RQAAoM2r1yueffv2OX3funWrgoODlZaWppEjRzq2W61W2WxVN4Hm5+dry5YtevHFFzV+/HhJ0vbt2xUeHq4DBw5o4sSJ9b0GAADQzjRqmHF+fr4kKSAgwGn74cOHFRwcrKioKM2dO1c5OTmOfWlpaSorK9OECRMc28LCwhQdHa3U1NQqz2O321VQUOD0AQAA7VeDA4phGFq0aJFuuOEGRUdHO7bHxsZqx44dOnjwoH7/+9/r2LFjGjt2rOx2uyQpOztbHh4e6tKli9PxQkJClJ2dXeW5EhMT5efn5/iEh4c3tNoAAKANaPAongULFujjjz/WkSNHnLZPnz7d8efo6GgNGTJEERERevPNNxUXF1ft8QzDkMViqXLfsmXLtGjRIsf3goICQgoAAO1YgwLKwoUL9dprr+mdd95Rt27daiwbGhqqiIgIpaenS5JsNptKS0uVm5vr1IqSk5OjESNGVHkMq9Uqq9XakKoCaKTaFtoDgOZQr4BiGIYWLlyoPXv26PDhw+rRo/Zhd+fPn1dmZqZCQ0MlSYMHD5a7u7uSk5M1bdo0SVJWVpZOnTql1atXN+ASgA4sL1MqPl9zGe9AyZ8WRwBtS70Cyvz587Vz50799a9/lY+Pj6PPiJ+fn7y8vFRUVKSEhATddtttCg0N1eeff67ly5crKChIt956q6PsnDlztHjxYgUGBiogIEBLlixR//79HaN6ANRBXqb07I+lsuKay7l7S/PfJ6QAaFPqFVA2bdokSRo9erTT9q1bt2r27NlydXXVyZMn9cILLygvL0+hoaEaM2aMXn75Zfn4+DjKr127Vm5ubpo2bZpKSko0btw4JSUlydXVtfFXBHQUxecrwkncZikoquoy356Vds+tKEtAAdCG1PsVT028vLz09ttv13ocT09PrV+/XuvXr6/P6QFUJShKCru2tWsBAE2qUfOgAAAANAcCCgAAMB0CCgAAMB0CCgAAMJ0GzyQLwJx+OLFacNE5zZS0471zyuns3XqVakcCijNqLVPi7i+pmpFVAOqEgAIAdVDi7q8yF0/Fpj9Sa9kyF0/phuMM7QYagYACAHVQaLVp26BX5FWWV2O5gOKMihDD3DNAoxBQAKCOCq02FVptrV0NoEOgkywAADAdAgoAADAdAgoAADAdAgoAADAdAgoAADAdRvEAHUBNk4vVZeIxAGhpBBSgHavr5GJlLp7/mf0UAMyBgAK0Y3WdXKzE3Z/5PZrYW4dS9J33uWr3X77nD8QwJT5QFQIK0M4xuVjLqk+r1bZBr4g1e4CqEVAAoAnVpdXq8nT4XmV5Tos71oSWFnQ0BBQAaGK0WgGNxzBjAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOgQUAABgOvUKKImJibr++uvl4+Oj4OBgTZ06VWfOnHEqYxiGEhISFBYWJi8vL40ePVqnT592KmO327Vw4UIFBQWpU6dOuvnmm/Xll182/moAAEC7UK+AkpKSovnz5+vo0aNKTk7W999/rwkTJujChQuOMqtXr9aaNWu0YcMGHTt2TDabTTExMSosLHSUiY+P1549e7Rr1y4dOXJERUVFmjx5si5dutR0VwYAANost/oU3rdvn9P3rVu3Kjg4WGlpaRo5cqQMw9C6dev00EMPKS4uTpK0bds2hYSEaOfOnZo3b57y8/O1ZcsWvfjiixo/frwkafv27QoPD9eBAwc0ceLEJro0AADQVjWqD0p+fr4kKSAgQJKUkZGh7OxsTZgwwVHGarVq1KhRSk1NlSSlpaWprKzMqUxYWJiio6MdZa5kt9tVUFDg9AEAAO1XgwOKYRhatGiRbrjhBkVHR0uSsrOzJUkhISFOZUNCQhz7srOz5eHhoS5dulRb5kqJiYny8/NzfMLDwxtabQAA0AY0OKAsWLBAH3/8sV566aVK+ywWi9N3wzAqbbtSTWWWLVum/Px8xyczM7Oh1QYAAG1AgwLKwoUL9dprr+nQoUPq1q2bY7vNZpOkSi0hOTk5jlYVm82m0tJS5ebmVlvmSlarVb6+vk4fAADQftUroBiGoQULFmj37t06ePCgevTo4bS/R48estlsSk5OdmwrLS1VSkqKRowYIUkaPHiw3N3dncpkZWXp1KlTjjIAAKBjq9convnz52vnzp3661//Kh8fH0dLiZ+fn7y8vGSxWBQfH69Vq1YpMjJSkZGRWrVqlby9vTVjxgxH2Tlz5mjx4sUKDAxUQECAlixZov79+ztG9QAAgI6tXgFl06ZNkqTRo0c7bd+6datmz54tSXrwwQdVUlKi++67T7m5uRo6dKj2798vHx8fR/m1a9fKzc1N06ZNU0lJicaNG6ekpCS5uro27moAAEC7UK+AYhhGrWUsFosSEhKUkJBQbRlPT0+tX79e69evr8/pAQBAB8FaPAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHQIKAAAwHTqNcwYQAvJy5SKz1favOO9c44/BxRnKPY/23I6e7dg5QCg+RFQALPJy5Se/bFUVlxp18wrvpe5eKrE3b9FqgUALYmAApjA2uSzjj8HF/1TM8uK9Vbko/rOu0cNv5JK3P1VaLU1d/UAoMURUACT+s67h3I692ntagBAq6CTLAAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB0CCgAAMB3mQQFaWhXT2AcXOU9hDwAdHQEFaEnVTGPPFPYA4IyAArSk4vMV4SRusxQU5dj8w0UAJaawBwACCtAagqKksGsdX1mNGACc0UkWAACYDgEFAACYDgEFAACYDn1QAKANWJt8tk7lHoiJqr0Q0AbQggIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyHgAIAAEyn3gHlnXfe0ZQpUxQWFiaLxaJXX33Vaf/s2bNlsVicPsOGDXMqY7fbtXDhQgUFBalTp066+eab9eWXXzbqQgAAQPtR75lkL1y4oIEDB+ruu+/WbbfdVmWZSZMmaevWrY7vHh4eTvvj4+P1+uuva9euXQoMDNTixYs1efJkpaWlydXVtb5VAoA2KaA4o8b9Je7+KrTaWqg2gLnUO6DExsYqNja2xjJWq1U2W9X/UuXn52vLli168cUXNX78eEnS9u3bFR4ergMHDmjixIn1rRIAtCkl7v4qc/FUbPojNZYrc/HUtkGvEFLQITXLWjyHDx9WcHCw/P39NWrUKD3++OMKDg6WJKWlpamsrEwTJkxwlA8LC1N0dLRSU1MJKADavUKrTdsGvSKvsrxqywQUZyg2/RH9KP9Dfefdo8bj0dKC9qjJA0psbKxuv/12RUREKCMjQw8//LDGjh2rtLQ0Wa1WZWdny8PDQ126dHH6XUhIiLKzs6s8pt1ul91ud3wvKCho6moDQIsqtNpqDBV1bWWRaGlB+9TkAWX69OmOP0dHR2vIkCGKiIjQm2++qbi4uGp/ZxiGLBZLlfsSExO1YsWKpq4qAJhWXVpZpP+2tHiV5RFQ0K40yyueHwoNDVVERITS09MlSTabTaWlpcrNzXVqRcnJydGIESOqPMayZcu0aNEix/eCggKFh4c3b8UBoJXV1soCtGfNPg/K+fPnlZmZqdDQUEnS4MGD5e7uruTkZEeZrKwsnTp1qtqAYrVa5evr6/QBAADtV71bUIqKivTpp586vmdkZOjEiRMKCAhQQECAEhISdNtttyk0NFSff/65li9frqCgIN16662SJD8/P82ZM0eLFy9WYGCgAgICtGTJEvXv398xqgcA0LzWJp+ttcwDMVEtUBOgavUOKMePH9eYMWMc3y+/epk1a5Y2bdqkkydP6oUXXlBeXp5CQ0M1ZswYvfzyy/Lx8XH8Zu3atXJzc9O0adNUUlKicePGKSkpiTlQAACApAYElNGjR8swjGr3v/3227Uew9PTU+vXr9f69evre3oAANABsBYPAAAwHQIKAAAwHQIKAAAwnWafBwXoyK4cKRFcdE4zJe1475xyOnu3TqUAoA2gBQUAAJgOAQUAAJgOr3iAppSXKRWfd3wNLjrntDugOKOlawQAbRIBBWgqeZnSsz+Wyoodm2ZWUazMxVMl7v4tVi0AaIsIKEBTKT5fEU7iNktBFVOE73jvXKViJe7+LAAHALUgoABNLShKCrtWkhipAwANRCdZAABgOrSgAEA74OiA/VVx1QW8AyX/8JarENBIBBQAaMNK3P1V5uKp2PRHKjZ8VE1Bd29p/vuEFLQZBBQAaMMKrTZtG/SKvMryJEkzh3avXOjbs9LuuRUduQkoaCMIKADQxhVabf8dGRYW1bqVAZoInWQBAIDpEFAAAIDpEFAAAIDp0AcFANqRtclnK20LLjqnmaqY2ZjJA9FW0IICAABMh4ACAABMh4ACAABMh4ACAABMh4ACAABMh4ACAABMh4ACAABMh4ACAABMh4ACAABMh4ACAABMh4ACAABMh4ACAABMh8UCgbrKy5SKz0uqWHTtSgHFGYoVC7IBQFMgoAB1kZcpPftjqaxYkjSzmmJlLp4qcfdvsWoBQHtFQAHqovh8RTiJ2ywFRVXZgiJJJe7+KrTaWrhyQN0EFGfUuJ/nF2ZS74Dyzjvv6KmnnlJaWpqysrK0Z88eTZ061bHfMAytWLFCzz//vHJzczV06FA9++yz6tevn6OM3W7XkiVL9NJLL6mkpETjxo3Txo0b1a1btya5KKDZBEVJYdfyCgdtSom7v8pcPBWb/kiN5cpcPLVt0CuEFJhCvQPKhQsXNHDgQN1999267bbbKu1fvXq11qxZo6SkJEVFRWnlypWKiYnRmTNn5OPjI0mKj4/X66+/rl27dikwMFCLFy/W5MmTlZaWJldX18ZfFQDAodBq07ZBr8irLK/aMgHFGYpNf0ReZXkEFJhCvQNKbGysYmNjq9xnGIbWrVunhx56SHFxcZKkbdu2KSQkRDt37tS8efOUn5+vLVu26MUXX9T48eMlSdu3b1d4eLgOHDigiRMnNuJyAABVKbTaCB5oU5p0mHFGRoays7M1YcIExzar1apRo0YpNTVVkpSWlqaysjKnMmFhYYqOjnaUuZLdbldBQYHTBwAAtF9NGlCys7MlSSEhIU7bQ0JCHPuys7Pl4eGhLl26VFvmSomJifLz83N8wsPDm7LaAADAZJplojaLxeL03TCMStuuVFOZZcuWKT8/3/HJzMxssroCAADzadKAYrNVvN+8siUkJyfH0apis9lUWlqq3NzcastcyWq1ytfX1+kDAADaryYNKD169JDNZlNycrJjW2lpqVJSUjRixAhJ0uDBg+Xu7u5UJisrS6dOnXKUAQAAHVu9R/EUFRXp008/dXzPyMjQiRMnFBAQoO7duys+Pl6rVq1SZGSkIiMjtWrVKnl7e2vGjBmSJD8/P82ZM0eLFy9WYGCgAgICtGTJEvXv398xqgcAAHRs9Q4ox48f15gxYxzfFy1aJEmaNWuWkpKS9OCDD6qkpET33XefY6K2/fv3O+ZAkaS1a9fKzc1N06ZNc0zUlpSUxBwoAABAUgMCyujRo2UYRrX7LRaLEhISlJCQUG0ZT09PrV+/XuvXr6/v6QEAQAfAWjzAD6xNPlvl9uCic5opVioGgJbSLMOMAQAAGoOAAgAATIeAAgAATIeAAgAATIeAAgAATIeAAgAATIeAAgAATIeAAgAATIeAAgAATIeZZIG8TKn4vKSKGWOrElCc0ZI1AoAOj4CCji0vU3r2x1JZsSRpZg1Fy1w8VeLu3yLVAoCOjoCCjq34fEU4idssBUVpx3tVt6BIUom7vwqtthasHAB0XAQUQJKCoqSwa1kIEPiB6hbPvNIDMVHNXBN0RHSSBQAApkNAAQAApkNAAQAApkNAAQAApkNAAQAApkNAAQAApsMwYwCAQ22zJjMfEFoKAQUAoBJ3f5W5eCo2/ZEay5W5eGrboFcIKWh2BBQAgAqtNm0b9Iq8yvKqLRNQnKHY9EfkVZZHQEGzI6AAACRVhBSCB8yCTrIAAMB0CCgAAMB0eMWD9i0vs2LF4up8W7fF0AAALYuAgvYrL1N69sdSWXHN5dy9Je/AlqkTAKBOCChov4rPV4STuM1SUA3LwXsHSv7hLVcvAECtCCho/4KipLBrW7sWQLtRaTK3r6popST4o5EIKACAOql2MrePqijs7i3Nf5+QggYjoAAA6qS6ydxmDu3uXPDbs9LuuRWvWQkoaCACCgCgzqqczC2shj5eQAMxDwoAADCdJg8oCQkJslgsTh+b7b9p2zAMJSQkKCwsTF5eXho9erROnz7d1NUAAABtWLO84unXr58OHDjg+O7q6ur48+rVq7VmzRolJSUpKipKK1euVExMjM6cOSMfH5/mqA4gSVqbzKRsANBWNEtAcXNzc2o1ucwwDK1bt04PPfSQ4uLiJEnbtm1TSEiIdu7cqXnz5jVHdQAAJlDX/0l4IIY+LWimgJKenq6wsDBZrVYNHTpUq1at0tVXX62MjAxlZ2drwoQJjrJWq1WjRo1SampqtQHFbrfLbrc7vhcUFDRHtdHWMI09YApXBo/gonOaKWnHe+eU09m7dSqFNq/JA8rQoUP1wgsvKCoqSl9//bVWrlypESNG6PTp08rOzpYkhYSEOP0mJCREX3zxRbXHTExM1IoVK5q6qmjLmMYeML1KE7pdocTdv/KIIOA/mjygxMbGOv7cv39/DR8+XD179tS2bds0bNgwSZLFYnH6jWEYlbb90LJly7Ro0SLH94KCAoWHM7a+Q2Mae8C0qp3Q7QplLp7aNugVQgqq1OzzoHTq1En9+/dXenq6pk6dKknKzs5WaGioo0xOTk6lVpUfslqtslqtzV1VtEVMYw+YTnUTuv1QQHGGYtMfkVdZHgEFVWr2gGK32/XJJ5/oxhtvVI8ePWSz2ZScnKzrrrtOklRaWqqUlBQ9+eSTzV0VAEALqXJCN6AemjygLFmyRFOmTFH37t2Vk5OjlStXqqCgQLNmzZLFYlF8fLxWrVqlyMhIRUZGatWqVfL29taMGTOauioAAKCNavKA8uWXX+pnP/uZvv32W3Xt2lXDhg3T0aNHFRERIUl68MEHVVJSovvuu0+5ubkaOnSo9u/fzxwoAADAockDyq5du2rcb7FYlJCQoISEhKY+NTooJmADgPaHtXgAAIDpsJoxzKe2CdgkJmEDgHaOgAJzqesEbNIPJmErafZqAQBaFgEF5lLXCdikH0zCRmsKALQ3BBSYExOwAUCHRkBBi6tp1M0PFxmbeeu1LVYnAIC5EFAAAK2mygUFv/pBHzTW0+qwCCgwLeY3AdqvGhcU/Oi/f6xpQcEHYmrpp4Y2jYACAGhxLCiI2hBQAACtggUFURNmkgUAAKZDQAEAAKZDQAEAAKZDHxS0rLxMBRf9s9rdVQ45BIAq1GWkHyN92i4CClrOf9bZmVnLOjtlLp4qcfdvmToBAEyJgIKW8591dt6KfFTfefeotliJuz89+wGggyOgoMV9591DOZ37tHY1AAAmRidZAABgOrSgAABMrS6d53k13P4QUAAAplTjej1XqG7Nnrqu6cVoH/MhoAAATKku6/VIrNnTXhFQAACmxXo9HRcBBU2mtqbU4KJzmtlCdQEAtG2M4gEAAKZDCwoAoF2obbQPI33aFgIKmoyPPbvGzmysswOgOdR1tE91I30kRvuYEQEFdVLbv7w+9mzN+uB2uZdfrLEc6+wAaGp1Ge3DSJ+2h4CCJuFVlif38ousswOgVdR1tA+vgdoOAgqaFOvsADCjpngNhJZFQEGd0L8EQFvGa6C2h4CC2uVl0r8EQJvHpG9tCwGlo8vLlIrP11zm27P0LwEAMdqnJRFQOrK8TOnZH0tlxbUWLXPx1L/9riOAAOjwanvlLUn6qljyDpT8w1ukTu0RAaUjKz5fEU7iNktBNaf9bR8UEE4AdAg19anzKsvVlH8+WOsrb30kyd1bmv8+IaWBWjWgbNy4UU899ZSysrLUr18/rVu3TjfeeGNrVqlD2vEvT+V87V1zIWst+wGgjavPSJ/dfZ9RiXuXastc7nC749CHyulcUm05XgVVr9UCyssvv6z4+Hht3LhRP/nJT/Tcc88pNjZW//jHP9S9e/fWqlarufxes7amw5lD/3Nv6tJ0WFv/km/r9i4VADqCuoz0kerX367WEY55XrX+t7yj9nuxGIZhtMaJhw4dqkGDBmnTpk2Obddcc42mTp2qxMTEGn9bUFAgPz8/5efny9fXt8nr1qCHoS6dTWuw471zdW86lCqaDqe/KHkHVbl795GP6nQsxvwDQNOrz+zaTfXf4LoGlNYMPPX5+7tVWlBKS0uVlpamX//6107bJ0yYoNTU1Erl7Xa77Ha743t+fr6kigttDhcvFNWpnOP8eV9Km8dI31ffjFebKZePabHqragnaxyu61WWp9izD8t9S1y1ZcbX8VgX3f1U9H1n6fu6XTMAoHYX1Vl/6JMkz7L8assEFH+hmM8el5H/lS526tzocya++kGjj/FDzfF37OVj1qVtpFUCyrfffqtLly4pJCTEaXtISIiys7MrlU9MTNSKFSsqbQ8Pb92OR8ub5aiFkv7PhMcCADSPn7V2BarUPH/HVSgsLJSfn1+NZVq1k6zFYnH6bhhGpW2StGzZMi1atMjxvby8XN99950CAwOrLF8XBQUFCg8PV2ZmZrO8JkJl3POWxf1uedzzlsX9bllNcb8Nw1BhYaHCwsJqLdsqASUoKEiurq6VWktycnIqtapIktVqldVqddrm7+/fJHXx9fXlwW5h3POWxf1uedzzlsX9blmNvd+1tZxc5tLgMzSCh4eHBg8erOTkZKftycnJGjFiRGtUCQAAmEirveJZtGiR7rzzTg0ZMkTDhw/X888/r3Pnzunee+9trSoBAACTaLWAMn36dJ0/f16PPvqosrKyFB0drb179yoiIqJFzm+1WvXb3/620qsjNB/uecvifrc87nnL4n63rJa+3602DwoAAEB1WqUPCgAAQE0IKAAAwHQIKAAAwHQIKAAAwHTafUBJSEiQxWJx+ths/12UyTAMJSQkKCwsTF5eXho9erROnz7dijVuW9555x1NmTJFYWFhslgsevXVV5321+X+2u12LVy4UEFBQerUqZNuvvlmffnlly14FW1Hbfd79uzZlZ73YcOGOZXhftddYmKirr/+evn4+Cg4OFhTp07VmTNnnMrwjDetutxznvOms2nTJg0YMMAx+drw4cP11ltvOfa35vPd7gOKJPXr109ZWVmOz8mTJx37Vq9erTVr1mjDhg06duyYbDabYmJiVFhY2Io1bjsuXLiggQMHasOGDVXur8v9jY+P1549e7Rr1y4dOXJERUVFmjx5si5dutRSl9Fm1Ha/JWnSpElOz/vevXud9nO/6y4lJUXz58/X0aNHlZycrO+//14TJkzQhQsXHGV4xptWXe65xHPeVLp166YnnnhCx48f1/HjxzV27FjdcsstjhDSqs+30c799re/NQYOHFjlvvLycsNmsxlPPPGEY9vFixcNPz8/4w9/+EML1bD9kGTs2bPH8b0u9zcvL89wd3c3du3a5Sjz73//23BxcTH27dvXYnVvi66834ZhGLNmzTJuueWWan/D/W6cnJwcQ5KRkpJiGAbPeEu48p4bBs95c+vSpYvxxz/+sdWf7w7RgpKenq6wsDD16NFDP/3pT/XZZ59JkjIyMpSdna0JEyY4ylqtVo0aNUqpqamtVd12oy73Ny0tTWVlZU5lwsLCFB0dzT+DBjp8+LCCg4MVFRWluXPnKicnx7GP+904+fn5kqSAgABJPOMt4cp7fhnPedO7dOmSdu3apQsXLmj48OGt/ny3+4AydOhQvfDCC3r77be1efNmZWdna8SIETp//rxjscIrFygMCQmptJAh6q8u9zc7O1seHh7q0qVLtWVQd7GxsdqxY4cOHjyo3//+9zp27JjGjh0ru90uifvdGIZhaNGiRbrhhhsUHR0tiWe8uVV1zyWe86Z28uRJde7cWVarVffee6/27Nmjvn37tvrz3WpT3beU2NhYx5/79++v4cOHq2fPntq2bZujU5XFYnH6jWEYlbah4Rpyf/ln0DDTp093/Dk6OlpDhgxRRESE3nzzTcXFxVX7O+537RYsWKCPP/5YR44cqbSPZ7x5VHfPec6bVu/evXXixAnl5eXpL3/5i2bNmqWUlBTH/tZ6vtt9C8qVOnXqpP79+ys9Pd0xmufKlJeTk1MpMaL+6nJ/bTabSktLlZubW20ZNFxoaKgiIiKUnp4uifvdUAsXLtRrr72mQ4cOqVu3bo7tPOPNp7p7XhWe88bx8PBQr169NGTIECUmJmrgwIF6+umnW/357nABxW6365NPPlFoaKh69Oghm82m5ORkx/7S0lKlpKRoxIgRrVjL9qEu93fw4MFyd3d3KpOVlaVTp07xz6AJnD9/XpmZmQoNDZXE/a4vwzC0YMEC7d69WwcPHlSPHj2c9vOMN73a7nlVeM6blmEYstvtrf98N6qLbRuwePFi4/Dhw8Znn31mHD161Jg8ebLh4+NjfP7554ZhGMYTTzxh+Pn5Gbt37zZOnjxp/OxnPzNCQ0ONgoKCVq5521BYWGh8+OGHxocffmhIMtasWWN8+OGHxhdffGEYRt3u77333mt069bNOHDggPHBBx8YY8eONQYOHGh8//33rXVZplXT/S4sLDQWL15spKamGhkZGcahQ4eM4cOHGz/60Y+43w30f//3f4afn59x+PBhIysry/EpLi52lOEZb1q13XOe86a1bNky45133jEyMjKMjz/+2Fi+fLnh4uJi7N+/3zCM1n2+231AmT59uhEaGmq4u7sbYWFhRlxcnHH69GnH/vLycuO3v/2tYbPZDKvVaowcOdI4efJkK9a4bTl06JAhqdJn1qxZhmHU7f6WlJQYCxYsMAICAgwvLy9j8uTJxrlz51rhasyvpvtdXFxsTJgwwejatavh7u5udO/e3Zg1a1ale8n9rruq7rUkY+vWrY4yPONNq7Z7znPetH7+858bERERhoeHh9G1a1dj3LhxjnBiGK37fFsMwzAa1wYDAADQtDpcHxQAAGB+BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6BBQAAGA6/x+At/Pa5ZWtgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.hist(control_sele_sum, bins = bins_control_sele, alpha=0.5, label='control_sele_sum')\n",
    "pyplot.hist(BCG_sele_sum, bins = bins_BCG_sele, alpha=1, label='BCG_sele_sum', histtype=u'step')\n",
    "pyplot.legend(loc='upper right')\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "d1b1f5fa-4c2e-4a15-9a53-2fecb71c7b07",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.3417217995921577e-07"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U1, p = mannwhitneyu(control_sele_sum, BCG_sele_sum)\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71f3a765-1625-4226-aaf0-399e1f911cc0",
   "metadata": {},
   "source": [
    "Similar as in full set, H0 is < 0.01 so data doesn't comes from same distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6ded118-1381-4ccf-9ed0-0ab979fde688",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}