ANOVA/ANOVA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import f as f_test\n",
    "from scipy.stats import ncf\n",
    "from scipy.stats import t as t_test\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set_style('darkgrid')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ANOVA\n",
    "\n",
    "* Assumptions\n",
    "\n",
    "1. Independence of the observations \n",
    "\n",
    "2. No significant outliers\n",
    "\n",
    "3. Normality\n",
    "\n",
    "4. Homogeneity of variances : The variance among the group should be approximately equal\n",
    "\n",
    "It is important to note that ANOVA is not robust to violations to the assumption of independence. This is to say, that even if you violate the assumptions of homogeneity or normality, you can conduct the test and basically trust the findings. However, the results of the ANOVA are invalid if the independence assumption is violated. In general, with violations of homogeneity the analysis is considered robust if you have equal sized groups. With violations of normality, continuing with the ANOVA is generally ok if you have a large sample size.\n",
    "\n",
    "* When to use ANOVA\n",
    "\n",
    "1. No matter what distribution the population is, if the sample size is over 30, the sameple means approach normal distribution. If we want to compare multiple sample means, we use ANOVA.\n",
    "\n",
    "2. Given the population is normal distribution, the sample means is thus normal distribution. If we want to compare multiple sample means, we use ANOVA.\n",
    "\n",
    "---\n",
    "\n",
    "### One-Way ANOVA\n",
    "\n",
    "![](./resources/1_way_SS.PNG)\n",
    "\n",
    "![](./resources/1_way_Fi.PNG)\n",
    "\n",
    "[REFERENCE](https://www.alexanderdemos.org/ANOVA5.html)\n",
    "\n",
    "---\n",
    "\n",
    "### Two-Way ANOVA\n",
    "\n",
    "We now have two factors, and we will test the main effect and interaction\n",
    "\n",
    "* Main effect (Factor A: Rows): The independent effects of Length of Exposure on Perceived Ability\n",
    "\n",
    "* Main effect (Factor B: Columns): The independent effects of Type of Exposure on Perceived Ability\n",
    "\n",
    "* Interaction: The joint effects of Length and Type of Exposure on Perceived Ability\n",
    "\n",
    "![](./resources/2_way_SS.PNG)\n",
    "\n",
    "![](./resources/2_way_F.PNG)\n",
    "\n",
    "[REFERENCE](https://www.alexanderdemos.org/ANOVA8.html)\n",
    "\n",
    "---\n",
    "\n",
    "### Main Effects vs Interaction\n",
    "\n",
    "If the interaction is telling a conflicting story with the main effects, you may either not want to follow up the main effects or carefully explain them.\n",
    "\n",
    "![](./resources/Interaction.PNG)\n",
    "\n",
    "The above graph shows there is an interaction between (Watching, Reading, Thinking) and (Low, High). It means that even if we conclude there is significant difference in group (Low, High) from 2-way ANOVA. This conclusion does not always correct in different group (Watching, Reading, Thinking).\n",
    "\n",
    "[REFERENCE](https://www.alexanderdemos.org/ANOVA9.html#Basic_Rules)\n",
    "\n",
    "---\n",
    "\n",
    "### Effect Size\n",
    "\n",
    "The effect size tells us how much do the data differ from the hypothosis. SSB is the variation caused by treatment, and SST is the variation of Treatment + Noise.\n",
    "\n",
    "![](./resources/effect_size.png)\n",
    "\n",
    "| d | Effect Size |\n",
    "| --- | --- |\n",
    "| 0.14 | large |\n",
    "| 0.06 | medium |\n",
    "| 0.01 | small |\n",
    "\n",
    "---\n",
    "\n",
    "### Power Analysis\n",
    "\n",
    "Power analysis answers questions like \"how much statistical power does my study have?\" and \"how big a sample size do I need?\".\n",
    "\n",
    "![](./resources/power_analysis.PNG)\n",
    "\n",
    "If H0 is not rejected, F score is in f-distribution.\n",
    "\n",
    "If H0 is rejected, F score is in non-central f-distribution\n",
    "\n",
    "![](./resources/f_distribution.PNG)\n",
    "\n",
    "![](./resources/non_central_parameter.PNG) \n",
    "![](./resources/parameter.PNG)\n",
    "\n",
    "We use MSW to estimate the variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "### Post Hoc Test\n",
    "\n",
    "ANOVA can only know whether the mean of three or more groups are the same. However, it can not tell us how one group differs from a certain group. We can not use t-test in this situation because the type I error will increase if we use t-test on the same group multiple times. Instead, we use post hoc test which is a more rigorous test than t-test.\n",
    "\n",
    "![](./resources/FWER.PNG)\n",
    "\n",
    "Tukey HSD(Honestly Sinificant Difference): This test determines whether there is a difference between the mean of all possible pairs. \n",
    "\n",
    "![](./resources/choose_post_hoc.jpg)\n",
    "\n",
    "[REFERENCE 1](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6193594/)\n",
    "\n",
    "[REFERENCE 2](https://www.youtube.com/watch?v=u9h8qjxMV-Y&t=185s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "### Question 1\n",
    "\n",
    "From our survey we can calculate Caffeiene consumption as a function of how much you\n",
    "said you liked math. Here is a table of statistics based on our survey:\n",
    "\n",
    "|| Not at all | Just a little | A fair amount | Very much |\n",
    "|---|---|---|---|---|\n",
    "| n | 14 | 56 | 75 | 6 |\n",
    "|mean | 1.21 | 1.01 | 1.41 | 3 |\n",
    "|SS | 7.8574 | 41.7456 | 354.1363 | 108 |\n",
    "| ntotal | 151 |\n",
    "|grand mean | 1.3065 |\n",
    "|SStotal | 534.834 |\n",
    "\n",
    "Calculate the standard errors of the mean for each of the 4 groups.\n",
    "\n",
    "Make a bar graph of the means for each of the 4 groups with error bars as the standard\n",
    "error of the means.\n",
    "\n",
    "Using an alpha value of α = 0.05, is there difference in Caffeine consumption across the 4 groups of students who vary their preference for math?\n",
    "\n",
    "* from [one_factor_ANOVA_tutorial.pdf](./resources/one_factor_ANOVA_tutorial.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = np.array([14, 56, 75, 6])\n",
    "n_total = np.sum(n)\n",
    "mean = np.array([1.21, 1.01, 1.41, 3])\n",
    "SSW_g = np.array([7.8574, 41.7456, 354.1363, 108])\n",
    "mu = np.average(mean, weights=n)\n",
    "SST = 534.834\n",
    "alpha = 0.05\n",
    "\n",
    "SSW = np.sum(SSW_g)\n",
    "SSB = np.sum( np.square(mean - mu) * n )\n",
    "\n",
    "k = 4\n",
    "df_SSW = n_total - k\n",
    "df_SSB = k - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "H<sub>0 : There is no difference among group means\n",
    "\n",
    "H<sub>1 : At lease one group differs significantly from the overrall mean of the dependent variable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F:  2.208471453651477\n",
      "F_crit:  2.6661487758521765\n",
      "p_value:  0.08956999737639892\n",
      "eta_squared:  0.04312463900725823\n"
     ]
    }
   ],
   "source": [
    "MSW = SSW / df_SSW\n",
    "MSB = SSB / df_SSB\n",
    "f = MSB / MSW\n",
    "print(\"F: \", f)\n",
    "\n",
    "# dfn: numerator, dfd: denominator\n",
    "f_crit = f_test.ppf(1 - alpha, dfn=df_SSB, dfd=df_SSW)\n",
    "print(\"F_crit: \", f_crit)\n",
    "\n",
    "p_value = 1 - f_test.cdf(f, dfn=df_SSB, dfd=df_SSW)\n",
    "print(\"p_value: \", p_value)\n",
    "\n",
    "eta_squared = SSB / SST \n",
    "print(\"eta_squared: \", eta_squared)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x267d5374970>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAAD7CAYAAACoomWyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUdUlEQVR4nO3de3BU9d3H8c+ym6sElktoRQpNcCgWFUE60El9oKJIcKINEG5xqdLakbZPxqLcQoaJT0EMCFTocAk1BUJqIBYslCn1xkDLKMOleCsSCQRNSBEIISTkttl9/qCuBDbZEDa7P9z36581m7PnfPec5Z3DYTda3G63WwAAI3UI9gAAgOYRaQAwGJEGAIMRaQAwGJEGAIMRaQAwmK01C61du1bvvvuuGhoaNHnyZKWkpLT3XAAAtSLS+/fv17/+9S+99tprqqmpUU5OTiDmAgBIsvj6MMvSpUtlsVj02WefqaqqSrNmzdI999wTqPkAIKT5PJO+cOGCTp8+rTVr1qikpETTp0/Xrl27ZLFYvC7vdrvFZxgB4MZ06OC9qT4jbbfbFR8fr/DwcMXHxysiIkLl5eXq1q2b1+WdTpcqKi7f3LQAEGJiY2O83u/z3R3333+//vGPf8jtduvMmTOqqamR3W7393wAAC98nkn/+Mc/1oEDBzR+/Hi53W7Nnz9fVqs1ELMBQMjz+Q+HN6qhoZHLHQBwg9p8uQMAEDxEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMRqQBwGBEGgAMZmvNQj/5yU8UExMjSerVq5cWLVrUrkMBAK7wGem6ujpJUm5ubrsPAwBoyufljk8//VQ1NTWaNm2apk6dqiNHjgRgLACQDh8+qBdemKfDhw8Ge5Sg8XkmHRkZqZ/97GdKSUlRcXGxnn76ae3atUs2W6uulABAmxUU/EknT55QbW2NBg8eEuxxgsJnaePi4tSnTx9ZLBbFxcXJbrfr7Nmzuv32270ub7VaZLdH+31QAKGnvr7OcxuqXfEZ6ddff12FhYXKzMzUmTNnVFVVpdjY2GaXb2x0q6Lisl+HBBCaGhvdnttveldiY2O83u8z0uPHj9fcuXM1efJkWSwWvfjii1zqAIAA8Vnb8PBwLV26NBCzAACuwYdZAMBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgrYr0+fPnNXz4cBUVFbX3PACAq/iMdENDg+bPn6/IyMhAzAMAuIrPSGdlZWnSpEnq0aNHIOYBAFzF1tI3t27dqq5du+qBBx5QdnZ2q1ZotVpkt0f7ZTgAoc1qtXhuQ7UrLUb6z3/+sywWi9577z0dPXpUs2fP1urVqxUbG9vsYxob3aqouOz3QQGEnsZGt+f2m96V2NgYr/e3GOm8vDzPfzscDmVmZrYYaACAf/EWPAAwWItn0lfLzc1tzzkAAF5wJg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwIg0ABiPSAGAwm68FGhsblZGRoZMnT8pqtWrRokXq3bt3IGYDgJDn80x69+7dkqT8/HylpaVp0aJF7T4UAOAKn2fSDz30kEaMGCFJOn36tLp3797eMwEwVJeOYbJFRQZse1arxXMbGxsTsO1KkrOmVheqGgK6TW98RlqSbDabZs+erbfeeksrVqxocVmr1SK7PdovwwEwiy3Mqj3/Mzxg26uxWSWLRTUlJQHdriQN37tHdltYQLfpTasiLUlZWVl6/vnnNWHCBO3cuVPR0d5D3NjoVkXFZb8NCMAcgT6bDbZAtqy5fevzmvQbb7yhtWvXSpKioqJksVhktVr9Ox0AwCufZ9KjRo3S3LlzlZqaKqfTqfT0dEVERARiNgAIeT4jHR0drVdeeSUQswAArsGHWQDAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxGpAHAYEQaAAxma+mbDQ0NSk9PV2lpqerr6zV9+nSNHDkyULMBQMhrMdLbt2+X3W7XkiVLdOHCBSUnJxNpAAigFiM9evRoPfLII56vrVZruw8EAPhai5G+7bbbJElVVVVKS0vTs88+63OFVqtFdnu0X4YDgGAyoWUtRlqSysrK9Ktf/UpTpkxRUlKSzxU2NrpVUXHZL8MBMEtsbEywRwioQLasuX3bYqTPnTunadOmaf78+frhD3/YLoMBAJrX4lvw1qxZo8rKSq1atUoOh0MOh0O1tbWBmg0AQl6LZ9IZGRnKyMgI1CwAgGvwYRYAMBiRBgCDEWkAMBiRBgCDEWkgSA4fPqgXXpinw4cPBnsUGMznh1kAtI+Cgj/p5MkTqq2t0eDBQ4I9DgzFmTQQJDU1tU1uAW+INAAYjEgDgMGINAAYjEgDgMGINAAYjEgDgMGINAAYjA+zAP/VsXOYosIjA7Y9q9XiuQ3k//Gkpr5WVRcbArY93BwiDfxXVHikElYmBGx74RXh6qAO+qLii4Bud9//7lOViPStgssdN4jftwAgkDiTvkH8vgUAgcSZ9A3i9y0ACCQiDQAGI9IAYDAiDQAGI9IAYDAiDQSL7ZpbwAsiDQSJ8y6nGrs3ynmXM9ijwGC3/M/wjp0iFRURFrDtBe2jvHUNqqrkbX/fJK5vu+T6tivYY8Bwt3ykoyLCdP/MjQHbXsy5S7JK+vzcpYBu99CSqarSrRPpw4cPaseObUpKSuZDP8BNuOUjDTPxyUzAP1p1TfqDDz6Qw+Fo71nwDcInMwH/8HkmvW7dOm3fvl1RUVGBmAcAcBWfZ9K9e/fWypUrAzELADQRcc1tKPIZ6UceeUQ2G5euAQTe8EaX+rhcGt4Yuu+C8Xt9rVaL7PZof68WUpv3q1VOdQgL7LlIsN6q6GqoUyP/Hu7TrfJntJ/brX6N7qBt34T95PdXc2OjWxUVl/292mYFMgDB1tb9Ghsbo8//7x4/T9MyZ3lXSTY5y08FdNu953+k82cvtemxvJZ8C6V9JLV9P7VFc/uWTxzeIHcHW5NbAGhPrYp0r169tGXLlvae5ZZQ23OQGjp+W7U9BwV7FAAhgNPBG+Ts3EvOzr2CPYbxIq3uJrcA2obLHWgXyd+tVv/O9Ur+bnWwRwFuaZxJo10M7Favgd3qgz0GcMvjTBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgRBoADEakAcBgNl8LuFwuZWZm6tixYwoPD9eCBQvUp0+fQMwGACHP55n022+/rfr6em3evFnPPfecXnrppUDMBQBQKyJ96NAhPfDAA5Kk++67Tx9//HG7DwUAuMLidrvdLS0wb948jRo1SsOHD5ckjRgxQm+//bZsNp9XSgAAN8nnmXTHjh1VXV3t+drlchFoAAgQn5EePHiw9u7dK0k6cuSI+vXr1+5DAQCu8Hm546t3dxQWFsrtduvFF19U3759AzUfAIQ0n5EGAAQPH2YBAIMRaQAwWEhEev/+/RoyZIjKyso897388svaunVrs4+pqKjQjh07Wr2NTZs2tXrZhIQESZLD4VBRUVGrH9ee2jLLsWPHdODAgTZva+vWrXrnnXckfb3/6urqVFBQIElauXKlXnvttRte/83Kzs7Wj370I9XV1TW5v7KyUhMnTtS0adOafeyvf/3r9h6vzdp6vFojNTVV7733XpP7FixY4DmWppgzZ47njRC3ipCItCSFhYVp7ty5au0l+GPHjundd99t9fpXr17d1tFuWW+++aaOHz/e5sePHTtWI0eOlPT1/jt79mzQ/2Dv2LFDY8aM0c6dO5vcX1hYqB49eignJ6fZx/7+979v7/Ha7GaPV0smTJigv/zlL56v6+vrtXv3bj366KPtsr1QEjJveB42bJhcLpfy8vL0xBNPNPleTk6Odu7cKZvNpiFDhmjmzJlas2aNPv30U23evFkTJ070LLtr1y7l5eV5vn7llVe0efNmXbx4UZmZmcrMzPR8r7CwUC+99JJcLpcqKyuVkZGhwYMHt/tzvRkrV65U9+7dNXnyZBUVFSkzM1O5ublavny53n//fblcLj366KNKTEzUtm3bFBYWpgEDBujee+/1rGPp0qX6+OOPVV1drb59+2rRokUtbquiosKz/5xOp44fP35d7JYuXaoDBw7I7XbrySefVGJiYrs8//3796t3796aNGmSZs6cqbFjx0q6Ep3f/va3+vLLL7VixQqNHj3a67FNSEjQvn375HA41KVLF1VWVurVV1+V1WqV1Pxr4uGHH9agQYN06tQpDRs2TJcuXdKHH36ouLg4LVmyRCUlJZo3b56cTqcsFosyMjLUv39/z/Yk6Te/+Y0mTZqk0tJS7dmzR7W1tfr888/19NNPKyEhodnj5Q+jR4/W7373O9XU1CgqKkrvvPOOEhISFB0d7fXYXb1/unbtqscee0wjRoxQUVGRsrKylJ2d7Vn3nDlzZLPZdPr0adXX12vMmDHavXu3ysrKtGrVKpWVlSk/P1/Lly+XJM8+KS4uVkZGhhoaGhQZGen5/ubNm/WHP/xBVVVVyszM9Pu+8LeQOZOWpMzMTK1fv17FxcWe+44dO6a//e1vys/PV35+vk6dOqXdu3frmWee0bBhw5oEWpKKi4uVnZ2t3NxcxcXF6Z///KemT5+uzp07Nwm0JB0/flyzZ8/W+vXr9dRTT7V4ecV0b7zxhl5++WXl5eUpMjJS3/rWt5ScnKwnn3yyyYu8qqpKnTp10h//+Efl5+fryJEjOnPmTIvrvnr/PfPMM7rzzjubXDbYs2ePSkpKlJ+fr40bN2rNmjWqrKxsl+dZUFCglJQUxcfHKzw8XB988IEkKTw8XOnp6Ro2bJjS0tJadWyTkpK0fv16T6Cl5l8TpaWlevbZZ7Vp0yZt3LhRU6ZMUUFBgQ4dOqTKykotXrxYDodDeXl5mjdvntLT01t8HlVVVVq7dq1Wr16t7OzsZo+Xv0RERGjkyJF66623JElbt27VxIkTWzx2X+2fCRMmaNu2bZKk119/XePHj79u/XfccYdycnIUHx+vkpISrVu3TqNGjWrxb7tZWVn6xS9+4TnR+ve//y1JGjBggDZu3KgnnnjilvgzGTJn0pLUpUsXpaena86cOZ4z2hMnTmjgwIEKCwuTJA0ZMkSfffaZBg4c6HUd3bp10+zZs3XbbbfpxIkTuu+++5rdXo8ePbRq1SpFRkaqurpaHTt29PtzuhnV1dUKDw/3PHeLxdLsssuWLdOyZct07tw5z+9y8SYiIkLl5eWaMWOGoqOjdfnyZTU0NNzUnIWFhfrkk0/kcDgkSU6nU6dPn1anTp1uar3Xunjxovbu3avy8nLl5uaqqqpKmzZt8vpaaM2xjYuLa/Xj7Ha7evbsKUmKjo7WnXfeKUmKiYlRXV2dioqK9IMf/ECSdNddd+k///nPdeu++lJe//79JUm333676uvrb3RXtElKSooWL16soUOHqrKyUgMGDNC6deu8Hjvp6/0zdOhQLVy4UOfPn9e+ffs0Y8aM69b9/e9/X5LUqVMnxcfHe/7b23P7aj+cPHlSgwYNkiSNGTNGkvTXv/5VAwYMkCR1795dtbW1fnv+7SWkzqQl6cEHH1RcXJznJ3d8fLw+/PBDOZ1Oud1uHThwQHFxcerQoYNcLleTx166dEkrVqzQ8uXLtWDBAkVERHheEN6udS9cuFBpaWnKyspSv379Wn09PFDmzJmjQ4cOyeVy6fz58+ratasiIiJ09uxZSdInn3wi6cpf9Xft2qVly5Zpw4YN2rZtm0pLS2WxWK7bR3v37lVZWZmWLVumGTNmqLa2tlXP+6tlvO33+Ph4DR06VLm5udqwYYMSExPVq1cvf+yCJrZv365x48YpJydHr776qrZs2aJ9+/apvLz8umVbc2y9/dBr7nEt/YCUpL59++rgwYOSpKNHj6p79+6SrkSvurpa9fX1Ta43e1uft+PlT9/73vdUXV2tjRs3aty4cZJaPnZfzWixWJSUlKSFCxcqISHBc9Jw7ezNufo1W1paqosXL0q6ss8++ugjSVeObW5urs91mSjkIi1d+aVRkZGRkq68sBITEzV58mSNHz9ed9xxhx566CH17t1bhYWFWr9+vedxHTt21ODBg5WcnKzU1FRFRkbqyy+/lHTlBfH888832c5jjz2mX/7yl5oyZYqKi4s9y5riqaee0uLFi5Wamqrk5GTZ7XYlJiZqz549cjgcOnr0qKQrf9Xv3LmzHn/8cf30pz9VQkKCevbsqbvvvlt5eXl6//33Peu899579cUXX2jChAlKS0vTd77znVY976/2X7du3dTQ0KAlS5Z4vvfggw8qOjpaU6ZM8Vwjbo+/lRQUFOjxxx/3fB0VFaVRo0Zpy5Yt1y3b1mPb1sfNmjVLmzZtUmpqqjIzM7Vw4UJJ0tSpUzVx4kSlpaV5zsSb4+14+du4ceNUUFDg+QfD1h67sWPH6s033/R6qcOXu+++WzExMUpJSdHKlSs9PwRmzZqltWvXyuFwaMeOHUpKSrqJZxY8fOIQQNCdOXNGs2bN0oYNG4I9inFC8kwagDn+/ve/6+c//7mee+65YI9iJM6kAcBgnEkDgMGINAAYjEgDgMGINAAYjEgDgMGINAAY7P8BAtR4zCL2SJsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD3CAYAAAAALt/WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6B0lEQVR4nO3deXwU9fnA8c/s7J1NsgkJNwmHRECEEPBAiICK1nofbRDFo1qr1WorP+vRipRSxN6tt+2vHqglP9GqaL0QFUUFOYKACghy3+TcTbLHzPz+2GQBIeTaa5Ln/XrtK7szOzPPN5BnZ5/5zverGIZhIIQQwrQsyQ5ACCFE+0giF0IIk5NELoQQJieJXAghTE4SuRBCmJw10QfUdR1Na1tHGVVV2rytWUmbOwdpc+fQnjbbbGqT6xKeyDXNoLKytk3ber3uNm9rVtLmzkHa3Dm0p825uelNrpPSihBCmJwkciGEMDlJ5EIIYXIJr5ELIURraFqYiop9hMPBZIfSbnv2KDQ3KorVaicrKxdVbXl6lkQuhEhpFRX7cDrdpKV1R1GUZIfTLqpqQdP0JtcbhoHfX01FxT5ycnq0eL9SWhFCpLRwOEhaWobpk3hLKIpCWlpGq799SCIXQqS8zpDEG7WlrR0qkRuGwRtr97Ctoi7ZoQghRMJ0qBr5qh3VTH9rHXZV4dpT8rjmpD7YrR3qs0oIkQQrVizj1Vdf4je/eSC67LHHHiI/vy+7d+/i008/RlWt3HbbHQwZMjTh8XWoLLdww37sqsLpA7rw5CdbuOLZ5Ww+0LnuHBNCJE51dRVlZSt48slnmD59Fn/+8++TEkeHOSM3DIOFG/ZzSn4WD1wwhAs3l3Pv61/xj0+38LvzByc7PCFEDLyxdg+vrdkd031eOLQ7553QrU3bqqqVk046FUVR6N69e0NXyQqysrJiGmNzOkwi/3J3DXtqAtw8pi8Ao/tmc96Qbvzni11U1YXIdNmSG6AQwtSWL1/GrbfeGH29c+cOLrzwEjIzvdFlbncafr8v9RK5rutMnz6ddevWYbfbmTlzJvn5+QDs27ePO+64I/rer776iqlTp3LFFVfEL+ImLNywH9WiUDwgO7rswqHdKV25k7e+2ktJUa+ExySEiK3zTujW5rPn9ho5ctQRNfK0NA+1tf7ostpaPx5P04NbxUuzNfIFCxYQDAYpLS1l6tSpzJ49O7ouNzeXOXPmMGfOHO644w6GDBnCD3/4w7gGfDSNZZWT87xkOA+eeRd09TC4myfmX8WEEAIgHA6zdOln6LrO7t270XUDr9eb8DiaPSNfvnw5xcXFABQWFrJmzZoj3mMYBr/97W/54x//iKo2PWYuRMbj9XrdbQpWVS1H3fbLXdVsr6zn5nEDjlhfclIe01//kh21IU7omdmm4yZTU21u8/7OOgMAbcHCmO0z1mLdZjOQNjdtzx4FVU1uvwxVtaAoh8dhsShkZXkpLBzBTTf9CMPQufPOu5uNtSVtUZTW5clmE7nP58Pj8RwShEo4HMZqPbjpwoULGThwIP3792/2gPEYj/zVFdtRFTi5Z8YR60/P9+KwWnjuk83cddbANh03mWI9ZnNmOHJ7cFUKjwMt41R3Di1ts2EYx7ytPRGGDy9i+PCiw+L4yU9ujT6/7rqDtfNjxdrcLfqNDOPIPNmu8cg9Hg9+/8EakK7rhyVxgNdeey0pJZVGC9fvY0QfL173kRc0051WJgzM4e2v91Ef0pIQnRBCxFezibyoqIhFixYBUFZWRkFBwRHvWbt2LUVFRbGPrgW2V9axubyOCcflNPmeC4d2oyYQ5oNvDiQwMiGESIxmSysTJ05k8eLFTJo0CcMwmDVrFvPnz6e2tpaSkhLKy8tJS0tL2lgIO6rqATgut+l60sg+Xrp67HzwzX6+N7hrokITQoiEaDaRWywWZsyYcdiyAQMGRJ9nZ2fz6quvxj6yFjrgj4wS1sVtb/I9FkXh1L5ZvL/hAJpuoFo6zwA8QoiOz/S36EcTeVrTiRzglPwsagJhvtpTk4iwhBAiYTpAIg/hsFpIsx+72+PJeVkowGebKxITmBBCJIjpE3l5bZAuafZma/Ret41B3Tws2SKJXAjROitWLOP+++85bNljjz3Ef/87H4Dt27cxZUryeu6ZfqyVA/7gMevjhzq1bxbPLt2GLxDG4zB904UQKeCtt97gxRfnUlVVlbQYTJ/NDtQG6eN1tei9p+Rn8dSSbSzfVsm4Y3RXFEKkJsfX83B+NTem+6wfPInAoMvbvH16egYPP/wkJSUXxy6oVjJ/IveHKOzVslvvh/XMwGWz8NnmCknkQohWOdrohzfccBNjxhQnMaoIUyfysKZTWRdqcWnFploY2ccrdXIhTCow6PJ2nT23x9FGP0wVpr7YWV4bAqBLWsvHGj81P4ttlfVsr5R5PYUQHYOpE/mB2pb1IT/UKX0jA74vlbNyIUQHYerSSktvBjpUfpaLbukOlm6t5NLhPeMVmhCiAykqGkVR0ajDlt18888Oe/3aa28nMqTDmPuMvA2JXFEUinpnsnJ7FYZhxCs0IYRIGFMn8sYaeXYLL3Y2KuydSXltiG2V9fEISwghEsrUifyAP0i6w4rD2rpmjGjorrhye2UcohJCiMQyfSLPPspkEs3pm+3C67Kxckd1HKISQojEMn0ib019vJGiKBT2yqBse/JuqRVCiFgxdyKvDbUpkQOM6J3Jjqp69tYEYhyVEEIklum7H7Y1kTfe1l+2o4qzB8msQUKIpu3atZNrrrmCgoLjo8tGjjyJ6677cRKjOsi0ibwupOEPanRpQ40coKCrB7dNZeV2SeRCiOb17duPhx9+MtlhHJVpE3lb+pAfympRGNYzg5U7pE4uhJlcfHHLRjttqVdeMf9wHc0mcl3XmT59OuvWrcNutzNz5kzy8/Oj67/44gtmz56NYRjk5ubyhz/8AYfDEdegof2JHKCwdwaPL95CVV2ITFfbzuyFEJ3D5s3fHjb64f33zyQ3NzW+zTebyBcsWEAwGKS0tJSysjJmz57NY489BoBhGNx33338/e9/Jz8/nxdffJEdO3bQv3//uAd+cMCsdiTyhjr5qp3VnD6gS0ziEkLEV7LOoFO5tNJsr5Xly5dTXBwZb7ewsJA1a9ZE13377bd4vV6eeeYZrrrqKiorKxOSxOGQM/I21sgBTuiejk1VWCndEIUQJtbsGbnP58Pj8URfq6pKOBzGarVSUVHBypUrue+++8jPz+emm25i6NChjB49usn9qaqC1+tuU7Cqaolu69cMFAX69sjEqra9F+WwXpms2VPT5pji7dA2x2R/DXfBpmp7IfZtNgNpc9P27FFQ2/E3HguqakFRiEkcLdmHorQuTzabyD0eD36/P/pa13Ws1shmXq+X/Px8jjvuOACKi4tZs2bNMRO5phlUVta2OMBDeb3u6LY7yv1kuWz4ato3XsqQrh7+vWIHe/b7Wn2rfyIc2uZYyAzrAFTFcJ+xFus2m4G0uWmGYaBpegIialrXrt154omn2x2HqlpatA/DODJP5uamN/n+ZjNXUVERixYtAqCsrIyCgoLouj59+uD3+9myZQsAy5YtY+DAgc0GGQsH/G2/GehQw3tlENYNvtpdE4OohBAi8Zo9I584cSKLFy9m0qRJGIbBrFmzmD9/PrW1tZSUlPC73/2OqVOnYhgGI0aMYPz48QkIu+FmoFaOeng0w3pmAJELnoW9Wzb3pxBCpJJmE7nFYmHGjBmHLRswYED0+ejRo5k3b17sI2vGAX+Qvtnt70+a5baTl+Vi1Y4qoE/7AxNCiARLvaJwCxiGwYHatt+e/13Dembwxc5qmWhCCGFKpkzkNYEwIc1o9YQSTRneM4Oq+jBbKsx/h5cQovMxZSI/4G//zUCHGt5wY9AXMj65EMKETJnIy2sjNwNlteNmoEPlZ7vIcFr5YqckciGE+Zhy0KxAQ19oZ4z6fVuUyABaq3bKHZ5CiCPFchjbDRvW8fHHi7juuh/z4Yfvc8IJQ8nJyW1XfKZM5GE9clGyPXd0ftewnhl8vKmcyroQXhlAS4iUlXnx92O6v6pX/tui98VqrJWBA49n4MDIB8KLL/6bvn3v7eSJ3KLEbJ+N/cm/kAG0hBDtsG3bVh58cCahUAin08n06bN49NG/UVVVRXV1FVdcMYWFC9/hnHPO45tv1jNz5jQeffR/sdnafgJpzkTecItrLBP5Cd3TUS2KJHIhUlxLz6BjraXD2D7yyF+56qprOfXU03jvvXfZsGEdACNHjmLy5Cl8/vlSAE47bSzHHVfAnXfe264kDmZN5HE4I3faVAZ19fCFTDQhhDiKlpZWtm7dwtChwwA488yJALz77lvk5eUfa7N2MWWvlcZEbovxiGjDembw5R4foSQP0COEMK/8/H589dVaAN55503mzZsLgKIcma8sFgu63v58Y+pEHsszcogMoBUI63y9xxfT/QohOo9bbrmd5557mltvvZF33nmTs88+t8n3Dh06jJkz76e6un2VAHOWVrTGXisxTuSHDKB1YsNzIYTo0aMnTz75dIve27t3H/72t8cOW/arX02PPi8qGkVR0SgAbrzxp9x440/bHZ85E7ke+4udADkeBz0znazaUcVVo3rHdN9CiI4jFArxi1/ccsTyvLx8fvnLXyU8HlMmci1aWol9ZWh4zwyWbKnAMAwUJbYfFEKItkm1v0ebzRa3+TvbMnif1Mi/Y3ivDMprQ+yoat/MQ0KI2LBa7fj9nWN0UsMw8PursVpbN46UKc/I41UjBxjeMzKA1qod1fT2tn+8cyFE+2Rl5VJRsQ+frzLZobSboijNfiBZrXayslp3p6c5E7muY1EiY6TEWv8cNx6HyqqdVZx3QreY718I0TqqaiUnp0eyw4iJeM3NatrSSjzKKhD5cDixRwarZEhbIYRJmDKRhzQjLhc6Gw3vlcGmA7VU14fidgwhhIiVZksruq4zffp01q1bh91uZ+bMmeTnH7zV9KmnnmLevHlkZ2cD8Jvf/Ib+/fvHL2IazsjjUB9v1DiA1uqdNYzpnx234wghRCw0m8gXLFhAMBiktLSUsrIyZs+ezWOPHezsvnbtWh588EGGDh0a10APFdb1uJVWAIb2yEBVYNXOKknkQoiU12wiX758OcXFxQAUFhayZs2aw9avXbuWJ598kn379jF+/Hh+8pOfxCfSQ4S1+NXIAVw2lYKuHqmTCyFModlE7vP58Hg80deqqhIOh7FaI5ued955TJ48GY/Hw6233sr777/PhAkTmtyfqip4ve42BauqFrxeNxarit2qtnk/LXFy/y7M/Xwbbo8Te4xmImqLxjbHbH8NbYnn7669Yt1mM5A2dw7xanOzidzj8eD3+6OvdV2PJnHDMLjmmmtIT08HYNy4cXz55ZfHTOSaZrS5+01j153a+hAWhbh042k0JMdNIKzz2fq90Zp5MsS6u1JmwzR5VXH83bVXvLpopTJpc+fQnjbn5qY3ua7ZU82ioiIWLVoEQFlZGQUFBdF1Pp+P888/H7/fj2EYLFmyJCG18nh2P2w0vFfkxqCy7TI+uRAitTV7Rj5x4kQWL17MpEmTMAyDWbNmMX/+fGpraykpKeEXv/gFV199NXa7ndGjRzNu3Li4Bx3W4nuxE6BLmp28LBcrd1RxNX3ieiwhhGiPZhO5xWJhxowZhy0bMGBA9PnFF1/MxRdfHPPAjiXS/TD+desRvTJ5/5v96IYRl7tIhRAiFkx5Q1AiSisAhb0zqK4Ps2l/56rjCSHMRRL5MRQ21MlXyjyeQogUZs5EHud+5I16ZTrJ9djlgqcQIqWZM5Hrelxv0W+kKAqFvTIp21HVKcZCFkKYk0kTeXwHzTpUYa9M9vqC7KyWiSaEEKnJxIk8Mb1IRvSO3AxUtl1u1xdCpCZzJvIE9CNvNCAnjXSHVS54CiFSljkTeZyHsT2URVEY3itDLngKIVKWeRN5gs7IAYp6Z7Kloo79vkDCjimEEC1lzkQe5xmCvmtkHy8AK+SsXAiRgsyZyBN8Rl7Q1UOaXWX5NknkQojUY9JEnph+5I2sFoURvTNZtq0yYccUQoiWMmkiT2xpBSLlla0VdeyTOrkQIsWYM5Frieu10mhkn8i4K1JeEUKkGtMlck03MCChNXKAglwPHofKcimvCCFSjOkSeViPjHmS6ESuWhRG9MqUnitCiJRjwkQemXMy0YkcYFRepE6+t0bq5EKI1GG+RK41nJEnYIag7xrZ2wvA8u2VCT+2EEI0xXyJPEmlFYCBXdPIcFpZvlXKK0KI1CGJvBUsSqROLv3JhRCppNlErus606ZNo6SkhClTprBly5ajvu++++7jj3/8Y8wD/K5k1sghUiffUVXPzioZn1wIkRqaTeQLFiwgGAxSWlrK1KlTmT179hHvmTt3LuvXr49LgN8VitbIk5PIT873ArBkS0VSji+EEN9lbe4Ny5cvp7i4GIDCwkLWrFlz2PqVK1eyatUqSkpK2LRpU7MHVFUFr9fdpmBV1YI7zQFApsfZ5v20x4hMF93SHazcWcN1pw+I+/FU1RLTdqrWyGd3Mn53LRXrNpuBtLlziFebm03kPp8Pj8dzSCAq4XAYq9XK3r17efjhh3n44Yd58803W3RATTOorKxtU7Ber5uKyjoAAvWhNu+nvU7K8/LRxv0cKPejxrnE4/W6Y9rOzHCkNFWVpN9dS8S6zWYgbe4c2tPm3Nz0Jtc1m8g9Hg9+vz/6Wtd1rNbIZm+99RYVFRXceOON7Nu3j/r6evr378+ll17apkBbItk1coBT8rN4fe0e1u31MaR7079cIYRIhGYTeVFREe+//z7f//73KSsro6CgILru6quv5uqrrwbg5ZdfZtOmTXFN4nCw14otSTVyOLxOLolcCJFszV7snDhxIna7nUmTJvHAAw9wzz33MH/+fEpLSxMR3xEOdj9MXs/JbLedgblpLJULnkKIFNDsGbnFYmHGjBmHLRsw4MiLfPE+E28UvbMziaUViJRXSlfuoD6k4bSpSY1FCNG5mfeGoCSWVgBOyfcS0gwZREsIkXQmTOTJv9gJUNgrE7uqSH9yIUTSmTCRJ79GDuC0qQzvlcnSLZVJjUMIIcyXyFOkRg6ROvk3+/0y/ZsQIqnMl8hTpEYOMKZfNgCffFue5EiEEJ2ZCRN5atTIAQbkuOnqsfPxJknkQojkMWEiT53SiqIojO3fhaVbKglperLDEUJ0UuZL5FpqXOxsNKZ/NrUhjZXSDVEIkSSpkQ1bIZRCNXKIDKBlVxUWS51cCJEkpkvkYS11auQALptKUR+v1MmFEEljvkSeQjXyRmP7ZbO1oo5tFXXJDkUI0QmZNpHHexzw1hjTP9INUcorQohkMGUit1oUFCV1Enlvr4u+2S4WS3lFCJEE5kvkmpFSZZVGp/XLZvn2SmqDWrJDEUJ0MuZL5LqeMj1WDnX6gC6ENIPPNstZuRAisUyYyI2U6UN+qOG9MvG6bCzcsD/ZoQghOpnUy4jNaKyRpxqrRWHcgC58vKmcYFju8hRCJI4k8hiaMDAHf1Dj822VyQ5FCNGJmC+Ra6lZI4fIXZ5pdpUPpLwihEigZhO5rutMmzaNkpISpkyZwpYtWw5b//bbb3PZZZdx+eWX8+KLL8Yt0EZaCp+R260WxvTL5sNvDqA19HcXQoh4azaRL1iwgGAwSGlpKVOnTmX27NnRdZqm8ac//Ymnn36a0tJS/vnPf1JeHt9eG6l6sbPRhIE5VNSFWLVTBtESQiRGsxlx+fLlFBcXA1BYWMiaNWui61RV5b///S/p6elUVlYCkJaWFp9IG6RyjRwi/cntqsL7Gw4kOxQhRCdhbe4NPp8Pj8cTfa2qKuFwGKs1sqnVauWdd95hxowZjBs3Lrq8Kaqq4PW62xSsqlrAouCwq23eR7x5gbHH5bBo4wFmXDy03Xegqqolpm1VrZHP7lT9/UHs22wG0ubOIV5tbjaRezwe/H5/9LWu60ck67PPPpuzzjqLu+++m1deeYXLLrusyf1pmkFlZW2bgvV63dQFwihG2/eRCGP7ZrFw3T4+WbeXE7qnt2tfXq87pm3NbOgaWZXCv79Yt9kMpM2dQ3vanJvbdC5ptrRSVFTEokWLACgrK6OgoCC6zufzcdVVVxEMBrFYLLhcLixxrl+n6i36hxp3XBesFoV3vt6b7FCEEJ1As2fkEydOZPHixUyaNAnDMJg1axbz58+ntraWkpISLrjgAq688kqsVivHH388F154YVwDDusGLpsa12O0V4bTxph+2by7bh+3nd4/pUZqFEJ0PM0mcovFwowZMw5bNmDAgOjzkpISSkpKYh9ZE8K6kbL9yA91zuCufLjxACu3VzEqz5vscIQQHVjq9uNrQljXU760AlDcPxu3TeUtKa8IIeLMfIlcS+1+5I2cNpXxA7uwcP1+GXtFCBFXqZ8Rv8MspRWAcwZ1pSYQ5tPNFckORQjRgZkzkZugtAJwcp4Xr8vG21JeEULEkSTyOLKqFs4qiNwc5A+Gkx2OEKKDMl8i18xxsbPR9wZ3JRDWWbheRkQUQsSH+RK5bmBVzRP2sJ4Z5GW5eG3N7mSHIoTooMyTERuYqbQCoCgKFw7tTtmOajaXd67bkYUQiSGJPAHOO6EbqgLz5axcCBEHpkrkhmGk9MQSTclJszO2fxdeX7uHsCZ9yoUQsWWqRB5umHXHLP3ID3Xhid0prw2x+Nv4TrwhhOh8TJXIQw1ns2a4s/O7TuuXTU6anVdXS3lFCBFbpsqIYa3hjNxkpRWIxHzeCd1Y/G05+3yBZIcjhOhAmh39MJWEGkorNhOWVgAuGtqdZ5Zu45XVu/nx6Pym36iHUau3olZsxBLYStr+7VjqK1DqK1C0gx8ChmrHcGSiO7wY7ly0jDy0jD5o3v4YzqwEtEgIkQpMlcjD0dKKORN5nywXp/XL4qVVu7j25D7YGvvD62FsO5dg2/Eptl1Lse1ZgRKuj27nsrrRnVnoziywuaLLlUAVlopvIj8Dh0/2rHl6Es45gXDXYYR6nEyoW9Fh2wohOg5zJfLGi50mrJE3KhnRi9tfXsPCdXs437MOxzfzcXz7Dpb6CgzFQrjLEOoGX0E490S0rAF48odSGXA0v+NQHWrNtsiZfPkGrPvXYt2/FvvmBSgYGBYb4W6FWGq2Rz4QDAPaOZ+oECI1mCuRN56Rm7S0AjC6m8Fdnrc468Nf4NX3oNs8BPueRWDA9wn1HovhyDh8A5cbAi24kcjmQssuQMsugL5nRRcrgWpsuz7HtvMzbNsXo1ZvQ63eRvYzJxHofy7B484n1OMkUMz74ShEZ2eqRB4y8cVOpXY/7hWP4lrzDDdrAZbog9h/6i/pVnQxqC04424jw5FBsO+ZBPueCYB3zlkogUq0rsNxffkC7tVPoaV1IzDwYuoH/QCty6C4xSKEiA9TJfKDpRXzJHIlWIN7+SO4vvgXaPUECi7lwNAf86MXKxm7N5uZcUziR2Oodgx3V6q//78oQR/2ze/i2PAari/+F3fZE4RyT6R+yGQCBZdg2D0JjU0I0TbmSuQNpRXVDDVyw8DxzeukLZ6Oxb+XwMALqT3pDrSsAdiBC4du5P/KdnL7uAC5nsQm82iIdg+BgksIFFyCUrsf54ZXcH71f6R/eA9pn8wkMPBi6k68Bi1nSFLiE0K0TLOJXNd1pk+fzrp167Db7cycOZP8/INd515//XWeeeYZVFWloKCA6dOnY4lTog2Z5M5OS/U20j+4C/u2RYRyT6T63H8S7jbisPf8oLAnc1fs4P9W7uSW4n5JivQgw51D3fAbqBt2PdY9K3GufR7n+pdwffk8wV6jqRt2PcG+E8GiJjtUIcR3NJtxFyxYQDAYpLS0lKlTpzJ79uzouvr6ev7617/y7LPPMnfuXHw+H++//37cgjVD90PHupfJKj0b6+4V1Jw+k8rLXz8iiUOkK+IZBTm8WLaTmvoUmnRCUQh3L8J35p84cM3n+Eb/CrVqK5lv3kD286fjXP0MhOqSHaUQ4hDNnpEvX76c4uJiAAoLC1mzZk10nd1uZ+7cubhckf7J4XAYh+PYZQJVVfB63W0KVq+M9K32ZrjavI+4CVSjvjkVy9qX0HufgnbREzi9eTiPscnPzizgvcc+Yf7X+/jp+AFHfY+qWmLaVtUa+exu2T7d0H0q+vjbMda9geWzh0hf9Cs8y/6MPurH6KNuBJc3ZrFFY4xxm81A2tw5xKvNzSZyn8+Hx3PwopeqqoTDYaxWKxaLhZycHADmzJlDbW0tY8aMOeb+NM2gsrJt43IHQxoA9bWBNu8jHtSKjWS8eT1K5bf4T7mT2qJbACs0E2Mvt5Ux/bJ56pPNXHJCV1y2I8sWXq87pm3NDEe+1VS1dp89JsLFZ2HbtRTXysdwLJqN8ulD1A+9mtrCGzHcuTGLMdZtNgNpc+fQnjbn5qY3ua7Z0orH48Hv90df67qO1Wo97PWDDz7I4sWLeeihh1DieJPJwRp56lzstG9+D++887HUlVN10VxqR90OlpZfQ77ulD5U1oX4zxe74hhljCgKoZ6nUH3e05SXvEOw71m4yp6gy7OnkvbxdCx+GRBMiGRoNiMWFRWxaNEiAMrKyigoKDhs/bRp0wgEAjz66KPREku8hMKpVSN3rfonGW9ci5aRT8UP/kuo1+hW72N4r0xG9snkuWXbCYbNM1a5ljOEmrMfoWLyB9QPvBjXF0+RPWcMaYvuk4QuRII1m8gnTpyI3W5n0qRJPPDAA9xzzz3Mnz+f0tJS1q5dy7x581i/fj3XXHMNU6ZM4d13341bsCnTj9wwSPvkd3g+nk5wwLlUXvof9Izebd7ddSfnsc8X5PW15kuAmrc/vjP/RPmVi6g//lJca+dEEvpH96P49yY7PCE6hWZrABaLhRkzZhy2bMCAgxfmvv7669hH1YSUGMZWC5H+wS9xfv0idUOvwVc8o91d8k7O9zKsZwb//Gwr3x/SDedRauWpTs/MxzfhD9QW3Yp72d9xrX4a15fPU3fitdSO+CmGKzvZIQrRYaVOsbkFQnqk9GBLVo1cC5Lxzs04v34R/8n/g+/0mTHpV60oCrcW92OfL0jpyp0xCDR59Mx8fGf+iYrJ7xPo/31cK58ge85o3Ev+gBKoTnZ4QnRIpkrkST0j14JkvPNTHJvewjf2N9Se9POYjh44oncmY/tn88zSbVTXh2K232TRvP2pmfh3Kq54j1DeeNKW/Y3sOaNxrXhE+qELEWPmSuR6kkY/PCSJ1xTPoG749XE5zC1j++ELhHlm6ba47D8ZtOwCqr/3BBU/fJNQ95F4Pn2A7OfG4Fz9NGjBZIcnRIdgqkSelNEPdY30BbdHk3j9sB/F7VDH5aZx7pCulK7cyZ6ajjUdXDj3RKrPf5aKS/+D5u1H+qJfk/38OBxfzwNdS3Z4QpiaqRJ5wieWMAw8i36N85v5+E77dVyTeKOfnNYX3TB4bPHmuB8rGcI9TqLq4nlUXvAcutNLxns/J6v0bOyb3opMdiGEaDVzJfIEj7XiXvonXGvnUDviZupG3JSQY/bMdHJFUW/eWLuHVTuqmt/AjBSFUN54Kn/wBlXnPA56iMw3b8A77wJs2z5OdnRCmI7JEnniRj90rn6GtGV/pW5wCf7R98b9eIe6/tQ8unrs/P69b6IfXh2SYiF43PlUXLGQmgl/wFK7F+9rk1CfvwTr7hXJjk4I0zBVIg/pOhYFLHGea9K++T08H91HoO9Z+MY/mPC5Ld12lTsmDGD9Pj///rzjXPhsksVK/ZArKL9yEb6x01H2riXrpQvJ+O/1qAe+SnZ0QqQ8UyXysGbEvayiHviK9HduIdxlMNUTH2nVuCmxdMbAHE7O8/KX9zZwwN9JendYndQNv4HwLSvwn3Inth2fkjX3bNLfuQW1clOyoxMiZZkrket6XC90Kv69ZL5+LYY9jerzngJ7WtyO1WwsisKdZxxHfUjjrx92siRm91A76nbKpyymruinOL59h6wXJuBZ+D9YqrcnOzohUo65ErlmxK8+rgXIfOvHWOrLqf7+U+ienvE5Tiv07eLmJ8X9eeurvXz4zf5kh5NwhjML/+h7OHDVYupOvBbnupfJfr4Yz4e/wuIzwWiRQiSIqRJ5SI9facXz0f3Ydi+n+sy/EO46LC7HaIubxw1gYG4as97dQGWd+e/4bAsjrSv+4t9QPuVj6gdPwvnlC2Q/N5a0j+7H4t+T7PCESDpzJXJNj0sid375Aq61z1FbdAvB486P+f7bw261MP17x1NdH+YP732T7HCSSvf0xDf+gchIiwWX4Fr9NNlzTiPt49/ISIuiUzNVIo/HxU7r7hV4Pvw1wT7j8J/yy5juO1YKunq4YXQe76zbx3vr9yU7nKTTM/rgO+OPlE/+gMDAi3B98S+6PBdJ6HKGLjojkyVyPaazAyl15WS8/RN0T3eqz344pWeIv+akPgzu5mHWuxvYWVWf7HBSgu7tR82Zf44k9OMuwPXFvyJn6Ivuw+Iz9yiSQrSGqRJ5TGvkhk7Ggtuw1B6g+ntPYDizYrPfOLGqFmadPxhNN7jn9a8IdeQbhVopktD/QvmVH0ZKLg2TW3je/yWWqs3JDk+IuDNVIg/HsEbuWvEo9q0f4Cv+DeHcE2Oyz3jr7XUx7XvH8+XuGv7W2boktoCe2TdScrnyY+qHTMa57iWynz+d9Hd/hnogcROgCJFo5krkuhGT0optx6ekLfk99QMvov6Eq2IQWeKcMTCHK4p6UbpyJwvWSb38aPSM3vjG/S7SD33YDTg2vU323LPIeONarLs+T3Z4QsScuRJ5DC52KnXlpL97K1pm36Tcfh8LPzu9HzXzxnD9lZms2+tLdjgpS0/rjn/sNA5cswT/yVOx7V5O1suX4H3poobRFqU8JTqGZhO5rutMmzaNkpISpkyZwpYtW454T11dHZMmTWLjxo1xCbJRSG9nacUwSF84FUtdBdVnP4Zh98QuuASyqRYG5KRhtSjc8Z817O1gY5fHmuHMovakX3Dg6iXUFP8WS+0+Mt+8gaznx+Fc86zMWCRMr9lEvmDBAoLBIKWlpUydOpXZs2cftn716tVceeWVbNsW/8Gd2ntG7lz9NI7N7+I/7V603BNiGFni2VSFATlp+AIad7yyltqgTM7QLJub+mHXUX7lIqrPfhTDkUH6h/fS5ZmTSPt0tvR0EabVbCJfvnw5xcXFABQWFrJmzZrD1geDQR555BH69+8fnwgPEdb1Nt+ir+7/Es8nMwnkn0HdsPhM1ZZobrvKrPMHs2Gfj1+/8VXHHvI2lixWAgMvpPLy16m45GVCvU7FteIRsuecRvrbP8W6e7lMciFMpdmh/Xw+Hx7PwRKEqqqEw2Gs1simI0eObNUBVVXB63W3MsyIsG7gcthav32oDuv/3QbOTCyXPIY3LXmDYbWWqlqO2l6rNfIZfF5RbyrDOtPnf8nM977hT5cPRz3Gtxa1Ybu2/hskQlNtjous8TBkPOHKLViW/RNH2Ryc37yG0W0Y2qgbME64FGzxjyWhbU4R0ubYaTaRezwe/H5/9LWu69Ek3haaZlBZWdumbUNhA0PTW7192kfTsO37msoLnicUSoM2Hj8ZvF73UdsbDrsAqKys47yCHMpP78ffF32Lohvcd05Bk2O2Z4YjZ+1VKfw7aKrN8ZULo+6BYbfhXP8fXKufwvrGbegL7qN+0A+oP+EqtKzj4nb05LQ5uaTNrZObm97kumZLK0VFRSxatAiAsrIyCgoK2hRELITbcLHTtvVD3F/8i9phPyKUNy5OkSXflJP6cOPofF5fu4ffv/cNupQG2saeRv3Qq6iYtIDKS+YRzBuPa/UzZL8wnsz/XIZj3UsQloujIrU0e2o9ceJEFi9ezKRJkzAMg1mzZjF//nxqa2spKSlJRIxRrR3GVqmvIP29OwhnFeAffU8cI0sNN4zOI6jpPL10G/6gxv3nFMR0SINORVEI9TyVUM9T8dXux/l1Ka61L5Cx4Hb0j6YRKLiY+kElkZvJTNiFVXQszSZyi8XCjBkzDls2YMCAI943Z86c2EXVhFZ1PzQM0j+4G0t9ORXnPwtWV3yDSwGKovDTsX1Js6s88vFmqutDPHjBEJy21B1DxgwMdw51RbdQN+JmbDs+xfnlCzi/nItr9TOEuwyiftAPqR94MUZa12SHKjopU52uhTSjxTMEOda/jGPjG/hP+R/TdzVsDUVRuPaUPO6dOJDPNlfw0xdXU17bSaaKizfFQqj3GGrOfoQD1y6nZtwDGKoTz+IZdHnmJDLnX4Vj/X8g1LnqviL5kjMhZRu1dKwVS81OPIvuI9TjJOoKb0pAZKnnkmE9yHRamfbmOq5+biV/vGgIg7o1fbFEtI7h9FI/dAr1Q6eglm/Asf5lnOteJuPdn2FYXQT6nU1g4EUE88aB6kh2uKKDM1ci11tQIzd00hdORdHDVJ/5l5QemjbezijIpWemk/959UtumLuK+84u4IfJDqoD0rIHUnvqXdSecie2XUtxrH8Vx8bXcW54Fd2eQbDf2QSOu4Bgn7GS1EVcmCqRR2YIOnZpxbnmWezbP6Jm3APomX0TE1gKG9QtnWevGsFdr33Jr//7NadW1NLb2/GvFySFYjl4gbR4BvZti3BsfAP7prdwrpuHbvMQ7Hsmgf7nEsybkNTJvUXHYqpE3twZuaXyWzyf/I5g3jjTjWoYT9luO4/+YBiPL97M/heC+AIavr0+ju9qzrFmTEG1Eex7JsG+Z8L42di3fYR905s4vn0H54ZXMVQHwd5jCPY9O/Ie75EdCIRoKdMkck03MAyarpEbOhkLp2JYbNRM+IN0CfsOm2rhZ6f3x5LrYXN5Ldc+v5LrTunDtSfnYbea6pq3+aj2aFL36WFsu5Zi//YdHN++g2PLQvgQjG7DcPeZQDD/DMJdh4PFNH+aIgWY5n9LWI/c4NJUInd98S9su5ZSfeZf0D09ExmaqaQ7rQzuls6ZBTn849OtvLtuH/dOLGBE78xkh9Y5WKyEep1GqNdp+Mfcj1q+DvuW93Bvfx/38odIW/Y3dEcmwd7FhPJOJ9j7dPSM3smOWqQ4EyXyyK3lR0vkauUm0j59gEDfswgcf3miQzMdq6ow87zBnDukGw8u2MCNpas4Z1Autxb3o3uGM9nhdR6KgtZlEHVdBuE4406qdu/Avu1jbNs+wL71Q5wbXwcgnNmXUO+xhHqNIdhrNIY7J7lxi5RjnkSuNZyRf/dORV0j/b07MKxOfONnS0mlFcb0y6b02lE8tWQrLyzfwfsb9nPFyN5cc1If0p2m+a/RYRjOLAIDLyAw8AIwDNSKDdi3fYRt+0c41r+Ca+1zAISzjyfU85TIo8fJ6J4eSY5cJJtp/lqbKq24Vj+Fbfcyqs/6K3pa92SEZmoum8pPx/bj0mE9ePTjzTyzdBvzynYyqagXVxT1ItNlS3aInZOioGUXUJddQN3w60EPY923Gtv2xdh3foZj3cu41jwLgJbem1D3UYR6nES4exHhLoOlxt7JmOZf+2iJ3FK9jbTPfk8gbwKBgsuSFVqH0D3DyYzvD+KqUb3515Kt/O9nW5m7YgcXndidHxT2lC6LyWaxEu42gnC3EdSNvDWS2Pd/iW3XUmy7Po8MHbDhFQAMq4tQ12GEuxYS7jqcUNfh6Bl58m21AzNRIv9OjdwwSP/wbgB84x6Q/6QxUtDVw+wLhvDNfj9PL9lK6cqd/Hv5Dk4f0IXLCntwcl7WMcc7FwlisRLuOoxw12HUDb8BDANLzXZse1Zg3b0C254VuFY/jaJFpgHUHZmEc09seAwlnHMCWma/Tn3DXEdinkQerZFHkohj/X+wb/2QmuIZclU/Do7LSWPmeYO57fQAL63ayctf7ObDjQfo6rFz/tDunDu4K32zO9ekAClNUdAz+hDI6ENg4EWRZVoQ64Gvse5dhXXfaqz71uBa9U8UPQSAYXUSzj6ecPYgtJzBhLML0LIL0N3d5MTIZMyTyKOlFQtKoBrPx/cT6lZE/dBrkhxZx9Y13cHNY/tx/an5fLTpAK+t2c3TS7byr8+2MjA3jbMKcpkwMIe+2S4U+eNPLao9etYepQVRK77Bun9t5HHgaxxbFmD5ujT6Ft2RiZY1kHDWcWiND29/tIw8qb2nKNP8qxxWI9eChHNOwFc8Q74aJojdauHMglzOLMhlny/Ae+v3s2DdPh5bvJnHFm+mV6aTsf2zGd03mxG9M3Hb5d8lJal2tJwhaDlDCPCD6GKldh/W8vWo5esiPys24Nj8Lpav5kbfY1hsaBl90DL7Rh96Rj5aZj5aem+wStfVZDFPItcO1sgNdxeqLprbzBYiXnI9DiYV9WJSUS/21ARYvOkAH20q55XVuylduRPVojC0ezoj+2QyrGcmJ/ZMJ8MpvV9SmeHOJeTOJdR7zGHLlfoK1MpNqBUbsVZuRK3ajFr5LfYdn6GEDx+uV3N3Q8/ojZbeGz29V+Snpxeapwe6pyeGI1NKNnFinkTecEZua8UMQSL+uqU7uHR4Ty4d3pP6kMaqndUs21rJ51sreWbpNjRjGwD9st0M6uZhcPd0BnX1cFxOmvRVNwHDmUW4+0jC3UcSOGyFgVK3H7VqC2r1ZtTq7VhqtqFWb8O2pwzLxv9Ga/HRTaxutLRu6J7u6GndsWT3xqVmo6d1RXc3PnIx7OmS8FvJNH9Jh9bIRWpy2lROyc/ilPwsAGqDGl/urmHVzirW7Krh862VvPnV3uj7u6U7GJDjpm+2m/xsN/lZLnp7XWRkSFfHlKcoGO5cwu5cwj1GHble17DU7sXi24nFtwu14afFvxvVvxvbrs+xbHwDj3bkpCeG6kB3dUF356K7umC4uqA7s9Fd2RjO7MhzZxaGMyvy05EJauf+xmeeRK4de6wVkXrcdpVReV5G5Xmjy/b5Aqzf6+eb/ZHHxv1+lm+rIhDWo++xWy30zHDQPcNJjwwHPTKc5Hrs5KY5yE2308VtJ8NplYurqcyiont6HPOuU2+mi6o9O7H492Cp3R9J/LV7sdTtx1J3AEvtvsjyA19jqTsQ7Up5NLotDcPhxXBkojsyMByZkef2dAx7OoYjA8OefvC13YNh80Se29wYdo+pL+Q2G7mu60yfPp1169Zht9uZOXMm+fn50fULFy7kkUcewWq1ctlll/HDH8Zn6oLoGbmUVkwt1+Mg1+NgTP/s6DLdMNhbE2BLeR3bq+rYVxdm014fu6vrWbfHR0Vd6Ij9WC0K2W4bWW47WS4bXreNTKeVTKeNTJcVj8NKesPD47CS5lBJs6u47VY5GUgVioLhzEJzZqF1aea9hgHhOix15Vjqy1HqK7DUV0R+BqpQApWRn/VVKIEq1KrNKMEalEA1lpCvReEYqgPDltbwcGNYXZGfNjeGNfIamyuyPPpwgtWJ0fBAdUSeqw4MtWGdagfVHtme+HTZbTaRL1iwgGAwSGlpKWVlZcyePZvHHnsMgFAoxAMPPMC8efNwuVxcccUVTJgwgdzc3JgHeqxBs4S5WRSF7hlOumc4OYUsvF43lZUHL6TVhzT2+4Ps9QXYVxPkQG2QA/4QB2qDVNaGqKwLsbWyjur6EL6A1uzx7KqCy6ZGH06bBadNxWm14Gh4OK0qdqsFu2rBblWwqZHnNlWJ/rRaGn8qqBYL1obnhz5Ui4JFifxUG35aLKAqkeUWi4KqgGG3UlMfiixTFCxK5PeiKJF5WC0KKNB5v4UoCtjc6DZ36+8b0TWUkA8lUIMSrEYJ+rCEfChBP0qopuGnDyXoQwnXoYT8DY86lHAtim9Xw/LayM9w/TG/HRxL+OJ/QK9z27TtsTSbyJcvX05xcTEAhYWFrFmzJrpu48aN5OXlkZkZGQJ15MiRLFu2jHPPjX2gmtTIOy2nTaW319WiYQLCukFNfYiagEZNIExNfYjaoIYvqOELhKkLadQGNfxBjfqQRn1Ypy6kURfSqaoPEwhr1Id0QppOIBx5hDSdhspe0imAJZLRo8+VhueKAgpK9Drhoa8b0/+h72183bjf6LLvHvMo7znifd/5gFGaWHXoclW1oOst/8W25COsZZ9zCpDe8Gh6fCblaEdUIw+LoWEniIMgDiOIgwB2gtiNEA6ChzwPYCOMzQhhVTROVk4kHlOINJvIfT4fHs/BmWRUVSUcDmO1WvH5fKSnH5zQNy0tDZ/v2F9jVFXB623914szhvbgpuoAhf27HDkCYgemqpaj/r4++KDxWSt/lw0betsRU7w11eaWiscgr2FNJ6jphDSDkKYTDOuEdINQWCesG4S1yOuwpqPpRvS5rke+TWq6gWYY0de6QXSZYRgYKIQ1HcMwousan+tGZFIVg4bXuoEB0fWNzxvfY0SXRV7zndfGIds0anx6xPJD1hscknQPffqdXHzoy6Ptq3EbRTnatk0k9hbk+1h+1n43rhZvB9Q3PL7LYoEzc3q26/92U5pN5B6PB7/fH32t6zpWq/Wo6/x+/2GJ/Wg0zTjsa3NrTD2roM3bmtV3ywydQaq32drwwALYG08q2ncDVKq3OR6kza2Tm9t0bm321LaoqIhFixYBUFZWRkFBQXTdgAED2LJlC5WVlQSDQZYtW8aIESPaFKQQQoi2afaMfOLEiSxevJhJkyZhGAazZs1i/vz51NbWUlJSwt13383111+PYRhcdtlldOvWLRFxCyGEaKAYRlurQW0TCmlt/mohX8U6B2lz5yBtbp12lVaEEEKkNknkQghhcpLIhRDC5CSRCyGEyUkiF0IIk0t4rxUhhBCxJWfkQghhcpLIhRDC5CSRCyGEyUkiF0IIk5NELoQQJieJXAghTE4SuRBCmJwpErmu60ybNo2SkhKmTJnCli1bkh1S3IVCIe68804mT57M5ZdfznvvvZfskBLiwIEDjBs3jo0bNyY7lIR54oknKCkp4dJLL+XFF19MdjhxFQqFmDp1KpMmTWLy5Mkd/t951apVTJkyBYAtW7ZwxRVXMHnyZO6//370hnmIY8EUifzQCaCnTp3K7Nmzkx1S3L322mt4vV5eeOEF/vGPf/Db3/422SHFXSgUYtq0aTidzmSHkjBLlixh5cqV/Pvf/2bOnDns3r072SHF1Ycffkg4HGbu3Lnccsst/PWvf012SHHzj3/8g1//+tcEApGJmh944AF+/vOf88ILL2AYRkxPzkyRyI81AXRH9b3vfY/bb789+lpV2zeVmBk8+OCDTJo0ia5duyY7lIT5+OOPKSgo4JZbbuGmm25i/PjxyQ4prvr164emaei6js/ni04b2RHl5eXx0EMPRV+vXbuWk08+GYDTTz+dTz75JGbHMsVv8VgTQHdUaWlpQKTtt912Gz//+c+TG1Ccvfzyy2RnZ1NcXMyTTz6Z7HASpqKigp07d/L444+zfft2br75Zt56663ozPUdjdvtZseOHZx77rlUVFTw+OOPJzukuDnnnHPYvn179LVhGNF/17S0NGpqamJ2LFOckR9rAuiObNeuXVx99dVcdNFFXHDBBckOJ65eeuklPvnkE6ZMmcJXX33FXXfdxb59+5IdVtx5vV7Gjh2L3W6nf//+OBwOysvLkx1W3Dz99NOMHTuWt99+m1dffZW77747Wnro6CyWg+nW7/eTkZERu33HbE9xdKwJoDuq/fv386Mf/Yg777yTyy+/PNnhxN3zzz/Pc889x5w5cxg8eDAPPvggubm5yQ4r7kaOHMlHH32EYRjs2bOHuro6vF5vssOKm4yMDNLTI1OWZWZmEg6H0TQtyVElxpAhQ1iyZAkAixYtYtSoUTHbtylOa482AXRH9/jjj1NdXc2jjz7Ko48+CkQunnSmC4GdwYQJE/j888+5/PLLMQyDadOmdejrIddeey333nsvkydPJhQK8Ytf/AK3253ssBLirrvu4r777uPPf/4z/fv355xzzonZvmUYWyGEMDlTlFaEEEI0TRK5EEKYnCRyIYQwOUnkQghhcpLIhRDC5CSRCyGEyUkiF0IIk/t/eRo8WuQA5REAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### standard deviation of the sample mean of each group\n",
    "SW_g = np.sqrt(SSW_g / (n - 1)) / np.sqrt(n)\n",
    "data = [0 for i in range(k)]\n",
    "for i in range(k):\n",
    "    data[i] = t_test.rvs(df=n[i]-1, loc=mean[i], scale=SW_g[i], size=100)\n",
    "\n",
    "plt.figure(1)\n",
    "plot = sns.barplot(data=data, ci='sd')\n",
    "plot.set_xticklabels([\"Not at all\",  \"Just a little\", \"A fair amount\",  \"Very much\"])\n",
    "plot.margins(0.1)\n",
    "plot.set_ylim(0, 6)\n",
    "\n",
    "### non-central parameter\n",
    "nc = SSB / MSW\n",
    "\n",
    "x = np.linspace(0, 10, num=100)\n",
    "y0 = f_test.pdf(x, dfn=df_SSB, dfd=df_SSW)\n",
    "y1 = ncf.pdf(x, dfn=df_SSB, dfd=df_SSW, nc=nc)\n",
    "plt.figure(2)\n",
    "plot = sns.lineplot(x=x, y=y0, label=\"H0\")\n",
    "plot = sns.lineplot(x=x, y=y1, label=\"H1\")\n",
    "plt.vlines(f, 0, f_test.pdf(f, dfn=df_SSB, dfd=df_SSW), 'b', label='F')\n",
    "plt.vlines(f_crit, 0, 0.7, 'r', label='F_crit')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ans:\n",
    "\n",
    "Since p_value > α, we failed to reject H<sub>0. \n",
    "\n",
    "There is not a significant difference in mean Caffeine consumption across the 4 groups of students who vary their preference.\n",
    "\n",
    "The effect size calculated by eta-squared is 0.04 which is medium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "### Question 2\n",
    "\n",
    "The Wine dataset consists of 3 different type of wine. \n",
    "\n",
    "Using an alpha value of α = 0.05, is there difference in the malic_acid across the 3 type of wine?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean_g: [2.01067797 1.93267606 3.33375   ]       mean: 2.425701340813241\n",
      "SSW_g:  [27.49777288 72.19659155 55.626325  ]\n",
      "SSW: 155.32068943065167     SSB: 66.99918243986482     SST: 222.31987187051658\n"
     ]
    }
   ],
   "source": [
    "from sklearn.datasets import load_wine\n",
    "wine = load_wine()\n",
    "data_x = wine.data[:, 1]\n",
    "data_y = wine.target\n",
    "n_g = np.array([0, 0, 0])\n",
    "for y in data_y:\n",
    "    n_g[y] += 1\n",
    "\n",
    "x = {0:data_x[:n_g[0]], 1:data_x[n_g[0]:n_g[0]+n_g[1]], 2:data_x[n_g[0]+n_g[1]:]}\n",
    "\n",
    "alpha = 0.05\n",
    "n = np.sum(n_g)\n",
    "\n",
    "k = 3\n",
    "df_SSB = k - 1\n",
    "df_SSW = n - k\n",
    "df_SST = n - 1\n",
    "\n",
    "mean_g = np.array([np.mean(x[i]) for i in range(k)])\n",
    "mean = np.mean(mean_g)\n",
    "print(\"mean_g: {}       mean: {}\".format(mean_g, mean))\n",
    "\n",
    "SSW_g = np.array([np.sum(np.square(x[i] - mean_g[i])) for i in range(k)])\n",
    "print(\"SSW_g: \", SSW_g)\n",
    "SSW = np.sum(SSW_g)\n",
    "SSB = np.sum(n_g * np.square(mean_g - mean))\n",
    "SST = np.sum([np.sum(np.square(x[i] - mean)) for i in range(k)])\n",
    "print(\"SSW: {}     SSB: {}     SST: {}\".format(SSW, SSB, SST))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "H<sub>0 : There is no difference among group means\n",
    "\n",
    "H<sub>1 : At lease one group differs significantly from the overrall mean of the dependent variable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F:  37.7440280813694\n",
      "F_crit:  3.0476049210928564\n",
      "p_value:  2.353672812205332e-14\n",
      "eta_squared:  0.3013638946269565\n",
      "Beta:  1.8242912179970812e-10\n",
      "Power:  0.9999999998175709\n"
     ]
    }
   ],
   "source": [
    "MSB = SSB / df_SSB\n",
    "MSW = SSW / df_SSW\n",
    "f = MSB / MSW\n",
    "print(\"F: \", f)\n",
    "\n",
    "# dfn: numerator, dfd: denominator\n",
    "f_crit = f_test.ppf(1 - alpha, dfn=df_SSB, dfd=df_SSW)\n",
    "print(\"F_crit: \", f_crit)\n",
    "\n",
    "p_value = 1 - f_test.cdf(f, dfn=df_SSB, dfd=df_SSW)\n",
    "print(\"p_value: \", p_value)\n",
    "\n",
    "eta_squared = SSB / SST\n",
    "print(\"eta_squared: \", eta_squared)\n",
    "\n",
    "nc = SSB / MSW\n",
    "beta = ncf.cdf(f_crit, dfn=df_SSB, dfd=df_SSW, nc=nc)\n",
    "print(\"Beta: \", beta)\n",
    "print(\"Power: \", 1 - beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.5, 4.5)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD7CAYAAAB37B+tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUxElEQVR4nO3df2xV9f3H8dflliu9vZfWzivLRjpHM/PdxEzo/uk2dmllQtKylBS5UCw1dAkaGAMD6Y8vIeB+dFeUBY1Fqlm3EQhgx5fQuCzqaGais0kxGr01M2kNOiGmCFe497b29nK+fzDv1rW9t7S9Pz7wfPxzufeUnvfl0zw5ObfnXptlWZYAAMaalekBAADTQ8gBwHCEHAAMR8gBwHCEHAAMR8gBwHCTCvlnn30mr9ervr6+UY+3t7eroqJCtbW1qq2tVX9/f0qGBABMLCfZF0SjUe3evVtz5swZsy0QCMjv92vhwoUpGQ4AkFzSI3K/36+1a9fqzjvvHLMtEAiora1N69at06FDh1IyIAAgsYRH5CdPnlRhYaGWLFmitra2MdsrKipUU1Mjl8ulLVu2qKurS2VlZQl3aFmWuJYUAG7MrFm2CbfZEl2iv379etlsNtlsNr3//vu66667dPDgQXk8HlmWpVAoJLfbLUk6cuSIgsGgNm/enHCYaDSmYDAyxacCALcmj8c94baER+RHjhyJ/7m2tlZ79uyRx+ORJIVCIVVWVurPf/6znE6nuru7VV1dPUMjAwAmK+mLnf+ts7NTkUhEPp9P27dv14YNG+RwOFRaWiqv15uKGQEACSQ8tZIKnFoBkMxbb/Wos/P/tHLlKi1e/L1Mj5MVpnxqBQAy4cUXj+rDD/s1NDRIyCeBKzsBZJ3BwaFRt0iMkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4SYV8s8++0xer1d9fX2jHj9z5oyqq6vl8/l04sSJlAwIAEgs6Wd2RqNR7d69W3PmzBnzeEtLizo6OpSbm6t169aprKxMHo8nZcMCAMZKekTu9/u1du1a3XnnnaMe7+vrU1FRkfLz8+VwOFRSUqKenp6UDQoAGF/CkJ88eVKFhYVasmTJmG2hUEhutzt+Py8vT6FQaOYnBAAklPDUyp/+9CfZbDb9/e9/1/vvv6+GhgYdPHhQHo9HLpdL4XA4/rXhcHhU2Cdit9tUUOCc/uQAblp2uy1+Sy+SSxjyI0eOxP9cW1urPXv2xM+BFxcX69y5cwoGg3I6nerp6VF9fX3SHcZiloLByDTHBnAzi8Ws+C29uM7jmfhAOemLnf+ts7NTkUhEPp9PjY2Nqq+vl2VZqq6u1rx586Y1KADgxk065IcPH5Z0/Uj8S+Xl5SovL5/5qQAAk8YFQQBgOEIOAIYj5ABgOEIOAIYj5ABgOEIOAIYj5ABgOEIOAIYj5ABgOEIOAIYj5ABgOEIOAIYj5ABgOEIOAIYj5ABgOEIOAIYj5ABgOEIOAIYj5ABgOEIOAIZL+uHLsVhMu3bt0ocffii73a6WlhYVFRXFt7e3t6ujo0OFhYWSpL1792rBggWpmxgAMErSkHd1dUmSjh07pu7ubrW0tOjgwYPx7YFAQH6/XwsXLkzdlACACSUN+bJly7R06VJJ0vnz53XHHXeM2h4IBNTW1qaBgQEtXbpUmzZtSsmgAIDxJQ25JOXk5KihoUGvvPKKnn766VHbKioqVFNTI5fLpS1btqirq0tlZWUTfi+73aaCAuf0pgZwU7PbbfFbepHcpEIuSX6/Xzt27NCaNWv00ksvyel0yrIs1dXVye12S5K8Xq96e3sThjwWsxQMRqY/OYCbVixmxW/pxXUej3vCbUl/a+XUqVM6dOiQJCk3N1c2m012u12SFAqFVFlZqXA4LMuy1N3dzblyAEizpEfkDzzwgJqamrR+/XqNjIyoublZL7/8siKRiHw+n7Zv364NGzbI4XCotLRUXq83HXMDAP4lacidTqcOHDgw4faqqipVVVXN5EwAgBvABUEAYDhCDgCGI+QAYDhCDgCGm/TvkQO4ubjyZyvXMSfTY4zrPy8ISvT705k0ODyk0OfRTI8hiZADt6xcxxz94JkfZHqMcTmCDs3SLH0c/DhrZ3z9Z68rpOwIOadWAMBwhBwADEfIAcBwhBwADEfIAcBwhBwADEfIAcBwhBwADEfIAcBwhBwADEfIAcBwhBwADEfIAcBwSUMei8XU1NSktWvXav369froo49GbT9z5oyqq6vl8/l04sSJlA0KABhf0pB3dXVJko4dO6atW7eqpaUlvi0ajaqlpUW/+93vdPjwYR0/flwDAwOpmxYAMEbSkC9btky/+MUvJEnnz5/XHXfcEd/W19enoqIi5efny+FwqKSkRD09PambFgAwxqQ+WCInJ0cNDQ165ZVX9PTTT8cfD4VCcrv//ekdeXl5CoVCCb+X3W5TQYFziuMCQPbIlpZN+hOC/H6/duzYoTVr1uill16S0+mUy+VSOByOf004HB4V9vHEYpaCwcjUJwYwI7L1I9RMks6WJVqvpKdWTp06pUOHDkmScnNzZbPZZLfbJUnFxcU6d+6cgsGghoeH1dPTo0WLFs3Q2EB6vfVWj/bu/V+99RanB2GWpEfkDzzwgJqamrR+/XqNjIyoublZL7/8siKRiHw+nxobG1VfXy/LslRdXa158+alY25gxr344lF9+GG/hoYGtXjx9zI9DjBpSUPudDp14MCBCbeXl5ervLx8RocCMmFwcGjULWAKLggCAMMRcgAwHCEHAMMRcgAwHCEHAMMRcgAwHCEHAMMRcgAwHCEHAMMRcgAwHCEHAMMR8gzi3fYAzIRJvx85Zh7vtgdgJnBEnkG82x6AmUDIAcBwhBwADEfIAcBwhBwADEfIAcBwCX/9MBqNqrm5WZ988omGh4f16KOP6v77749vb29vV0dHhwoLCyVJe/fu1YIFC1I7MQBglIQhP336tAoKCrRv3z5dvnxZq1atGhXyQCAgv9+vhQsXpnxQAMD4EoZ8xYoVWr58efy+3W4ftT0QCKitrU0DAwNaunSpNm3alJopAQATShjyvLw8SVIoFNLWrVu1bdu2UdsrKipUU1Mjl8ulLVu2qKurS2VlZQl3aLfbVFDgnN7UNwm73Ra/5d8k81gP3Khs+TlJeon+hQsXtHnzZtXU1GjlypXxxy3LUl1dndxutyTJ6/Wqt7c3achjMUvBYGSaY98cYjErfsu/Sebdauvh8bgzPYLx0vlzkmi9Ev7WysWLF7Vx40bt3LlTq1evHrUtFAqpsrJS4XBYlmWpu7ubc+UAkAEJj8ife+45XblyRa2trWptbZUkPfjggxocHJTP59P27du1YcMGORwOlZaWyuv1pmVoAMC/JQz5rl27tGvXrgm3V1VVqaqqaqZnAgDcAC4IAgDDEXIAMBwhBwDDEXIAMBwhBwDDEXIAMBwhBwDDEXIAMBwhB5B9cv7rFgkRcgBZZ+TbI4rdEdPIt0cyPYoR+P8OQNa59tVruvbVa5kewxgckQOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABgu4QVB0WhUzc3N+uSTTzQ8PKxHH31U999/f3z7mTNn9OyzzyonJ0fV1dVas2ZNygcGAIyWMOSnT59WQUGB9u3bp8uXL2vVqlXxkEejUbW0tKijo0O5ublat26dysrK5PF40jI4AOC6hKdWVqxYoZ///Ofx+3a7Pf7nvr4+FRUVKT8/Xw6HQyUlJerp6UndpACAcSU8Is/Ly5MkhUIhbd26Vdu2bYtvC4VCcrvdo742FAol3aHdblNBgXOK4964a7J02+zsfEsZu90Wv/V43Em+OjO+iI5olmyZHiMt/nM90vkzCnNly89J0sJduHBBmzdvVk1NjVauXBl/3OVyKRwOx++Hw+FRYZ9ILGYpGIxMcdwb5/G4VbLzj2nb341wX7wqu6SPLl7N2hnP7tuggYGrmR4jLWIxK36bzp/RTMnWgweTpLtlE0l4auXixYvauHGjdu7cqdWrV4/aVlxcrHPnzikYDGp4eFg9PT1atGjRzEwMAJi0hEfkzz33nK5cuaLW1la1trZKkh588EENDg7K5/OpsbFR9fX1sixL1dXVmjdvXlqGBgD8W8KQ79q1S7t27Zpwe3l5ucrLy2d8KNy8CvNny+6Yk+kxxpXtr1nEhod06fNopsdAFsrOVwFx07I75uijx+/N9BjjGrlUKClHI5fOZeWMRbvflUTIMRZXdgKA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qg5ABhuUiF/5513VFtbO+bx9vZ2VVRUqLa2VrW1terv75/xAQEAiSX9zM7nn39ep0+fVm5u7phtgUBAfr9fCxcuTMlwAIDkkh6RFxUV6Zlnnhl3WyAQUFtbm9atW6dDhw7N+HAAgOSSHpEvX75c//znP8fdVlFRoZqaGrlcLm3ZskVdXV0qKytL+P3sdpsKCpxTmxYZwXplD9Yiu2TLeiQN+UQsy1JdXZ3cbrckyev1qre3N2nIYzFLwWBkqru9YR6PO237ulnN5HqxHtPDWmSXbGnZlH9rJRQKqbKyUuFwWJZlqbu7m3PlAJABN3xE3tnZqUgkIp/Pp+3bt2vDhg1yOBwqLS2V1+tNxYwAgAQmFfL58+frxIkTkqSVK1fGH6+qqlJVVVVKBgMATA4XBAGA4Qg5ABiOkAOA4Qg5ABiOkAOA4Qh5BlmzckbdAsBUEPIMGvraIkVdX9XQ1xZlehQABuNQMING8udrJH9+pscAYDiOyAHAcIQcAAxHyAHAcIQcAAxHyAHAcIQcAAxHyAHAcIQcAAxHyAHAcIQcAAxHyAHAcJMK+TvvvKPa2toxj585c0bV1dXy+Xzxz/QEAKRX0jfNev7553X69Gnl5uaOejwajaqlpUUdHR3Kzc3VunXrVFZWJo/Hk7JhAQBjJT0iLyoq0jPPPDPm8b6+PhUVFSk/P18Oh0MlJSXq6elJyZAAgIklDfny5cuVkzP2wD0UCsntdsfv5+XlKRQKzex0AICkpvx+5C6XS+FwOH4/HA6PCvtE7HabCgqcU90tMoD1yh6sRXbJlvWYcsiLi4t17tw5BYNBOZ1O9fT0qL6+Punfi8UsBYORqe72hnk8yf9zQWIzuV6sx/SwFtklW1p2wyHv7OxUJBKRz+dTY2Oj6uvrZVmWqqurNW/evGkNCgC4cZMK+fz58+O/Xrhy5cr44+Xl5SovL0/NZACASeGCIOBf5titUbeAKQg58C+r7grrf/KHtequcPIvBrLIlF/sBG423/3KsL77leFMjwHcMI7IAcBwhBwADEfIAcBwhBwADEfIAcBwhBwADEfIAcBwhBwADEfIAcBwhBwADEfIAcBwhBwADEfIAcBwhBwADEfIAcBwhBwADEfIAcBwST8h6Nq1a9qzZ4/+8Y9/yOFw6Je//KW+8Y1vxLe3t7ero6NDhYWFkqS9e/dqwYIFqZsYADBK0pC/+uqrGh4e1vHjx/X222/rN7/5jQ4ePBjfHggE5Pf7tXDhwpQOCgAYX9KQnz17VkuWLJEk3XfffXrvvfdGbQ8EAmpra9PAwICWLl2qTZs2pWZSAMC4koY8FArJ5XLF79vtdo2MjCgn5/pfraioUE1NjVwul7Zs2aKuri6VlZVN+P1mz7bL43HPwOiTd3bfhrTu72Yz0+tVtPvdGf1+t5KZXovXf/b6jH6/W026WzaRpC92ulwuhcPh+P1r167FI25Zlurq6lRYWCiHwyGv16ve3t7UTQsAGCNpyBcvXqzXXntNkvT222/r7rvvjm8LhUKqrKxUOByWZVnq7u7mXDkApJnNsiwr0Rd8+VsrH3zwgSzL0q9//Wv19vYqEonI5/Pp1KlTOnz4sBwOh0pLS7V169Z0zQ4A0CRCDgDIblwQBACGI+QAYDhCDgCGI+Qz6OTJk3ryySdT9v0vXbqkjRs3qqamRtu2bdPg4GDK9nUzSPV6fOn3v/99WvZjslSvxfnz5/Xwww+rtrZWDz30kPr7+1O2r2xEyA3S2tqqyspKHT16VN/5znd0/PjxTI90SxsaGtKOHTt09OjRTI9yyztw4IAeeughHT58WJs2bdL+/fszPVJaJb2yExMbGhpSU1OTzp8/r2g0quXLl8e3PfXUU3rvvfcUDodVXFyslpYWnT17Vn6/Xzk5OZo7d66efPJJDQwMqKmpSTk5ObLb7XriiSc0b968cfd39uzZ+Fsg/OhHP9L+/fv18MMPp+OpGiHd6/HFF1+oqqpK3//+92+5I8Bk0r0WDQ0NcruvX2UZi8V02223peV5ZgtCPg3Hjh3T17/+df32t7/VBx98oDfeeENXr15VKBTS3Llz1d7ermvXrqmiokKffvqpXn31Vf34xz9WfX29zpw5oytXruiNN97QPffco8bGRvX09Ojzzz+f8Ic1FArFf1jz8vJ09erVdD7drJfu9cjPz9cPf/hDnTx5Ms3PNPuley2+fPfV/v5++f1+Pfvss+l8uhnHqZVp6O/v13333SdJuvvuuzV37lxJ0m233aZLly7pscce0+7duxWJRBSNRvXII4/o0qVLqqur01/+8hfl5ORo9erVuv322/XTn/5UR44ckd1un3B///l2CeFwOL4/XJfu9cDEMrEWb775pjZv3qwnnnjilnsrbUI+DcXFxXr33etvAPXxxx/Hz8u99tprunDhgvbv36/HHntMQ0NDsixLnZ2dWrVqlQ4fPqxvfetbOnHihP7617+qpKREf/jDH7RixQq98MILE+5v8eLF+tvf/hbfR0lJSeqfpEHSvR6YWLrX4s0339SvfvUrvfDCC7r33nvT8hyzCVd2TsMXX3yh5uZmffrpp4rFYlq2bJkuX76suro6PfLII7Lb7XI4HPHzhTk5Odq7d6+cTqdmz56txx9/XJZlaefOnbLb7Zo1a5aampp0zz33jLu/ixcvqqGhQeFwWLfffrueeuopOZ3OND/r7JXu9fjSyZMn1d/frx07dqTpmWa/dK/FT37yEw0PD8vj8UiSvvnNb+rxxx9P51POKEIOAIbjxc4sMzw8rPr6+jGP32pHGNmC9cgerMXEOCIHAMPxYicAGI6QA4DhCDkAGI6QA4DhCDkAGO7/AfSQKuzgfisvAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(1)\n",
    "plot = sns.barplot(data=[x[i] for i in range(k)], ci='sd')\n",
    "plot.set_xticklabels(wine.target_names)\n",
    "plot.margins(0.1)\n",
    "plot.set_ylim(0.5, 4.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ans:\n",
    "\n",
    "Since p_value < α, we rejected H<sub>0. \n",
    "\n",
    "There is significant difference in the mean of malic_acid across the 3 type of irises.\n",
    "\n",
    "The effect size calculated by eta-squared is 0.3 which is large.\n",
    "\n",
    "β (Type II Error) = 0 which means if we rejected H1, there is a 0% chance that we are wrong.\n",
    "\n",
    "The power of this analysis is 1 which mean that we have a 100% chance that we correctly reject H0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Post Hoc Test:\n",
    "\n",
    "Since we knew there was difference in group means, we wanted to find which group differs from other groups. To compare all group pairs, we employed a Tukey HSD test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Multiple Comparison of Means - Tukey HSD, FWER=0.05\n",
      "===================================================\n",
      "group1 group2 meandiff p-adj   lower  upper  reject\n",
      "---------------------------------------------------\n",
      "     0      1   -0.078 0.8855 -0.4703 0.3143  False\n",
      "     0      2   1.3231    0.0  0.8902 1.7559   True\n",
      "     1      2   1.4011    0.0  0.9849 1.8172   True\n",
      "---------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd\n",
    "\n",
    "tukey = pairwise_tukeyhsd(endog=data_x, groups=data_y, alpha=alpha)\n",
    "print(tukey)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "### Question 3\n",
    "\n",
    "They hypothesized that the more people merely watch others, the more they believe they can perform the skill themselves.\n",
    "\n",
    "They compared extensively watching the “tablecloth trick” with extensively reading or thinking about it. Participants assessed their own abilities to perform the “tablecloth trick.” \n",
    "\n",
    "Participants were asked to rate from 1 to 7. \n",
    "\n",
    "The study was a 3 (type of exposure: watch, read, think) ? 2 (amount of exposure: low, high) between-subjects design.\n",
    "\n",
    "* From [REFERENCE](https://www.alexanderdemos.org/ANOVA9.html#Basic_Rules)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[1 3 3 4 2]\n",
      "  [3 2 3 1 4]\n",
      "  [3 3 2 3 1]]\n",
      "\n",
      " [[6 4 7 5 5]\n",
      "  [5 3 2 3 2]\n",
      "  [3 2 1 4 4]]]\n"
     ]
    }
   ],
   "source": [
    "### data\n",
    "num_cols = 3\n",
    "num_rows = 2 \n",
    "alpha = 0.05\n",
    "\n",
    "data = {\"rate\": [1, 3, 3, 4, 2, \n",
    "                 3, 2, 3, 1, 4,\n",
    "                 3, 3, 2, 3, 1,\n",
    "                 6, 4, 7, 5, 5, \n",
    "                 5, 3, 2, 3, 2, \n",
    "                 3, 2, 1, 4, 4],\n",
    "        \"type\": ([\"Watching\"]*5 + [\"Reading\"]*5 + [\"Thinking\"]*5) * 2,\n",
    "        \"exposure\": [\"Low\"]*15 + [\"High\"]*15}\n",
    "df = pd.DataFrame(data)\n",
    "data_array = np.array(data[\"rate\"]).reshape(2, 3, 5)\n",
    "print(data_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean_cell: \n",
      "[[2.6 2.6 2.4]\n",
      " [5.4 3.  2.8]]\n",
      "SSW_cell: \n",
      "[[5.2 5.2 3.2]\n",
      " [5.2 6.  6.8]]\n",
      "\n",
      "n_grand: 30\n",
      "n_rows: [15 15]\n",
      "n_cols: [10 10 10]\n",
      "\n",
      "mean_grand: 3.1333333333333333\n",
      "mean_rows: [2.53333333 3.73333333]\n",
      "mean_cols: [4.  2.8 2.6]\n",
      "\n",
      "SSW_c: 31.6000\tSSB_col: 11.4667\tSSB_row: 10.8000\n",
      "SST: 1386.8456\tSS_rxc: 1332.9789\n",
      "\n",
      "df_SSW_c: 24\tdf_SSB_col: 2\tdf_SSB_row: 1\n",
      "df_SST: 29\tdf_SS_rxc: 2\n"
     ]
    }
   ],
   "source": [
    "n_cell = 5 * np.ones((2, 3), dtype='int')\n",
    "mean_cell = np.mean(data_array, axis=2)\n",
    "SSW_cell = np.sum(np.square(data_array - mean_cell.reshape(2, 3, 1)), axis=2)\n",
    "print(\"mean_cell: \")\n",
    "print(mean_cell)\n",
    "print(\"SSW_cell: \")\n",
    "print(SSW_cell)\n",
    "print(\"\")\n",
    "\n",
    "n_grand = np.sum(n_cell)\n",
    "n_cols = np.sum(n_cell, axis=0)\n",
    "n_rows = np.sum(n_cell, axis=1)\n",
    "print(\"n_grand: {}\".format(n_grand))\n",
    "print(\"n_rows: {}\".format(n_rows))\n",
    "print(\"n_cols: {}\".format(n_cols))\n",
    "print(\"\")\n",
    "\n",
    "mean_grand = np.average(mean_cell, weights=n_cell)\n",
    "mean_rows = np.average(mean_cell, axis=1, weights=n_cell)\n",
    "mean_cols = np.average(mean_cell, axis=0, weights=n_cell)\n",
    "print(\"mean_grand: {}\".format(mean_grand))\n",
    "print(\"mean_rows: {}\".format(mean_rows))\n",
    "print(\"mean_cols: {}\".format(mean_cols))\n",
    "print(\"\")\n",
    "\n",
    "SSW_c = np.sum(SSW_cell)\n",
    "SSB_col = np.sum( n_cols * np.square(mean_cols - mean_grand) )\n",
    "SSB_row = np.sum( n_rows * np.square(mean_rows - mean_grand) )\n",
    "SST = 1386.8456\n",
    "SS_rxc = SST - SSW_c - SSB_col - SSB_row\n",
    "print(\"SSW_c: {:.4f}\\tSSB_col: {:.4f}\\tSSB_row: {:.4f}\".format(SSW_c, SSB_col, SSB_row))\n",
    "print(\"SST: {:.4f}\\tSS_rxc: {:.4f}\".format(SST, SS_rxc))\n",
    "print(\"\")\n",
    "\n",
    "df_SSW_c = n_grand - num_cols * num_rows\n",
    "df_SSB_col = num_cols - 1\n",
    "df_SSB_row = num_rows - 1\n",
    "df_SS_rxc = (num_cols - 1) * (num_rows - 1)\n",
    "df_SST = n_grand - 1\n",
    "print(\"df_SSW_c: {}\\tdf_SSB_col: {}\\tdf_SSB_row: {}\".format(df_SSW_c, df_SSB_col, df_SSB_row))\n",
    "print(\"df_SST: {}\\tdf_SS_rxc: {}\".format(df_SST, df_SS_rxc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Main Effects:\n",
    "\n",
    "H<sub>0 : All groups have equal means\n",
    "\n",
    "H<sub>1 : At least 1 group is different from the other\n",
    "\n",
    "Interaction:\n",
    "\n",
    "H<sub>0 : There is no interaction between factors A and B. All the mean differences between treatment conditions are explained by the main effects.\n",
    "\n",
    "H<sub>1 : There is an interaction between factors A and B. The mean differences between treatment conditions are not what would be predicted from the overall main effects of the two factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSW_c: 1.3167\tMSB_col: 5.7333\tMSB_row: 10.8000\n",
      "MS_rxc: 666.4895\n",
      "\n",
      "f_col: 4.354\tf_col_crit: 3.403\n",
      "p_value_col: 0.02435\n",
      "eta_2_col: 0.0083\n",
      "Beta_col: 0.3018\tPower_col: 0.6982\n",
      "\n",
      "f_row: 8.203\tf_row_crit: 4.260\n",
      "p_value_row: 0.00855\n",
      "eta_2_row: 0.0078\n",
      "Beta_row: 0.2154\tPower_row: 0.7846\n",
      "\n",
      "f_rxc: 506.195\tf_rxc_crit: 3.403\n",
      "p_value_rxc: 0.00000\n",
      "eta_2_rxc: 0.9612\n",
      "Beta_rxc: 0.0000\tPower_row: 1.0000\n"
     ]
    }
   ],
   "source": [
    "MSW_c = SSW_c / df_SSW_c\n",
    "MSB_col = SSB_col / df_SSB_col\n",
    "MSB_row = SSB_row / df_SSB_row\n",
    "MS_rxc = SS_rxc / df_SS_rxc\n",
    "print(\"MSW_c: {:.4f}\\tMSB_col: {:.4f}\\tMSB_row: {:.4f}\".format(MSW_c, MSB_col, MSB_row))\n",
    "print(\"MS_rxc: {:.4f}\".format(MS_rxc))\n",
    "print(\"\")\n",
    "\n",
    "\n",
    "f_col = MSB_col / MSW_c\n",
    "f_col_crit = f_test.ppf(1 - alpha, dfn=df_SSB_col, dfd=df_SSW_c)\n",
    "p_value_col = 1 - f_test.cdf(f_col, dfn=df_SSB_col, dfd=df_SSW_c)\n",
    "eta_2_col = SSB_col / SST\n",
    "### non-central parameter\n",
    "nc_col = SSB_col / MSW_c\n",
    "beta_col = ncf.cdf(f_col_crit, dfn=df_SSB_col, dfd=df_SSW_c, nc=nc_col)\n",
    "print(\"f_col: {:.3f}\\tf_col_crit: {:.3f}\".format(f_col, f_col_crit))\n",
    "print(\"p_value_col: {:.5f}\".format(p_value_col))\n",
    "print(\"eta_2_col: {:.4f}\".format(eta_2_col))\n",
    "print(\"Beta_col: {:.4f}\\tPower_col: {:.4f}\".format(beta_col, 1 - beta_col))\n",
    "print(\"\")\n",
    "\n",
    "\n",
    "f_row = MSB_row / MSW_c\n",
    "f_row_crit = f_test.ppf(1 - alpha, dfn=df_SSB_row, dfd=df_SSW_c)\n",
    "p_value_row = 1 - f_test.cdf(f_row, dfn=df_SSB_row, dfd=df_SSW_c)\n",
    "eta_2_row = SSB_row / SST\n",
    "### non-central parameter\n",
    "nc_row = SSB_row / MSW_c\n",
    "beta_row = ncf.cdf(f_row_crit, dfn=df_SSB_row, dfd=df_SSW_c, nc=nc_row)\n",
    "print(\"f_row: {:.3f}\\tf_row_crit: {:.3f}\".format(f_row, f_row_crit))\n",
    "print(\"p_value_row: {:.5f}\".format(p_value_row))\n",
    "print(\"eta_2_row: {:.4f}\".format(eta_2_row))\n",
    "print(\"Beta_row: {:.4f}\\tPower_row: {:.4f}\".format(beta_row, 1 - beta_row))\n",
    "print(\"\")\n",
    "\n",
    "\n",
    "f_rxc = MS_rxc / MSW_c\n",
    "f_rxc_crit = f_test.ppf(1 - alpha, dfn=df_SS_rxc, dfd=df_SSW_c)\n",
    "p_value_rxc = 1 - f_test.cdf(f_rxc, dfn=df_SS_rxc, dfd=df_SSW_c)\n",
    "eta_2_rxc = SS_rxc / SST\n",
    "### non-central parameter\n",
    "nc_rxc = SS_rxc / MSW_c\n",
    "beta_rxc = ncf.cdf(f_rxc_crit, dfn=df_SS_rxc, dfd=df_SSW_c, nc=nc_rxc)\n",
    "print(\"f_rxc: {:.3f}\\tf_rxc_crit: {:.3f}\".format(f_rxc, f_rxc_crit))\n",
    "print(\"p_value_rxc: {:.5f}\".format(p_value_rxc))\n",
    "print(\"eta_2_rxc: {:.4f}\".format(eta_2_rxc))\n",
    "print(\"Beta_rxc: {:.4f}\\tPower_row: {:.4f}\".format(beta_rxc, 1 - beta_rxc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x267d54892e0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAESCAYAAADpO/4pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdgklEQVR4nO3de1xUZf4H8M9cGEBuw12URDJdb3nDLfF+6eeqr4osdDUbM1KTSFPSUiREdwEVzRJD7OJ6XVlBbHW3zEqzMjUUU9dwDfCGqAhCijBX5veH22yugCjnOPDwef8zzsyZ5/meOc7Hx+eceUZhtVqtICIiYSjtXQAREUmLwU5EJBgGOxGRYBjsRESCYbATEQmGwU5EJBgGO9E9SElJwaJFiwAAQ4cOxYkTJ+xcEdGdGOxERIJR27sAInvLzMzEX/7yFyiVSnh6emLJkiXIzc3F6tWrYTKZ4OTkhLfeegs9e/a0d6lE9cJgp2bt1KlTWLZsGbZv346AgACsW7cOc+fOxbVr17BhwwZ4enri559/xksvvYTdu3fbu1yiemGwU7N24MAB9O/fHwEBAQCASZMmwcHBAStXrsSkSZNs2ykUCpw/f95OVRLdGwY7NWsqlQoKhcJ2X6/Xo7q6GqGhoXj33Xdtj1+6dAl+fn744osv7FAl0b3hyVNq1h5//HEcOHAAxcXFAID09HTs3bsX+/fvR35+PgBg3759ePrpp6HX6+1ZKlG9ccROzdrvfvc7zJkzB5MnTwYA+Pr6IikpCTk5OYiOjobVaoVarcbq1avh4uJi52qJ6kfBZXuJiMTCqRgiIsEw2ImIBCPLHHtWVha2b98OADAYDMjNzcX+/fvh7u4uR3dERPQbss+xL1y4EB07dsQf//hHObshIqL/kPWqmBMnTiAvLw8LFiyoczur1QqewiUiujdKpaLGx2UN9jVr1iAqKuqu25nN1Sgvr5SzFCIi4fj6utX4uGwnT69fv46CggL06dNHri6IiKgGsgV7dnY2+vbtK1fzRERUC9mC/cyZMwgMDJSreSIiqkWj+OapyWThHDsR3TeLxYyysqswm432LkUWarUGnp6+UKluPy1a2xw714ohoiavrOwqnJxawMWl5W2rdYrAarXi5s3rKCu7Ch+fgHq9ht88JaImz2w2wsXFXbhQB279FoCLi/s9/W+EwU5EQhAx1H91r/vGYCciEgyDnYhIMAx2IiLB8KoYIhKe2WxGcnIiCgsvoLq6GhMmTERa2iosXJgEpVKJ+PgYpKZ+hFdeeQnduvXAmTMFcHd3R3x8IhwcHJCUtBAXL16ExWLBuHETMGzYcGRlZeCzz/4BpVKJbt16ICrqdSQkxGPYsOHo06cvDh78Hl99tRvz58fjueeeRFBQWwQFBWPcuAlYujQRRqMBGo0j3nwzBv7+LSXdXwY7EQlv585P4OGhxbx5cfjll3JERU1FTMwCLFnyZ1itVsTGLoSLiyv0ej2GDx+JHj16ITX1Pfz979vg4OAADw8t3n77T6isvImIiBcQEvIYPv10J2bOnIOuXR/F9u2ZMJvNtfZfXHwFa9dugoeHFnFx8xAe/keEhvbD4cM/IC1tFRYs+LOk+8tgJyLh5efn4fjxo/jpp38BuPWFptatA+Hq6gYHBzXat/8dAECtVqNHj14AgK5du+Pgwf1QqdTo3fsxAECLFi5o2zYYFy8WIiYmDlu2bEJaWgq6dHn0jj5/+91PDw8tPDy0AICCgjxs3PgXbN683tan1BjsRCS8oKC28PPzw8SJETAY9Fi/fi0OH85GixbOqK62Yu/eLzFkyBMwm834+efTaN++A06cOIbg4HZQKhU4fvwoBg0agsrKm8jPz0erVq2wfv1azJ49D46OjoiOfg0nThyDRqNBaWkJAOD06VO2/pXK/57ObNOmLcaPfwGPPtod586dxdGjRyTfXwY7EQkvLOxZLFnyZ7z22lTcvFmBAQMGY+3aNXj//Y9gtVbj1Vcno1OnLgCAzZvX48qVy/D3b4kpUyKhUCiwZMmfERn5MgwGAyIipsDT0wvt2j2CKVMmQqv1hK+vLzp37gpnZ2ckJS3C7t278NBDbWqsJSrqdSxfvhhGoxEGgx6vvz5b8v3lWjFE1ORdvnwOLVsGNbid8PCnsHlzJhwdHSWoSlo17eMDX4+diIjsg1MxRET/kZm5094lSIIjdiIiwTDYiYgEw2AnIhIM59iJSDiu7k5wdnSQrL0qgwkV1/WStSc3BjsRCcfZ0QEhczZI1t6R5ImoQN3BnpNzGH//+zYsXJgkWb/3i1MxRESCYbATEckkO/sgpkx5Ea+9NhUxMXNw48YNzJv3Bk6d+gkAMH78s9i3by8AYNasKFy9WixJv5yKISKSgdVqxdKliUhN/Qi+vn7YunUL1q//GAMHDsHBg9/D3d0DGo0jsrMPISTk9zAajfD19ZOkb47YBZeTcxgLF85HTs5he5dC1KyUl5ejRQsXW1j36NETZ84UoF+/gcjOPoRDhw5gwoQXkZt7EgcP7ke/fgMk65vBLriMjL8iN/ckMjL+au9SiJoVrVaLysqbKCm5tdrjjz/m4KGH2sDd3R2Ojk746qvd6NMnFP7+LbF16xYMGjRUsr45FSO4qir9bbdEzUGVwYQjyRMlba8+fvjhEF5+WWe7/8ILkzB//hwolQq4ubkjJiYeADBgwCB8+ukOuLt74LHH+mD79ky0bh0oWb1c3VFwM2e+isuXi9CyZSu8+26qvcshkoVUqzs2ZlzdkYioGWOwExEJRrY59jVr1mDPnj0wmUwYP348xowZI1dXRET0G7IE+6FDh3D06FFs2bIFVVVVWLt2rRzdEBFRDWQJ9u+++w4dOnRAVFQUKioq8Oabb8rRDRER1UCWYC8rK0NRURHS0tJQWFiIyMhI7Nq1CwqFosbtVSoFtNoWcpTS7KlUCtst32MS1ZUrCqhU/z1l6OGqgkrjJFn7FqMev1RYJGvvfigU9f8MyxLsWq0WDz/8MDQaDR5++GE4Ojri2rVr8Pb2rnF7i8XKyx1lYrFYbbd8j0lUVqsVFku17b5K44Lzix6VrP02cSdgsdyoc5uaVndcvToF3t7euHnzJl56aUqNr/v00504d+4sIiOn19m+1XrnZ/iBXu4YEhKCb7/9FlarFVeuXEFVVRW0Wq0cXRERNWqurm61hrpcZBmxDxkyBNnZ2QgPD4fVakVcXBxUKpUcXRERNXoLFszDwoVJ+Mc/PsG2bVvh7u4BtdoBw4b9HwDg5MkTmDUrCuXlZXjmmXCEhT3boP5ku9yRJ0yJqLk5cuQwXnttqu1+UdFFTJ48DcCtRcE2bdqAdev+CgcHB8yYMc22nVqtxjvvrMLly5cwZ87rjTfYiYiam5CQ3nfMsf+qsPACgoOD4eR066Ru167dbM916NARCoUCXl7e0Osbvq4Tv3lKRPQABAY+hHPnzsJg0KO6uhq5uSdtz9V2xeD94oidiIRjMerRJu6EpO01lFarxYQJL+LVV6fA3d0dBoMBarUaZrNZggpvx2AnIuFc+8UEoH5L7UqlV6/e6NWr922P/XoJ46hRT8FsNqOk5Co+/ngjACAqagr8/PzRo0cv2/aOjo7IzNzZ4FoY7ERED4BarYZer0dExASo1Q7o3LkrunfvKU9fsrRKRER3eOWVKLzySpTs/fDkKREJoRH8ZpBs7nXfGOxE1OSp1RrcvHldyHC3Wq24efM61GpNvV/DqRgiavI8PX1RVnYVFRXl9i5FFmq1Bp6evvXfXsZaiIgeCJVKDR+fAHuX0WhwKoaISDAMdiIiwTDYiYgEw2AnIhIMg52ISDAMdiIiwfByx0bIy8NBsh/i/e2PWdf2+4j3wmLU/2eBJSJqrBjsjZBK4yTZD/Gar3kBUMN87Zwkbd5aCpXBTtSYcSqGiEgwDHYiIsEw2ImIBMNgJyISDIOdiEgwDHYiIsEw2ImIBMNgJyISDIOdiEgwDHYiIsHItqTAM888Aze3W2uTBAYGIikpSa6uiIjoN2QJdoPBAADYuHGjHM0TEVEdZJmKOXXqFKqqqhAREYGJEyfixx9/lKMbIiKqgSwjdicnJ7z88ssYM2YMzp49iylTpmDXrl1Qq2vuTqVSQKttIUcpJAMeK6LGTZZgDw4ORlBQEBQKBYKDg6HVanH16lUEBATUuL3FYkV5eaUcpTRJUqybLiceK6LGobaskGUqJjMzE4sXLwYAXLlyBRUVFfD19ZWjKyIi+h+yjNjDw8Mxb948jB8/HgqFAomJibVOwxARkbRkSVuNRoPly5fL0TQREd0Fv6BERCQYBjsRkWAY7EREgmGwExEJhsFORCQYBjsRkWAY7EREgmGwExEJhsFORCQYBrvgnFTW226JSHwMdsGNbnsTHT2MGN32pr1LIaIHhCtzCa67txHdvY32LoOIHiCO2ImIBMNgJyISDIOdiEgwDHYiIsEw2ImIBMNgJyISDIOdiEgwDHYiIsEw2ImIBMNgJyISDIOdiEgwDHYiIsEw2ImIBMNgJyISTL2W7T19+jTi4+Nx48YNPPXUU2jfvj2GDBkid21ERHQf6jViT0hIQFJSErRaLcLDw5GSkiJ3XUREdJ/qPRUTFBQEhUIBLy8vuLi43HX70tJSDBo0CPn5+Q0qkIiI7k29gt3DwwPp6emoqqrCP//5T7i7u9e5vclkQlxcHJycnCQpkoiI6q9ewZ6YmIjCwkJ4enriX//6FxITE+vcfsmSJRg3bhz8/PwkKZJIBDk5h7Fw4Xzk5By2dykkuHqdPN2wYQNmz55tu798+XK88cYbNW6blZUFLy8vDBgwAB988EG9ilCpFNBqW9RrW7I/Hqv7k5WVjry8PJhMBgwdOtDe5ZDA6gz2jIwMZGZmIj8/H9988w0AoLq6GiaTqdZg37ZtGxQKBQ4cOIDc3Fy89dZbWL16NXx9fWvtx2Kxory8sgG7IRZfXzd7l1AnHqv7U1FRabvle0hSqC0r6gz2sLAwhIaGYs2aNZg2bRoAQKlUwtvbu9bXbN682fZnnU6H+Pj4OkOdiIikVeccu0ajQWBgIOLi4lBcXIyioiJcuHABu3fvflD1ERHRParXHPv06dNhMplQXFwMi8UCPz8/PPnkk3d93caNGxtcIBER3Zt6XRVTUVGBjz/+GN26dUNWVhYMBoPcdRER0X2qV7CrVCoAQFVVFZycnGAymWQtiohIbiJfflqvqZhhw4Zh1apV6NixI8aOHQtXV1e56yIiklVGxl9x5kwB9Poq9OrV297lSKpewf7555/brnYZNGgQ2rZtK2dNRESyq6rS33YrknoFu0KhQFRUFIKDg6FU3pq9iY6OlrUwIiK6P/UK9ueee07uOoiISCL1CvbRo0fLXQcREUmEv6BERCQYBjsRkWAY7EREgmGwExEJhsFORCQYBjsRkWAY7EREgmGwExEJhsFORCQYBnsDibz0Z1PDY0F0S72WFKDaibz0Z1PDY0F0C0fsDSTy0p9NDY8F0S0MdiIiwTDYiYgEw2AnIhIMg52ISDAMdiIiwTDYiYgEw2AnIhIMg52ISDAMdiIiwTDYiYgEI8taMRaLBbGxsThz5gxUKhWSkpLQpk0bOboiIqL/IcuIfe/evQCA9PR0zJgxA0lJSXJ0Q0SNGFfbtB9ZRuxPPPEEBg8eDAAoKiqCj49PndurVApotS3kKEV2KpXCdttU9+FeNdb9lONYqGCG0sFRmrZ+U5+vr5skbVabDLA00kVas7LSkZeXB5PJgKFDB9q7nDuI/NmV7W+EWq3GW2+9hS+++AIrV66sc1uLxYry8kq5SpGVxWK13Uq1D1J96OXSWI+VXMfi/KJHJWnLfM0LgBrma+cka7NN3AmUXr0hSVtSq6iotN02xr8zcvx9edBqywpZT54uWbIEn3/+Od5++21UVjbNN46IqKmRJdg/+eQTrFmzBgDg7OwMhUIBlUolR1dERPQ/ZJmKGT58OObNm4cJEybAbDYjJiYGjo7SzFMSEVHdZAn2Fi1a4L333pOjaSIiugt+QYmISDAMdiIiwTTOC2CJiGrg5eEAlcZJkrbk+F6BxajHtV9MkrTVEAx2ImoyVBqnRv+9AsD+wc6pGCIiwTDYiYgEw2AnIhIMg52ISDAMdiIiwTTLq2Jc3Z3g7OggSVtyXDLVnPBYEEmvWQa7s6MDQuZskKQtt5IbUAE4X3JDsjaPJE+UpJ2mgMeCSHqciiEiEgyDnYhIMAx2IiLBMNiJiATDYCciEkyzvCqGiGrGy0/FwGAnIhtefioGTsUQEQmGwU5EJBgGOxGRYBjsRESCYbATEQmGwU70gDiprLfdEsmFwU70gIxuexMdPYwY3famvUshwfE6dqIHpLu3Ed29jfYug5oBjtiJiATDYCciEozkUzEmkwkxMTG4ePEijEYjIiMjMWzYMKm7ISKiWkge7Dt27IBWq0VycjLKysowevRoBjsR0QMkebCPGDECf/jDH2z3VSqV1F0QEVEdJA92FxcXAEBFRQVmzJiBmTNn3vU1KpUCWm0LqUshmfBYNS48Ho1LYzgeslzueOnSJURFReH555/HU089ddftLRYryssr5SilRlwbumGkPFY8Fg3H49G4NIYskzzYS0pKEBERgbi4OISGhkrdPBER3YXklzumpaXh+vXrSE1NhU6ng06ng16vl7obIiKqheQj9tjYWMTGxkrdLBER1RO/oEREJBgGOxGRYBjsJAyrUn3bLVFdRF5GmcFOwtC36gmTa0voW/W0dynUBIi8jDKHNiQMs0cgzB6B9i6DmgiRl1HmiJ2ISDAMdiIiwTDYiYgEw2AnIhIMg72BeIkdUc342bAfBnsD8RI7oprxs2E//Ke0gXiJHVHN+NmwH47YiYgEw2AnIhIMg52ISDAMdiIiwTDYiYgEw2AnIhIMg52ISDAMdiIiwTDYiYgEw2AnIhIMg52ISDAMdiIiwTDYiYgEw2AnIhIMg52ISDAMdiIiwcgW7MeOHYNOp5OreSIiqoUsv6D04YcfYseOHXB2dpajeSIiqoMsI/Y2bdogJSVFjqaJiOguFFar1SpHw4WFhYiOjsbWrVvlaJ6IiGrBk6dERIJhsBMRCYbBTkQkGNnm2ImIyD44YiciEgyDnYhIMAx2IiLBMNgbaOLEiTh+/DgAwGg0IiQkBB9//LHt+RdeeAGnTp2643VFRUXYs2dPre3qdDrk5+ff9lhubi5WrVolUeViO3ToEEJDQ6HT6aDT6fDss89ixowZMBqN993m2LFjUVhYiKysLHz11VcSViu2xYsXQ6fTYcSIERg8eDB0Oh369OmDWbNm3bFtQkICioqKam1r6NChMBgMtz02a9asBh1XEcmypEBz0r9/fxw+fBjdunXDkSNH0L9/f3z99dd4+eWXYTAYcOnSJXTs2PGO1x08eBAFBQUYOnRovfvq1KkTOnXqJGX5QuvTpw9WrFhhu//GG29gz549GDFiRIPaffbZZxtaWrMyd+5cAEBWVhYKCgowe/ZsHDp0COnp6XdsO3/+/Htu/7fHmG5hsDdQ3759kZqaioiICOzbtw9jxozBsmXLcOPGDZw8eRKPPfYY5s+fj8uXL6OsrAwDBw7E9OnT8cEHH0Cv16Nnz57w8fFBQkICrFYr/P39sWzZMgDA+++/j5KSElRVVeGdd95BUVER0tPTsWLFCgwfPhy9evXCmTNn4O3tjZSUFJhMJrz55psoLi5GQEAAsrOz8d1339n5HWocjEYjiouL4eHhgeXLlyM7OxtWqxWTJk3CyJEj8cMPP9j+N6TX67FkyRIEBwdjxYoV+Pbbb9GyZUuUlZUBAFJSUuDj44OHH34YH374IRwcHFBYWIhRo0YhMjIS586dw9y5c6FWq9G6dWtcvHgRGzdutOfuN0rnzp3D5MmTce3aNQwZMgTTp0+HTqdDfHw8Pv30UxQWFqK0tBRFRUWYN28eBgwYYHvtli1bsH//frzzzjsYMWIEPvvsMyxYsAAajQYXL15EcXExFi9ejC5duiAjIwObN2+Gh4cHHBwcMGrUKOH/cWawN1Dnzp1RUFAAq9WK7OxsREdHIzQ0FN9//z3+/e9/Izg4GJ6enhgzZgwMBgMGDhyImTNnYurUqSgoKMCwYcPw9NNPY8WKFWjXrh02b95sm4IZNGgQwsLCkJKSgl27dqFbt262fi9cuID169cjICAA48aNw4kTJ3Ds2DEEBgZi5cqVyM/Px5NPPmmvt6VROHjwIHQ6HUpLS6FUKjF27FgYjUYUFhYiPT0dBoMBY8eORb9+/fDzzz8jOTkZ/v7+SEtLw65duzBs2DBkZ2cjMzMTlZWVGD58+B19FBUVYceOHTAajRgwYAAiIyOxdOlSTJs2DYMGDcLWrVtx8eJFO+x942cwGJCamgqLxYLBgwdj+vTptz2v0Wjw0UcfYf/+/Vi7dq0t2Ddu3Ijc3Fy89957UKlUt72mVatWWLRoEbZu3Yq//e1vmDlzJj766CN88skn0Gg0mDhx4gPbP3tisDeQUqlEx44d8c0338DX1xcajQYDBw7E119/jVOnTiEpKQlpaWk4ePAgXF1da5wLLC0tRbt27QAAEyZMsD3etWtXAICPjw9KSkpue42npycCAgIAAAEBATAYDMjPz8fAgQMBAO3atYOXl5cs+9xU/DoVU1ZWhoiICAQGBuL06dM4efKkbUlps9mMoqIi+Pv7IyEhAS1atMCVK1fQq1cv5OXloWvXrlAqlXB1dUWHDh3u6KNDhw5Qq9VQq9VwcnICAOTn56Nnz54AgJCQEOzcufPB7XQT0r59e2g0GgCAWn1nFP067diyZcvbPjcHDhyASqW6I9T/9zU5OTk4f/482rVrZ1tp9tfjIjqePJVAv379sGbNGtuIIiQkBD/99BMA4Msvv4SbmxuWL1+OiIgI6PV6WK1WKJVKVFdXAwD8/Pxw9uxZAMAHH3yAL7744q59KhSKOx7r0KEDjh49CgA4f/68beqgufP09ERycjJiY2Ph4+ODxx9/HBs3bsT69esxcuRIBAYGIjY2FomJiVi8eDH8/PxgtVoRHByM48ePo7q6GpWVlcjLy7uj7bsdh2PHjsm+f01VTe9dfZ5PTU2Fu7s7tmzZctfXtGnTBgUFBdDr9aiurrZd6CA6jtgl0LdvX8TGxmLp0qUAbv0X0s3NDZ07d0ZoaCiio6Nx5MgRODs7IygoCMXFxejQoQNWr16NLl26YOHChYiJiYFSqYSvry8mTZqEDRs23HMd4eHhmDt3LiZMmIBWrVrB0dFR6l1tsh555BHodDrs3bsXAQEBeP7551FZWYknnngCrq6uCAsLw9ixY+Hu7g4fHx8UFxejU6dOGDFiBMLDw+Hn5wdvb+969TV79mzExMRg7dq1cHNzq3E0Sg0TGxuLMWPGIDQ0tM7tvLy8MGXKFDz//PPQarUwGAzN4nhwSQGB5OTkoLKyEv3798fZs2cxefJkfPnll/Yuq9nZsWMHunfvjqCgIGRkZCAnJwdJSUn2LqtZMpvN+PDDDxEZGQng1lTnzJkz8fvf/97OlclL/H+6mpGHHnoI0dHRWLVqFcxmM+Li4uxdUrMUEBCAWbNmwdnZGUqlEomJifYuqdlSq9WoqqrC6NGj4eDggG7duqF37972Lkt2HLETEQmGJ0+JiATDYCciEgyDnYhIMAx2arYMBgMyMjLsXQaR5Bjs1GxdvXqVwU5C4uWO1GylpaUhLy8PHTt2xM6dO9G+fXvs27cPX3/9Nby8vFBQUIDS0lJcv34dsbGx6N27Nz777DOsW7cOSqUSISEhmD17tr13g+gOHLFTszVt2jQ88sgjSEhIwPbt2wEA27ZtQ3h4OADAyckJGzZsQHJyMhYtWoTy8nKkpKRg3bp12LJlC65cuYL9+/fbcxeIasRgp2Zv1KhR2LNnD0pLS3H58mV06dIFwK1FxIBbi1WVlJTg/PnzuHbtGqZOnWr7IZQLFy7Ys3SiGnEqhpqtXxdic3Z2xuOPP46EhASEhYXZnj958iTCwsJw+vRp+Pv7IzAwEAEBAVi7di0cHByQlZXFHz6hRonBTs2Wt7c3TCYTkpOTMXbsWIwfPx7x8fG253Nzc/Hiiy+iqqoKf/rTn+Dl5YVJkyZBp9PBYrGgdevWGDlypP12gKgWXFKACMDx48exadMm2wqdv/5K0vjx4+1cGdG944idmr1NmzZh27ZtWLlypb1LIZIER+xERILhVTFERIJhsBMRCYbBTkQkGAY7EZFgGOxERIL5f2L3KcUMwxWBAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEFCAYAAAD+A2xwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuTklEQVR4nO3deXxTdb7/8dc5J3vTNi1tWYSyFAGXq2yzoBaVkavooCJqC3NxGea6XL3jCIPL/BQVF3C5Dvei4siM43YVFLkKzKgj44Iig4Kg1hFE0Cp7gZY2adMs5/z+SJq2FJqWNk1P8nk+7CNn68nnm+A7335zFsUwDAMhhBCmpSa7ACGEEB0jQS6EECYnQS6EECYnQS6EECYnQS6EECYnQS6EECZnSXYBQrTV0KFDGTJkCKravP/x+OOP07dv3yRVJUTySZALU3n22WfJzc1NdhlCdCsS5CIl/N///R+PP/44r7/+OoqiMHnyZK699lp69+7NI488Qp8+fdi+fTsOh4N58+ZRVFRETU0N99xzD5s3b0ZRFIqLi5kxYwYWi4X/+Z//4e2338ZqtZKTk8PcuXMpKChg6NChrF27NvZh0jC/detW7r//flwuFz6fj1dffZUPP/yQhQsXEgwGcTgc3HrrrYwYMSLJr5RIRRLkwlSuvPLKZkMrffv25fHHH2fSpEl8+OGHPPzwwwQCAUaPHs3FF1/MunXrKCsr49Zbb2X06NG89NJLzJo1i2XLlnHffffh8XhYsWIFwWCQ66+/nqeffpqJEyfy7LPPsnbtWmw2G08//TSff/4555xzTqu1bd26lVWrVnHcccfx3Xff8fvf/57nnnuOnJwctm7dytVXX83f/vY3XC5Xol8mkWYkyIWptDa0cs8993DRRRfhcDhYtmxZbPmwYcMYPXo0AJMnT2bOnDlUVlayevVqXnrpJRRFwWazUVpayrPPPsuvfvUrhg0bxqRJkxg7dixjx45lzJgxcWvr3bs3xx13HABr1qxh3759XHXVVbH1iqLw/fffM2zYsA68AkK0JEEuUsaBAweor68nEAiwb98++vXrB4CmaS221TQNXddRFCW2TNd1QqEQqqrywgsv8MUXX7B27VoeeOABiouLueWWW5rtIxAINJtv2tPWdZ0xY8Ywf/782LLdu3dTUFDQGU0Vohk5/FCkhGAwyIwZM7jpppu48cYbufnmmwkGgwBs3ryZzZs3A7BkyRJGjBhBVlYWZ5xxBi+88AKGYRAIBHj55Zc57bTT2Lx5Mz//+c8pKiri2muv5aqrruKLL74AIDc3Nza9cuXKo9YzZswY1qxZw7Zt2wB4//33ufDCC/H7/Yl8GUSakh65MJXDx8gBZsyYwT/+8Q/y8vK47LLLAFi1ahW///3vOfPMM8nLy2P+/Pns3LmT3NxcHnroIQDuuOMO7rvvPiZOnEgwGKS4uJjrrrsOm83GhAkTmDx5Mi6XC4fDwR133BH7nTlz5pCVlcVpp51Gfn7+EescPHgwc+bMYcaMGRiGgcViYeHChWRkZCTw1RHpSpHL2IpUtm7dOu69995We89CmJ0MrQghhMlJj1wIIUxOeuRCCGFyEuRCCGFyXX7Uiq7rhMNtH83RNKVd26eKdGx3OrYZ0rPd6dhm6Fi7rdaW50M06PIgD4cNqqpq27y9x+Nq1/apIh3bnY5thvRsdzq2GTrW7vz8zKOuk6EVIYQwOQlyIYQwuTYF+Weffca0adNaLF+5ciWXXXYZpaWlzJ49G13XO71AIYQQrYs7Rr5o0SKWL1+O0+lsttzv9zN//nxWrFiB0+lkxowZvPvuu/zsZz9LWLFCiO4rHA5RWVlBKBSIu+3evQrpeApLW9ptsdjIyclH09r+FWbcLQsLC1mwYEGLK7/ZbDYWL14cC/hQKITdbm/zEwshUktlZQUOh4uMjF7Nrip5JJqmEg6n31/w8dptGAY+XzWVlRXk5fVu837jBvm5557Ljh07WixXVZW8vDwAnn/+eWprazn99NPjPqGmKXg8bb+wvqap7do+VaRju9OxzZA67d63L0RWliduiDfQtPT8ii5eu7OyPNTWVrfr30SHDj/UdZ2HH36Yb7/9lgULFrTpDZTDD9smHdudjm2G1Gm3ruvougHEHzKRHnnrdF1v8W8iYYcfzp49m/r6ep544okWY+idzR8M85cv96bluJoQQrSm3T3yFStWUFtby8knn8zSpUsZPXo0V155JQBXXHEF48eP7/QiAT7+voq739zCsJ5uivLkms5CiJY+/XQ9r7/+KvfcMze2bOHCBfTvP4A9e3azdu2HaJqFX/96BieeeHISK+1cbQryvn378vLLLwMwceLE2PKGu650BasWGbbx1oe67DmFEKmhuvoQmzZ9ylNPPcvevXu5445b+OMfn0t2WZ3GNHcIckWvM1AbDCe5EiFEPH/5ci/Ly/Ycdb2iQHtHSS88uRcXnNTzmOrRNAs/+tFPURSFXr16RQ+VrCQnJ+eY9tfdmCfIbZEgrwtIkAshjm7DhvXceOM1sfldu3Zy4YWTyM72xJa5XBn4fF4J8q7mlB65EKZxwUk9W+09J/KolVGjRrcYI8/IcFNb64stq6314XYf/SgQszHNgZwNPfLaQPodsiSE6JhQKMTHH/8DXdfZs2cPum7g8XiSXVanMU2PPDZGHpAvO4UQ7ZOdnc0ppwzn2muvxjAMZsy4NdkldSrTBLndoqIqUCdDK0KIoxg5cjQjR45utuz66/8zNj19+rVdXVKXMM3QiqIoOK0atUEZWhFCiKZME+QQGSeXo1aEEKI5UwW506rhkyAXQohmTBXkGTZNxsiFEOIwpgryyBi5BLkQQjRlqiCXMXIhhGjJVEEuPXIhRGs+/XQ9d911e7NlCxcu4K9/XQHAjh0/MG3a5ckoLaFMcxw5RHrktdIjF0Icgzff/AuvvLKYQ4cOJbuUTmeuILfKl51CmIF981IcXy0+6npFaf/Nl/0nlFI/7NJjrikzM4vHHnuKkpKLj3kf3ZWpgtwZ7ZEbhtHm+wIKIdLLka5++KtfXcfppxcnsarEMlWQu6waBlAf0nFEr70ihOh+6odd2mrvuauvfpjqTPVlZ8MVEOWkICGEaGSuII/2wmWcXAghGplqaMUZuya5BLkQoqV4Vz8EWL78ra4sqUuYrEceKVd65EII0chcQW6L/AEhY+RCCNHIXEEuY+RCCNGCqYLcaYuUK2PkQgjRyFRBLj1yIYRoyVxBHh0jlx65EEI0alOQf/bZZ0ybNq3F8nfeeYfJkydTUlLCyy+/3OnFHc6mKWgKcgVEIYRoIu5x5IsWLWL58uU4nc5my4PBIHPnzmXp0qU4nU6mTJnC2WefTX5+fsKKVRQldr0VIYQ43O7du7jyyikMGTI0tmzUqB9x9dX/nsSqEi9ukBcWFrJgwQJuueWWZsu3bdtGYWEh2dnZAIwaNYr169czYcKEVvenaQoej6vNBWqa2mz7DLuFsNK+fZjR4e1OB+nYZkiddu/dq6BpjX/kX3iho1P3v3y5P+42mqYycOBAFi78Y6c+d2dq+hodjdLOjIsb5Oeeey47duxosdzr9ZKZmRmbz8jIwOv1xn3CcNigqqq2zQV6PK5m2zs0lSpvoF37MKPD250O0rHNkDrtNgyj2YWwWrtM7bFcxrYtF9kKh3UMo23b3n//3Rw6dIjq6kM89NB8nn32T3z++SYAxo8/j3PPncBNN/0HzzzzImVlnzNr1m9YufJtDhzYz7x59/Loo4+1q35o+8XCDKNlTubnZx5l6w6cou92u/H5fLF5n8/XLNgTxWXTqA2GEv48QoiOee21uqOuS+TVD7/77ttml7G96677yM8vOOK2o0aNpqTkF6xZ8wG7d+/iqaeeIRwOc/310xk16kdkZ2ezd+8e1q1bS0FBT7Zs+YrNm79i7NizElL7sTrmIC8qKqK8vJyqqipcLhfr169n+vTpnVnbEcl9O4UQrRkwYCCPPfZUm7YtLOwPQHn5t5x66nAURcFisXDSSf/Cd99tZ+zYs1i7dg1ffPEZv/jFlXzyyTrKyj7nttvuTGQT2q3dhx+uWLGCJUuWYLVaue2225g+fTqlpaVMnjyZnj17JqLGZiL37UzMJ7kQIr0oSiQC+/cfGBtWCYVClJV9Tt++hRQXn8WqVW+RkZHBT396Gh988D6BQJAePfKSWHVLbeqR9+3bN3Z44cSJE2PLx40bx7hx4xJT2VHI7d6EEJ3t9NOL2bhxA9deezXBYJBx485h6NBhANTX1zNq1I/JyspC0zROO+30JFfbkmK09xuHDgoGwx36svP+v33NB9sP8uZ1P01Eed1GqnwB1h7p2GZInXbv2VNOr17927RtIsfIu7O2tvtIr2VCvuxMFhkjF0K0VTAY5Oabb2ixvLCwP7fc8v+SUFFimC7IndGhFd0wUOUGzEKIVlit1jZ/8WlmprrWCjS/AbMQQggzBrnc7k0IIZqRIBdCCJMzXZA7o9cklysgCiFEhOm+7IzdXEJ65EKIw5jh6of/+MdH7N27h4suuoTXX1/GBRdciMXSsSg2X5DbpEcuhBlkX3z+UdcpCrT3DJZDr/21Tdu15xT9ZPjpT0+LTT///J8577wL0i/InTJGLoToBIm6+uGXX5bx3//9CIZhkJ9fwF133cvMmb/G48nB663hnHP+lR9++IF+/fpx8OAB7r77d8yd+18daovpgtwlY+RCmEJrPehUvvrhQw/dzz33PMCAAQNZtuwVvvvuOyDy4TBu3M9YseJ1AH7+84t55pk/cffdD3S4zaYNchkjF0IcSbKvflhZeZABAwYCcMkll7V4rkQw3VErMkYuhOgsibj6YV5eHj/88D0AL7zwDO+//y4AqtoybhVFbfcNNo7EdEFu1RQ0VZExciFEpzn99GJ69z6Oa6+9mmuuuYqzzhrH0KHDKCjo2e6rH86a9Tvmzp3DjTdew9atWxgz5ujbn3rqcH772193OMxNd/VDgHGPfcT5Jxbw23GDO7u8biNVrojXHunYZkiddsvVD+OTqx824bSq0iMXQsSVqKsf7tmzh/vum91i+YgRo5g+/dpj3u+xMmWQZ9gscnMJIbohwzBQutFVSRN19cNevXol7Fj1YxkkMd0YOUSOJZcvO4XoXiwWGz5fdad8eZeuDMPA56vGYrG16/dM2SN3ydCKEN1OTk4+lZUVeL1VcbdVFCUtA78t7bZYbOTk5Ldrv6YMcqdVY29NfbLLEEI0oWkW8vJ6t2nbVPmCt70S1W5TDq24bHIDZiGEaGDaIK8Npt+hS0IIcSSmDHKnVaM2EEp2GUII0S2YMshdVo26oI6ehl+WCCHE4cwZ5NHrrfhleEUIIeIHua7rzJ49m5KSEqZNm0Z5eXmz9cuXL2fSpElMnjyZF198MWGFNiUXzhJCiEZxDz9ctWoVgUCAJUuWsGnTJubNm8fChQtj6x966CFWrlyJy+Xiggsu4IILLiA7OzuhRcfu2xkIQ0ZCn0oIIbq9uEG+YcMGiouLARg+fDhlZWXN1g8dOpSamhosFkuXnZ4r1yQXQohGcYPc6/Xidrtj85qmEQqFYveYO/7445k8eTJOp5Px48eTlZXV6v40TcHjcbW5QE1TW2yfnxuZV+2Wdu3LTI7U7lSXjm2G9Gx3OrYZEtfuuEHudrvx+XyxeV3XYyG+efNm3nvvPf7+97/jcrmYNWsWb7zxBhMmTDjq/sJho8OXsVWiY+O79nup8jjavC8zSccz39KxzZCe7U7HNkPH2t3aZWzjftk5cuRIVq9eDcCmTZsYMmRIbF1mZiYOhwO73Y6maeTm5lJdXX1MRbZHvjtyQZn9vkDCn0sIIbq7uD3y8ePHs2bNGkpLSzEMgwceeIAVK1ZQW1tLSUkJJSUlTJ06FavVSmFhIZMmTUp40TkuG6oC+7wS5EIIETfIVVVlzpw5zZYVFRXFpqdMmcKUKVM6v7JWWFSFHhk2KuTCWUIIYc4TggDy3XYqZGhFCCFMHOQZNiq80iMXQgjzBrnbRoWMkQshhHmDvCDTTrU/hF9O0xdCpDnTBnlehhyCKIQQYOIgL3DbAdgn4+RCiDRn2iDPazgpSMbJhRBpzrRB3tgjlyAXQqQ30wa5267hsKhyCKIQIu2ZNsgVRZFDEIUQAhMHOUTP7pQeuRAizZk8yKVHLoQQJg/ySI/cMIxklyISLPvi88m++PxklyES7OKLnVx8sTPZZZiOyYPcRiBscMgfSnYpQgiRNKYO8oZDEOVYciFEOjN1kDfcKUjO7hRCpDOTB3mkRy5Hrggh0pmpg7zhwlly5IoQIp2ZOshtFhWP0ypBLoRIa6YOcoiMk8sYuRAinaVEkMtRK0KIdJYCQW6XHrkQIq2ZPsgL3DYqa4OEwnqySxFCiKQwfZDnue0YyC3fhBDpy/RBXuCWQxCFEOnN9EEuJwUJIdKdJd4Guq5z9913s2XLFmw2G/fddx/9+/ePrf/888+ZN28ehmGQn5/Pww8/jN1uT2jRTfXOijzXjip/lz2nEEJ0J3F75KtWrSIQCLBkyRJmzpzJvHnzYusMw+DOO+9k7ty5vPTSSxQXF7Nz586EFny4LIeVfLeN7Qd8Xfq8QgjRXcTtkW/YsIHi4mIAhg8fTllZWWzdt99+i8fj4dlnn+Xrr7/mzDPPZNCgQa3uT9MUPB5XmwvUNDXu9kN6ZvJdlb9d++3u2tLuVNNamzVLpM+Riq+JvNeNLCn8PkPi3uu4Qe71enG73U0K0QiFQlgsFiorK9m4cSN33nkn/fv357rrruPkk09mzJgxR91fOGxQVVXb5gI9Hlfc7ft7HKwvr+TAQR+aqrR5391ZW9qdalprc3YocnjpoRR8TeS9bhQKRW4qUVVV19UldYmOvNf5+ZlHXRd3aMXtduPzNQ5b6LqOxWKJFuWhf//+DB48GKvVSnFxcbMee1cZ1MNFfUhn1yEZJxdCpJ+4QT5y5EhWr14NwKZNmxgyZEhsXb9+/fD5fJSXlwOwfv16jj/++ASVenRFeRkAMk4uhEhLcYdWxo8fz5o1aygtLcUwDB544AFWrFhBbW0tJSUl3H///cycORPDMBgxYgRnnXVWF5Td3MAekTGn7QdqOXNwlz+9EEIkVdwgV1WVOXPmNFtWVFQUmx4zZgxLly7t/MraIcNmoXeWnW37pUcuhEg/pj8hqMGgHhlsP5BeXxgJIQSkUJAX5bn47mAtId1IdilCCNGlUibIB/XIIBg22FGZmoctCSHE0aRMkBflRb7w3CZHrggh0kzKBPmAXBcKsH2/jJMLIdJLygS5w6rR1+OQHrkQIu2kTJBD9MgV6ZELIdJMSgV5UZ6L7ytrCYTktm9CiPSRUkE+qEcGYQO+lyNXhBBpJKWCvOGaK3KGpxAinaRUkBfmOLGoClv2eZNdihBCdJmUCnKbReWkXpls3Hko2aUIIUSXSakgBxjZL5uv9nqpDYSTXYoQQnSJlAvyEX2zCesGX+yqTnYpQgjRJVIuyE/pk4WmwKcyvCKESBMpF+QZNgtDe2ay8YeqZJcihBBdIuWCHGBk32zK9tRQLycGCSHSQEoG+Yi+2QTDBmW7ZZxcCJH6UjLIhx+XhQJs3CHj5EKI1JeSQZ7lsDI4P0OCXAiRFlIyyCEyTv75rmpCYRknF0KktpQOcn9I56u9crq+ECK1pWyQj+ibDcCnMrwihEhxKRvkOS4bx+dnsGb7gWSXIoQQCZWyQQ5w9uA8Nu2sZr8vkOxShBAiYeIGua7rzJ49m5KSEqZNm0Z5efkRt7vzzjt55JFHOr3Ajhg3JA8DeG/r/mSXIoQQCRM3yFetWkUgEGDJkiXMnDmTefPmtdhm8eLFfP311wkpsCMG9XDRP8fJOxLkQogUFjfIN2zYQHFxMQDDhw+nrKys2fqNGzfy2WefUVJSkpgKO0BRFMYNyePTH6qoqg0muxwhhEgIS7wNvF4vbrc7Nq9pGqFQCIvFwr59+3jsscd47LHHeOONN9r0hJqm4PG42lygpqnt2v5wF43sy5/X/cAnu2u4bFTfY95PV+tou82otTZrlkifIxVfE3mvG1lS+H2GxL3XcYPc7Xbj8zXeA1PXdSyWyK+9+eabVFZWcs0111BRUYHf72fQoEFccsklR91fOGxQVVXb5gI9Hle7tj9cH6eFPtkOVn62k/FFuce8n67W0XabUWttzo5eAO1QCr4m8l43CoWcAFRVpeYN1DvyXufnZx51XdwgHzlyJO+++y7nn38+mzZtYsiQIbF1V1xxBVdccQUAy5YtY/v27a2GeDIoisK44/NY/OlOavwhMh1xmyyEEKYSd4x8/Pjx2Gw2SktLmTt3LrfffjsrVqxgyZIlXVFfpxh3fB4h3eADOaZcCJGC4nZPVVVlzpw5zZYVFRW12K679cSbOql3JgVuG3/bXMH5J/ZMdjlCCNGpUvqEoAaqovDzk3ry0bcH2XXIn+xyhBCiU6VFkANMOqU3igLLPt+d7FKEEKJTpU2Q98pyMLaoB69/sYeA3AJOCJFC0ibIAS49tQ9VdUH+vrUi2aUIIUSnSasg/1F/D4U5TpZukuEVIUTqSKsgVxWFyaf25vNd1WzZJzecEEKkhrQKcoCfn9QTu0Vl6aZdyS5FCCE6RdoFeZbDyoQTCvjLP/eyt6Y+2eUIIUSHpV2QA1z9k0IMA/687vtklyKEEB2WlkHeJ9vBRf/Si9e/2MPuajlBSAhhbmkZ5BDplasK/Gmt9MqFEOaWtkHeM9POpFN6s/LLPexI0UtmCiHSQ9oGOcBVP+6HRVP549oj34dUCCHMIK2DPM9t5/LhffjrP/dRtrs62eUIIcQxSesgB5g+ppB8t425b28lpBvJLkcIIdot7YM8w2ZhxtlFfF3h4xU5SUgIYUJpH+QQuYPQmAE5/GHNd1R45SQhIYS5SJATua/nLT8bTEg3ePTdbckuRwgh2kWCPKqvx8kvf1LIqq/389ZX+5JdjhBCtJkEeRNX/Lgfp/TJYu6qrXJsuRDCNCTIm7CoCveePwxFgTv/uplQWO4kJITo/iTID9Mn28Hvxg+hbHcNf/hIThQSQnR/EuRHMH5oPhed3ItnP/6B97/Zn+xyhBCiVRLkR/HbcUUM6+nmjr9sZsteuZuQEKL7kiA/CodV49GLTyLLYWHGa2VyfLkQotuSIG9FntvOo5NOpqY+xMzXvqQuGE52SUII0ULcINd1ndmzZ1NSUsK0adMoL2/+BeDKlSu57LLLKC0tZfbs2eh6ah3pMbTAzX0XnMCWfV5mvPYlfglzIUQ3EzfIV61aRSAQYMmSJcycOZN58+bF1vn9fubPn89zzz3H4sWL8Xq9vPvuuwktOBnGFvXgrvOGsuH7KmYt/yf1odT6sBJCmJsl3gYbNmyguLgYgOHDh1NWVhZbZ7PZWLx4MU6nE4BQKITdbm91f5qm4PG42lygpqnt2j5Rpp42EKvdyu9eK+PON7fwWOkIbJbEjUx1l3Z3pdbarEVf61R8TeS9bmRJ4fcZEvdexw1yr9eL2+1uUohGKBTCYrGgqip5eXkAPP/889TW1nL66ae3ur9w2KCqqrbNBXo8rnZtn0jji3KpOWcwc1d9w9V//piHLjqRDFvcl/CYdKd2d5XW2pwd/SvoUAq+JvJeNwqFIp3CqhQ9s7oj73V+fuZR18XtUrrdbnw+X2xe13UsFkuz+QcffJA1a9awYMECFEU5piLN4pJT+zD73CFs+KGK65Z8zgFfINklCSHSXNwgHzlyJKtXrwZg06ZNDBkypNn62bNnU19fzxNPPBEbYkl1E0/uxX9NOpnvDtYy/aVNlB9Mr96UEKJ7iRvk48ePx2azUVpayty5c7n99ttZsWIFS5Ys4csvv2Tp0qV8/fXXXHnllUybNo233367K+pOutMH5rLw8lPwBcJc9eJGPth2INklCSHSVNwBXlVVmTNnTrNlRUVFsenNmzd3flUmcXLvLJ77txHc8vo/mfHal1xzWn+m/7QQNcWHl4QQ3YucENRBvbMcLCo9lQtOLOCpj8q5aVkZ++UsUCFEF5Ig7wQOq8Zd5w3ltnMGs3HHIaY896lcbEsI0WUkyDuJoihMPrUPL/zbSHpl2vnt6//knje3UFUXTHZpQogUJ0HeyQb0cPH01OFc/ZN+vPHPvVz+5/W8+dU+DMNIdmlCiBQlQZ4AVk3lP84YyPPTRtIn28Gdf93MDUu/4JsKX/xfFkKIdpIgT6Dj8938acpwZo0bzJZ9Xn7x/Abmvr2Vg7VyEpEQovMk5vxyEaOpCpeP6MO5w/JZtLacpZ/t5s2v9lE6sg+/GN2XLIc12SUKIUxOeuRdJNtp5bfjBrP4ylGcNjCXp9f9wIWLPuapj76TL0SFEB0iQd7FBuS6mDvxBF68YiQ/KvSwaO33THxqHf/17jZ2V/uTXZ4QwoRkaCVJjs938/BFJ/FNhY8X1v/AK5t28crGnZw5OI/LhvfhZ9npcd0aIUTHSZAn2eD8DO6eMIzrTh/AK5t28foXe3hn634Gv7ediScVMOGEAnJctmSXKYToxiTIu4leWQ7+c+wg/n1Mf97eUsFrX+7l9+9tZ8Hqbyku6sH5JxRw2sDchN7MQghhThLk3YzDqjHx5F5MO2MQ67+pYEXZHt745z7e3boft11j3PF5/GxIPj8q9GDVJNSFEBLk3drgvAxuPquI/xw7iE++r+Str/axast+lpftxW3XOGNQD84s6sFPB+TgtstbKUS6kv/7TcCiKowZkMuYAbncPl7n4/JK3t26n9XbDvDmV/vQVIXhx2Vx2oBcftI/h+MLMuRSukKkEQlyk7FbVIqLelBc1IOQblC2q5oPvz3Ih9sPsOCDb1nwwbfkOK2MLvQwql82o/p66J/rTO4t+AwDwvUo4XoI1aOEA5HpcGSacCCybL+CrdobWaYHUMJB0CM/mncXGAauT+aDHkTRQxD9UYwQ6GEwwih6OLLcCIMeRjHCYOhN5vXovA7ooOsoGE2WRacb6saIPhKZbpPoa60okenYoxqZVlQM1MhmiopmtZIdMiLLFTW6XeOPoWigatF5LTqvgmIBVcVQLBDbRsNQNVAtoFiaTGsYqgVUS+wxNq1EpzVrdDsraFZQrdFto4+aNbJOtTbbBuk0JJ0EuYlZVIXhfbMZ3jebG4sHUuGt5+PyKtaVV7L+hyre3lIBQK7Lyil9smI/QwvcOKxa447CAZSgDyXgizwGfSjBWpRQbeN0bL4OQrUooTqUYF3kMeRHCdVByB+ZDvubTEcDvI2yj7biUOQ6NRkfPwIQDZSGINKi01rjdJNgi4Sf2iQMGwJRiQQhapOQjYSu0RC80BjEkZk4LWgS+tEPAiX2AUHzD4uGn3AAJRQEQ2/2QdM43fRDKByd15t8gOmRDzaj4cMs1ObXuzM0hr21Mew1G4ZmiyyPPVpBtWFoVjSHk8ywGtlOtUV+T7OjVV+HgYpzw0uNv6vZMDR79LFhH/Ym62xw2HzD9qjpEXHp0cpUpYdQAjUo9YdQ66vpE6hhkuMQlwysQelTTU11JfsP7qemuhL/rkNo5T7cSh0+6rBp9bgVPw6jFs1o+//4hmrBsLgwLE4MqxMsDgzNgWF1YrjyGuctkR80e2Ras0fXNf2fsvn/gO7sLGpq9cg20Z5fQzBkv1WKoagcun5lNJxTpxfYkTurH5WhN/mLJdw4HXsMRj8EosvCwSbLG7Y7bFk40DgfbvirqPlfTpFtgpFl4UD0d6PzeiDy7zUcQPGFsATqm2wfmVa9F4Jh4P7HvM55GRQ19kHQ+O/OfoRltubLLPYm/z6bTFvskQ+S2HpH4/oWv9OwH3vC/3KRIE82Q0epr0apr0L1V6L4q1Drq1CVWlxVFZHl9YdQ/IeigV2FEqhGqa9GDbZ+NcVMoLfFhW5zY2RlErRkUK3ncDDUk80BK7vqLFSF7XgNJ3WKE5c7i+wsD3k5OeTl5NC7Ry45Wdlgy8CwRsIbLXHHtBseF6GjBJrR0LNKkx5Wh0UDDM3WbECou1xM+WgfXsE3IifCVVz7TbOQbxh+IxwNfT3QZNlhQ3TNtm86X99sXWw6VIdSXxVZFvJHPnBCTbbTO+cid4bFRXjSIig4s1P215T8X9GZDD0Stv5KlLoDqHUHUf0HI9P+ymhQH/ZYXxX907ulDEC3ZmDYszHs2ej2bMLZA9DtHgx7FoY9C92ejWHLwrBlRuZtWRj2zMi8LbNF8KlAXvRnmGGws8rP5n1eNu/1sqHCyzf7fVR83/AP14fb7mdAriv2U5jjpDDHSV+PE7sc0y4SpeEvOrrBh0/D8NcRv+OJfoiE/E0+VOobp0MN034UPYy14MSElChB3pqGYK47gFq3PxrOB1Br96PWHYgGdDSw6w6g+Csjf8YeaVeaHd2Zi2HPQXfkEO5xAoYjMt346Ik82j1kFvSiyh8ZO0wUVVHol+OkX46T8UPzY8ur6oJsP+Bj2/5atu/38V1lHevKK1n55d7YNgpQkGmnr8fBcdkOjst20ifbQe8sO72zHOS5bXLkjEgNitr4wWLv2AeLx+OCzh5GI92C3DAgWBsJ32gYNwvouuiy2orYsqMFs27PRnfmYTh7EPYMJNhrNLqzB4YzF92RGwlkZw90Zw90Ry5Y23ntlAwXBDv/DW8Lj9PKyL4eRvb1NFvurQ/xQ1Ud3x+s4/vKOnYcqmNnlZ8131ZywLe32bYWVaEg006vTDsFmXYK3HYK3DbyM+3kZ9jId9vokWGTk5qE6ASmCXLFX4lrw2MYNjfhrEL0zOMwUKPf3IciR1fU16AGaiLDFf7ImLNaX4lSVxkJaP+ByBjYEejWDAxnHrozl3BmP/SewzEcPdBdeejOvEggO3tEwtuRk9Cecnfltls4oWcmJ/TMbLHOHwyzu7qeXYf87Knxs6e6nj019eyt9vP5rmr21dQT0lv2ZbIdFvLcNnJdNnp5nGRaVXKcVnJdNnJcVjzOyI9bN9BU6eELcSTmCfL6auzb30SrLm/T9rotKzpUkYvuyiOcOyQ63SPWk44EdCSk291jFs04rBoDe7gY2MN1xPW6YVBZG2S/L0CFt54Kb4D9vgAHYj9BNv1QRUVNPf5Qy+8MFu88BMD1T6wl22kh22Ely2Ehy2kly24h02Eh0x75cdstZDo03LbItNuukWGzyAeBSFmmCXI9uz8Hp62BkB+tZgeqd3dkhaphKBYMawaGzR370k+ObuheVEWhR0ZkOGVogfuI2zQcyeAPhqmsC3KwNkhVXZBDdUH6/tVJSDc4+/g8DvmDHPKH2FtTz9cVPmr8IWqDRx4Ca8ppVcmwWXDZNDJsWvTRgtOq4rJpuKyN0w6rhtOq4rQ2TjssWnRexW6JzNstqnxAiKSLm3a6rnP33XezZcsWbDYb9913H/3794+tf+edd3j88cexWCxMnjyZyy+/PKEFY3EQzhlMOGdwYp9HJI3DqtHbqtE7yxFblp1pB+D28ccf8XdCuoHXH6K6PoS3PkRN9NFXH8YbiDz6Ao3TtcEQtYEwu6v9+IORdbWB8BH/GojHoirYLQ3hrmK3aNgsKjZNxW5VsWlKZNqiYm3yaNMUrJpKVoadcDCETVOxagoWTW2cViOPVk3BqqpYNAWLGtnGoirNttEUJbo+sk4+YNJH3CBftWoVgUCAJUuWsGnTJubNm8fChQsBCAaDzJ07l6VLl+J0OpkyZQpnn302+fn5cfYqROeyqAoelxWPq2PfXeiGQX1IpzYQpi4Yxh/UI4+hptM69SEdf5PphvlAuHE+ENIJhHV89ZHHYNjAHwwTDBsEwo3LEkUBLNGA11QlFu4Nj5rauK7ZckVBjU1H/ppqtjy67EjbKETuU6sqSnQZKEqT31EUVBVcThv1/mDkdxQFNbp+n7cQBXj1s12R5USWK0rzR1Uh9nsN2ylNl6HETtRVm05Hj6Rqupwmz6EokTY07CP6X2RddCa2/+iLHFmvNG7X9Pkb3ovo6+PxHHnosaPiBvmGDRsoLi4GYPjw4ZSVlcXWbdu2jcLCQrKzIydWjxo1ivXr1zNhwoSEFCtEoqmKgtMaGULpCoZhkJHppOKgl2DYIBgN+FDYIKhH5iPTOsGQQUg3COmRD4CgHlkX0g2CYYOwYRAK69Ftoj9hA90wCEfnw3rDtE7YIPr7OrpBbF04un0wpFNvGIQN0Jss1w2j2fZ6k20a1jU8pwHNfieePZWRTuC8Vd8k9oVPkgcv+RfGDczp9P3GDXKv14vb3TimqWkaoVAIi8WC1+slM7PxCIaMjAy8Xm+r+9O09n0qaZqasE+x7iwd291qm997DwBPl1XTdTRN5biCrGSX0SUaPhAUVSEU0mNh37DcuC3y4RY2zkKPXq4mMm80Tjf8TvTDIzJtYES3afggMYhsrxvNHyPLG7ehYd1h2xtNljX8Dk23oXEbmszrR9jeINJJ+NeTeuG2dX4nIW6Qu91ufL7GU8F1XcdisRxxnc/naxbsRxIOG+26rkRCrkNhAunY7nRsM6Rnuz0eF15f84upNb00WYuoi61UQDPv2L/bph3ze52ff/RsjXs2xsiRI1m9ejUAmzZtYsiQIbF1RUVFlJeXU1VVRSAQYP369YwYMeKYihRCCHFs4vbIx48fz5o1aygtLcUwDB544AFWrFhBbW0tJSUl3HbbbUyfPh3DMJg8eTI9e/bsirqFEEJEKYZhtOEriM4TDIZlaKUN0rHd6dhmSM92p2OboWPt7tDQihBCiO5NglwIIUxOglwIIUxOglwIIUxOglwIIUyuy49aEUII0bmkRy6EECYnQS6EECYnQS6EECYnQS6EECYnQS6EECYnQS6EECYnQS6EECbXbYNc13Vmz55NSUkJ06ZNo7y8PNklJVwwGGTWrFlMnTqVSy+9lL///e/JLqlLHThwgDPPPJNt27Ylu5Qu8Yc//IGSkhIuueQSXnnllWSX0yWCwSAzZ86ktLSUqVOnpvx7/dlnnzFt2jQAysvLmTJlClOnTuWuu+5C19t/o++j6bZB3vSmzzNnzmTevHnJLinhli9fjsfj4cUXX2TRokXce++9yS6pywSDQWbPno3D4Uh2KV1i3bp1bNy4kZdeeonnn3+ePXv2JLukLvH+++8TCoVYvHgxN9xwA/Pnz092SQmzaNEi7rjjDurrI3dCmjt3Lr/5zW948cUXMQyjUztq3TbIW7vpc6o677zzuOmmm2LzmtY1NwDuDh588EFKS0spKChIdild4sMPP2TIkCHccMMNXHfddZx11lnJLqlLDBw4kHA4jK7reL3e2G0jU1FhYSELFiyIzX/55Zf8+Mc/BmDs2LF89NFHnfZc3fZVbO2mz6kqIyMDiLT917/+Nb/5zW+SW1AXWbZsGbm5uRQXF/PUU08lu5wuUVlZya5du3jyySfZsWMH119/PW+++SaKYt77UbaFy+Vi586dTJgwgcrKSp588slkl5Qw5557Ljt27IjNG4YRe38zMjKoqanptOfqtj3y1m76nMp2797NFVdcwUUXXcTEiROTXU6XePXVV/noo4+YNm0aX331FbfeeisVFRXJLiuhPB4PZ5xxBjabjUGDBmG32zl48GCyy0q4Z555hjPOOIO33nqL119/ndtuuy029JDqVLUxbn0+H1lZWZ23707bUydr7abPqWr//v388pe/ZNasWVx66aXJLqfL/O///i8vvPACzz//PCeccAIPPvgg+fn5yS4roUaNGsUHH3yAYRjs3buXuro6PB5PsstKuKysLDIzI7csy87OJhQKEQ6Hk1xV1zjxxBNZt24dAKtXr2b06NGdtu9u28U90k2fU92TTz5JdXU1TzzxBE888QQQ+cIkXb4ATCdnn302n3zyCZdeeimGYTB79uy0+E7kqquu4ne/+x1Tp04lGAxy880343K5kl1Wl7j11lu58847efTRRxk0aBDnnntup+1bLmMrhBAm122HVoQQQrSNBLkQQpicBLkQQpicBLkQQpicBLkQQpicBLkQQpicBLkQQpjc/wcU92B4j0OP2AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEFCAYAAAD+A2xwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxBklEQVR4nO3de3wU9b3/8dfM7P2SbO5cE0ggKCog4BWpoqXealFAAxyxPXpOq8dzqj3Uqj1KPVYF6/m19mjRamtb9agoYgGrVpEqivUCChTkGiCRWxJCAtndJHuZ+f0xYUkkkJBsstndz/PxyGN3Z2dmP5PAO9985zvfUQzDMBBCCJG01EQXIIQQonskyIUQIslJkAshRJKTIBdCiCQnQS6EEElOglwIIZKcJdEFCNFTHnjgAT777DMAysvLGThwIA6HA4CFCxfGnguR7BQZRy7SwcUXX8yvf/1rzjjjjESXIkTcSYtcpJVgMMiFF17Iyy+/zNChQwH43ve+x/XXX8/y5cux2+1s3ryZ2tpaJkyYwD333IPVaqW8vJwHH3yQ+vp6otEos2fPZvr06Qk+GiFM0kcu0orL5eLqq6/mlVdeAaCyspJdu3YxadIkANavX88zzzzDG2+8QXl5OQsXLiQSifDDH/6QOXPmsHjxYp5//nmeeeYZ1q5dm8AjEeIoCXKRdmbNmsWSJUsIh8MsXLiQ6dOno2kaANdccw1utxubzcaUKVP48MMP2bVrF5WVlfz0pz9lypQpXH/99TQ1NfHll18m+EiEMEnXikg7Q4cOZcSIEbz77ru8/vrrvPzyy7H3jgQ6gGEYqKpKNBrF6/WyZMmS2HsHDhzA6/X2at1CHI+0yEVamjVrFr/4xS8YNWoUBQUFseVvvvkmoVCI5uZmXnvtNSZNmsTQoUNxOByxIN+3bx/f/va32bBhQ6LKF6INCXKRliZNmkQwGGTGjBltljscDmbNmsVVV13F+PHjmTZtGjabjQULFrBo0SKuuuoqbrzxRm677TbGjRuXoOqFaEu6VkRaWLFiRZvX69evJycnh/PPP7/N8vPOO4+bbrrpmO1POeUUnnvuuR6tUYiukiAXaefOO+/k008/5Ve/+hWKoiS6HCG6TS4IEkKIJCd95EIIkeQkyIUQIsn1eh+5rutEo13rzdE0pcvbJis55vQgx5weunPMVqt23Pd6PcijUYP6+mCXtvX5XF3eNlnJMacHOeb00J1jzss7/gVo0rUihBBJToJcCCGSnAS5EEIkObkgSAgRd9FohLq6GiKR0HHXqapSSLfLWDpzzBaLjaysPDSt8/EsQS6EiLu6uhocDhdud7/jXj2raSrRqN7LlSVWR8dsGAaBwGHq6mrIze3f6f1K14oQIu4ikRBud4ZMgXCSFEXB7c444V8y7elUkK9bt47Zs2cfs3zFihVMmzaNsrKyNnM6CyGEhHjXdOX71mHXytNPP83SpUtxOp1tlofDYebNm8eiRYtwOp3MnDmTSZMmkZeXd9JFdEZjOMr7a/fyjaJM+QcihBCtdBjkhYWFPPbYY/zkJz9ps7y8vJzCwkIyMzMBGDduHKtXr+byyy8/4f40TcHnc510of/YfoAfv7qeRT84l9GDfCe9fbLSNLVL369kJsec/KqqFDSt4z/4O7NOV33++Wpee20RP//5/NiyBQv+l6KiIezfv49Vqz5E0zRuv/3HnHba6T1Wx9d15pgV5eRyssMgv/TSS9m9e/cxy/1+f5tbXbndbvx+f4cf2NUrO70tjfANFXUUeWwnvX2ykqvf0kOqHbNhGB2eyOzpk53RqH5MHbpuUF9fz+efr+Gpp/5IVVUV99zzE373u2d7rI7WOnvMhnFsTp7oys4uj1rxeDwEAoHY60Ag0KP3MBzoc6CpCpV1qfOPXYh08JeNVSzdsP+Y5YoCXR19+J3T+3HlaQUdr9gOTbNw1lnnoigK/fr1axkqWUdWVlbXiukDuhzkJSUlVFRUUF9fj8vlYvXq1e3eWSVerJrKIJ+TirrGHvsMIURqWbNmNf/+79+Pvd67dw/f+c41ZGb6YstcLjeBgD+9gnzZsmUEg0HKysq46667uOmmmzAMg2nTprW5iW1PGJrrlha5EEnmytMK2m0998Y48nHjxvPf/z0v9vqJJx7D7fYQDB7tTQgGA3g8Pdeb0Bs6FeSDBg2KDS+86qqrYssvvvhiLr744p6prB3FuW7+vqMW3TBQZeSKEKILIpEIn376MTNnzqa6uhpdN/D5fIkuq1uS6srOIbkumiM61Q3N9MtwJLocIUQSyszMZNSoMfzgB/+MYRj853/emeiSui2pgrw41w1AxcFGCXIhxAmNHTuesWPHt1l2yy3/EXt+000/6O2SekxSXaI/JKclyKWfXAghYpIqyPO9dlxWjUoZuSKEEDFJFeSKolCY5aTioAS5EEIckVRBDlCU7ZQhiEII0UrSBXlhlpN9h5tpCkcTXYoQQvQJSRfkRVkuDGD3oaZElyKEEH1C8gV5tjmdbuVB6V4RQhzf55+v5mc/u7vNsieeeIw33lgGwO7dXzF79nWJKC3ukmocOcDgLDPIZc4VIURXvfXWX3jllZc4dOhQokuJi6QLcrfNQp7HJkEuRJKwb16EY9NLxyxXlK7ffLnp1Bk0nzK9yzV5vRk8/vhTlJVd3eV99CVJF+RgnvCUrhUhREfam/3wX/7lZiZMmJjAquIvKYO8KMvFu1trEl2GEKITmk+Z3m7rOVGzH6aipDvZCWaL/FBThPpgONGlCCFEwiVlkB8ZubJLuleEECI5u1aG5pg3Jd1xMMiYQZkJrkYI0Rd1NPshwNKlf+3NknpMUrbI+2c4cFhUdhwIdLyyEEKkuKQMclVRKM51s6NWulaEECIpgxygOMdFubTIhRAiuYP8YDBMfaOMXBFCpLekDfKSltu+7aiVVrkQIr0lbZAXt4xcKT8g/eRCiPSWtEFe4LXjtmkyckUIkfaSchw5mBPuFOfIyBUhRPv27dvLd787k9LSEbFl48adxT//8792a7/f+c6lfW78edIGOUBxrov3th3AMAwURUl0OUKI47j6aucxy7oz++Gf/9y52U+HDBnK448/1aXPSCbJHeQ5Lpb8I8LBYJgcty3R5QghktRXX1Xy8MMPEA6HcTgc3HffQzQ1NTJ//s+JRCIoisJtt/2Y4cNLE11qu5I6yFuPXJEgF6Lvaq8F3RuzH+7atbPNNLY/+9kD5OXlH7Peb37zKNdf/z3OPfd83n33HbZt28LSpYuZPr2MiRMvYtu2Lcyf/3N+//vnerTerkruIG81cuWswqwEVyOE6Gs627VSWVnB6aePAuCSSyYD8L//+0tGjx4LwPDhI6iuruq5QrspaUetAOS4bWQ6LDKWXAjRLUVFQ9m0aSMAb7/9JosWvcSQIUNYv/4LALZt20J2dk4iSzyhpG6RmyNXXOyQseRCiG649dbbeOSRh/jTn36Pw+Fg7tyfM2HCN3j44Qd48cXniUQi3H33vYku87gUo6unjbsoHI5SX9+14PX5XMdsO3/5Nv66uZoVt56fkiNX2jvmVCfHnPz276+gX7+iE67TG33kfU1nj7m9719enve463fYItd1nfvuu48tW7Zgs9l44IEHKCo6+gFLly7lD3/4A6qqMm3aNGbNmtVhkfFUnOPG3xylxh8i32vv1c8WQiSPcDjMj3506zHLCwuL+MlP/isBFcVPh0G+fPlyQqEQCxcuZO3atcyfP58nnngi9v4vfvELXn/9dVwuF1deeSVXXnklmZm9d7OHktyWE561AQlyIcRxWa3WlB1T3uHJzjVr1jBxonnH6TFjxrBhw4Y2748YMYKGhgZCoVBCLswZ1jIEcVu1nPAUQqSnDlvkfr8fj8cTe61pGpFIBIvF3HT48OFMmzYNp9PJ5MmTycjIOOH+NE3B53N1qVhNU4/Z1ueDAZkOdtQ3dnm/fVl7x5zq5JiTX1WVgqZ1PCiuM+ukms4cs6KcXE52GOQej4dA4GhrV9f1WIhv3ryZ9957j3fffReXy8Udd9zBm2++yeWXX37c/UWjRlxPdgIMz3Wzcc+hlDpZdESqnQTrDDnm5GcYRocn9eRk5/EZxrE5eaKTnR3+ahg7diwrV64EYO3atZSWHr1E1ev14nA4sNvtaJpGdnY2hw8f7rDIeCvNd1NxsJHGcLTXP1sIIRKtwxb55MmTWbVqFTNmzMAwDB566CGWLVtGMBikrKyMsrIyZs2ahdVqpbCwkGuuuaY36m5jRL4HA9hWE2DUgBN37Qgh0kNfnf3w17/+f5SV/RMul5OPP/473/rWZd2qBzoR5Kqqcv/997dZVlJSEns+c+ZMZs6c2e1CumNEvtmHv7XaL0EuRB+UefUVxyxTFOjqVSyH/vxGp9bri7Mf3nbbHAA+/3w1q1a93ztBngwKvHYyHBa2VPsTXYoQIgnFc/bDP/7xd3zwwftEo1Guvnoa55xzHnfe+SMyMjKZMOECVq36kDvu+CnPPvsM27dvY8mSxUyZMrVb9adEkCuKQmm+h601MgRRiL6ovRZ0Ks5+uHXrZj755COeeuqPhMNhnnzycc4++1wOHqzl979/HofDzqpVHwJwww03smTJq90OcUiRIAcYkedh0bq9RHQDi5p6l+oLIU5eb89+WFlZwamnnoamaWiaxu23/5h9+/bSv/8ArFZrN47kxFJmEGdpvpvmiM6ug6kzhEsI0TviNfthUdEQtm7dgq7rRCIRbr/93wiHQyjKsVGrqiq6Hp+prlKnRd7qhOeRqz2FEKIz4jX74fDhIzjnnPO45Zab0HWda66ZjtXa/k1vBg4cxI4d23n55Re47rruzVGV9LMfHhHRDS56bBXTRvfnRxeVtLtOMkq1C0U6Q445+cnsh+1L2OyHycKiKgzLdbNVRq4IIdoRz9kPlyxZzDvvvHXM8ptv/vdYX3tvSpkgB7Of/N2tBxIyeZcQoq2+9v8wnrMfTpkyNS6jTdrTlU6SlDnZCWY/+eGmCPsbmhNdihBpzWKxEQgc7lIopTPDMAgEDmOxnNzN5FOqRX7khOeWKj/9MxwJrkaI9JWVlUddXQ1+f/1x11EUJe2CvjPHbLHYyMrKO6n9plSQD8t1oymwqaqBi4bnJrocIdKWplnIze1/wnVS7QRvZ/TUMadU14rDqjEsz8OGfQ2JLkUIIXpNSgU5wOn9vWzc34CeZn+yCSHSV8oF+ch+XgKhKJUHGxNdihBC9IqUC/LT+5uD5jfs7/0bXAghRCKkXJAPyXbhtmnST97LMq++ot05p3vS1Vc7ufpqZ69+phB9UcoFuaoonNrPy5f7JciFEOkh5YIc4PR+XrbWBGiSe3gKIdJASgb5af28RHVD7hgkhEgLKRnkR054bpTuFSFEGkjJIM/12Cnw2tkoJzyFEGkgJYMczO6VDdIiF0KkgZQN8tP7e9l7qIm6YCjRpQghRI9K2SAf2U/6yYUQ6SFlg/zUAi+aAv/YK1d4CiFSW8oGucumUZrv4Ys9EuRCiNSWskEOcOagTDbuO0wokl43eBVCpJeUDvIxAzMJRQ02VUk/uRAidaV4kGcA8MXuQwmuRAghek5KB3mWy8aQbCdrpZ9cCJHCUjrIwewnX7f3EFFd7hgkhEhNKR/kYwZm4m+OUn4gkOhShBCiR1g6WkHXde677z62bNmCzWbjgQceoKioKPb++vXrmT9/PoZhkJeXxyOPPILdbu/Rok/GmYMyAbOfvDTfk+BqhBAi/jpskS9fvpxQKMTChQuZM2cO8+fPj71nGAb33nsv8+bN48UXX2TixIns2bOnRws+Wf0zHBR47azdIyc8hRCpqcMW+Zo1a5g4cSIAY8aMYcOGDbH3du7cic/n409/+hNbt27lwgsvpLi4+IT70zQFn8/VpWI1Te3StmcPzebjHQfJzHSiKEqXPjtRunrMvU2zmG2CeNTa2WO2xPEzEy1Zfs7xJMccPx0Gud/vx+M52iWhaRqRSASLxUJdXR1ffPEF9957L0VFRdx8882cfvrpnHfeecfdXzRqUF8f7FKxPp+rS9uenu9m2fp9bNh1kMFZyXWPx64ec2/LbLno6lAcau3sMUci5s+yvr6x25+ZaMnyc44nOeaTk5fnPe57HXateDweAoGjJwp1XcdisbQU5aOoqIhhw4ZhtVqZOHFimxZ7XzHmSD+5dK8IIVJQh0E+duxYVq5cCcDatWspLS2NvTd48GACgQAVFRUArF69muHDh/dQqV03NNuFz2llzVf1iS5FCCHirsOulcmTJ7Nq1SpmzJiBYRg89NBDLFu2jGAwSFlZGQ8++CBz5szBMAzOPPNMLrrool4o++QoisJZhT4+rajHMIyk6ycXQogT6TDIVVXl/vvvb7OspKQk9vy8885j0aJF8a8szs4p8vHOlhp2HgxSnONOdDlCCBE3KX9B0BFnF2UB8ElFfWILEUKIOEubIO+f4aAwy8mnFXWJLkUIIeIqbYIc4KxCH59/dYhIVOYnF0KkjrQK8nOKsgiGo/xjn8xPLoRIHWkV5OMH+1AVpHtFCJFS0irIvQ4LI/t55YSnECKlpFWQgzl65cv9h/E3RxJdihBCxEX6BXmhj6iBXOUphEgZaRfkowZk4LSqfLxL+smFEKkh7YLcqqmcVZjFqp0HMQy5/ZsQIvmlXZADXFCczb7DzZTXptcUmkKI1JS2QQ7wQXltgisRQojuS8sgz/PYObXAwwflBxNdihBCdFtaBjnAxOIcNuw7TF0wlOhShBCiW9I3yEuyMYBVO6VVLoRIbmkb5CPyPeR5bHy4Q4JcCJHc0jbIFUXhguJsPt5VR1hmQxRCJLG0DXKAC4pzCISifP6V3JRZCJG80jrIzy70YbeofLBDhiEKIZJXWge5w6pxTlEWf9t2AF2u8hRCJKm0DnKAb47Ipdof4h97Dye6FCGE6JK0D/KJxTnYNIXlWw8kuhQhhOiStA9yj93CeUOyWbG1RrpXhBBJKe2DHOCbI/Kke0UIkbQkyDGv8rRpCu9sqUl0KUIIcdIkyAG3zcL5Q7NZIaNXhBBJSIK8xTdL86jxh1i/R7pXhBDJRYK8xQUl2dgtqnSvCCGSjgR5C7fNwoSh2SzfWkNE5l4RQiQRCfJWrjytgIPBMB/JjZmFEElEgryV84dkke2y8vrGqkSXIoQQndZhkOu6zty5cykrK2P27NlUVFS0u969997L//zP/8S9wN5k0VQuOzWfD8prqW8MJ7ocIYTolA6DfPny5YRCIRYuXMicOXOYP3/+Meu89NJLbN26tUcK7G3fPq2AiG7w9ubqRJcihBCd0mGQr1mzhokTJwIwZswYNmzY0Ob9L774gnXr1lFWVtYzFfay4XkeRuR7pHtFCJE0LB2t4Pf78Xg8sdeaphGJRLBYLFRXV/P444/z+OOP8+abb3bqAzVNwedzdalYTVO7vO3JuHb8IB54YzPVzVFKC7w9/nkn0lvH3F2axWwTxKPWzh6zJY6fmWjJ8nOOJznm+OkwyD0eD4FAIPZa13UsFnOzt956i7q6Or7//e9TU1NDU1MTxcXFTJ069bj7i0YN6uuDXSrW53N1eduT8Y0iH5qq8OLHFdx2YXGPf96J9NYxd1dmxByyeSgOtXb2mCMRJwD19Y3d/sxES5afczzJMZ+cvLzjNyo7DPKxY8fyt7/9jSuuuIK1a9dSWloae++GG27ghhtuAGDx4sXs2LHjhCGeLLJcNiYWZ/OXjVXcMmEINosM7hFC9F0dJtTkyZOx2WzMmDGDefPmcffdd7Ns2TIWLlzYG/UlzNTR/alrDLNim8xTLoTo2zpskauqyv33399mWUlJyTHrpUJLvLVzirIY7HPwytq9XHZqfqLLEUKI45I+g+NQFYXpYwawfu9htlT7E12OEEIclwT5CXz7tALsFpVFa/cmuhQhhDguCfITyHBYuezUfN7aVE1DUyTR5QghRLskyDtw7egBNEV0lm3cn+hShBCiXRLkHRhR4OGM/hksWruXqC53DxJC9D0S5J0wc9xAvqpv4v3y2kSXIoQQx5Ag74SLh+cyyOfg2U+/wpB7egoh+hgJ8k7QVIXrxw9i4/4GPt99KNHlCCFEGxLknXTlyAKyXVb+9OlXiS5FCCHakCDvJIdVY8bYgfx9Vx1b5QIhIUQfIkF+EqaN7o/LqvHsZ9IqF0L0HRLkJyHDYeWaUf15Z0sNFQfTa/pNIUTfJUF+kmafNQibpvL039u/d6kQQvQ2CfKTlOO2MWPsQN7eXMP2mkDHGwghRA+TIO+C68cPwmXT+O1HuxJdihBCSJB3RabTyvXjB/He9lo27m9IdDlCiDQnQd5FM8cNJNNh4ckPdyW6FCFEmpMg7yK3zcL3zink44o6Pt51MNHlCCHSmAR5N1w3ZgCDfA5++d4OIjIzohAiQSTIu8FmUbn9wmJ21gZZvE7uIiSESAwJ8m76RkkOZxf6+O1HFdQ3hhNdjhAiDUmQd5OiKPxoUgn+5ghPfyQXCQkhep8EeRwMy3UzbfQAXl23ly0yoZYQopdJkMfJzROKyHRaefDtrXLiUwjRqyTI4yTDYeWOi4exqcrPS5/vSXQ5Qog0IkEeR5eU5vKNkhyeXLWL3fWNiS5HCJEmJMjjSFEUfnLJMCyqwrx3tsn9PYUQvUKCPM4KvHb+feJQPq2sZ/H6fYkuRwiRBiTIe8DU0f05tyiLX723g121cgMKIUTPkiDvAaqi8LPLSnFYVO55YzOhiJ7okoQQKUyCvIfkeuzce+kItlT7eXLVrkSXI4RIYRLkPejCYTlMG92f51bv5u8yQ6IQood0GOS6rjN37lzKysqYPXs2FRVtL0N//fXXufbaa5kxYwZz585F16UbobXbLyxmeJ6be/6ymT2HZEiiECL+Ogzy5cuXEwqFWLhwIXPmzGH+/Pmx95qamnj00Ud59tlneemll/D7/fztb3/r0YKTjcOq8YvvjMQw4CdLvqQpHE10SUKIFGPpaIU1a9YwceJEAMaMGcOGDRti79lsNl566SWcTicAkUgEu91+wv1pmoLP5+pSsZqmdnnbRPL5XPzqutH86/Nr+J/3d/LItDNQFKVT2ybLMWsWs00Qj1o7e8yWOH5moiXLzzme5Jjjp8Mg9/v9eDyeVoVoRCIRLBYLqqqSm5sLwHPPPUcwGGTChAkn3F80alBf37UheT6fq8vbJtrofDc/OL+IJ1dVUJRpZ/ZZgzu1XbIcc2bLyJxDcai1s8cciZgNiPoUuIo2WX7O8STHfHLy8rzHfa/DIPd4PAQCgdhrXdexWCxtXj/yyCPs3LmTxx57rNMtzXT0z+cUUn4gyP+u3EmB1863TslPdElCiBTQYR/52LFjWblyJQBr166ltLS0zftz586lubmZBQsWxLpYRPvM8eUjOHNgBve9tYXPd9cnuiQhRAroMMgnT56MzWZjxowZzJs3j7vvvptly5axcOFCNm7cyKJFi9i6dSvf/e53mT17Nu+8805v1J207BaVR6acxsBMBz/+85eUHwh0vJEQQpxAh10rqqpy//33t1lWUlISe7558+b4V5XiMp1Wfj31DG58cS23LvoHv71uFEXZ6XXSRwgRP3JBUIIMyHSw4Noz0HWDf3tlvUx7K4ToMgnyBCrOcfOba8+gOaLzb6+sZ9/hpkSXJIRIQhLkCTY8z8Pj08/A3xzlX15cK7MlCiFOmgR5H3BKgZcnrxtFRDf414Xr2FTVkOiShBBJRIK8jyjN9/C7GWNwWVVueXk9n1XWJbokIUSSkCDvQwZnOXl6xhgKvHb+49UNvCZ3GBJCdIIEeR+T77Xz+5ljOLvQx0PvbOOhNzcR1eXen0KI4+twHLnofR67hV9eczqPvlfOHz6qYNOew9x/xQiyXLae+UA9ihI6jNpUh9JUhxJqQG1uQAkdRgkHUSJBlHAAIk0okWaUaDPoYRQ9DNEIGFEstZvAgMyl/2TuUwEDFVQLKCqGagXNiqHaQLNhWOygOTAsDgyLE8PqwrA6MawelKxsLM0WDHsGhs2LbssAiwNk+gch2iVB3kdZVIUfXzyMMwqzuP/1L7n+uc956NunMnpgZud3YuiowRrUhj2o/r1o/n2ogf2ogSrUYBVqsBa18YAZ3py41W+oNjNoNbsZxKoVVCuoFgzVAnoEACXUAC37Ugzd/CVhRECPoETDoIdQoiGINKNEm1Batvu6rPY+356Bbs/EcPjQHVlodQ+CasW55mUMRxa6MwfdmYvhzDYfbV4Jf5EWJMj7uLLxgxnitXPX61/yg4XruOWCoVw/fhCa2hJQhoHq34dWvx2tfof5dWgX2uFKtMNfma3nVgyLA91VgO4uIJo9jLDzXHRHNrozG8Puw7BnmmFp82LYvRhWN4bFBZr1hHVmvnAFAIemLz25A9QjKJFGlHAQwkHUcACvPULg4AGUUANKyI/SfAg1dBil6RBqc7352r8fpfkwih7B8/H8dndtaHZ0Zy66Kw/dld/ymIfuLjj66C5Ad+Z1eHxC9GUS5ElgRIGH564fy6NvrmbNqtfJ2ljN1f0OkuHfjla3HTXcanZKq5to5hCiWcMJFV1CNGMwuncQUc8AdM8ADHtm32qlqhbzl4bNnKJTBwyfi1BGJ6axfdGcpK3mB9tQG+tQmw6iNB5AbTzQ8tdGDWrwAGpjDVrDbqxVn6M01h7z14eBguHMJerph+5u+fL0J+rpj+7uH3uOVaZREH2TBHkfpTQfwlK9HnXDBjK+WkN2zQZ+3bAbbEAQasszqPKVkn3KdUSzS4lmlRD1laC78vtWUPcGixPd60T3Duh4XT1iBn2gGjVYbXY1+febXU3+/Wbg7/sMtbn+2E3tmeieAeYvRe/A2C9H8/lAdHeBtOxFQkiQ9wWGgXZoJ9a9n2LZvxrr/tVY6rYffTtzKOF+42g8/btEckeyz17CfR8c5JOKes5Qvdx1ynBK8z0n+AARo1pire4TijSiBqrQ/HtR/fvM7iv/XvNcQ8NerPvXHBP2hqKaXTWegUS9A83Q9w4y/yLyDkD3Dor95SFEPEmQJ4JhoNWXY929Cuvej7Hu/QQtWA2AbvcR7jeO5tKphAvOxD38XOqb2rbycoHHpg3kzU3V/Oq9Hdzw/OeUjR3Iv5xbhNchP9K4sDjRM4egZw45/jrhoBnuDXtaPe5BbdiDtWotavkb5sieVsxW/UCi3kFm2Lc8KgOGoZCL4cxJv7+oRLfJ//peogSqse1eia1yJdY9q9ACVQBEPQMID7qA4IBzCA84h6ivGJRWw/sdLmg6tr9YURSuGFnAhKHZPPbBTl5cs4e/bKzi++cPYeqoflg0uUSgx1ldRLOGEc0aRri99w3d7L5p2IPWsAe1YffRx8MVWPd8hBr2x1bPxTwZfbQVP6jNc93b0n2jyM9WtCVB3lP0CNb9a7BVrMBWscIcZw3ozhxCgy4gOPB8QoMmoGcUdasFlum0cs+3Srl29AAefb+cR1Zs5+Uv9vD984v45og8VGndJY6ixrpxIv3GHfu+YZgjcBr2kKFX07S/vCX0d6M27MZevQ61qe1UDYZqbemyOdqaPxr2ZncOWg9dbyD6LAnyOFKa6rFVvodt1zvYKt9DbT6EoVoI9xuP/9y7CBdeRCR3ZI+0qEYUeFhw7ShWlh9kwYc7+a+/bOaZTyr5/nlFXDQ8VwK9L1IUDIePqMOH4TuLxoJ2RuqEAmj+o+FuPpqvrV+9jz1Q3WYUjoGC7s5vGak00DwR29KaP3KS1rCfxLUIIilIkHeT6t+Lbefb2He8hXXP31GMqNnqHnopzUMuITxoIoY9o1dqURSFC4flMLEkm+VbanjqowruXLaJwiwns8cP4oqRBdgs8md5UrG5zVFJ2aXtvx8NtZyA3dMq6PeaQV+9DnXHm8f209u8bUbfmH32rUbiuPtJqz7JSJB3gXqoAnv5G9h3vIG16gsAIlnDaDzzZpqHfotI/hhQtcTVpyh865R8LinNY8W2Azz76Vc8+M42nli1i6mj+jN1dH/yPPaE1SfiSLOd+KRsrJ/+yIib3a2Cf48Z9k0H226CYl5A5emP7h0QG2ZpPprj6nVXvjn9gugT5CfRSerhSuzbX8e+/XWsNesBCOeNwn/uXYSKLyOaNSzBFR5LUxUmj8jjm6W5fFZZzwtr9vD7jyv5w6dfMWlYDlPO6MdZhVlHrxIVqadVPz2MbX+dcGOroZV7UP1Hp3TQardgq3gPJdK228dQ1KNh7+lPtOXCKfNiqn7ma3eBOUeO6HES5Ceg+vea4b1tKdbqtQCE88fgP/8emkuuQM8oTGyBnaQoCmcXZXF2URa76xt5dd0+lm3Yz/KtByjw2rnytAIuPyWfITly5WJasjrNC8qySo4z+qblpGxsvp4j4+pbHmu3YK18v80VxkfojqyWqRD6EXX3iz3X3f1QCgpRjUx0Z6607rtJvntfozTVYd/+F+zbXsO29xMAwnln4D/vpzQP+3bShPfxDPI5ue3CYm6ZMISV5bUs3bCfP35SyTMfV1Ka5+bSU/K5uDSXQT5noksVfUWrk7LR3JHHXy3UYF481XK1rBbY3/K6CjWwH9uBL81J2gw9tk0OLa17Z67ZwncXmCdrjzyPzZFjPkoLv30S5ADhRuy73sa+9TVsle+h6BEiWcMInD2H5uFTzLHdKcZmUfnmiDy+OSKPGn8zy7ce4O3N1Tz2wU4e+2Anw3LdXDQshwtKcji1wCOjXkSHDJuXaLb3+CdmwZwiIViDGtiPl0M01lS2BH1VLPCt1etRGg+0OyOnbs80g/1I8LtyMZx56K7cVsvz0J05aRX66RvkegTr7lU4ti7GtuMt1HCAqKc/jaNuoql0qtnySJPwyvPYmTl2IDPHDmTPoUbe317Le9treeaTSn73cSVZTivnDsni3CFZnFXokxOloutUS6xf3fC5aGpvyCUcnRMnWNMyL06VOQFasBo1WIMSPICl5h/m+60uqmqzC5sX3ZmD4cxtmeL4yDTHOa1e55jTHjuyk3qkTnoFuWFgqVmPfctiHNuWojbWoNszaR7+HZpLryE84Ny0v2puYKaTWeMGMWvcIOqDYT6uqGPVzoN8tPMgb24ypxEYku1k3GAfYwZmMmZgBv0y0qflI3pJ6zlx8jpYN9LYEvI1Zvg3HkANHmiZCbMWNViDdmgX1v2fozTVtunaaU23eTFapnTWHdmxgNcdWS1z4Geb897HXmeB1jcaNWkR5OqhXTi2/hn71sVY6ndgqDZCQy6hacRUQkUX95kfRl/jc1m57NR8Ljs1H90w2FYd4NPKOj6rrOetTdW8us68p2iB186fagO4bRa+3F3PiHwPblta/NMSfYHFiZ4xGD1jcMfr6lHzxG1jLWpTLUpjbcsUyEeeH0RtqjO7eWo3m1MjRxqPuzvD4jJP6Dp85nz+Dh967DHz6KO9Jfgz27nCNw5S9n+bEqjGsX0Z9m1/jo31Dg04l4Yzb6a5+AoMhy+xBSYZVVEYUeBhRIGH2WcNJqIbbK/x88Wew2zYe5hgKEp9MMwPFq5HAQqznIzI9zAsz01pnoeSXBcFXjtKmnRXiT5K1TCc2USd2UQZ3rltIo3mbRAb68yQb6pDaa6P3RpRba5veTyEenAr1qY6lOZDx1yIBRC54lEYOj2+x0SKBbnSVI99x1vYty3BumcViqETyRmJ/7z/onn4lM7NVy06xaIqnFLg5ZQCL4wdSObvMohEdX51zWlsqvKzucrPP/Yd5u0tNbFt3DaNoTkuhmS7KMpyUpTtojDLycBMBw5r4i6gEuKELE50jxM8A4h2dhvDMO94deSOVs31KOEgrpGXQMf3TDn5EuO/y96lNB/Gtutt7NuWYftqJYoeJppRRHDcf5gjTk50Bl3ElUVTuaA4hwuKc2LLGpoibDvgZ8eBIDtrg+w4GOTjXXW8vrGqzbb5HhsDfU4GZNgpLvCSZdPol2GnwOugwGvHLlMLiGSiKGBzo9vc4B0Y+wXgsrkgGP8kT84gNwzsW1/Dvu3P2L76wAxvz0AaR91I87CriOSPTpsRJ32d12Fh7CAfYwf52iz3N0eorGuksq6R3fWN7D7UxJ76Rj6rrOeNTdUYXxt5lumwkO+1k++xk+u2keOxUe0vwaoqfLH7ENkuK9kuGx67Jt03Iu0kZZCrDbvJWP5Dot5BZniXXEGk4My0H3GSTDx2CyP7eRnZ79g75jg9DrZ+VUdVQzP7G5qoamimxh+KPW6u9lMXDLG3zuwq+/7CdbFtLaqCz2nF57SS6bTgc1rJcFjIcFjJsFvwOix4Wx49dgsem2Y+2i3YNEV+CYik1GGQ67rOfffdx5YtW7DZbDzwwAMUFRXF3l+xYgW/+c1vsFgsTJs2jeuuu65HCwbQMwZz4KYNfe9GwiIu7BaVwVlOBmcd/+rSqG7wnb85CUd17p92OgeDYeqCYeoaw9QHwxxqClPfGKb8QIDDTREON0WI6MdeYNKapip4bBquI19WCy6bitOqxb4cVvO1w6LGXtst5mu7VcVuUbFr5jKbRcWmKdgtastzFYsqvyxE/HUY5MuXLycUCrFw4ULWrl3L/PnzeeKJJwAIh8PMmzePRYsW4XQ6mTlzJpMmTSIvr6OBn90no07Sm6YqWDUFq6Zx7pDsDtc3DIPGsE5Dc8T8aorgb47gD0XwN0fxN0cIhqIEQlGCoQiBUJTGcJRgKEptIEwwHKUxFKUpEqUprLdzzWHn2TQFq2YGu7Xlud2qoSnmXxQW1VxuURUsmvnaXH7ktblMa1mmtfqyKEefq4r5fdJalqmqgqaYI5DM94+uc+S5qpjrxZ7HHtsuU1o9KpjLFeVrj9Bmmfm6ZTsgYtE4HAzFlrden5bnCke35WvrKS0rtX6drr8kOwzyNWvWMHHiRADGjBnDhg0bYu+Vl5dTWFhIZqY5Uf24ceNYvXo1l19+eQ+VK0TXKIoSa2kXeLt33YBhGDRHdJoiOk1hM9jN11GaIzqhqE6o5f1wVKc5YhCKms9DLe+Ho0eXhaMGiqYSaAwT0Q1zmW7+4ok060R0g4huENUNItGvvW55jBoGkajRrV8wqeRIuB95fuTF0V8IbX85tL+d0uYP/tbLj132tc//2j7B/IU5f+oZjM53n/wBdaDDIPf7/Xg8R+/QrmkakUgEi8WC3+/H6z3ax+l2u/H7279c9uj2Cj5f12bZ0zS1y9smq6Q55vfeA8AXh1119phbPhJIgu9PBzRNJRpt/4rDk6G3CvWoYaC3hL1utAS+fnS5+fzoNoZhENVp9dzczjCIbWMAumF2bWFg7rdlHSD2OQbmLzzdaPtoLjfXUxSFqK5jmLuKfRatnh/Zz5HlrfdhtHqNAQatt2lZn5YFxB6ObtNqP0eW02a50c46R7/XRqsdt/4F2u46gKYoDMlz98j/5w6D3OPxEAgcnZ5S13UsFku77wUCgTbB3p5o1KC+vmvDb3w+V5e3TVZyzOmhJ49ZbfmK/WdXAE0xvxJIfs4nJy/v+Nna4TCPsWPHsnLlSgDWrl1LaenRcdklJSVUVFRQX19PKBRi9erVnHnmmV0qUgghRNd02CKfPHkyq1atYsaMGRiGwUMPPcSyZcsIBoOUlZVx1113cdNNN2EYBtOmTaOgoKA36hZCCNFCMVp3DPWCcDgqXSsnQY45Pcgxp4eEda0IIYTo2yTIhRAiyUmQCyFEkpMgF0KIJCdBLoQQSa7XR60IIYSIL2mRCyFEkpMgF0KIJCdBLoQQSU6CXAghkpwEuRBCJDkJciGESHIS5EIIkeSSIsh1XWfu3LmUlZUxe/ZsKioqEl1SjwuHw9xxxx3MmjWL6dOn8+677ya6pF5RW1vLhRdeSHl5eaJL6TW//e1vKSsrY+rUqbzyyiuJLqdHhcNh5syZw4wZM5g1a1bK/5zXrVvH7NmzAaioqGDmzJnMmjWLn/3sZ+h69+8IdURSBHnrG0DPmTOH+fPnJ7qkHrd06VJ8Ph8vvPACTz/9ND//+c8TXVKPC4fDzJ07F4fDkehSes0nn3zCF198wYsvvshzzz3H/v37E11Sj3r//feJRCK89NJL3HrrrTz66KOJLqnHPP3009xzzz00NzcDMG/ePG6//XZeeOEFDMOIa+MsKYL8RDeATlWXXXYZt912W+y1pmkJrKZ3PPzww8yYMYP8/PxEl9JrPvzwQ0pLS7n11lu5+eabueiiixJdUo8aOnQo0WgUXdfx+/2x20amosLCQh577LHY640bN3L22WcD8I1vfIOPPvoobp+VFN/FE90AOlW53eadtv1+Pz/84Q+5/fbbE1tQD1u8eDHZ2dlMnDiRp556KtHl9Jq6ujr27t3Lk08+ye7du7nlllt466232tyFPZW4XC727NnD5ZdfTl1dHU8++WSiS+oxl156Kbt37469NlpuOA3m/++Ghoa4fVZStMhPdAPoVLZv3z5uuOEGpkyZwlVXXZXocnrUq6++ykcffcTs2bPZtGkTd955JzU1NYkuq8f5fD4uuOACbDYbxcXF2O12Dh48mOiyeswf//hHLrjgAv7617+yZMkS7rrrrljXQ6pT1aNxGwgEyMjIiN++47anHnSiG0CnqgMHDnDjjTdyxx13MH369ESX0+P+7//+j+eff57nnnuOU089lYcffpi8vLxEl9Xjxo0bxwcffIBhGFRVVdHY2IjP50t0WT0mIyMDr9e8ZVlmZiaRSIRoNJrgqnrHyJEj+eSTTwBYuXIl48ePj9u+k6JZ294NoFPdk08+yeHDh1mwYAELFiwAzJMn6XQiMB1MmjSJzz77jOnTp2MYBnPnzk3p8yHf+973+OlPf8qsWbMIh8P86Ec/wuVyJbqsXnHnnXdy77338stf/pLi4mIuvfTSuO1bprEVQogklxRdK0IIIY5PglwIIZKcBLkQQiQ5CXIhhEhyEuRCCJHkJMiFECLJSZALIUSS+/9Shnp1ffoAigAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEFCAYAAAD+A2xwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAePElEQVR4nO3df1xUZb4H8M+Zc0B+DDG60o/XFqgk1Uq+EN19bT/I1ULXH6hJXdAWusnNH7tmpZbSppEiYnX31U2vvaytrumuQmStdNcsNlvu4ktbUFLKH+UaaRppijgzwjAzz/0DOTLM4MjIMDzD5/2PzDkzZ77PqB+eec5znqMIIQSIiEhahkAXQEREV4dBTkQkOQY5EZHkGORERJJjkBMRSY5BTkQkOQY59XijR4/G/v37O9x//vx5ZGdnd0st+/btw9KlSwEA+/fvx7x587rlfYkuh0FO0jt37txlg74rff3116irqwMA3H777XjllVe65X2JLodBTtK4/fbbsXr1amRmZmL06NH485//DADIzc1FY2MjJk+eDIfDgSNHjmDGjBmYOnUqJk+ejJKSEgDA7t27MWnSJGRmZiItLQ02mw35+fl48MEHMX78eIwbNw5VVVUAAIvFgtzcXIwdOxbjx4/HH/7wB5w8eRKvvPIKKisrkZubi927d2PixIkAWr4VLFy4EBMnTkRaWhpeeOEF2O32y9ZN1GUEUQ83atQosW/fPpGQkCA2bNgghBBi//79IjExUTQ2Nopjx46JpKQkIYQQzc3NYvz48aKmpkYIIURDQ4MYN26c2Lt3r9i1a5e49dZbxfHjx4UQQuzZs0c89thjwuFwCCGEWLdunZg1a5YQQoiCggLx5JNPCrvdLpqamsRDDz0kdu3aJd59910xc+ZMIYQQu3btEhMmTBBCCPH000+L5cuXC6fTKZqamsSMGTPEunXrhBCiw7qJuooW6F8kRJ1x7733AgCGDBkCm80Gq9Xqsv+bb77Bt99+i2eeeUbf1tjYiC+//BLx8fG44YYb8NOf/hQAMGzYMERHR2Pz5s04duwYdu/ejcjISADAzp07kZubC1VVoaoqNm7cCADYsmWLx7rKy8uxadMmKIqC0NBQZGZmYv369Zg5c2aHdffp06cLPxnqzRjkJJXW8FMUBQAg2i0V5HA4EBUVhb/85S/6ttOnTyMqKgrV1dWIiIjQt3/66adYsWIFHnnkEdx7770YNGgQtm7dCgDQNE1/DwA4efIkwsLCOqzL6XS6PN/pdOpDK1dSN9HV4Bg5SU/TNDgcDgghMHDgQISFhelBfvLkSUycOBE1NTVur6uoqMCoUaMwffp0JCYmoqysDA6HAwBwxx134L333oPT6YTNZsO8efPwz3/+E6qqugR0q7vvvhsbN26EEAI2mw3FxcW48847/dtwoosY5CS9mJgYDB06FBMmTIDFYsHatWtRUlKCtLQ0zJgxA48//jiGDx/u9rrMzEx89tlnSEtLw/3334+bbroJx48fh9PpxNy5cxESEoLJkydjypQpGDlyJMaMGYOkpCQcO3YMc+fOdTnWs88+izNnziAtLQ1paWkYOHAgZs+e3V0fAfVyiuB3PCIiqbFHTkQkOQY5EZHkGORERJJjkBMRSa7b55E7nU44HL6dX1VVxefXyqq3tZntDW5sr+9CQtQO93V7kDscAvX1Vu9P9MBkivD5tbLqbW1me4Mb2+u7mJioDvdxaIWISHIMciIiyTHIiYgkx0WziKjLORx2nD17Cna7zWV7XZ3SqxYM86W9mhaKvn1joKpXHs8MciLqcmfPnkJYWAQiI693WRVSVQ1wOJwBrKx7dba9QghYLA04e/YU+ve/4Ypfx6EVIupydrsNkZHXuIQ4eacoCiIjr3H7JuPNFQX5559/jqysLLftn3zyCdLT05GRkYHi4uJOvTERBTeGuG98+dy8Dq28/vrr2Lp1K8LDw122Nzc3Y+XKlSgpKUF4eDimTZuGUaNGISYmptNFXInGZge27P0OowaY+A+EiKgNr0EeGxuL1atX4+mnn3bZfuTIEcTGxiI6OhoAMHz4cFRWVmLcuHGXPZ6qKjCZIi77HE8++rIOi7bsx/Z5d2NQTGSnXy8rVTX49HnJiu0NDnV1ClTV8xf+jrZ3pT17KvHeeyVYvrxQ37Z27SuIixuA778/iYqKf0BVVTzxxEIMGZLo11p8aa+idC4nvQb52LFjcfz4cbftZrMZUVGXrjSKjIyE2Wz2+oa+XtnZcL4RAHDm3AX0C+k9Q/u8Ei64BWt7hRAeT/J118lOh8PpVoPTKVBfX489e6rw2mv/g7q6Ojz77NP44x/f9lsdvrZXCPecvNyVnT7PWjEajbBYLPpji8XiEuxdTR9N6T0zl4iCwv9+UYetNd8DaPl/3BWzDyclXo8JQ67r9OtUVcPPf/5LKIqC66+//uI0ybPo27fv1RcVQD4HeXx8PGpra1FfX4+IiAhUVlYiJyenK2tzwVFxIuqMqqpKzJ07U3984sR3mDTpfkRHm/RtERGRsFjMvS/IS0tLYbVakZGRgcWLFyMnJwdCCKSnp+O66zr/G7KzBLvkRFKZMOQ6vffcnfPIhw8fgeefX6k/fvXV1YiMNMJqvTSSYLVaYDT6byShu1xRkN9444369MK0tDR9++jRozF69Gj/VNbexbGVXnRRGBF1Mbvdjs8+24Vp07Lwww8/wOkUMJlMgS7rqklzZSeHyInoakVHR2Po0CTMmvUIhBCYP39RoEvqEtIFOZOciLxJTh6B5OQRLtvmzHlM/zknZ1Z3l+RX0szj4zVARESeSRPkrXiyk4jIlURBfvFkZ4CrICLqaaQJ8tahFc5aISJyJU+QX/yTOU5E5EqeIOfJTiIij6QJch3HVojIiz17KvHcc7ku2159dTX++tdSAMDx48eQlfVvgSjNLySaR86TnUR09T788H/xzjubce7cuUCX0mWkCXLwZCeRlPocLEHYgc0AWtbZ7oqbLzfelommWx/w6bVRUddgzZrXkJEx5arr6CmkCXIOkRNRZ3ha/fA//mM27rorJYBV+Yc0Qd6KHXIiuTTd+oDeew706ofBSpqTnZfmkTPKiYjakifIA10AEVEPJc3Qij5rhR1yIvLC2+qHALB16/buLMmvpOmRs0tOROSZPEF+ETvkRESupAnyS2utMMqJiNqSJ8h5QRARkUfyBDkHyYmIPJInyNkjJyLySJogJyIiz6SZR96KJzuJ6EqcPHkCDz88DQkJt+jbhg//OR555NEAVuUf0gQ5h1aI5DVlSjiArlv98P33L1zR8wYMGIg1a1676vfr6eQJcq5HTkR+smJFHs6dO4eGhnOYMiUdW7e+hzVrXsNbb72OpqZGzJ79GF5++UUcOPAFmpvtyMmZiZSUXwW6bJ1EQX4Rk5xIOq096O5c/RAAvvnmqMtSts89l4+YmGs9Pnf48BHIyHgIAHDw4AGsWPEc6urq8PLLa/F///d3nDtXj9dffxs//nga775bzCD3hT60wiQnoivUmaGV2Ng4/eeHHsrGAw+kYdmyldA0Dd9+W4shQ4YCAH7yk/6YOfO3fqnXV5y1QkQEQFEuxeFLL63EvHkL8MYb69DQ0IABAwbg4MEvAQBmsxnz588NVJkeSdMjb8X+OBH5U3HxJvTt2w/p6f+G8PBwrFq1HPn5L6Cy8jPMmZMDh8PR42a+KKKb79TQ3OxAfb2106/bd6IBOZuq8V9TE3HnwH5+qKxnMpkifPq8ZMX2Bofvv6/F9dfHuW3v7jHyQPO1vZ4+v5iYqA6f77VH7nQ6kZeXh0OHDiE0NBT5+fmIi7v0Blu3bsVbb70Fg8GA9PR0TJ8+vdNFX4lLi2YREXVec3Mznnzyd27bY2Pj8PTTvw9ARV3Ha5CXlZXBZrOhqKgI1dXVKCwsxKuvvqrvf+GFF/DBBx8gIiICEyZMwIQJExAdHd3lhSpMciK6CiEhIUE7p9xrkFdVVSElpeWu00lJSaipqXHZf8stt+D8+fPQNA1CCCiKfxa34jK2RESeeQ1ys9kMo9GoP1ZVFXa7HZrW8tLBgwcjPT0d4eHhSE1NxTXXXHPZ46mqApMpotOFRpmbAQARkX18er2sVNXA9gaxYG1vXZ0CVfU8Ka6j7cHKl/YqSudy0muQG41GWCwW/bHT6dRD/ODBg/j000/xt7/9DREREXjqqaewbds2jBs3rsPjORzCp5M7ZnPjxT+bgvLkUEeC9WRYR9je4CCE8HiSjyc7r4wQ7jl5uZOdXn9VJCcno7y8HABQXV2NhIQEfV9UVBTCwsLQp08fqKqKfv36oaGhodNFXwmutUJE5JnXHnlqaioqKiqQmZkJIQQKCgpQWloKq9WKjIwMZGRkYPr06QgJCUFsbCzuv/9+vxR6aeSdSU5E3vXU1Q+/+uoQ/vGPcjzyyKP4+993YMiQRPTvH3NVx/Qa5AaDAcuWLXPZFh8fr/88bdo0TJs27aqKuBL6olnMcSLpRE8ZD6Dlm3VX/B8+9/5fr+h5PXH1w8GDb8HgwS2/XN55ZxMGDHjG/0HeY+hrrRARda2uXP3w2LFvsWpVPpqbmxEWFoa8vAKsXftf+vGnTcvCJ598hLFjJ+Drrw8jP38p1q59AyEhIT7XL02Qcxo5kbxae9C9YfXD//7vl/Gb3/w7fvnLO7FjRxm++uqQy/H37KkEANx55924+eYEPPXUM1cV4oBMQa4nOaOciK5MIFY//PbbWiQmtjz3vvvGwOFw4uOPP3Q5fleTZkInbyxBRP7UVasfxsUNxIEDXwAAtm//K0pKNrsdv5XBYIDTefXfUKQJcnD6IRF1g7arH2Zm/garVi3H3XePRFRUFObMycH8+XPx4IMdT/D43e8ex8aN/4O5c2di+/ZtGDOm4+tqEhOHIj//OTQ0nLuqmqVZ/fDIaQsy11ehYOJtSL3l6s7wyiRYLxjpCNsbHLj6YYses/phT3HpgiB2yYmo87py9cOetpKiPEEO/yzGRUT+4c9F9HzRlasf+nMlRV86q9KMkXPSCpE8NC0UFksDv0F3khACFksDNC20U6+TpkfOC4KI5NG3bwzOnj0Fs7neZbuiKL0q3H1pr6aFom/fzp0HlCbIuR45kTxUVUP//je4bQ/Wk7sd6a72yjO0onCtFSIiT+QJ8kAXQETUQ8kT5LwgiIjII2mCvBXHyImIXEkT5K098l8/8Rt9bWMiIllMmRKOKVPC/XJseYKco+RERB7JE+TMcSIij+QJ8kAXQETUQ0kT5K04a4WIyJU0Qd6TFt8hIupJ5AnyQBdARNRDyRPkepJzbIWIqC15gjzQBRAR9VDSBHlrl5z9cSIiV9IEOXvkRESeSRPkBiY5EZFH0gS5wlsEERF5JE2QM8eJiDyTJsg5skJE5Jk0QW7glZ1ERB5JE+TMcSIizzRvT3A6ncjLy8OhQ4cQGhqK/Px8xMXF6fv37duHwsJCCCEQExODF198EX369PFbwRwjJyJy5bVHXlZWBpvNhqKiIixYsACFhYX6PiEElixZgpUrV2LTpk1ISUnBd99955dC2SMnIvLMa4+8qqoKKSkpAICkpCTU1NTo+44ePQqTyYT169fj8OHDGDlyJAYNGnTZ46mqApMpotOFhtrsAACDQYGmGXw6hoxUtfe0FWB7g11vbq+mtfSb/dF+r0FuNpthNBrbFKbCbrdD0zScPXsWe/fuxZIlSxAXF4fZs2cjMTERd9xxR4fHczgE6uutnS60sdkBAHA6BOx2J875cAwZmUwRPn1esmJ7g1tvbq/d3nK/zvr6Cz4dKyYmqsN9XodWjEYjLBaL/tjpdELTtItFmhAXF4ebb74ZISEhSElJcemx+4PgKDkRkQuvQZ6cnIzy8nIAQHV1NRISEvR9N910EywWC2prawEAlZWVGDx4sF8K5Y0liIg88zq0kpqaioqKCmRmZkIIgYKCApSWlsJqtSIjIwMrVqzAggULIITAsGHD8Ktf/covhXKtFSIiz7wGucFgwLJly1y2xcfH6z/fcccdKCkp6frK2mGOExF5Js0FQVyPnIjIM2mCnEMrRESeSRPkvGUnEZFn8gQ5Z60QEXkkTZADLcPk7JATEbmSK8gDXQARUQ8kV5BzeIWIyI1UQc6ZK0RE7qQKcgAQHCQnInIhVZC3DK0wyYmI2pIqyDm0QkTkTrIgV9gfJyJqR6ogZ4eciMidXEHO6YdERG4kC3LOWiEiak+qIDewR05E5EaqIGeMExG5kyvImeRERG4kC3JOPyQiak+qIG+5IIhRTkTUlmRBzgXJiYjakyrIudIKEZE7uYKcZzuJiNxIFeRcNIuIyJ1UQc4hciIid5IFObvkRETtSRXknLVCROROqiDnrBUiIndSBTlPdhIRuZMqyDlGTkTkTqogb7nVGwdXiIja8hrkTqcTS5cuRUZGBrKyslBbW+vxeUuWLMFLL73U5QW2ZeAgORGRG69BXlZWBpvNhqKiIixYsACFhYVuz9m8eTMOHz7slwLb4o0liIjceQ3yqqoqpKSkAACSkpJQU1Pjsn/v3r34/PPPkZGR4Z8K2+DsQyIid5q3J5jNZhiNRv2xqqqw2+3QNA0//PAD1qxZgzVr1mDbtm1X9IaqqsBkivCpWNWgQFEUaJrB52PIRlV7T1sBtjfY9eb2alpLv9kf7fca5EajERaLRX/sdDqhaS0v+/DDD3H27FnMnDkTp06dQmNjIwYNGoSpU6d2eDyHQ6C+3upTsYoCOJ0CdrsT53w8hmxMpgifPy8Zsb3BrTe3124PBwDU11/w6VgxMVEd7vMa5MnJydixYwfGjx+P6upqJCQk6Puys7ORnZ0NANiyZQv+9a9/XTbErxanHxIRufMa5KmpqaioqEBmZiaEECgoKEBpaSmsVmu3jIu3xQuCiIjceQ1yg8GAZcuWuWyLj493e54/e+J6LZxHTkTkRroLgoiIyJVUQa4ogGCHnIjIhVRBzh45EZE7qYJcVRSOkBMRtSNVkCtca4WIyI1UQc6hFSIid5IFOTj9kIioHbmCnFcEERG5kSvIFYXTD4mI2pEsyANdARFRzyNVkCucfkhE5EaqIFc5a4WIyI1UQW7gJfpERG7kCnLefZmIyI1cQc4xciIiN3IFuQHskBMRtSNVkHPRLCIid3IFuYFBTkTUnlRBzpsvExG5kyrIVUWB4PxDIiIXUgW5QapqiYi6h1TRyJOdRETu5ArylgXJiYioDamCnBcEERG5kyrIOf2QiMiddEHOVbOIiFzJFeScR05E5EaqIG+9ZyfnkhMRXSJVkGu81xsRkRupgvxSjzzAhRAR9SBSBXlrj5w5TkR0iVRB3jqywiAnIrpE8/YEp9OJvLw8HDp0CKGhocjPz0dcXJy+/4MPPsD69euhqioSEhKQl5cHg58WRdFaj8uxFSIindfELSsrg81mQ1FRERYsWIDCwkJ9X2NjI15++WW8/fbb2Lx5M8xmM3bs2OG3YlUOrRARufHaI6+qqkJKSgoAICkpCTU1Nfq+0NBQbN68GeHh4QAAu92OPn36XPZ4qqrAZIrwqdjQkLMXj2Hw+Riy6U1tBdjeYNeb26tpLf1mf7Tfa5CbzWYYjcY2hamw2+3QNA0GgwH9+/cHAGzYsAFWqxV33XXXZY/ncAjU11t9Krb160Oz3YnzPh5DNiZThM+fl4zY3uDWm9trt7d0eOvrL/h0rJiYqA73eQ1yo9EIi8WiP3Y6ndA0zeXxiy++iKNHj2L16tV+vYuPyguCiIjceB0jT05ORnl5OQCguroaCQkJLvuXLl2KpqYmrF27Vh9i8ZfWk52McSKiS7z2yFNTU1FRUYHMzEwIIVBQUIDS0lJYrVYkJiaipKQEI0aMwMMPPwwAyM7ORmpqql+KVTn/kIjIjdcgNxgMWLZsmcu2+Ph4/eeDBw92fVUd0FTOWiEiak+qC4JCOEZORORGqiDXVI6RExG1J1eQc9EsIiI3UgV5iN4jZ5ITEbWSLMjZIyciak+yIL/YI2eQExHppAxyJ4dWiIh0kgU5h1aIiNqTLMhbh1aY5EREraQK8lCNY+RERO3JFeT6GDkREbWSK8g1Dq0QEbUnVZC3XtnpZI4TEemkCnJFUWBQFDjZIyci0kkV5ACgKDzZSUTUlnRBzh45EZErSYM80FUQEfUcEgY52CMnImpDviA3cGiFiKgt6YJcVRQ4eUUQEZFOuiA3KICDPXIiIp10Qa4aFDh4tpOISCddkBsUhT1yIqI2pAtyzaDA6RTslRMRXSRdkKsX11s532gPcCVERD2DdEGuXbxLUH1jc4ArISLqGaQL8hBDS8k/WmwBroSIqGeQL8gv9shPmRnkRESAhEHeenOJkw2NAa6EiKhnkC7IDYqCENWAb85YA10KEVGPIF2QA0BEqIoD35sDXQYRUY/gNcidTieWLl2KjIwMZGVloba21mX/J598gvT0dGRkZKC4uNhvhbZl7KPh6Bkrjtdf6Jb3IyLqybwGeVlZGWw2G4qKirBgwQIUFhbq+5qbm7Fy5Uq8+eab2LBhA4qKinDq1Cm/FgwA/SJCoBkU/OeOIzh3gdMQiah307w9oaqqCikpKQCApKQk1NTU6PuOHDmC2NhYREdHAwCGDx+OyspKjBs3rssL7XOwBOpX7wCnv4AG4O/9X8Sx4xdw+o/AO2FjMP3h+YgIVbv8fYmIejqvQW42m2E0GvXHqqrCbrdD0zSYzWZERUXp+yIjI2E2X37sWlUVmEwRnS5UiQhtuV/nol8CAG4AEB0RinqrDUnXmhDzk0iEqFIO+V+Wqhp8+rxkxfYGt97c3k8/bd3a9e33GuRGoxEWi0V/7HQ6oWmax30Wi8Ul2D1xOATq632YcRI7CaahmW6v7QNgBADL+eCcjmgyRfj2eUmK7Q1ubK/vYmI6zlavXdjk5GSUl5cDAKqrq5GQkKDvi4+PR21tLerr62Gz2VBZWYlhw4Z1QclERHSlvPbIU1NTUVFRgczMTAghUFBQgNLSUlitVmRkZGDx4sXIycmBEALp6em47rrruqNuIiK6SBGiexf3bm52+PxVo7d9LQN6X5vZ3uDG9vruqoZWiIioZ2OQExFJjkFORCQ5BjkRkeQY5EREkuv2WStERNS12CMnIpIcg5yISHIMciIiyTHIiYgkxyAnIpIcg5yISHIMciIiyXldxrYncDqdyMvLw6FDhxAaGor8/HzExcUFuqwu8/nnn+Oll17Chg0bUFtbi8WLF0NRFAwePBjPPfccDAYDiouLsXnzZmiahjlz5mDUqFGBLtsnzc3NeOaZZ/Ddd9/BZrNhzpw5uPnmm4O2zQ6HA88++yyOHj0KVVWxcuVKCCGCtr2tfvzxR0ydOhVvvvkmNE0L6vZOmTJFv6HOjTfeiNmzZ3d/e4UEtm/fLhYtWiSEEGLv3r1i9uzZAa6o67z22mti4sSJ4sEHHxRCCDFr1iyxa9cuIYQQS5YsER999JH44YcfxMSJE0VTU5NoaGjQf5ZRSUmJyM/PF0IIcebMGTFy5MigbvPHH38sFi9eLIQQYteuXWL27NlB3V4hhLDZbOK3v/2tGDNmjPj666+Dur2NjY1i8uTJLtsC0V4phlYudwNo2cXGxmL16tX64y+++AK/+MUvAAD33HMPdu7ciX379mHYsGEIDQ1FVFQUYmNjcfDgwUCVfFV+/etf4/HHH9cfq6oa1G2+7777sHz5cgDAiRMn0L9//6BuLwCsWrUKmZmZuPbaawEE97/pgwcP4sKFC5gxYways7NRXV0dkPZKEeQd3QA6GIwdO1a/ByoACCGgKAqAlptZnz9/3qebXPdUkZGRMBqNMJvNmDdvHp544omgb7OmaVi0aBGWL1+OsWPHBnV7t2zZgn79+ukdLyC4/02HhYUhJycHb7zxBp5//nksXLgwIO2VIsgvdwPoYGMwXPorsVgsuOaaa3y6yXVPdvLkSWRnZ2Py5MlIS0vrFW1etWoVtm/fjiVLlqCpqUnfHmztfffdd7Fz505kZWXhwIEDWLRoEc6cOaPvD7b2Dhw4EJMmTYKiKBg4cCBMJhN+/PFHfX93tVeKIL/cDaCDzc9+9jPs3r0bAFBeXo4RI0Zg6NChqKqqQlNTE86fP48jR45I+xmcPn0aM2bMwFNPPYUHHngAQHC3+f3338e6desAAOHh4VAUBYmJiUHb3j/96U/YuHEjNmzYgNtuuw2rVq3CPffcE7TtLSkpQWFhIQCgrq4OZrMZd911V7e3V4rVD1tnrRw+fFi/AXR8fHygy+oyx48fx/z581FcXIyjR49iyZIlaG5uxqBBg5Cfnw9VVVFcXIyioiIIITBr1iyMHTs20GX7JD8/H9u2bcOgQYP0bb///e+Rn58flG22Wq3Izc3F6dOnYbfb8eijjyI+Pj6o/45bZWVlIS8vDwaDIWjba7PZkJubixMnTkBRFCxcuBB9+/bt9vZKEeRERNQxKYZWiIioYwxyIiLJMciJiCTHICcikhyDnIhIcgxyIiLJMciJiCT3/yhmxfRzRSvXAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(1)\n",
    "plot = sns.barplot(x=\"type\", y=\"rate\", hue='exposure', data=df)\n",
    "plot.set(title=\"cell\")\n",
    "plot.margins(0.1)\n",
    "plot.set_ylim(1, 7)\n",
    "\n",
    "x = np.linspace(0, 10, num=100)\n",
    "y0 = f_test.pdf(x, dfn=df_SSB_row, dfd=df_SSW_c)\n",
    "y1 = ncf.pdf(x, dfn=df_SSB_row, dfd=df_SSW_c, nc=nc_row)\n",
    "plt.figure(2)\n",
    "plot = sns.lineplot(x=x, y=y0, label=\"H0\")\n",
    "plot = sns.lineplot(x=x, y=y1, label=\"H1\")\n",
    "plot.set(title=\"Exposure\")\n",
    "plt.vlines(f_row, 0, 0.7, 'b', label='F_row')\n",
    "plt.vlines(f_row_crit, 0, 0.7, 'r', label='F_row_crit')\n",
    "plt.legend()\n",
    "\n",
    "y0 = f_test.pdf(x, dfn=df_SSB_col, dfd=df_SSW_c)\n",
    "y1 = ncf.pdf(x, dfn=df_SSB_col, dfd=df_SSW_c, nc=nc_col)\n",
    "plt.figure(3)\n",
    "plot = sns.lineplot(x=x, y=y0, label=\"H0\")\n",
    "plot = sns.lineplot(x=x, y=y1, label=\"H1\")\n",
    "plot.set(title=\"Type\")\n",
    "plt.vlines(f_col, 0, 0.7, 'b', label='F_col')\n",
    "plt.vlines(f_col_crit, 0, 0.7, 'r', label='F_col_crit')\n",
    "plt.legend()\n",
    "\n",
    "y0 = f_test.pdf(x, dfn=df_SS_rxc, dfd=df_SSW_c)\n",
    "y1 = ncf.pdf(x, dfn=df_SS_rxc, dfd=df_SSW_c, nc=nc_rxc)\n",
    "plt.figure(4)\n",
    "plot = sns.lineplot(x=x, y=y0, label=\"H0\")\n",
    "plot = sns.lineplot(x=x, y=y1, label=\"H1\")\n",
    "plot.set(title=\"Interaction\")\n",
    "plt.vlines(f_rxc, 0, 0.7, 'b', label='F_rxc')\n",
    "plt.vlines(f_rxc_crit, 0, 0.7, 'r', label='F_rxc_crit')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ans:\n",
    "\n",
    "Since p_value > α in the main effect, we rejected H<sub>0. \n",
    "\n",
    "There is a significant main effect for type (columns) on the rating.\n",
    "\n",
    "There is a significant main effect for exposure (rows) on the rating. \n",
    "\n",
    "Since p_value > α in the interaction, we rejected H<sub>0.\n",
    "\n",
    "There is a significant interaction between exposure and type.\n",
    "\n",
    "#### Main Effects vs Interaction\n",
    "\n",
    "They hypothesized that spending more time watching YouTube videos would cause someone to think the can do the tablecloth trick. Therefore, from their hypothesis, we would predict that the degree of exposure would impact only the watched condition. \n",
    "\n",
    "The high exposure group should have a higher mean on watching over the low exposure group, but there is no differences between low and high exposure on reading or thinking about the trick.\n",
    "\n",
    "The interaction is telling us \"when we test the exposure, there is differences in different type (watching, reading, thinking). As a result, we need to test the simple effects (one factor ANOVA) of level of exposure at each type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type = Watching:\n",
      "dfn: 1\tdfd: 24\n",
      "f: 14.886075949367093\tp_value: 0.0007532957173216337\n",
      "\n",
      "type = Reading:\n",
      "dfn: 1\tdfd: 24\n",
      "f: 0.3037974683544302\tp_value: 0.5866067187582621\n",
      "\n",
      "type = Thinking:\n",
      "dfn: 1\tdfd: 24\n",
      "f: 0.3037974683544302\tp_value: 0.5866067187582621\n",
      "\n"
     ]
    }
   ],
   "source": [
    "### Subset the data to only have type = Watching\n",
    "data_w = data_array[:, 0, :]\n",
    "mean_grand_w = np.mean(data_w)\n",
    "mean_b_w = np.mean(data_w, axis=1)\n",
    "SSB_w = np.sum(5 * np.square(mean_b_w - mean_grand_w))\n",
    "df_MSB_w = num_rows - 1\n",
    "MSB_w = SSB_w / df_MSB_w\n",
    "\n",
    "### To calculate the error term, we can not use\n",
    "### SS Within in the watching subset. Instead,\n",
    "### we use the SS Within calculated from the\n",
    "### 2-way ANOVA.\n",
    "MSW_w = MSW_c\n",
    "df_MSW_w = df_SSW_c\n",
    "f_w = MSB_w / MSW_w\n",
    "p_value_w = 1 - f_test.cdf(f_w, dfn=df_MSB_w, dfd=df_MSW_w)\n",
    "print(\"type = Watching:\")\n",
    "print(\"dfn: {}\\tdfd: {}\".format(df_MSB_w, df_MSW_w))\n",
    "print(\"f: {}\\tp_value: {}\".format(f_w, p_value_w))\n",
    "print(\"\")\n",
    "\n",
    "\n",
    "### Subset the data to only have type = Reading\n",
    "data_r = data_array[:, 1, :]\n",
    "mean_grand_r = np.mean(data_r)\n",
    "mean_b_r = np.mean(data_r, axis=1)\n",
    "SSB_r = np.sum(5 * np.square(mean_b_r - mean_grand_r))\n",
    "df_MSB_r = num_rows - 1\n",
    "MSB_r = SSB_r / df_MSB_r\n",
    "\n",
    "### To calculate the error term, we can not use\n",
    "### SS Within in the watching subset. Instead,\n",
    "### we use the SS Within calculated from the\n",
    "### 2-way ANOVA.\n",
    "MSW_r = MSW_c\n",
    "df_MSW_r = df_SSW_c\n",
    "f_r = MSB_r / MSW_r\n",
    "p_value_r = 1 - f_test.cdf(f_r, dfn=df_MSB_r, dfd=df_MSW_r)\n",
    "print(\"type = Reading:\")\n",
    "print(\"dfn: {}\\tdfd: {}\".format(df_MSB_r, df_MSW_r))\n",
    "print(\"f: {}\\tp_value: {}\".format(f_r, p_value_r))\n",
    "print(\"\")\n",
    "\n",
    "\n",
    "### Subset the data to only have type = Thinking\n",
    "data_t = data_array[:, 2, :]\n",
    "mean_grand_t = np.mean(data_t)\n",
    "mean_b_t = np.mean(data_t, axis=1)\n",
    "SSB_t = np.sum(5 * np.square(mean_b_t - mean_grand_t))\n",
    "df_MSB_t = num_rows - 1\n",
    "MSB_t = SSB_t / df_MSB_t\n",
    "\n",
    "### To calculate the error term, we can not use\n",
    "### SS Within in the watching subset. Instead,\n",
    "### we use the SS Within calculated from the\n",
    "### 2-way ANOVA.\n",
    "MSW_t = MSW_c\n",
    "df_MSW_t = df_SSW_c\n",
    "f_t = MSB_t / MSW_c\n",
    "p_value_t = 1 - f_test.cdf(f_t, dfn=df_MSB_t, dfd=df_MSW_t)\n",
    "print(\"type = Thinking:\")\n",
    "print(\"dfn: {}\\tdfd: {}\".format(df_MSB_t, df_MSW_t))\n",
    "print(\"f: {}\\tp_value: {}\".format(f_t, p_value_t))\n",
    "print(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Conclusion\n",
    "\n",
    "Only type watching has a p_value less than 0.05. We can conclude that the means of ratings are significantly diffent in the exposure group only on type watching.  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('py39tf')",
   "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.9.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "df238622a495e1563c034796ccb7979ee1cb944f3608c1d858728f86cec43c42"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}