Anova_13.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.stats import t\n",
    "from scipy.stats import f_oneway\n",
    "from sklearn.metrics import mean_squared_error\n",
    "import itertools\n",
    "from scipy.stats import f\n",
    "\n",
    "df = pd.read_excel(r\"C:\\Users\\ljwil\\Desktop\\Intro STATS\\Project Stats 2\\Chapter 13\\Practice Portfolio 13 Data-3.xlsx\",sheet_name=\"ANOVA\",header=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This functions handles the pairwise calculation for the number of mnm data set\n",
    "def pairwise_caluation_num(Num_Reg,Num_Peanut,Num_PB,desired_alpha,\n",
    "RegularvsPeanut,RegularvsPB,PeanutvsPB):\n",
    "     Number_array = [Num_Reg, Num_Peanut,Num_PB]\n",
    "     g1_mean = np.mean(Num_Reg)\n",
    "     sse_group1 = np.sum((Num_Reg - g1_mean)**2)\n",
    "     g2_mean = np.mean(Num_Peanut)\n",
    "     sse_group2 = np.sum((Num_Peanut - g2_mean)**2)\n",
    "     g3_mean = np.mean(Num_PB)\n",
    "     sse_group3 = np.sum((Num_PB - g3_mean)**2)\n",
    "     sse = sse_group1+sse_group2+sse_group3\n",
    "     mse = sse/dfw\n",
    "     Bonferroni_Correction = desired_alpha/3  \n",
    "     count = 0 \n",
    "     for pair in itertools.combinations(Number_array, 2): \n",
    "             element1, element2 = pair \n",
    "             # num does the caluation for the mean difference between arrays \n",
    "             num  = round(element1.mean()- element2.mean(),2)  #rounding for the sake of getting consistant results takes \n",
    "             #Takes the square root of the MSE* 1/n + 1/n1 which  n represents the observed sample size \n",
    "             denom = round(np.sqrt(mse*((1/25)+(1/25))),9)#rounding to get consistant results of excel sheet \n",
    "             t_value = round(num/denom,9)     \n",
    "             p_value = 2 * (1 - t.cdf(abs(t_value), dfw)) #two-tailed Student's t-distributio\n",
    "             if count == 0:\n",
    "              count+=1\n",
    "              RegularvsPeanut[\"numerator\"].append(num)\n",
    "              RegularvsPeanut[\"denominator\"].append(denom)\n",
    "              RegularvsPeanut[\"t_value\"].append(t_value)\n",
    "              RegularvsPeanut[\"p-value\"].append(p_value)\n",
    "              if p_value < Bonferroni_Correction:           #check p value against the bonferroni_correction  to see signifances\n",
    "                  RegularvsPeanut[\"Yes/No\"].append(\"Yes\")\n",
    "                  continue\n",
    "              else:\n",
    "                  RegularvsPeanut[\"Yes/No\"].append(\"No\")\n",
    "                  continue\n",
    "\n",
    "             if count ==1: \n",
    "                count+=1\n",
    "                RegularvsPB[\"numerator\"].append(num)\n",
    "                RegularvsPB[\"denominator\"].append(denom)\n",
    "                RegularvsPB[\"t_value\"].append(t_value)\n",
    "                RegularvsPB[\"p-value\"].append(p_value)\n",
    "                if p_value < Bonferroni_Correction:   #check p value against the bonferroni_correction  to see signifances\n",
    "                  RegularvsPB[\"Yes/No\"].append(\"Yes\")  \n",
    "                  continue\n",
    "                else:\n",
    "                  RegularvsPB[\"Yes/No\"].append(\"No\")\n",
    "                  continue\n",
    "             if count == 2:\n",
    "                count+=1\n",
    "                PeanutvsPB[\"numerator\"].append(num)\n",
    "                PeanutvsPB[\"denominator\"].append(denom)\n",
    "                PeanutvsPB[\"t_value\"].append(t_value)\n",
    "                PeanutvsPB[\"p-value\"].append(p_value)\n",
    "                if p_value < Bonferroni_Correction:     #check p value against the bonferroni_correction  to see signifances\n",
    "                  PeanutvsPB[\"Yes/No\"].append(\"Yes\")\n",
    "                  continue\n",
    "                else:\n",
    "                  PeanutvsPB[\"Yes/No\"].append(\"No\")\n",
    "                  continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This functions handles the pairwise calculation for the number of blue MnMs data set\n",
    "def pairwise_caluation_blue(Blue_Reg,Blue_Peanuts,Blue_PB,desired_alpha,\n",
    "BlueRegularvsPeanut,BlueRegularvsPB,BluePeanutvsPB):\n",
    "     Number_array = [Blue_Reg, Blue_Peanuts,Blue_PB]\n",
    "     g1_mean = np.mean(Blue_Reg)\n",
    "     sse_group1 = np.sum((Blue_Reg - g1_mean)**2)\n",
    "     g2_mean = np.mean(Blue_Peanuts)\n",
    "     sse_group2 = np.sum((Blue_Peanuts- g2_mean)**2)\n",
    "     g3_mean = np.mean(Blue_PB)\n",
    "     sse_group3 = np.sum((Blue_PB - g3_mean)**2)\n",
    "     sse = sse_group1+sse_group2+sse_group3\n",
    "     mse = sse/dfw\n",
    "     Bonferroni_Correction = desired_alpha/3  \n",
    "     count = 0 \n",
    "     for pair in itertools.combinations(Number_array, 2): \n",
    "             element1, element2 = pair \n",
    "             # num does the caluation for the mean difference between arrays \n",
    "             num  = round(element1.mean()- element2.mean(),2)  #rounding for the sake of getting consistant results takes \n",
    "             #Takes the square root of the MSE* 1/n + 1/n1 which  n represents the observed sample size \n",
    "             denom = round(np.sqrt(mse*((1/25)+(1/25))),9)#rounding to get consistant results of excel sheet \n",
    "             t_value = round(num/denom,9)     \n",
    "             p_value = 2 * (1 - t.cdf(abs(t_value), dfw)) #two-tailed Student's t-distributio\n",
    "             if count == 0:\n",
    "              count+=1\n",
    "              BlueRegularvsPeanut[\"numerator\"].append(num)\n",
    "              BlueRegularvsPeanut[\"denominator\"].append(denom)\n",
    "              BlueRegularvsPeanut[\"t_value\"].append(t_value)\n",
    "              BlueRegularvsPeanut[\"p-value\"].append(p_value)\n",
    "              if p_value < Bonferroni_Correction:           #check p value against the bonferroni_correction  to see signifances\n",
    "                  BlueRegularvsPeanut[\"Yes/No\"].append(\"Yes\")\n",
    "                  continue\n",
    "              else:\n",
    "                  BlueRegularvsPeanut[\"Yes/No\"].append(\"No\")\n",
    "                  continue\n",
    "\n",
    "             if count ==1: \n",
    "                count+=1\n",
    "                BlueRegularvsPB[\"numerator\"].append(num)\n",
    "                BlueRegularvsPB[\"denominator\"].append(denom)\n",
    "                BlueRegularvsPB[\"t_value\"].append(t_value)\n",
    "                BlueRegularvsPB[\"p-value\"].append(p_value)\n",
    "                if p_value < Bonferroni_Correction:   #check p value against the bonferroni_correction  to see signifances\n",
    "                  BlueRegularvsPB[\"Yes/No\"].append(\"Yes\")  \n",
    "                  continue\n",
    "                else:\n",
    "                  BlueRegularvsPB[\"Yes/No\"].append(\"No\")\n",
    "                  continue\n",
    "             if count == 2:\n",
    "                count+=1\n",
    "                BluePeanutvsPB[\"numerator\"].append(num)\n",
    "                BluePeanutvsPB[\"denominator\"].append(denom)\n",
    "                BluePeanutvsPB[\"t_value\"].append(t_value)\n",
    "                BluePeanutvsPB[\"p-value\"].append(p_value)\n",
    "                if p_value < Bonferroni_Correction:     #check p value against the bonferroni_correction  to see signifances\n",
    "                  BluePeanutvsPB[\"Yes/No\"].append(\"Yes\")\n",
    "                  continue\n",
    "                else:\n",
    "                  BluePeanutvsPB[\"Yes/No\"].append(\"No\")\n",
    "                  continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Also setting appart the two tables in the dataframe into 2 other data frames will make programming this easier\n",
    "\n",
    "\n",
    "Number_of_MnMs = {\n",
    "    \"Regular\": [],\n",
    "    \"Peanut\": [],\n",
    "    \"PB\": []\n",
    "}\n",
    "Blue_MnMs = {\n",
    "    \"Regular\": [],\n",
    "    \"Peanut\": [],\n",
    "    \"PB\": []\n",
    "}\n",
    "BlueRegularvsPeanut = {\n",
    "     \"numerator\":[],\n",
    "     \"denominator\":[],\n",
    "     \"t_value\": [],\n",
    "     \"p-value\": [],\n",
    "     \"Yes/No\": []\n",
    "\n",
    "}\n",
    "BlueRegularvsPB = {\n",
    "     \"numerator\":[],\n",
    "     \"denominator\":[],\n",
    "     \"t_value\": [],\n",
    "     \"p-value\": [],\n",
    "     \"Yes/No\": []\n",
    "\n",
    "}\n",
    "BluePeanutvsPB = {\n",
    "     \"numerator\":[],\n",
    "     \"denominator\":[],\n",
    "     \"t_value\": [],\n",
    "     \"p-value\": [],\n",
    "     \"Yes/No\": []\n",
    "}\n",
    "RegularvsPeanut = {\n",
    "     \"numerator\":[],\n",
    "     \"denominator\":[],\n",
    "     \"t_value\": [],\n",
    "     \"p-value\": [],\n",
    "     \"Yes/No\": []\n",
    "\n",
    "}\n",
    "RegularvsPB =  {\n",
    "     \"numerator\":[],\n",
    "     \"denominator\":[],\n",
    "     \"t_value\": [],\n",
    "     \"p-value\": [],\n",
    "     \"Yes/No\": []\n",
    "\n",
    "}\n",
    "PeanutvsPB =  {\n",
    "     \"numerator\":[],\n",
    "     \"denominator\":[],\n",
    "     \"t_value\": [],\n",
    "     \"p-value\": [],\n",
    "     \"Yes/No\": []\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This iterates through the index to access the row data for each table in the data set and places them in the dictionaries used later for dataframes\n",
    "#Then creates the dataframes \n",
    "for index, row in df.iterrows():\n",
    "    Number_of_MnMs[\"Regular\"].append(row[\"Regular\"])\n",
    "    Number_of_MnMs[\"Peanut\"].append(row['Peanut'])\n",
    "    Number_of_MnMs[\"PB\"].append(row['PB'])\n",
    "    Blue_MnMs[\"Regular\"].append(row[\"Regular.1\"])\n",
    "    Blue_MnMs[\"Peanut\"].append(row[\"Peanut.1\"])\n",
    "    Blue_MnMs[\"PB\"].append(row['PB.1'])\n",
    "\n",
    "Number_Anova = pd.DataFrame(Number_of_MnMs)\n",
    "Blue_Anova =  pd.DataFrame(Blue_MnMs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This makes a summary of the data from the two tables and saves them for later to be used in dictionaires\n",
    "Num_Reg = Number_Anova.loc[:, 'Regular']\n",
    "Num_Reg_dev = np.std(Num_Reg)           #caluating standard deviation \n",
    "Num_Reg_var = np.var(Num_Reg)           #caluating variance\n",
    "Num_Reg_mean = np.mean(Num_Reg)\n",
    "Num_Reg_sum = np.sum(Num_Reg)\n",
    "Num_Reg_count = len(Num_Reg)\n",
    "\n",
    "Num_Peanut = Number_Anova.loc[:,'Peanut']\n",
    "Num_Peanut_dev = np.std(Num_Peanut)    #caluating standard deviation \n",
    "Num_Peanut_var = np.var(Num_Peanut)           #caluating variance\n",
    "Num_Peanut_mean = np.mean(Num_Peanut)\n",
    "Num_Peanut_sum = np.sum(Num_Peanut)\n",
    "Num_Peanut_count = len(Num_Peanut)\n",
    "\n",
    "Num_PB =  Number_Anova.loc[:,'PB']\n",
    "Num_PB_dev = np.std(Num_PB)           #caluating standard deviation \n",
    "Num_PB_var = np.var(Num_PB)           #caluating variance\n",
    "Num_PB_mean = np.mean(Num_PB)\n",
    "Num_PB_sum = np.sum(Num_PB)\n",
    "Num_PB_count = len(Num_PB)\n",
    "\n",
    "Blue_Reg = Blue_Anova.loc[:,'Regular']\n",
    "Blue_Reg_dev = np.std(Blue_Reg)         #caluating standard deviation \n",
    "Blue_Reg_var = np.var(Blue_Reg)           #caluating variance\n",
    "Blue_Reg_mean = np.mean(Blue_Reg)\n",
    "Blue_Reg_sum = np.sum(Blue_Reg)\n",
    "Blue_Reg_count= len(Blue_Reg)\n",
    "\n",
    "Blue_Peanuts = Blue_Anova.loc[:,'Peanut']\n",
    "Blue_Peanuts_dev = np.std(Blue_Peanuts)    #caluating standard deviation \n",
    "Blue_Peanuts_var = np.var(Blue_Peanuts)           #caluating variance\n",
    "Blue_Peanuts_mean = np.mean(Blue_Peanuts)\n",
    "Blue_Peanuts_sum = np.sum(Blue_Peanuts)\n",
    "Blue_Peanuts_count =len(Blue_Peanuts)\n",
    "\n",
    "Blue_PB = Blue_Anova.loc[:,'PB'] \n",
    "Blue_PB_dev = np.std(Blue_PB)            #caluating standard deviation \n",
    "Blue_PB_var = np.var(Blue_PB)           #caluating variance\n",
    "Blue_PB_mean = np.mean(Blue_PB)\n",
    "Blue_PB_sum = np.sum(Blue_PB)\n",
    "Blue_PB_count =len(Blue_PB)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [],
   "source": [
    "desired_alpha = 0.05 \n",
    "Bonferroni_Correction =desired_alpha/3 #3 representing the number of categories \n",
    "dfw = len(Num_Reg) + len(Num_Peanut) + len(Num_PB) - 3  #dfw equals 72\n",
    "pairwise_caluation_num(Num_Reg,Num_Peanut,Num_PB,desired_alpha,RegularvsPeanut,RegularvsPB,PeanutvsPB)\n",
    "pairwise_caluation_blue(Blue_Reg,Blue_Peanuts,Blue_PB,desired_alpha,BlueRegularvsPeanut,BlueRegularvsPB,BluePeanutvsPB)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_numMnMs = {\"Regular vs Peanuts\":RegularvsPeanut,\n",
    "             \"Regular vs PB\":RegularvsPB,\n",
    "             \"Peanut vs PB\":PeanutvsPB}\n",
    "data_blueMnMs = {\"Regular vs Peanuts\":BlueRegularvsPeanut,\n",
    "             \"Regular vs PB\":BlueRegularvsPB,\n",
    "             \"Peanut vs PB\":BluePeanutvsPB}\n",
    "\n",
    "Numb_Anova_table =  pd.DataFrame(data_numMnMs).transpose()  # contains the first table for Number of Mnms for anova caluations\n",
    "Blue_Anova_table =  pd.DataFrame(data_blueMnMs).transpose()  # contains the second table for Blue Mnms for anova caluations\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "Number_Summary = {\n",
    "   \"Regular\": [Num_Reg_count,Num_Reg_sum,Num_Reg_mean,Num_Reg_var,Num_Reg_dev],\n",
    "   \"Peanut\":  [Num_Peanut_count,Num_Peanut_sum,Num_Peanut_mean,Num_Peanut_var,Num_Peanut_dev],\n",
    "   \"PB\": [Num_PB_count,Num_PB_sum,Num_PB_mean,Num_PB_var,Num_PB_dev]\n",
    "}\n",
    "Anova_Number_Summary = {\n",
    "   \"Between Groups\":[],\n",
    "   \"Within Groups\": [],\n",
    "   \"Total\":[]\n",
    "\n",
    "}\n",
    "Blue_Summary = {\n",
    "   \"Regular\": [Blue_Reg_count,Blue_Reg_sum,Blue_Reg_mean,Blue_Reg_var,Blue_Reg_dev],\n",
    "     \"Peanut\":  [Blue_Peanuts_count,Blue_Peanuts_sum,Blue_Peanuts_mean,Blue_Peanuts_var,Blue_Peanuts_dev],\n",
    "     \"PB\": [Blue_PB_count,Blue_PB_sum,Blue_PB_mean,Blue_PB_var,Blue_PB_dev]\n",
    "}\n",
    "Anova_Blue_Summary = {\n",
    "   \"Between Groups\":[],\n",
    "   \"Within Groups\":[],\n",
    "   \"Total\":[]\n",
    "   \n",
    "}\n",
    "labels = ['Count', 'Sum', 'Average','Variance',\"Standard Devivation\"]\n",
    "Number_Summary_df = pd.DataFrame(Number_Summary).transpose()\n",
    "Number_Summary_df.columns = labels\n",
    "Blue_Summary_df = pd.DataFrame(Blue_Summary).transpose()\n",
    "Blue_Summary_df.columns = labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Anova Summary table creation  for number of MnMs\n",
    "\n",
    "# Combine all data into one array\n",
    "all_data = np.concatenate([Num_Reg, Num_Peanut, Num_PB])\n",
    "\n",
    "# Calculate overall mean\n",
    "overall_mean = np.mean(all_data)\n",
    "\n",
    "# Calculate group means\n",
    "group1_mean = np.mean(Num_Reg)\n",
    "group2_mean = np.mean(Num_Peanut)\n",
    "group3_mean = np.mean(Num_PB)\n",
    "\n",
    "# Calculate SS between groups\n",
    "ssb_group1 = len(Num_Reg) * (group1_mean - overall_mean)**2\n",
    "ssb_group2 = len(Num_Peanut) * (group2_mean - overall_mean)**2\n",
    "ssb_group3 = len(Num_PB) * (group3_mean - overall_mean)**2\n",
    "\n",
    "ssb_total = round(ssb_group1 + ssb_group2 + ssb_group3,7) #rounding to get consistent result\n",
    "df = len(['Regular','Peanut','PB']) -1 #because this df is based off the number of catecoriges \n",
    "#print(df)\n",
    "#print(ssb_total)\n",
    "#print(ssb_total/df)\n",
    "msb = ssb_total/df\n",
    "\n",
    "g1_mean = np.mean(Num_Reg)\n",
    "sse_group1 = np.sum((Num_Reg - g1_mean)**2)\n",
    "g2_mean = np.mean(Num_Peanut)\n",
    "sse_group2 = np.sum((Num_Peanut - g2_mean)**2)\n",
    "g3_mean = np.mean(Num_PB)\n",
    "sse_group3 = np.sum((Num_PB - g3_mean)**2)\n",
    "sse = sse_group1+sse_group2+sse_group3\n",
    "\n",
    "dfw = len(Num_Reg) + len(Num_Peanut) + len(Num_PB) - 3  #df equals 72\n",
    "mse = sse/dfw\n",
    "#print(mse)\n",
    "f_table = msb/mse\n",
    "p_value = 1 - f.cdf(f_table, df, dfw)\n",
    "alpha = 0.05\n",
    "f_critical = f.ppf(1 - alpha, df, dfw)\n",
    "\n",
    "Anova_Number_Summary['Between Groups'].append(ssb_total)\n",
    "Anova_Number_Summary[\"Between Groups\"].append(df)\n",
    "Anova_Number_Summary['Between Groups'].append(msb)\n",
    "Anova_Number_Summary[\"Between Groups\"].append(f_table)\n",
    "Anova_Number_Summary['Between Groups'].append(p_value)\n",
    "Anova_Number_Summary['Between Groups'].append(f_critical)\n",
    "\n",
    "#print(Anova_Number_Summary[\"Between Groups\"])\n",
    "Anova_Number_Summary[\"Within Groups\"].append(sse)\n",
    "Anova_Number_Summary['Within Groups'].append(dfw)\n",
    "Anova_Number_Summary['Within Groups'].append(mse)\n",
    "Anova_Number_Summary[\"Within Groups\"].append(np.nan)\n",
    "Anova_Number_Summary[\"Within Groups\"].append(np.nan)\n",
    "Anova_Number_Summary[\"Within Groups\"].append(np.nan)\n",
    "#print(Anova_Number_Summary[\"Within Groups\"])\n",
    "Anova_Number_Summary[\"Total\"].append(sse+ssb_total)\n",
    "Anova_Number_Summary[\"Total\"].append(df+dfw)\n",
    "Anova_Number_Summary[\"Total\"].append(np.nan)\n",
    "Anova_Number_Summary[\"Total\"].append(np.nan)\n",
    "Anova_Number_Summary[\"Total\"].append(np.nan)\n",
    "Anova_Number_Summary[\"Total\"].append(np.nan)\n",
    "\n",
    "\n",
    "Anova_num_sum_df = pd.DataFrame(Anova_Number_Summary).transpose()\n",
    "labels_Anovas = ['SS','df','MS','F','p-value','F Critical']\n",
    "Anova_num_sum_df.columns = labels_Anovas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Anova Summary table creation  for number of blue MnMs\n",
    "# Combine all data into one array\n",
    "all_data = np.concatenate([Blue_Reg, Blue_Peanuts, Blue_PB])\n",
    "\n",
    "# Calculate overall mean\n",
    "overall_mean = np.mean(all_data)\n",
    "\n",
    "# Calculate group means\n",
    "group1_mean = np.mean(Blue_Reg)\n",
    "group2_mean = np.mean(Blue_Peanuts)\n",
    "group3_mean = np.mean(Blue_PB)\n",
    "\n",
    "# Calculate SS between groups\n",
    "ssb_group1 = len(Blue_Reg) * (group1_mean - overall_mean)**2\n",
    "ssb_group2 = len(Blue_Peanuts) * (group2_mean - overall_mean)**2\n",
    "ssb_group3 = len(Blue_PB) * (group3_mean - overall_mean)**2\n",
    "\n",
    "ssb_total = round(ssb_group1 + ssb_group2 + ssb_group3,7) #rounding to get consistent result\n",
    "df = len(['Regular','Peanut','PB']) -1 #because this df is based off the number of catecoriges \n",
    "#print(df)\n",
    "#print(ssb_total)\n",
    "#print(ssb_total/df)\n",
    "msb = ssb_total/df\n",
    "\n",
    "g1_mean = np.mean(Blue_Reg)\n",
    "sse_group1 = np.sum((Blue_Reg - g1_mean)**2)\n",
    "g2_mean = np.mean(Blue_Peanuts)\n",
    "sse_group2 = np.sum((Blue_Peanuts - g2_mean)**2)\n",
    "g3_mean = np.mean(Blue_PB)\n",
    "sse_group3 = np.sum((Blue_PB - g3_mean)**2)\n",
    "sse = sse_group1+sse_group2+sse_group3\n",
    "\n",
    "dfw = len(Blue_Reg) + len(Blue_Peanuts) + len(Blue_PB) - 3  #df equals 72\n",
    "mse = sse/dfw\n",
    "#print(mse)\n",
    "f_table = msb/mse\n",
    "p_value = 1 - f.cdf(f_table, df, dfw)\n",
    "alpha = 0.05\n",
    "f_critical = f.ppf(1 - alpha, df, dfw)\n",
    "\n",
    "\n",
    "Anova_Blue_Summary['Between Groups'].append(ssb_total)\n",
    "Anova_Blue_Summary[\"Between Groups\"].append(df)\n",
    "Anova_Blue_Summary['Between Groups'].append(msb)\n",
    "Anova_Blue_Summary[\"Between Groups\"].append(f_table)\n",
    "Anova_Blue_Summary['Between Groups'].append(p_value)\n",
    "Anova_Blue_Summary['Between Groups'].append(f_critical)\n",
    "\n",
    "\n",
    "Anova_Blue_Summary[\"Within Groups\"].append(sse)\n",
    "Anova_Blue_Summary['Within Groups'].append(dfw)\n",
    "Anova_Blue_Summary['Within Groups'].append(mse)\n",
    "Anova_Blue_Summary[\"Within Groups\"].append(np.nan)\n",
    "Anova_Blue_Summary[\"Within Groups\"].append(np.nan)\n",
    "Anova_Blue_Summary[\"Within Groups\"].append(np.nan)\n",
    "\n",
    "Anova_Blue_Summary[\"Total\"].append(sse+ssb_total)\n",
    "Anova_Blue_Summary[\"Total\"].append(df+dfw)\n",
    "Anova_Blue_Summary[\"Total\"].append(np.nan)\n",
    "Anova_Blue_Summary[\"Total\"].append(np.nan)\n",
    "Anova_Blue_Summary[\"Total\"].append(np.nan)\n",
    "Anova_Blue_Summary[\"Total\"].append(np.nan)\n",
    "\n",
    "Anova_blue_sum_df = pd.DataFrame(Anova_Blue_Summary).transpose()\n",
    "Anova_blue_sum_df.columns = labels_Anovas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SECTION NUMBER OF MNM\n",
      "\n",
      "\n",
      "Summary of Number of MnMs\n",
      "\n",
      "         Count    Sum  Average  Variance  Standard Devivation\n",
      "Regular   25.0  384.0    15.36    2.1504             1.466424\n",
      "Peanut    25.0  321.0    12.84    3.1744             1.781685\n",
      "PB        25.0  385.0    15.40    7.2800             2.698148\n",
      "\n",
      "\n",
      "Anova Table for Number of MnMs\n",
      "\n",
      "                        SS    df         MS          F   p-value  F Critical\n",
      "Between Groups  107.546667   2.0  53.773333  12.286367  0.000026    3.123907\n",
      "Within Groups   315.120000  72.0   4.376667        NaN       NaN         NaN\n",
      "Total           422.666667  74.0        NaN        NaN       NaN         NaN\n",
      "\n",
      "\n",
      "Fisher's LSD multiple comparison test on the number of Mnms\n",
      "\n",
      "                   numerator    denominator         t_value  \\\n",
      "Regular vs Peanuts    [2.52]  [0.591720655]   [4.258766326]   \n",
      "Regular vs PB        [-0.04]  [0.591720655]  [-0.067599465]   \n",
      "Peanut vs PB         [-2.56]  [0.591720655]  [-4.326365792]   \n",
      "\n",
      "                                     p-value Yes/No  \n",
      "Regular vs Peanuts   [6.116975715353412e-05]  [Yes]  \n",
      "Regular vs PB            [0.946291854992467]   [No]  \n",
      "Peanut vs PB        [4.8010340664150775e-05]  [Yes]  \n",
      "\n",
      "\n",
      "An analysis of variance showed that the effect of type of MnM was significant,F(2,75) = 12.29 and p < 0.001.\n",
      "Post hoc analyses using Fisher's LSD indicated that the number of peanut MnMs was significantly lower\n",
      "(M = 12.84, SD = 1.782 and penut butter (M = 15.4), SD = 2.698)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"SECTION NUMBER OF MNM\\n\\n\")\n",
    "print(\"Summary of Number of MnMs\\n\")\n",
    "print(Number_Summary_df)\n",
    "print(\"\\n\")\n",
    "print(\"Anova Table for Number of MnMs\\n\")\n",
    "print(Anova_num_sum_df)\n",
    "print(\"\\n\")\n",
    "print(\"Fisher's LSD multiple comparison test on the number of Mnms\\n\")\n",
    "print(Numb_Anova_table)\n",
    "print(\"\\n\")\n",
    "print(f\"An analysis of variance showed that the effect of type of MnM was significant,F(2,75) = {round(Anova_num_sum_df.iat[0,3],2)} and p < 0.001.\\nPost hoc analyses using Fisher's LSD indicated that the number of peanut MnMs was significantly lower\\n(M = {Number_Summary_df.iat[1,2]}, SD = {round(Number_Summary_df.iat[1,4],3)} and penut butter (M = {Number_Summary_df.iat[2,2]}), SD = {round(Number_Summary_df.iat[2,4],3)})\\n\") \n",
    "f_statistic, p_value = f_oneway(Num_Reg, Num_Peanut, Num_PB)\n",
    "f_stat_7d = round(f_statistic,7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SECTION BLUE MNM\n",
      "\n",
      "\n",
      "Summary of Blue of MnMs\n",
      "\n",
      "         Count    Sum  Average  Variance  Standard Devivation\n",
      "Regular   25.0   60.0     2.40    1.8400             1.356466\n",
      "Peanut    25.0  127.0     5.08    2.5536             1.597999\n",
      "PB        25.0   43.0     1.72    1.3216             1.149609\n",
      "\n",
      "\n",
      "Anova Table for Blue MnMs\n",
      "\n",
      "                        SS    df         MS          F       p-value  \\\n",
      "Between Groups  157.786667   2.0  78.893333  39.755879  2.332912e-12   \n",
      "Within Groups   142.880000  72.0   1.984444        NaN           NaN   \n",
      "Total           300.666667  74.0        NaN        NaN           NaN   \n",
      "\n",
      "                F Critical  \n",
      "Between Groups    3.123907  \n",
      "Within Groups          NaN  \n",
      "Total                  NaN  \n",
      "\n",
      "\n",
      "Fisher's LSD multiple comparison test on the Blue Mnms\n",
      "\n",
      "                   numerator    denominator         t_value  \\\n",
      "Regular vs Peanuts   [-2.68]  [0.398441408]  [-6.726208537]   \n",
      "Regular vs PB         [0.68]  [0.398441408]   [1.706649927]   \n",
      "Peanut vs PB          [3.36]  [0.398441408]   [8.432858464]   \n",
      "\n",
      "                                    p-value Yes/No  \n",
      "Regular vs Peanuts  [3.520710833981866e-09]  [Yes]  \n",
      "Regular vs PB         [0.09219698837298851]   [No]  \n",
      "Peanut vs PB        [2.384314967684986e-12]  [Yes]  \n",
      "\n",
      "\n",
      "An analysis of variance showed that the effect of the number of blue MnMs was significant,\n",
      "F(2,75) = 39.756 and p < 0.001.\n",
      "Post hoc analyses using Fisher's LSD indicated that the amount of Blue MnMs in peanut MnMs were significantly\n",
      "higher (M = 5.08, SD = 1.598) compared to Regular MnMs (M = 2.4, SD = 1.356)\n",
      "and Peanut Butter(M = 1.72, SD =1.15 \n"
     ]
    }
   ],
   "source": [
    "print(\"SECTION BLUE MNM\\n\\n\")\n",
    "print(\"Summary of Blue of MnMs\\n\")\n",
    "print(Blue_Summary_df)\n",
    "print(\"\\n\")\n",
    "print(\"Anova Table for Blue MnMs\\n\")\n",
    "print(Anova_blue_sum_df)\n",
    "print(\"\\n\")\n",
    "print(\"Fisher's LSD multiple comparison test on the Blue Mnms\\n\")\n",
    "print(Blue_Anova_table)\n",
    "print(\"\\n\")\n",
    "print(f\"An analysis of variance showed that the effect of the number of blue MnMs was significant,\\nF(2,75) = {round(Anova_blue_sum_df.iat[0,3],3)} and p < 0.001.\\nPost hoc analyses using Fisher's LSD indicated that the amount of Blue MnMs in peanut MnMs were significantly\\nhigher (M = {Blue_Summary_df.iat[1,2]}, SD = {round(Blue_Summary_df.iat[1,4],3)}) compared to Regular MnMs (M = {Blue_Summary_df.iat[0,2]}, SD = {round(Blue_Summary_df.iat[0,4],3)})\\nand Peanut Butter(M = {Blue_Summary_df.iat[2,2]}, SD ={round(Blue_Summary_df.iat[2,4],3)} \")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGdCAYAAADnrPLBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtY0lEQVR4nO3deXQUZb7/8U+ngSZkkwCGBANmAKPsqOiMosAVBxFBZPQqA4rgKCLIEhWIyhIEAjhyo4CAXiE4Z0BcwEEYw0VkGdzYBhVlWCKCP1kCgrQJ2ECnf38w6TEmJF3QTfVDv1/neE7q6arqb3K+dn+oeqrK4fP5fAIAADBUlN0FAAAAnA/CDAAAMBphBgAAGI0wAwAAjEaYAQAARiPMAAAAoxFmAACA0QgzAADAaFXsLiDUiouLtW/fPsXFxcnhcNhdDgAACIDP59NPP/2klJQURUVVfOzlog8z+/btU2pqqt1lAACAc/Ddd9/psssuq3Cdiz7MxMXFSTrzx4iPj7e5GgAAEAi3263U1FT/93hFLvowU3JqKT4+njADAIBhApkiwgRgAABgNMIMAAAwGmEGAAAYjTADAACMRpgBAABGI8wAAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEa76J+aDQDA+dq/f7/2799vebvk5GQlJyeHoCL8EmEGAIBKzJ49W1lZWZa3GzNmjMaOHRv8glAKYQYAgEr0799f3bp1KzV24sQJtW3bVpK0bt06RUdHl9mOozIXBmEGAIBKlHe6qKioyP9zq1atFBMTc6HLwr8xARgAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDTCDAAAMBphBgAAGM3WMLN27Vp17dpVKSkpcjgcevfdd8uss23bNnXr1k0JCQmKiYlRmzZttHfv3gtfLAAACEu2Ps6gqKhILVu2VL9+/dSjR48yr+fn56tt27Z66KGHlJWVpfj4eH311VeqXr26DdVGLp4WCwAIZ7aGmc6dO6tz585nff2ZZ57R7bffrilTpvjHGjZseCFKwy/wtFgAQDgL2wdNFhcXa9myZRo+fLg6deqkf/7zn0pLS1NmZqa6d+9+1u08Ho88Ho9/2e12S5K8Xq+8Xm+oy74o/elPf1KXLl1KjZ04cULt2rWTJK1Zs+asT4vlbw7gYvXLzze+Y4LPyt8zbMNMQUGBCgsLNWnSJI0fP16TJ09WXl6eevTooVWrVvm/SH8tOzu73KMI+fn5io2NDXXZF624uLhSy06n0/9zbGysatSoUWabwsJC7dy5M+S1AYAdjh8/7v95165d5X4O4twVFhYGvK7D5/P5QlhLwBwOhxYvXuw/6rJv3z7Vq1dPPXv21Pz58/3rdevWTTExMVqwYEG5+ynvyExqaqqOHDmi+Pj4kP4OkaSoqEgJCQmSpGPHjikmJsbmigDgwuJzMLTcbrcSExN17NixSr+/w/bITO3atVWlShU1adKk1PhVV12ldevWnXU7l8sll8tVZtzpdJY6moDz88u/JX9bAJGIz8HQsvL3DNv7zFSrVk1t2rTR9u3bS43v2LFDDRo0sKkqAAAQbmw9MlNYWKhdu3b5l3fv3q0tW7YoMTFR9evX11NPPaV7771XN998szp06KC8vDy99957Wr16tX1FAwCAsGJrmNm4caM6dOjgX87IyJAk9enTR7m5ubrrrrs0a9YsZWdna/DgwUpPT9c777yjtm3b2lUyAAAIM2EzAThU3G63EhISAppAhMAVFRX5rw4rLCxk4huAiMPnYGhZ+f4O2zkzAAAAgSDMAAAAoxFmAACA0QgzAADAaIQZAABgNMIMAAAwWtg+zsAYh8baXYE9ik7+5+dDE6Tj1eyrxU51xtpdAWA/Pgf5HLQZR2YAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDTCDAAAMBphBgAAGI0wAwAAjFbF7gIQ/vYf/En7DxaWGjtx4pT/5y1bDyg6umqZ7ZKTYpWcFBfy+gAAkY0wg0rNnrdJWX9ec9bX23adW+74mCfbaezw9iGqCgCAMwgzqFT/Pteo223plrdLTooNQTUAAJRGmEGlkpPiOF0EAAhbTAAGAABGI8wAAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADCarWFm7dq16tq1q1JSUuRwOPTuu++edd1HH31UDodDOTk5F6w+AAAQ/mwNM0VFRWrZsqVmzJhR4XqLFy/Wp59+qpSUlAtUGQAAMIWtN83r3LmzOnfuXOE633//vR5//HEtX75cXbp0uUCVAQDwHzyjLryF9R2Ai4uLdf/99+upp55S06ZNA9rG4/HI4/H4l91utyTJ6/XK6/WGoEhH8PcJc4SipwDTRMDn4MzcTXruBevPqBv1RDuNeapDqMoKDyH6HLTynR3WYWby5MmqUqWKBg8eHPA22dnZysrKKjOen5+v2NgQPCuoKCn4+4Q5ju60uwLAfhHwOdjx1s5q2foGy9vVqZ2gnf8vIQQVhZEQfQ4WFhZWvtK/hW2Y2bRpk1588UVt3rxZDkfgqT8zM1MZGRn+ZbfbrdTUVDVs2FDx8fHBL/TwwuDvE+ao3djuCgD7RcDnYOPLJMl1Dlv+/O//LmIh+hwsObMSiLANM//4xz9UUFCg+vXr+8e8Xq+eeOIJ5eTk6Ntvvy13O5fLJZerbMM5nU45nc7gFxrlC/4+YY5Q9BRgGj4HI1uIPgetfGeHbZi5//771bFjx1JjnTp10v3336++ffvaVBUAAAg3toaZwsJC7dq1y7+8e/dubdmyRYmJiapfv75q1apVav2qVauqbt26Sk9Pv9ClAgCAMGVrmNm4caM6dPjPLO+SuS59+vRRbm6uTVUBCDf79+/X/v37LW+XnJys5OTkEFQEIJzYGmbat28vny/wc61nmycD4OI2e/bscq9SrMyYMWM0duzY4BcEIKyE7ZwZACjRv39/devWrdTYiRMn1LZtW0nSunXrFB0dXWY7jsoAkYEwAyDslXe6qKioyP9zq1atFBMTc6HLAhAmeGo2AAAwGmEGAAAYjTADAACMRpgBAABGI8wAAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGhV7C4AwHk6NNbuCuxRdPI/Px+aIB2vZl8tdqoz1u4KANtxZAYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDTCDAAAMBphBgAAGI0wAwAAjEaYAQAARiPMAAAAoxFmAACA0WwNM2vXrlXXrl2VkpIih8Ohd9991//aqVOnNGLECDVv3lwxMTFKSUnRAw88oH379tlXMAAACDu2hpmioiK1bNlSM2bMKPPa8ePHtXnzZo0aNUqbN2/WokWLtH37dnXr1s2GSgEAQLiqYuebd+7cWZ07dy73tYSEBK1YsaLU2PTp03Xddddp7969ql+//oUoEQAAhDlbw4xVx44dk8Ph0CWXXHLWdTwejzwej3/Z7XZLkrxer7xeb/CLKnYEf58wRyh6yqoI7UHvL35vb7Gj1HJEoQdhtxD1oJXvbGPCzM8//6wRI0aoZ8+eio+PP+t62dnZysrKKjOen5+v2NjY4BdWlBT8fcIcR3faXUFE9GDB4WM6dPhYqbGffz7p//m9VR5Vr+4rs12d2gm6tHZCyOuzFT0Iu4WoBwsLCwNe14gwc+rUKf33f/+3fD6fZs6cWeG6mZmZysjI8C+73W6lpqaqYcOGFYagc3Z4YfD3CXPUbmx3BRHRg/MXrNJzL6w56+u9//RCueOjnminMU91CFVZ4YEehN1C1IMlZ1YCEfZhpiTI7NmzRx9++GGlgcTlcsnlcpUZdzqdcjqdwS8wquy/BhFBQtFTVkVADw548Bp175xuebvkpFg5L/a/Dz0Iu4WoB618Z4d1mCkJMjt37tSqVatUq1Ytu0sCYIPkpDglJ8XZXQaAMGVrmCksLNSuXbv8y7t379aWLVuUmJio5ORk3X333dq8ebOWLl0qr9erAwcOSJISExNVrVo1u8oGAABhxNYws3HjRnXo8J/z2SVzXfr06aOxY8dqyZIlkqRWrVqV2m7VqlVq3779hSoTAACEMVvDTPv27eXznf1ca0WvAQAASDybCQAAGI4wAwAAjGY5zGzevFlffvmlf/lvf/ubunfvrqefflonT56sYEsAAIDgsxxm+vfvrx07dkiSvvnmG913332qUaOG3nrrLQ0fPjzoBQIAAFTEcpjZsWOH/+qit956SzfffLPmz5+v3NxcvfPOO8GuDwAAoEKWw4zP51NxcbEk6YMPPtDtt98uSUpNTdXhw4eDWx0AAEAlLIeZa6+9VuPHj9df/vIXrVmzRl26dJF05oZ3SUk8bAwAAFxYlsNMTk6ONm/erEGDBumZZ55Ro0aNJElvv/22brjhhqAXCAAAUBHLN81r0aJFqauZSjz//POheZAjAABABc7rDsCFhYX++TMlqlatel4FAQAAWGH5NNPu3bvVpUsXxcTEKCEhQTVr1lTNmjV1ySWXqGbNmqGoEQAA4KwsH5np3bu3fD6f5syZo6SkJDkcjlDUBQAAEBDLYebzzz/Xpk2blJ6eHop6AAAALLF8mqlNmzb67rvvQlELAACAZZaPzPzv//6vHn30UX3//fdq1qxZmQm/LVq0CFpxAAAAlbEcZg4dOqT8/Hz17dvXP+ZwOOTz+eRwOOT1eoNaIAAAQEUsh5l+/fqpdevWWrBgAROAAQCA7SyHmT179mjJkiX+O/8CAADYyfIE4P/6r//S559/HopaAAAALLN8ZKZr164aNmyYvvzySzVv3rzMBOBu3boFrTgAAIDKWA4zjz76qCRp3LhxZV5jAjAAALjQLIeZXz+LCQAAwE6W58wAAACEk3N6avbKlSu1cuVKFRQUlDlSM2fOnKAUBgAAEAjLYSYrK0vjxo3Ttddeq+TkZO4zAwAAbGU5zMyaNUu5ubm6//77Q1EPAACAJZbnzJw8eVI33HBDKGoBAACwzHKY+dOf/qT58+eHohYAAADLLJ9m+vnnn/XKK6/ogw8+UIsWLcrcNG/q1KlBKw4AAKAylsPMF198oVatWkmStm7dWuo1JgMDAIALzXKYWbVqVSjqAAAAOCe23jRv7dq16tq1q1JSUuRwOPTuu++Wet3n82n06NFKTk5WdHS0OnbsqJ07d9pTLAAACEsBH5np169fQOtZuWleUVGRWrZsqX79+qlHjx5lXp8yZYpeeuklzZs3T2lpaRo1apQ6deqkr7/+WtWrVw/4fQAAwMUr4DCTm5urBg0aqHXr1vL5fEF5886dO6tz587lvubz+ZSTk6Nnn31Wd955pyTp9ddfV1JSkt59913dd999QakBAACYLeAwM2DAAC1YsEC7d+9W37591bt3byUmJoassN27d+vAgQPq2LGjfywhIUHXX3+9Pvnkk7OGGY/HI4/H4192u92SJK/XG5onehcz6TmihcNT4unByEYPwm4h6kEr39kBh5kZM2Zo6tSpWrRokebMmaPMzEx16dJFDz30kH7/+98H/UqmAwcOSJKSkpJKjSclJflfK092draysrLKjOfn5ys2NjaoNUqSipIqXwcXr6NhMIeLHoxs9CDsFqIeLCwsDHhdS1czuVwu9ezZUz179tSePXuUm5urxx57TKdPn9ZXX30VmrBgUWZmpjIyMvzLbrdbqampatiwoeLj44P/hocXBn+fMEftxnZXQA9GOnoQdgtRD5acWQnEOT01W5KioqLkcDjk8/lCcvqmbt26kqSDBw8qOTnZP37w4EH/fW7K43K55HK5yow7nU45nc6g16mo4MwfgqFC0VNW0YORjR6E3ULUg1a+sy1dmu3xeLRgwQLdeuutuuKKK/Tll19q+vTp2rt3b9CPyqSlpalu3bpauXKlf8ztduuzzz7T7373u6C+FwAAMFfAR2Yee+wxvfHGG0pNTVW/fv20YMEC1a5d+7zevLCwULt27fIv7969W1u2bFFiYqLq16+voUOHavz48WrcuLH/0uyUlBR17979vN4XAABcPAIOM7NmzVL9+vX1m9/8RmvWrNGaNWvKXW/RokUBv/nGjRvVoUMH/3LJXJc+ffooNzdXw4cPV1FRkR555BH9+OOPatu2rfLy8rjHDAAA8As4zDzwwANBv2Kpffv2Fd6zxuFwaNy4cRo3blxQ3xcAAFw8LN00DwAAINzY+mwmAACA80WYAQAARiPMAAAAoxFmAACA0QIKM1dffbWOHj0qSRo3bpyOHz8e0qIAAAACFVCY2bZtm4qKiiRJWVlZlh7+BAAAEEoBXZrdqlUr9e3bV23btpXP59Of//znsz6+YPTo0UEtEAAAoCIBhZnc3FyNGTNGS5culcPh0Pvvv68qVcpu6nA4CDMAAOCCCijMpKen64033pB05mnZK1eu1KWXXhrSwgAAAAIR8B2ASxQXF4eiDgAAgHNiOcxIUn5+vnJycrRt2zZJUpMmTTRkyBA1bNgwqMUBAABUxvJ9ZpYvX64mTZpo/fr1atGihVq0aKHPPvtMTZs21YoVK0JRIwAAwFlZPjIzcuRIDRs2TJMmTSozPmLECN16661BKw4AAKAylo/MbNu2TQ899FCZ8X79+unrr78OSlEAAACBshxm6tSpoy1btpQZ37JlC1c4AQCAC87yaaaHH35YjzzyiL755hvdcMMNkqSPPvpIkydPVkZGRtALBAAAqIjlMDNq1CjFxcXphRdeUGZmpiQpJSVFY8eO1eDBg4NeIAAAQEUshxmHw6Fhw4Zp2LBh+umnnyRJcXFxQS8MAAAgEOd0n5kShBgAAGA3yxOAAQAAwglhBgAAGI0wAwAAjGYpzJw6dUq33HKLdu7cGap6AAAALLEUZqpWraovvvgiVLUAAABYZvk0U+/evfXaa6+FohYAAADLLF+affr0ac2ZM0cffPCBrrnmGsXExJR6ferUqUErDgAAoDKWw8zWrVt19dVXS5J27NhR6jWHwxGcqgAAAAJkOcysWrUqFHUAAACck3O+NHvXrl1avny5Tpw4IUny+XxBKwoAACBQlsPMDz/8oFtuuUVXXHGFbr/9du3fv1+S9NBDD+mJJ54IeoEAAAAVsRxmhg0bpqpVq2rv3r2qUaOGf/zee+9VXl5eUIvzer0aNWqU0tLSFB0drYYNG+q5557jKBAAAPCzPGfm//7v/7R8+XJddtllpcYbN26sPXv2BK0wSZo8ebJmzpypefPmqWnTptq4caP69u2rhIQEDR48OKjvBQAAzGQ5zBQVFZU6IlPiyJEjcrlcQSmqxMcff6w777xTXbp0kSRdfvnlWrBggdavXx/U9wEAAOayHGZuuukmvf7663ruuecknbkcu7i4WFOmTFGHDh2CWtwNN9ygV155RTt27NAVV1yhzz//XOvWravwXjYej0cej8e/7Ha7JZ05ZeX1eoNanySpmMvRI1ooesoqejCy0YOwW4h60Mp3tuUwM2XKFN1yyy3auHGjTp48qeHDh+urr77SkSNH9NFHH1ndXYVGjhwpt9utK6+8Uk6nU16vVxMmTFCvXr3Ouk12draysrLKjOfn5ys2Njao9UmSipKCv0+Y42gYPKeMHoxs9CDsFqIeLCwsDHhdy2GmWbNm2rFjh6ZPn664uDgVFhaqR48eGjhwoJKTk63urkJvvvmm/vrXv2r+/Plq2rSptmzZoqFDhyolJUV9+vQpd5vMzExlZGT4l91ut1JTU9WwYUPFx8cHtT5J0uGFwd8nzFG7sd0V0IORjh6E3ULUgyVnVgJhOcxIUkJCgp555plz2dSSp556SiNHjtR9990nSWrevLn27Nmj7Ozss4YZl8tV7twdp9Mpp9MZ/CKjuLIqooWip6yiByMbPQi7hagHrXxnn1OYOXr0qF577TVt27ZNktSkSRP17dtXiYmJ57K7szp+/LiiokpfPe50OlVcXBzU9wEAAOayfJ+ZtWvX6vLLL9dLL72ko0eP6ujRo3rppZeUlpamtWvXBrW4rl27asKECVq2bJm+/fZbLV68WFOnTtVdd90V1PcBAADmsnxkZuDAgbr33ns1c+ZM/yEgr9erxx57TAMHDtSXX34ZtOKmTZumUaNG6bHHHlNBQYFSUlLUv39/jR49OmjvAQAAzGY5zOzatUtvv/12qXNZTqdTGRkZev3114NaXFxcnHJycpSTkxPU/QIAgIuH5dNMV199tX+uzC9t27ZNLVu2DEpRAAAAgQroyMwXX3zh/3nw4MEaMmSIdu3apd/+9reSpE8//VQzZszQpEmTQlMlAADAWQQUZlq1aiWHw1HqAY/Dhw8vs94f//hH3XvvvcGrDgAAoBIBhZndu3eHug4AAIBzElCYadCgQajrAAAAOCfndNO8ffv2ad26dSooKChzA7vBgwcHpTAAAIBAWA4zubm56t+/v6pVq6ZatWrJ4fjP01IdDgdhBgAAXFCWw8yoUaM0evRoZWZmlnnUAAAAwIVmOY0cP35c9913H0EGAACEBcuJ5KGHHtJbb70ViloAAAAss3yaKTs7W3fccYfy8vLUvHlzVa1atdTrU6dODVpxAAAAlTmnMLN8+XKlp6dLUpkJwAAAABeS5TDzwgsvaM6cOXrwwQdDUA4AAIA1lufMuFwu3XjjjaGoBQAAwDLLYWbIkCGaNm1aKGoBAACwzPJppvXr1+vDDz/U0qVL1bRp0zITgBctWhS04gAAACpjOcxccskl6tGjRyhqAQAAsMxymJk7d24o6gAAADgn3MYXAAAYzfKRmbS0tArvJ/PNN9+cV0EAAABWWA4zQ4cOLbV86tQp/fOf/1ReXp6eeuqpYNUFAAAQEMthZsiQIeWOz5gxQxs3bjzvggAAAKwI2pyZzp0765133gnW7gAAAAIStDDz9ttvKzExMVi7AwAACIjl00ytW7cuNQHY5/PpwIEDOnTokF5++eWgFgcAAFAZy2Gme/fupZajoqJUp04dtW/fXldeeWWw6gIAAAiI5TAzZsyYUNQBAABwTrhpHgAAMFrAR2aioqIqvFmeJDkcDp0+ffq8iwIAAAhUwGFm8eLFZ33tk08+0UsvvaTi4uKgFAUAABCogMPMnXfeWWZs+/btGjlypN577z316tVL48aNC2pxAAAAlTmnOTP79u3Tww8/rObNm+v06dPasmWL5s2bpwYNGgS7Pn3//ffq3bu3atWqpejoaDVv3pw7DQMAAD9LVzMdO3ZMEydO1LRp09SqVSutXLlSN910U6hq09GjR3XjjTeqQ4cOev/991WnTh3t3LlTNWvWDNl7AgAAswQcZqZMmaLJkyerbt26WrBgQbmnnYJt8uTJSk1N1dy5c/1jaWlpIX9fAABgjoDDzMiRIxUdHa1GjRpp3rx5mjdvXrnrLVq0KGjFLVmyRJ06ddI999yjNWvWqF69enrsscf08MMPn3Ubj8cjj8fjX3a73ZIkr9crr9cbtNr8iiu+wgsXuVD0lFX0YGSjB2G3EPWgle/sgMPMAw88UOml2cH2zTffaObMmcrIyNDTTz+tDRs2aPDgwapWrZr69OlT7jbZ2dnKysoqM56fn6/Y2NjgF1mUFPx9whxHd9pdAT0Y6ehB2C1EPVhYWBjwug6fz+cLSRVBUK1aNV177bX6+OOP/WODBw/Whg0b9Mknn5S7TXlHZlJTU3XkyBHFx8cHv8jD44O/T5ij9rN2V0APRjp6EHYLUQ+63W4lJibq2LFjlX5/W36cwYWUnJysJk2alBq76qqr9M4775x1G5fLJZfLVWbc6XTK6XQGvUZFhW0WxIUQip6yih6MbPQg7BaiHrTynR3WjzO48cYbtX379lJjO3bsCMkl4AAAwExhHWaGDRumTz/9VBMnTtSuXbs0f/58vfLKKxo4cKDdpQEAgDAR1mGmTZs2Wrx4sRYsWKBmzZrpueeeU05Ojnr16mV3aQAAIEyE9ZwZSbrjjjt0xx132F0GAAAIU2F9ZAYAAKAyhBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDTCDAAAMBphBgAAGI0wAwAAjEaYAQAARiPMAAAAoxFmAACA0QgzAADAaIQZAABgNMIMAAAwGmEGAAAYjTADAACMRpgBAABGI8wAAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMZFWYmTZokh8OhoUOH2l0KAAAIE8aEmQ0bNmj27Nlq0aKF3aUAAIAwYkSYKSwsVK9evfTqq6+qZs2adpcDAADCSBW7CwjEwIED1aVLF3Xs2FHjx4+vcF2PxyOPx+NfdrvdkiSv1yuv1xv84oodwd8nzBGKnrKKHoxs9CDsFqIetPKdHfZh5o033tDmzZu1YcOGgNbPzs5WVlZWmfH8/HzFxsYGuzypKCn4+4Q5ju60uwJ6MNLRg7BbiHqwsLAw4HXDOsx89913GjJkiFasWKHq1asHtE1mZqYyMjL8y263W6mpqWrYsKHi4+ODX+ThhcHfJ8xRu7HdFdCDkY4ehN1C1IMlZ1YCEdZhZtOmTSooKNDVV1/tH/N6vVq7dq2mT58uj8cjp9NZahuXyyWXy1VmX06ns8y6QRHlC/4+YY5Q9JRV9GBkowdhtxD1oJXv7LAOM7fccou+/PLLUmN9+/bVlVdeqREjRoQmnAAAAKOEdZiJi4tTs2bNSo3FxMSoVq1aZcYBAEBkMuLSbAAAgLMJ6yMz5Vm9erXdJQAAgDDCkRkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDTCDAAAMBphBgAAGI0wAwAAjEaYAQAARiPMAAAAoxFmAACA0QgzAADAaIQZAABgNMIMAAAwGmEGAAAYjTADAACMRpgBAABGI8wAAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNHCPsxkZ2erTZs2iouL06WXXqru3btr+/btdpcFAADCRNiHmTVr1mjgwIH69NNPtWLFCp06dUq///3vVVRUZHdpAAAgDFSxu4DK5OXllVrOzc3VpZdeqk2bNunmm2+2qSoAABAuwj7M/NqxY8ckSYmJieW+7vF45PF4/Mtut1uS5PV65fV6g19QsSP4+4Q5QtFTVtGDkY0ehN1C1INWvrONCjPFxcUaOnSobrzxRjVr1qzcdbKzs5WVlVVmPD8/X7GxscEvqigp+PuEOY7utLsCejDS0YOwW4h6sLCwMOB1HT6fzxeSKkJgwIABev/997Vu3Tpddtll5a5T3pGZ1NRUHTlyRPHx8cEv6vD44O8T5qj9rN0V0IORjh6E3ULUg263W4mJiTp27Fil39/GHJkZNGiQli5dqrVr1541yEiSy+WSy+UqM+50OuV0OoNfWJQxWRChEIqesooejGz0IOwWoh608p0d9mHG5/Pp8ccf1+LFi7V69WqlpaXZXRIAAAgjYR9mBg4cqPnz5+tvf/ub4uLidODAAUlSQkKCoqOjba4OAADYLezvMzNz5kwdO3ZM7du3V3Jysv+/hQsX2l0aAAAIA2F/ZMag+ckAAMAGYX9kBgAAoCKEGQAAYDTCDAAAMBphBgAAGI0wAwAAjEaYAQAARiPMAAAAoxFmAACA0QgzAADAaIQZAABgNMIMAAAwGmEGAAAYjTADAACMRpgBAABGI8wAAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKMRZgAAgNEIMwAAwGiEGQAAYDTCDAAAMBphBgAAGI0wAwAAjEaYAQAARiPMAAAAoxkRZmbMmKHLL79c1atX1/XXX6/169fbXRIAAAgTYR9mFi5cqIyMDI0ZM0abN29Wy5Yt1alTJxUUFNhdGgAACANhH2amTp2qhx9+WH379lWTJk00a9Ys1ahRQ3PmzLG7NAAAEAaq2F1ARU6ePKlNmzYpMzPTPxYVFaWOHTvqk08+KXcbj8cjj8fjXz527Jgk6ejRo/J6vcEv8pin8nVw8apy1O4K6MFIRw/CbiHqQbfbLUny+XyVlxCSCoLk8OHD8nq9SkpKKjWelJSkf/3rX+Vuk52draysrDLjl19+eShKRMSbZHcBiHj0IOwW2h786aeflJCQUOE6YR1mzkVmZqYyMjL8y8XFxTpy5Ihq1aolh8NhY2UXH7fbrdTUVH333XeKj4+3uxxEIHoQdqMHQ8fn8+mnn35SSkpKpeuGdZipXbu2nE6nDh48WGr84MGDqlu3brnbuFwuuVyuUmOXXHJJqEqEpPj4eP4nhq3oQdiNHgyNyo7IlAjrCcDVqlXTNddco5UrV/rHiouLtXLlSv3ud7+zsTIAABAuwvrIjCRlZGSoT58+uvbaa3XdddcpJydHRUVF6tu3r92lAQCAMBD2Yebee+/VoUOHNHr0aB04cECtWrVSXl5emUnBuPBcLpfGjBlT5rQecKHQg7AbPRgeHL5ArnkCAAAIU2E9ZwYAAKAyhBkAAGA0wgwAADAaYQYh9+CDD6p79+52lwEAuEgRZi5iDz74oBwOhxwOh6pWraq0tDQNHz5cP//8s92lIcL8sherVaumRo0aady4cTp9+rTdpfmtXr1aDodDP/74o92lIMxU1L8lfVPyX3R0tJo2bapXXnnF7rIjSthfmo3zc9ttt2nu3Lk6deqUNm3apD59+sjhcGjy5Ml2l2bJqVOnVLVqVbvLwHko6UWPx6O///3vGjhwoKpWrVrqQbJAuDpb/5bcwHX79u2Kj4/XiRMn9N5772nAgAFq2LChbrnlFpsrjwwcmbnIuVwu1a1bV6mpqerevbs6duyoFStWSDpzN+Xs7GylpaUpOjpaLVu21Ntvv11q+yVLlqhx48aqXr26OnTooHnz5pX61+vYsWPVqlWrUtvk5ORU+GDPvLw8tW3bVpdccolq1aqlO+64Q/n5+f7Xv/32WzkcDi1cuFDt2rVT9erV9de//jUofw/Yp6QXGzRooAEDBqhjx45asmSJPB6PnnzySdWrV08xMTG6/vrrtXr1av92P/zwg3r27Kl69eqpRo0aat68uRYsWFBq3+3bt9fgwYM1fPhwJSYmqm7duho7dqz/9ZKe2rJli3/sxx9/lMPh0OrVq/Xtt9+qQ4cOkqSaNWvK4XDowQcfDOFfA6Y5W/+WuPTSS1W3bl2lpaVp8ODBSktL0+bNm22sOLIQZiLI1q1b9fHHH6tatWqSzjxh/PXXX9esWbP01VdfadiwYerdu7fWrFkjSdq9e7fuvvtude/eXZ9//rn69++vZ5555rzrKCoqUkZGhjZu3KiVK1cqKipKd911l4qLi0utN3LkSA0ZMkTbtm1Tp06dzvt9EV6io6N18uRJDRo0SJ988oneeOMNffHFF7rnnnt02223aefOnZKkn3/+Wddcc42WLVumrVu36pFHHtH999+v9evXl9rfvHnzFBMTo88++0xTpkzRuHHj/MG9MqmpqXrnnXcknfkX9v79+/Xiiy8G9xfGRaWkf3/N5/MpLy9Pe/fu1fXXX29DZZGJ00wXuaVLlyo2NlanT5+Wx+NRVFSUpk+fLo/Ho4kTJ+qDDz7wHyb9zW9+o3Xr1mn27Nlq166dZs+erfT0dD3//POSpPT0dG3dulUTJkw4r5r+8Ic/lFqeM2eO6tSpo6+//lrNmjXzjw8dOlQ9evQ4r/dC+PH5fFq5cqWWL1+unj17au7cudq7d6//ybhPPvmk8vLyNHfuXE2cOFH16tXTk08+6d/+8ccf1/Lly/Xmm2/quuuu84+3aNFCY8aMkSQ1btxY06dP18qVK3XrrbdWWpPT6VRiYqKkM//C5uG0OJtf9u/jjz/uH7/sssskSR6PR8XFxRo3bpxuvvlmu8qMOISZi1yHDh00c+ZMFRUV6X/+539UpUoV/eEPf9BXX32l48ePl/mgP3nypFq3bi3pzL9Q27RpU+r1X355nKudO3dq9OjR+uyzz3T48GH/EZm9e/eWCjPXXnvteb8XwkdJsD516pSKi4v1xz/+UXfffbdyc3N1xRVXlFrX4/GoVq1akiSv16uJEyfqzTff1Pfff6+TJ0/K4/GoRo0apbZp0aJFqeXk5GQVFBSE9pdCxCivf8eOHasNGzZIkv7xj38oLi5OHo9H69ev16BBg5SYmKgBAwbYXHlkIMxc5GJiYtSoUSNJZ46AtGzZUq+99po/NCxbtkz16tUrtY2VZ4xERUXp10/EOHXqVIXbdO3aVQ0aNNCrr76qlJQUFRcXq1mzZmUO2cbExARcB8JfSbCuVq2aUlJSVKVKFS1cuFBOp1ObNm2S0+kstX5sbKwk6fnnn9eLL76onJwcNW/eXDExMRo6dGiZfvn1BHGHw+EPylFRZ86o/7JXK+tT4JfK699fSktL8x/Ra9q0qT777DNNmDCBMHOBEGYiSFRUlJ5++mllZGRox44dcrlc2rt3r9q1a1fu+unp6fr73/9eaqzkXyEl6tSpowMHDsjn88nhcEhSqUmWv/bDDz9o+/btevXVV3XTTTdJktatW3cevxVM8ctgXaJ169byer0qKCjw98OvffTRR7rzzjvVu3dvSWcmru/YsUNNmjQJ+L3r1KkjSdq/f7//yOOv+7RkLpnX6w14v4gc5fVvRZxOp06cOBHCivBLTACOMPfcc4+cTqdmz56tJ598UsOGDdO8efOUn5+vzZs3a9q0aZo3b54kqX///vrXv/6lESNGaMeOHXrzzTeVm5srSf7g0r59ex06dEhTpkxRfn6+ZsyYoffff/+s71+zZk3VqlVLr7zyinbt2qUPP/xQGRkZIf+9EZ6uuOIK9erVSw888IAWLVqk3bt3a/369crOztayZcsknZn/smLFCn388cfatm2b+vfvr4MHD1p6n+joaP32t7/VpEmTtG3bNq1Zs0bPPvtsqXUaNGggh8OhpUuX6tChQyosLAza74mLX0FBgQ4cOKA9e/borbfe0l/+8hfdeeeddpcVMQgzEaZKlSoaNGiQpkyZoszMTI0aNUrZ2dm66qqrdNttt2nZsmVKS0uTdOaw6dtvv61FixapRYsWmjlzpv9qppJTUVdddZVefvllzZgxQy1bttT69etLTdb8taioKL3xxhvatGmTmjVrpmHDhvknGCMyzZ07Vw888ICeeOIJpaenq3v37tqwYYPq168vSXr22Wd19dVXq1OnTmrfvr3q1q17TneUnjNnjk6fPq1rrrlGQ4cO1fjx40u9Xq9ePWVlZWnkyJFKSkrSoEGDgvHrIUKkp6crOTlZjRo10ogRI9S/f39NmzbN7rIihsP36wkPQAUmTJigWbNm6bvvvrO7FAAAJDFnBpV4+eWX1aZNG9WqVUsfffSRnn/+ef7FCgAIK4QZVGjnzp0aP368jhw5ovr16+uJJ57g9vMAgLDCaSYAAGA0JgADAACjEWYAAIDRCDMAAMBohBkAAGA0wgwAADAaYQYAABiNMAMAAIxGmAEAAEYjzAAAAKP9f1IQbxiTsd5KAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqq0lEQVR4nO3de3yMd6LH8e9kSEQuiOOSaKRZl7QSl7rWYnHotqpF76fLQetlbUtd4ppuK4lF0LOqWlWcSvSc09Iq3aptnJaiirrmtEiRlOUQ1UUTSdMRkzl/OGZNEzUPM55H5vN+vfJ6ZZ55Lt/k9WO+eZ7fPGNzuVwuAQAAWFCQ2QEAAACuhqICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsq5rZAW5EeXm5Tp48qYiICNlsNrPjAAAAL7hcLp0/f14xMTEKCvrlcya3dFE5efKkYmNjzY4BAACuw/Hjx3Xbbbf94jq3dFGJiIiQdOkHjYyMNDkNAADwRlFRkWJjY92v47/kli4qly/3REZGUlQAALjFeDNtg8m0AADAsigqAADAsigqAADAskwvKidOnNCgQYNUt25dhYaGqmXLltq1a5fZsQAAgAWYOpn23Llz6tKli3r27KmPP/5Y9erV0+HDh1WnTh0zYwEAAIswtajMnj1bsbGxyszMdC+Lj483MREAALASU4vKhx9+qHvvvVePPfaYNm3apEaNGunZZ5/V8OHDK13f4XDI4XC4HxcVFUmSnE6nnE7nTckMAABujJHXbFOLyrfffquFCxcqOTlZzz//vHbu3KnRo0crODhYQ4YMqbB+RkaG0tPTKyzPz89XeHj4zYgMAABuUHFxsdfr2lwul8uPWX5RcHCw2rdvr61bt7qXjR49Wjt37tS2bdsqrF/ZGZXY2FidPXuWG74BAHCLKCoqUlRUlAoLC6/5+m3qGZXo6Gi1aNHCY9mdd96p999/v9L1Q0JCFBISUmG53W6X3W73S0YAAOBbRl6zTX17cpcuXXTw4EGPZYcOHVJcXJxJiQAAgJWYWlTGjRun7du3a+bMmcrLy9Pbb7+txYsXa+TIkWbGAgAAFmFqUenQoYNWr16td955R0lJSfrTn/6kefPmaeDAgWbGAgAAFmHqZNobVVRUpFq1ank1GQeANRUUFKigoMDwdtHR0YqOjvZDIgD+ZuT129TJtACwaNGiSm87cC2pqalKS0vzfSAAlkJRAWCqESNGqF+/fh7LSktL1bVrV0nSli1bFBoaWmE7zqYAgYGiAsBUlV3CKSkpcX/fpk0bhYWF3exYACzC9E9PBgAAuBqKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsCyKCgAAsKxqZgcA8AvS0sxOYI4LF/7x/YwZUnCweVnMFqhjAPh/pp5RSUtLk81m8/i64447zIwEAAAsxPQzKomJifr000/dj6tVMz0SAACwCNNbQbVq1dSwYUOzYwAAAAsyvagcPnxYMTExqlGjhjp37qyMjAw1bty40nUdDoccDof7cVFRkSTJ6XTK6XTelLzATWWzmZ3AFM4rfm6nzebxOODwfxuqICOv2aYWlU6dOikrK0sJCQkqKChQenq6unXrpn379ikiIqLC+hkZGUpPT6+wPD8/X+Hh4TcjMnBzNWhgdgJT/HjFHyR59eurZkiIiWlMdviw2QkAnysuLvZ6XZvL5XL5MYshP/zwg+Li4jR37lwNGzaswvOVnVGJjY3V2bNnFRkZeTOjAjfH9OlmJzBFyYULqjVzpiSp8PnnFRbI7/p54QWzEwA+V1RUpKioKBUWFl7z9dv0Sz9Xql27tpo3b668vLxKnw8JCVFIJX9Z2e122e12f8cDbj7r/B1xU9mv+LntLpfH44DD/22ogoy8Zlvqhm/FxcXKz89XdHS02VEAAIAFmFpUJkyYoE2bNuno0aPaunWrHnroIdntdj355JNmxgIAABZh6qWf//3f/9WTTz6pM2fOqF69euratau2b9+uevXqmRkLAABYhKlFZfny5WYeHgAAWJyl5qgAAABciaICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAsi6ICAAAs64aLSlFRkT744APl5ub6Ig8AAICb4aLy+OOP67XXXpMklZaWqn379nr88cfVqlUrvf/++z4PCAAAApfhorJ582Z169ZNkrR69Wq5XC798MMPmj9/vqZPn+7zgAAAIHBVM7pBYWGhoqKiJEnZ2dl65JFHVLNmTfXt21cTJ070eUAAVVvB+fMqKC72WFZaVub+PufUKYVWr15hu+jwcEVHRPg9HwBzGS4qsbGx2rZtm6KiopSdna3ly5dLks6dO6caNWr4PCCAqm3R7t1K37Tpqs93zcysdHlq9+5K69HDT6kAWIXhojJ27FgNHDhQ4eHhiouLU4///49i8+bNatmypa/zAajiRrRrp34JCYa3iw4P90MaAFZjuKg8++yz6tixo44fP6577rlHQUGXprn86le/Yo4KAMOiIyK4hAPgqgwXFUlq37692rdv77Gsb9++PgkEAABwmddFZdq0aV6tN3Xq1OsOAwAAcCWvi0paWppiYmJUv359uVyuStex2WwUFQAA4DNeF5U+ffpow4YNat++vZ5++mk98MAD7vkpAAAA/uB101i7dq3y8/PVqVMnTZw4UY0aNdLkyZN18OBBf+YDAAABzNApkZiYGKWkpOjgwYNasWKFTp8+rQ4dOqhLly4qLS31V0YAABCgrutdP5LUoUMHHT16VAcOHNDevXtVVlam0NBQX2YDAAABzvAkk23btmn48OFq2LChXn31VQ0ZMkQnT55UZGSkP/IBAIAA5vUZlTlz5igrK0t///vfNXDgQH3++edq1aqVP7MBAIAA53VRmTJliho3bqzHH39cNptNWVlZla43d+5cX2UDAAABzuui8pvf/EY2m0379++/6jo2m+26g8yaNUspKSkaM2aM5s2bd937AQAAVYfXRWXjxo1+C7Fz504tWrSIS0kAAMCD6XdsKy4u1sCBA7VkyRLVqVPH7DgAAMBCDL892el0KisrS+vXr9fp06dVXl7u8fyGDRsM7W/kyJHq27evevfufc1PX3Y4HHI4HO7HRUVF7kxOp9PQcYFbwg1cTkUVwf9tqIKMvGYbLipjxoxRVlaW+vbtq6SkpBual7J8+XLt2bNHO3fu9Gr9jIwMpaenV1ien5+v8PDw684BWFaDBmYngNkOHzY7AeBzxcXFXq9ruKgsX75c7777ru6//36jm3o4fvy4xowZo08++UQ1atTwapuUlBQlJye7HxcVFSk2NlZNmjThPi6omlasMDsBzNasmdkJAJ+7fEXEG4aLSnBwsJo2bWp0swp2796t06dPq23btu5lTqdTmzdv1muvvSaHwyG73e6xTUhIiEJCQirsy263V1gXqBKu8knlCCD834YqyMhrtuHJtOPHj9crr7wi1w3+B9qrVy99/fXXysnJcX+1b99eAwcOVE5ODsUDAAAYP6OyZcsWffbZZ/r444+VmJio6tWrezy/atUqr/YTERGhpKQkj2VhYWGqW7duheUAACAwGS4qtWvX1kMPPeSPLAAAAB4MF5XMzEx/5JDk35vKAQCAW4/pN3wDAAC4Gq/PqPzzP/+zV+sZveEbAADA1Rj6rJ+4uDj17du3wgRaAAAAf/C6qMyePVuZmZl67733NHDgQD399NO8OwcAAPiV13NUJk6cqAMHDuiDDz7Q+fPn1aVLF3Xs2FFvvPGGoTvMAQAAeMvwZNrOnTtryZIlKigo0MiRI7V06VLFxMRQVgAAgM9d97t+9uzZo02bNik3N1dJSUnMWwEAAD5nqKicPHlSM2fOVPPmzfXoo48qKipKX375pbZv367Q0FB/ZQQAAAHK68m0999/vz777DP99re/1UsvvaS+ffuqWjXD94sDAADwms3l5acLBgUFKTo6WvXr15fNZrvqenv27PFZuGspKipSrVq1VFhYqMjIyJt2XOCmSUszOwHMxhhAFWTk9dvrUyKpqak3HAwAAMAIigoAALAsPusHAABYFkUFAABYFkUFAABYFkUFAABY1g0VlZ9++slXOQAAACowXFTKy8v1pz/9SY0aNVJ4eLi+/fZbSdKLL76oN9980+cBAQBA4DJcVKZPn66srCzNmTNHwcHB7uVJSUn693//d5+GAwAAgc1wUXnrrbe0ePFiDRw4UHa73b28devW+uabb3waDgAABDbDReXEiRNq2rRpheXl5eUqKyvzSSgAAADpOopKixYt9Pnnn1dYvnLlSt11110+CQUAACAZuIX+ZVOnTtWQIUN04sQJlZeXa9WqVTp48KDeeustffTRR/7ICAAAApThMyr9+/fXmjVr9OmnnyosLExTp05Vbm6u1qxZo3vuuccfGQEAQIAyfEZFkrp166ZPPvnE11kAAAA8cGdaAABgWYbPqAQFBclms131eafTeUOBAAAALjNcVFavXu3xuKysTHv37tWyZcuUnp7us2AAAACGi0r//v0rLHv00UeVmJioFStWaNiwYT4JBgAA4LM5KnfffbfWr1/vq90BAAD4pqiUlpZq/vz5atSokS92BwAAIOk6Lv3UqVPHYzKty+XS+fPnVbNmTf3nf/6nT8MBAIDAZriovPzyyx5FJSgoSPXq1VOnTp1Up04dn4YDAACBzXBRGTp0qB9iAAAAVORVUfnqq6+83mGrVq2uOwwAAMCVvCoqbdq0kc1mk8vl+sX1bDYbN3wDAAA+41VROXLkiL9zAAAAVOBVUYmLi/N3DgAAgAoMT6Y9c+aM6tatK0k6fvy4lixZotLSUvXr10/dunXzeUAAABC4vL7h29dff63bb79d9evX1x133KGcnBx16NBBL7/8shYvXqyePXvqgw8+8GNUAAAQaLwuKpMmTVLLli21efNm9ejRQw888ID69u2rwsJCnTt3TiNGjNCsWbP8mRUAAAQYry/97Ny5Uxs2bFCrVq3UunVrLV68WM8++6yCgi51neeee053332334ICAIDA4/UZlbNnz6phw4aSpPDwcIWFhXncibZOnTo6f/687xMCAICAZehDCa+8dX5ljwEAAHzJ0Lt+hg4dqpCQEEnSTz/9pD/84Q8KCwuTJDkcDt+nAwAAAc3rojJkyBCPx4MGDaqwzuDBg288EQAAwP/zuqhkZmb6MwcAAEAFhm/4hqqloKBABQUFhreLjo5WdHS0HxIBAPAPphaVhQsXauHChTp69KgkKTExUVOnTlWfPn3MjBVQFi1apPT0dMPbpaamKi0tzfeBAAC4gqlF5bbbbtOsWbPUrFkzuVwuLVu2TP3799fevXuVmJhoZrSAMWLECPXr189jWWlpqbp27SpJ2rJli0JDQytsx9kUAMDNYGpRefDBBz0ez5gxQwsXLtT27dspKjdJZZdwSkpK3N+3adPG/c4uAABuNq/uo9K2bVudO3dOkjRt2jT9+OOPPg/idDq1fPlylZSUqHPnzj7fPwAAuPV4dUYlNzdXJSUlqlOnjtLT0/WHP/xBNWvW9EmAr7/+Wp07d9ZPP/2k8PBwrV69Wi1atKh0XYfD4XG/lqKiIkmXSo7T6fRJHsjjd8nv1mTcVBH8+0MVZOR1xaui0qZNGz311FPq2rWrXC6X/u3f/k3h4eGVrjt16lSvDy5JCQkJysnJUWFhoVauXKkhQ4Zo06ZNlZaVjIyMSid+5ufnXzUPjLvyjFleXp7PSimuQ4MGZieA2Q4fNjsB4HPFxcVer2tzuVyua6108OBBpaamKj8/X3v27FGLFi1UrVrFjmOz2bRnzx5jaX+md+/eatKkiRYtWlThucrOqMTGxurs2bOKjIy8oePiH0pKSlSrVi1JUmFhIXNUzDR9utkJYLYXXjA7AeBzRUVFioqKUmFh4TVfv706o5KQkKDly5dLkoKCgrR+/XrVr1//xpNWory8/Kq34w8JCXHfwv9KdrtddrvdL3kC0ZW/S363Jrv23xGo6vj3hyrIyOuK4Xf9lJeXG93kqlJSUtSnTx81btxY58+f19tvv62NGzdq3bp1PjsGAAC4dV3X25Pz8/M1b9485ebmSpJatGihMWPGqEmTJob2c/r0aQ0ePFgFBQWqVauWWrVqpXXr1umee+65nlgAAKCKMVxU1q1bp379+qlNmzbq0qWLJOmLL75QYmKi1qxZY6hkvPnmm0YPDwAAAojhojJlyhSNGzdOs2bNqrB88uTJnA0BAAA+49UN366Um5urYcOGVVj+9NNP68CBAz4JBQAAIF1HUalXr55ycnIqLM/JyfHbO4EAAEBgMnzpZ/jw4fr973+vb7/9Vr/+9a8lXZqjMnv2bCUnJ/s8IAAACFyGi8qLL76oiIgI/fnPf1ZKSookKSYmRmlpaRo9erTPAwIAgMBluKjYbDaNGzdO48aN0/nz5yVJERERPg8GAABwXfdRuYyCAgAA/MnwZFoAAICbhaICAAAsi6ICAAAsy1BRKSsrU69evXT48GF/5QEAAHAzVFSqV6+ur776yl9ZAAAAPBi+9DNo0CA+TBAAANwUht+efPHiRS1dulSffvqp2rVrp7CwMI/n586d67NwAAAgsBkuKvv27VPbtm0lSYcOHfJ4zmaz+SYVAACArqOofPbZZ/7IAQAAUMF1vz05Ly9P69atU2lpqSTJ5XL5LBQAAIB0HUXlzJkz6tWrl5o3b677779fBQUFkqRhw4Zp/PjxPg8IAAACl+GiMm7cOFWvXl3Hjh1TzZo13cufeOIJZWdn+zQcAAAIbIbnqPz3f/+31q1bp9tuu81jebNmzfS3v/3NZ8EAAAAMn1EpKSnxOJNy2dmzZxUSEuKTUAAAANJ1FJVu3brprbfecj+22WwqLy/XnDlz1LNnT5+GAwAAgc3wpZ85c+aoV69e2rVrly5cuKBJkyZp//79Onv2rL744gt/ZAQAAAHK8BmVpKQkHTp0SF27dlX//v1VUlKihx9+WHv37lWTJk38kREAAAQow2dUJKlWrVr64x//6OssAAAAHq6rqJw7d05vvvmmcnNzJUktWrTQU089paioKJ+GAwAAgc3wpZ/Nmzfr9ttv1/z583Xu3DmdO3dO8+fPV3x8vDZv3uyPjAAAIEAZPqMycuRIPfHEE1q4cKHsdrskyel06tlnn9XIkSP19ddf+zwkAAAITIbPqOTl5Wn8+PHukiJJdrtdycnJysvL82k4AAAQ2AwXlbZt27rnplwpNzdXrVu39kkoAAAAyctLP1999ZX7+9GjR2vMmDHKy8vT3XffLUnavn27FixYoFmzZvknpVnS0sxOYI4LF/7x/YwZUnCweVnMFqhjAAAswqui0qZNG9lsNrlcLveySZMmVVjvd7/7nZ544gnfpQMAAAHNq6Jy5MgRf+cAAACowKuiEhcX5+8cAAAAFVzXDd9OnjypLVu26PTp0yovL/d4bvTo0T4JBgAAYLioZGVlacSIEQoODlbdunVls9ncz9lsNooKAOCWUlBQoIKCAsPbRUdHKzo62g+JcCXDReXFF1/U1KlTlZKSoqAgw+9uBgDAUhYtWqT09HTD26WmpiqNdwb6neGi8uOPP+pf/uVfKCkAgCphxIgR6tevn8ey0tJSde3aVZK0ZcsWhYaGVtiOsyk3h+GiMmzYML333nuaMmWKP/IAAHBTVXYJp6SkxP19mzZtFBYWdrNj4f8ZLioZGRl64IEHlJ2drZYtW6p69eoez8+dO9dn4QAAQGC7rqKybt06JSQkSFKFybQAAAC+Yrio/PnPf9bSpUs1dOhQP8QBAAD4B8MzYkNCQtSlSxd/ZAEAAPBguKiMGTNGr776qj+yAAAAeDB86WfHjh3asGGDPvroIyUmJlaYTLtq1SqfhQMAAIHNcFGpXbu2Hn74YX9kAQAA8GC4qGRmZvojBwAAQAXcXhYAAFiW4TMq8fHxv3i/lG+//faGAgEAAFxmuKiMHTvW43FZWZn27t2r7OxsTZw40dC+MjIytGrVKn3zzTcKDQ3Vr3/9a82ePdt9MzkAABDYDBeVMWPGVLp8wYIF2rVrl6F9bdq0SSNHjlSHDh108eJFPf/88/rtb3+rAwcO8LkKAADAd3NU+vTpo/fff9/QNtnZ2Ro6dKgSExPVunVrZWVl6dixY9q9e7evYgEAgFuYz4rKypUrFRUVdUP7KCwslKQb3g8AAKgaDF/6ueuuuzwm07pcLp06dUrff/+9Xn/99esOUl5errFjx6pLly5KSkqqdB2HwyGHw+F+XFRUJElyOp1yOp3XfeyrCtAPWXRe8XM7bTaPxwHHH+PKiED+3eMSs8dggLryNcVvrzEBzMjv03BRGTBggMfjoKAg1atXTz169NAdd9xhdHduI0eO1L59+7Rly5arrpORkaH09PQKy/Pz8xUeHn7dx76qBg18v89bwI9XlMG8+vVVMyTExDQmO3zY3OMH6BjEFcwegwHqxx9/dH+fl5enmjVrmpim6ikuLvZ6XZvL5XL5MYtXRo0apb/85S/avHmz4uPjr7peZWdUYmNjdfbsWUVGRvo+2PTpvt/nLaDkwgXVmjlTklT4/PMKCw42OZGJXnjB3OMH6BjEFcwegwGqpKREtWrVknRpWgJv8PCtoqIiRUVFqbCw8Jqv34bPqPiSy+XSc889p9WrV2vjxo2/WFKkS5/cHFLJX/d2u112u90fAX2/z1uA/Yqf2+5yeTwOOP4YV0YE8u8el5g9BgPUla8pfnuNCWBGfp9eF5WgoKBfvNGbJNlsNl28eNHrg48cOVJvv/22/vKXvygiIkKnTp2SJNWqVUuhoaFe7wcAAFRNXheV1atXX/W5bdu2af78+SovLzd08IULF0qSevTo4bE8MzNTQ4cONbQvAABQ9XhdVPr3719h2cGDBzVlyhStWbNGAwcO1LRp0wwd3ALTYwAAgIVd131UTp48qeHDh6tly5a6ePGicnJytGzZMsXFxfk6HwAACGCGikphYaEmT56spk2bav/+/Vq/fr3WrFlz1fueAAAA3AivL/3MmTNHs2fPVsOGDfXOO+9UeikIAADAl7wuKlOmTFFoaKiaNm2qZcuWadmyZZWut2rVKp+FAwAAgc3rojJ48OBrvj0ZAADAl7wuKllZWX6MAQAAUJHPPj0ZAADA1ygqAADAsigqAADAskz9UEKYr+D8eRX87OO2S8vK3N/nnDql0OrVK2wXHR6u6IgIv+cDAAQ2ikqAW7R7t9I3bbrq810zMytdntq9u9J+9hlNAKqgtDSzE5jjwoV/fD9jhhQcbF4Ws5k8BigqAW5Eu3bql5BgeLvo8HA/pAEAwBNFJcBFR0RwCQcAYFlMpgUAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZFUQEAAJZVzewAAACYqeD8eRUUF3ssKy0rc3+fc+qUQqtXr7BddHi4oiMi/J4v0FFUAAABbdHu3UrftOmqz3fNzKx0eWr37krr0cNPqXAZRQUAENBGtGunfgkJhreLDg/3Qxr8HEUFABDQoiMiuIRjYUymBQAAlmVqUdm8ebMefPBBxcTEyGaz6YMPPjAzDgAAsBhTi0pJSYlat26tBQsWmBkDAABYlKlzVPr06aM+ffqYGQEAAFjYLTWZ1uFwyOFwuB8XFRVJkpxOp5xOp+8PaLP5fp+4tfhjXBnBGARjEGbzwxg08pp9SxWVjIwMpaenV1ien5+vcH+8TaxBA9/vE7eWw4fNPT5jEIxBmM0PY7D4ZzfY+yW3VFFJSUlRcnKy+3FRUZFiY2PVpEkTRUZG+v6AK1b4fp+4tTRrZu7xGYNgDMJsfhiDl6+IeOOWKiohISEKCQmpsNxut8tut/v+gC6X7/eJW4s/xpURjEEwBmE2P4xBI6/Z3EcFAABYlqlnVIqLi5WXl+d+fOTIEeXk5CgqKkqNGzc2MRkAALACU4vKrl271LNnT/fjy/NPhgwZoqysLJNSAQAAqzC1qPTo0UMurn8CAICrYI4KAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLEsUlQULFuj2229XjRo11KlTJ+3YscPsSAAAwAJMLyorVqxQcnKyUlNTtWfPHrVu3Vr33nuvTp8+bXY0AABgMtOLyty5czV8+HA99dRTatGihd544w3VrFlTS5cuNTsaAAAwWTUzD37hwgXt3r1bKSkp7mVBQUHq3bu3tm3bVmF9h8Mhh8PhflxYWChJOnfunJxOp+8DXnEsBKhz58w9PmMQjEGYzQ9jsKioSJLkcrmuvbLLRCdOnHBJcm3dutVj+cSJE10dO3assH5qaqpLEl988cUXX3zxVQW+jh8/fs2uYOoZFaNSUlKUnJzsflxeXq6zZ8+qbt26stlsJiareoqKihQbG6vjx48rMjLS7DgIQIxBmI0x6D8ul0vnz59XTEzMNdc1taj80z/9k+x2u7777juP5d99950aNmxYYf2QkBCFhIR4LKtdu7Y/Iwa8yMhI/oHCVIxBmI0x6B+1atXyaj1TJ9MGBwerXbt2Wr9+vXtZeXm51q9fr86dO5uYDAAAWIHpl36Sk5M1ZMgQtW/fXh07dtS8efNUUlKip556yuxoAADAZKYXlSeeeELff/+9pk6dqlOnTqlNmzbKzs5WgwYNzI4W0EJCQpSamlrhUhtwszAGYTbGoDXYXC5v3hsEAABw85l+wzcAAICroagAAADLoqgAAADLoqjghgwdOlQDBgwwOwYAoIqiqNyihg4dKpvNJpvNpurVqys+Pl6TJk3STz/9ZHY0BKArx2NwcLCaNm2qadOm6eLFi2ZHc9u4caNsNpt++OEHs6PAYn5p/F4eN5e/QkNDlZiYqMWLF5sdO2CY/vZkXL/77rtPmZmZKisr0+7duzVkyBDZbDbNnj3b7GiGlJWVqXr16mbHwA26PB4dDof++te/auTIkapevbrHh44CVnW18Xv55qMHDx5UZGSkSktLtWbNGj3zzDNq0qSJevXqZXLyqo8zKrewkJAQNWzYULGxsRowYIB69+6tTz75RNKlO/xmZGQoPj5eoaGhat26tVauXOmx/YcffqhmzZqpRo0a6tmzp5YtW+bxF2daWpratGnjsc28efN0++23XzVTdna2unbtqtq1a6tu3bp64IEHlJ+f737+6NGjstlsWrFihbp3764aNWrov/7rv3zy+4C5Lo/HuLg4PfPMM+rdu7c+/PBDORwOTZgwQY0aNVJYWJg6deqkjRs3urc7c+aMnnzySTVq1Eg1a9ZUy5Yt9c4773jsu0ePHho9erQmTZqkqKgoNWzYUGlpae7nL4+rnJwc97IffvhBNptNGzdu1NGjR9WzZ09JUp06dWSz2TR06FA//jZwq7na+L2sfv36atiwoeLj4zV69GjFx8drz549JiYOHBSVKmLfvn3aunWrgoODJUkZGRl666239MYbb2j//v0aN26cBg0apE2bNkmSjhw5okcffVQDBgzQ//zP/2jEiBH64x//eMM5SkpKlJycrF27dmn9+vUKCgrSQw89pPLyco/1pkyZojFjxig3N1f33nvvDR8X1hMaGqoLFy5o1KhR2rZtm5YvX66vvvpKjz32mO677z4dPnxYkvTTTz+pXbt2Wrt2rfbt26ff//73+td//Vft2LHDY3/Lli1TWFiYvvzyS82ZM0fTpk1zF/NriY2N1fvvvy/p0l/GBQUFeuWVV3z7A6NKuTx+f87lcik7O1vHjh1Tp06dTEgWeLj0cwv76KOPFB4erosXL8rhcCgoKEivvfaaHA6HZs6cqU8//dR92vJXv/qVtmzZokWLFql79+5atGiREhIS9NJLL0mSEhIStG/fPs2YMeOGMj3yyCMej5cuXap69erpwIEDSkpKci8fO3asHn744Rs6FqzJ5XJp/fr1WrdunZ588kllZmbq2LFj7k9JnTBhgrKzs5WZmamZM2eqUaNGmjBhgnv75557TuvWrdO7776rjh07upe3atVKqampkqRmzZrptdde0/r163XPPfdcM5PdbldUVJSkS38Z82GmuJorx+9zzz3nXn7bbbdJkhwOh8rLyzVt2jT95je/MStmQKGo3MJ69uyphQsXqqSkRC+//LKqVaumRx55RPv379ePP/5Y4T/wCxcu6K677pJ06a/KDh06eDx/5YvC9Tp8+LCmTp2qL7/8Un//+9/dZ1KOHTvmUVTat29/w8eCtVwuzmVlZSovL9fvfvc7Pfroo8rKylLz5s091nU4HKpbt64kyel0aubMmXr33Xd14sQJXbhwQQ6HQzVr1vTYplWrVh6Po6Ojdfr0af/+UAgYlY3ftLQ07dy5U5L0+eefKyIiQg6HQzt27NCoUaMUFRWlZ555xuTkVR9F5RYWFhampk2bSrp05qJ169Z688033YVg7dq1atSokcc2Rj6zIigoSD//hIWysrJf3ObBBx9UXFyclixZopiYGJWXlyspKanCKdSwsDCvc+DWcLk4BwcHKyYmRtWqVdOKFStkt9u1e/du2e12j/XDw8MlSS+99JJeeeUVzZs3Ty1btlRYWJjGjh1bYcz8fMK1zWZzF+GgoEtXsa8cr9caq8CVKhu/V4qPj3efiUtMTNSXX36pGTNmUFRuAopKFREUFKTnn39eycnJOnTokEJCQnTs2DF179690vUTEhL017/+1WPZ5b8cLqtXr55OnToll8slm80mSR6TFX/uzJkzOnjwoJYsWaJu3bpJkrZs2XIDPxVuJVcW58vuuusuOZ1OnT592j0mfu6LL75Q//79NWjQIEmXJoIfOnRILVq08PrY9erVkyQVFBS4zxr+fKxenr/ldDq93i8CR2Xj95fY7XaVlpb6MREuYzJtFfLYY4/Jbrdr0aJFmjBhgsaNG6dly5YpPz9fe/bs0auvvqply5ZJkkaMGKFvvvlGkydP1qFDh/Tuu+8qKytLktylpEePHvr+++81Z84c5efna8GCBfr444+vevw6deqobt26Wrx4sfLy8rRhwwYlJyf7/eeGdTVv3lwDBw7U4MGDtWrVKh05ckQ7duxQRkaG1q5dK+nSfJNPPvlEW7duVW5urkaMGKHvvvvO0HFCQ0N19913a9asWcrNzdWmTZv0wgsveKwTFxcnm82mjz76SN9//72Ki4t99nOi6jt9+rROnTqlv/3tb3rvvff0H//xH+rfv7/ZsQICRaUKqVatmkaNGqU5c+YoJSVFL774ojIyMnTnnXfqvvvu09q1axUfHy/p0mnMlStXatWqVWrVqpUWLlzoftfP5ctDd955p15//XUtWLBArVu31o4dOzwmPf5cUFCQli9frt27dyspKUnjxo1zT9ZF4MrMzNTgwYM1fvx4JSQkaMCAAdq5c6caN24sSXrhhRfUtm1b3XvvverRo4caNmx4XXc7Xrp0qS5evKh27dpp7Nixmj59usfzjRo1Unp6uqZMmaIGDRpo1KhRvvjxECASEhIUHR2tpk2bavLkyRoxYoReffVVs2MFBJvr55MQELBmzJihN954Q8ePHzc7CgAAkpijEtBef/11dejQQXXr1tUXX3yhl156ib8yAQCWQlEJYIcPH9b06dN19uxZNW7cWOPHj+d25wAAS+HSDwAAsCwm0wIAAMuiqAAAAMuiqAAAAMuiqAAAAMuiqAAAAMuiqAAAAMuiqAAAAMuiqAAAAMuiqAAAAMv6P9WR9p9GUkFwAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Bar Graphs with standard deviation error lines \n",
    "categories = ['Regular',\"Peanut\",\"PB\"]\n",
    "num_means = [Num_Reg_mean,Num_Peanut_mean,Num_PB_mean]\n",
    "num_dev =  [Num_Reg_dev/2,Num_Peanut_dev/2,Num_PB_dev/2]\n",
    "plt.ylabel('Number of Mnms')\n",
    "plt.bar(categories, num_means, yerr=num_dev, capsize=5, color='gold', alpha=0.5)\n",
    "plt.grid(axis='y', linestyle='-', alpha=0.5)\n",
    "plt.show()\n",
    "blue_means =  [Blue_Reg_mean,Blue_Peanuts_mean,Blue_PB_mean]\n",
    "blue_dev = [Blue_Reg_dev/2,Blue_Peanuts_dev/2,Blue_PB_dev/2]\n",
    "plt.ylabel(\"Number of Blue MnMs\")\n",
    "plt.bar(categories, blue_means, yerr=blue_dev, capsize=5, color='red', alpha=0.5)\n",
    "plt.grid(axis='y', linestyle='-', alpha=0.5)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The bar graphs are here to give visual representation to the interperation that was printed out."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}