Lab5 geog 172. .ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geog 172: ANOVA, Correlation, Regression (Lab 5)Due on December 9, 2020before 9.30 am\n",
    "\n",
    "By Jack Holbrook, Wednesday section, Greenbeg, 12/8/2020. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](img/p1lab5.PNG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question 1:\n",
    "\n",
    "    Assume that an analysis of variance is conducted for a study where there are N = 50 observations and k = 5 categories. Fill in the blanks in the following ANOVA table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Source:     \tSS  \tDF      \tMS      \tF\n",
    "# group: \tSST \tk−1 \tSST/(k−1) \tMST/MSE\n",
    "# Error:\tSSE \tN−k \tSSE/(N−k)\n",
    "# Total: SST + SSE, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "group SS = 465.2\n",
      "group SS = 116.3\n",
      "group DF = 4\n",
      "error sum of squares = 2000\n",
      "error DF = 45\n",
      "error DF = 44.44444444444444\n",
      "group F stat = 2.61675\n",
      "total SS = 2465.2\n",
      "totsl DF = 49\n"
     ]
    }
   ],
   "source": [
    "#given:\n",
    "N = 50\n",
    "k = 5\n",
    "\n",
    "#group:\n",
    "groupDF = k-1\n",
    "groupMS = 116.3\n",
    "groupSS = groupMS * (k-1)\n",
    "groupMS = groupSS / (k-1)  # is simply the inverse\n",
    "print(\"group SS = \" +str(groupSS))\n",
    "print(\"group SS = \" +str(groupMS))\n",
    "print(\"group DF = \" +str(groupDF))\n",
    "\n",
    "#error:\n",
    "errorDF = N-k\n",
    "ErrorSS = 2000\n",
    "errorMs = ErrorSS/ errorDF\n",
    "Fstat = groupMS / errorMs\n",
    "print(\"error sum of squares = \" + str(2000))\n",
    "print(\"error DF = \" +str(errorDF))\n",
    "print(\"error DF = \" +str(errorMs))\n",
    "print(\"group F stat = \" +str(Fstat))\n",
    "\n",
    "#total:\n",
    "print(\"total SS = \" +str(ErrorSS + groupSS))\n",
    "print(\"totsl DF = \" +str(groupDF + errorDF))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question 2: Caffeine Problem.\n",
    "\n",
    "Does the mean amount of nectar taken depend on the concentration of caffeine in the nectar? Carry out ANOVA step by step in the following subsections:\n",
    "\n",
    "a.State the null and alternate hypotheses appropriate to this question.\n",
    "b.Calculated the following summary statistics for each group: sample size, sample mean, sample standard deviation.\n",
    "c.What is the error mean square?\n",
    "d.How many degrees of freedom are associated with errormean square?e.Calculate the estimate of the grand mean.\n",
    "f.Calculate the group sum of squares.\n",
    "g.Calculate the group degrees of freedom and the group mean square.\n",
    "h.What is F for this example?\n",
    "i.Find the P-value for this test?\n",
    "j.Report your results in an ANOVA tab;e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data:\n",
    "import numpy as np\n",
    "        # this was just copied from the worksheet\n",
    "ppm50 =  [-0.4, 0.34, 0.19, 0.05, -0.14]\n",
    "ppm100 = [0.01, -0.39, -0.08, -0.09, -0.31]\n",
    "ppm150 = [0.65, 0.53, 0.39, -0.15, 0.46]\n",
    "ppm200 = [0.24, 0.44, 0.13, 1.03, 0.05]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A:\n",
    "    Null Hypothesis: there will be no relationship between the caffeine concentration in each feeder than the difference in consumption between the caffeinated feeders and the non caffeinated feeders. \n",
    "    \n",
    "    Alternate Hypothesis: There will be a correlation between the concentration in caffeine from each of the non control feeders and the the feeders.The difference in consumption between the two sets will change predictabley with the changing concentration in caffeine. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5\n",
      "5\n",
      "5\n",
      "5\n",
      "0.007999999999999997\n",
      "-0.17200000000000001\n",
      "0.37600000000000006\n",
      "0.378\n",
      "0.2582556872558667\n",
      "0.15157836257197133\n",
      "0.2766658634526493\n",
      "0.3513061343045407\n"
     ]
    }
   ],
   "source": [
    "# B:\n",
    "#cast to array\n",
    "ppm50 = np.array(ppm50) \n",
    "ppm100= np.array(ppm100) \n",
    "ppm150= np.array(ppm150) \n",
    "ppm200= np.array(ppm200) \n",
    "\n",
    "# sample size for each group\n",
    "print(len(ppm50))  \n",
    "print(len(ppm100))\n",
    "print(len(ppm150))\n",
    "print(len(ppm200))\n",
    "\n",
    "# sample means for each group:\n",
    "print(sum(ppm50)/len(ppm50))\n",
    "print(sum(ppm100)/len(ppm100))\n",
    "print(sum(ppm150)/len(ppm150))\n",
    "print(sum(ppm200)/len(ppm200))\n",
    "\n",
    "# Sample standard deviations for each group:\n",
    "print(ppm50.std())\n",
    "print(ppm100.std())\n",
    "print(ppm150.std())\n",
    "print(ppm200.std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.14750000000000002\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>features</th>\n",
       "      <th>concentration</th>\n",
       "      <th>grandmean</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ppm50</td>\n",
       "      <td>-0.40</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ppm50</td>\n",
       "      <td>0.34</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ppm50</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ppm50</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ppm50</td>\n",
       "      <td>-0.14</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>ppm100</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>ppm100</td>\n",
       "      <td>-0.39</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>ppm100</td>\n",
       "      <td>-0.08</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>ppm100</td>\n",
       "      <td>-0.09</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>ppm100</td>\n",
       "      <td>-0.31</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>ppm150</td>\n",
       "      <td>0.65</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>ppm150</td>\n",
       "      <td>0.53</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>ppm150</td>\n",
       "      <td>0.39</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>ppm150</td>\n",
       "      <td>-0.15</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>ppm150</td>\n",
       "      <td>0.46</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>ppm200</td>\n",
       "      <td>0.24</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>ppm200</td>\n",
       "      <td>0.44</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>ppm200</td>\n",
       "      <td>0.13</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>ppm200</td>\n",
       "      <td>1.03</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>ppm200</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   features  concentration  grandmean\n",
       "0     ppm50          -0.40     0.1475\n",
       "1     ppm50           0.34     0.1475\n",
       "2     ppm50           0.19     0.1475\n",
       "3     ppm50           0.05     0.1475\n",
       "4     ppm50          -0.14     0.1475\n",
       "5    ppm100           0.01     0.1475\n",
       "6    ppm100          -0.39     0.1475\n",
       "7    ppm100          -0.08     0.1475\n",
       "8    ppm100          -0.09     0.1475\n",
       "9    ppm100          -0.31     0.1475\n",
       "10   ppm150           0.65     0.1475\n",
       "11   ppm150           0.53     0.1475\n",
       "12   ppm150           0.39     0.1475\n",
       "13   ppm150          -0.15     0.1475\n",
       "14   ppm150           0.46     0.1475\n",
       "15   ppm200           0.24     0.1475\n",
       "16   ppm200           0.44     0.1475\n",
       "17   ppm200           0.13     0.1475\n",
       "18   ppm200           1.03     0.1475\n",
       "19   ppm200           0.05     0.1475"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "# begin by concating arrays to simplify calculations:\n",
    "ppmtot = np.concatenate((ppm50, ppm100, ppm150, ppm200), axis=0)\n",
    "values = (['ppm50'] * len(ppm50)) + (['ppm100']*len(ppm100)) + (['ppm150'] * len(ppm150)) + (['ppm200']*len(ppm200))\n",
    "caffeine = pd.DataFrame({'features': values, 'concentration': ppmtot})\n",
    "grandmean = caffeine['concentration'].mean()\n",
    "caffeine['grandmean'] = grandmean\n",
    "print(grandmean)\n",
    "caffeine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>features</th>\n",
       "      <th>concentration</th>\n",
       "      <th>grandmean_x</th>\n",
       "      <th>groupmean</th>\n",
       "      <th>grandmean_y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ppm50</td>\n",
       "      <td>-0.40</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.008</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ppm50</td>\n",
       "      <td>0.34</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.008</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ppm50</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.008</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ppm50</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.008</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ppm50</td>\n",
       "      <td>-0.14</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.008</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>ppm100</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>-0.172</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>ppm100</td>\n",
       "      <td>-0.39</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>-0.172</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>ppm100</td>\n",
       "      <td>-0.08</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>-0.172</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>ppm100</td>\n",
       "      <td>-0.09</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>-0.172</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>ppm100</td>\n",
       "      <td>-0.31</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>-0.172</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>ppm150</td>\n",
       "      <td>0.65</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.376</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>ppm150</td>\n",
       "      <td>0.53</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.376</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>ppm150</td>\n",
       "      <td>0.39</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.376</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>ppm150</td>\n",
       "      <td>-0.15</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.376</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>ppm150</td>\n",
       "      <td>0.46</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.376</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>ppm200</td>\n",
       "      <td>0.24</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.378</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>ppm200</td>\n",
       "      <td>0.44</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.378</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>ppm200</td>\n",
       "      <td>0.13</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.378</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>ppm200</td>\n",
       "      <td>1.03</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.378</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>ppm200</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.1475</td>\n",
       "      <td>0.378</td>\n",
       "      <td>0.1475</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   features  concentration  grandmean_x  groupmean  grandmean_y\n",
       "0     ppm50          -0.40       0.1475      0.008       0.1475\n",
       "1     ppm50           0.34       0.1475      0.008       0.1475\n",
       "2     ppm50           0.19       0.1475      0.008       0.1475\n",
       "3     ppm50           0.05       0.1475      0.008       0.1475\n",
       "4     ppm50          -0.14       0.1475      0.008       0.1475\n",
       "5    ppm100           0.01       0.1475     -0.172       0.1475\n",
       "6    ppm100          -0.39       0.1475     -0.172       0.1475\n",
       "7    ppm100          -0.08       0.1475     -0.172       0.1475\n",
       "8    ppm100          -0.09       0.1475     -0.172       0.1475\n",
       "9    ppm100          -0.31       0.1475     -0.172       0.1475\n",
       "10   ppm150           0.65       0.1475      0.376       0.1475\n",
       "11   ppm150           0.53       0.1475      0.376       0.1475\n",
       "12   ppm150           0.39       0.1475      0.376       0.1475\n",
       "13   ppm150          -0.15       0.1475      0.376       0.1475\n",
       "14   ppm150           0.46       0.1475      0.376       0.1475\n",
       "15   ppm200           0.24       0.1475      0.378       0.1475\n",
       "16   ppm200           0.44       0.1475      0.378       0.1475\n",
       "17   ppm200           0.13       0.1475      0.378       0.1475\n",
       "18   ppm200           1.03       0.1475      0.378       0.1475\n",
       "19   ppm200           0.05       0.1475      0.378       0.1475"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute the sum of squares error:\n",
    "sstotal = sum((caffeine['concentration'] -caffeine['grandmean'])**2)\n",
    "sstotal\n",
    "# now get the sum of squares error:\n",
    "groupmeans = caffeine.groupby('features').mean()\n",
    "groupmeans = groupmeans.rename(columns ={'concentration':'groupmean'})\n",
    "groupmeans\n",
    "\n",
    "# add data to original data frame\n",
    "\n",
    "\n",
    "# add group means and overall mean to the original data frame\n",
    "caffeine = caffeine.merge(groupmeans, left_on = 'features', right_index = True)\n",
    "caffeine\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.5825750000000003\n",
      "1.4481600000000001\n",
      "1.1344150000000002\n"
     ]
    }
   ],
   "source": [
    "caffeine\n",
    "sserror = sum((caffeine['concentration'] - caffeine['groupmean'])**2)\n",
    "print(sstotal)\n",
    "print(sserror)\n",
    "ssgroup = sstotal-sserror\n",
    "print(ssgroup)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.07621894736842105"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#C: \n",
    "#Error mean square:\n",
    "errorDF = 20-1\n",
    "k = 4\n",
    "errorMS = sserror/ errorDF\n",
    "errorMS  # gives the error mean square"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#D:\n",
    "The sample size is 20 in total, and the error degrees of freedom is the sample size N - 1, so 20 - 1 is 19 so the errorDF == 19. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#E:\n",
    "The grand mean is the the sum of all the means form all the groups. The calculation for the grand mean is done 5 cells above, it is the mean of each grouped mean. for this data the grand mean calcualtedf from each individual record is the same as the grand mean calculated from each aggregated group mean. The grand mean value is 0.147 in difference between the caffenaited stations and the non caffeinated stations. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#F:\n",
    "The group sum of squares is 1.1344, it is calculated 3 code cdells above. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#G:\n",
    "The group degrees of freedom is the total number of responding groups which in this scenario is four groups. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "group F stat = 4.961211698523185\n",
      "0.3781383333333334\n"
     ]
    }
   ],
   "source": [
    "#H:\n",
    "GroupMS = ssgroup / (k-1)\n",
    "Fstat = GroupMS / errorMS\n",
    "print(\"group F stat = \" +str(Fstat))\n",
    "print(GroupMS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#I:\n",
    "Using the table with the F statistic, the numerator at 4, and the denominator at 20, the calculated p value is significant at p < 0.02"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# J\n",
    "#### recap: full ANOVA table:\n",
    "# Source:     \tSS  \tDF      \tMS      \tF                                            \n",
    "# group: \t SS  \tDF      \tMS      \tF\n",
    "        1.134415      4       0.3781383    3.720908773892388\n",
    "# Error: \t SS  \tDF      \tMS     \n",
    "        1.44816,    19   .07621894736842105\n",
    "# Total: \n",
    "  2.5825750000000003                 \n",
    "    \n",
    "    # completed ANOVA table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question 3:\n",
    "Question considering correlation in hyena age and giggle fequency in hertz. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f1637d13588>"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATPklEQVR4nO3db4xc1XnH8e9TY9INjVj+WBQvtCYNctUWBTsrREqK0tDGQCPsopQSRcUhSG4k0iaN6sRWpCgvWgF1W5q0FRUNaSFKExJqDGpIHQqJ+grUNSaYhLgYCsGLwRvCQhq2inGevpizZrzsrmft2flz5vuRVnvn3DPeZ+7e+fnuOffOjcxEklSXn+l2AZKk9jPcJalChrskVchwl6QKGe6SVKHjul0AwKmnnporVqzodhmS1Fd27Njxg8xcNtu6ngj3FStWMDY21u0yJKmvRMTTc61zWEaSKmS4S1KFDHdJqpDhLkkVMtwlqUI9cbbM0di2c5wt23fz7OQUy4eH2LhmJetWjXS7LEnqCX0Z7tt2jrN56y6mDhwEYHxyis1bdwEY8JJEnw7LbNm++1CwT5s6cJAt23d3qSJJ6i19Ge7PTk4tqF2SBk1fhvvy4aEFtUvSoOnLcN+4ZiVDS5cc1ja0dAkb16zsUkWS1Fv6ckJ1etLUs2UkaXZ9Ge7QCHjDXJJm15fDMpKk+RnuklQhw12SKmS4S1KFDHdJqpDhLkkVMtwlqUKGuyRVyHCXpAoZ7pJUoZbCPSL+JCK+ExGPRsSXIuJnI+KsiHgwIvZExO0RcXzp+4byeE9Zv2JRX4Ek6XWOGO4RMQL8MTCamb8GLAGuBG4AbszMtwAvAteUp1wDvFjabyz9JEkd1OqwzHHAUEQcB7wR2Ae8C7ijrL8VWFeW15bHlPUXRUS0pVpJUkuOGO6ZOQ78JfB9GqH+ErADmMzMV0u3vcD0RzSOAM+U575a+p8y89+NiA0RMRYRYxMTE8f6OiRJTVoZljmJxtH4WcBy4ATg4mP9wZl5c2aOZubosmXLjvWfkyQ1aWVY5reA/8nMicw8AGwFLgCGyzANwBnAeFkeB84EKOtPBF5oa9WSpHm1Eu7fB86PiDeWsfOLgO8C3wTeW/qsB+4qy3eXx5T192dmtq9kSdKRtDLm/iCNidGHgF3lOTcDnwA+FhF7aIyp31KecgtwSmn/GLBpEeqWJM0jeuGgenR0NMfGxrpdhiT1lYjYkZmjs63zClVJqpDhLkkVMtwlqUKGuyRVyHCXpAoZ7pJUIcNdkipkuEtShQx3SaqQ4S5JFTLcJalChrskVchwl6QKGe6SVCHDXZIqZLhLUoUMd0mqkOEuSRUy3CWpQoa7JFXIcJekChnuklQhw12SKmS4S1KFDHdJqpDhLkkVMtwlqUKGuyRVyHCXpAoZ7pJUIcNdkipkuEtShQx3SaqQ4S5JFTLcJalChrskVeiI4R4RKyPi4aavlyPioxFxckTcGxGPl+8nlf4REZ+NiD0R8UhErF78lyFJanbEcM/M3Zl5bmaeC7wNeAW4E9gE3JeZZwP3lccAlwBnl68NwE2LULckaR4LHZa5CHgiM58G1gK3lvZbgXVleS1wWzY8AAxHxOntKFaS1JqFhvuVwJfK8mmZua8sPwecVpZHgGeanrO3tB0mIjZExFhEjE1MTCywDEnSfFoO94g4HrgM+OrMdZmZQC7kB2fmzZk5mpmjy5YtW8hTJUlHsJAj90uAhzLz+fL4+enhlvJ9f2kfB85set4ZpU2S1CELCff38dqQDMDdwPqyvB64q6n9qnLWzPnAS03DN5KkDjiulU4RcQLw28AfNjVfD3wlIq4BngauKO33AJcCe2icWXN126qVJLWkpXDPzB8Dp8xoe4HG2TMz+yZwbVuqkyQdFa9QlaQKGe6SVCHDXZIqZLhLUoUMd0mqUEtny/SibTvH2bJ9N89OTrF8eIiNa1aybtXrPuVAkgZSX4b7tp3jbN66i6kDBwEYn5xi89ZdAAa8JNGnwzJbtu8+FOzTpg4cZMv23V2qSJJ6S1+G+7OTUwtql6RB05fDMsuHhxifJciXDw91oRr1IudkNOj68sh945qVDC1dcljb0NIlbFyzsksVqZdMz8mMT06RvDYns22nH06qwdGX4b5u1QjXXX4OI8NDBDAyPMR1l5/jkZkA52Qk6NNhGWgEvGGu2TgnI/Xpkbs0n7nmXpyT0SAx3HXItp3jXHD9/Zy16WtccP39fTtG7ZyM1MfDMmqvmi4Mm67Xs2U0yAx3AfNPQvZjKDono0HnsIwAJyGl2hjuApyElGpjuAtwElKqjWPuApyElGpjuOsQJyGlejgsI0kVMtwlqUKGuyRVyHCXpAoZ7pJUIcNdkipkuEtShQx3SaqQ4S5JFTLcJalChrskVchwl6QKGe6SVCHDXZIq1FK4R8RwRNwREd+LiMci4u0RcXJE3BsRj5fvJ5W+ERGfjYg9EfFIRKxe3JcgSZqp1SP3zwD/npm/DLwVeAzYBNyXmWcD95XHAJcAZ5evDcBNba1YknRERwz3iDgRuBC4BSAzf5KZk8Ba4NbS7VZgXVleC9yWDQ8AwxFxepvrliTNo5Uj97OACeCfImJnRHwuIk4ATsvMfaXPc8BpZXkEeKbp+XtL22EiYkNEjEXE2MTExNG/AknS67QS7scBq4GbMnMV8GNeG4IBIDMTyIX84My8OTNHM3N02bJlC3mqJOkIWgn3vcDezHywPL6DRtg/Pz3cUr7vL+vHgTObnn9GaZMkdcgRwz0znwOeiYiVpeki4LvA3cD60rYeuKss3w1cVc6aOR94qWn4RpLUAce12O+PgC9GxPHAk8DVNP5j+EpEXAM8DVxR+t4DXArsAV4pfSVJHdRSuGfmw8DoLKsumqVvAtceW1mSpGPR6pF7z9m2c5wt23fz7OQUy4eH2LhmJetWve6kHGng+V4ZTH0Z7tt2jrN56y6mDhwEYHxyis1bdwG400pNfK8Mrr78bJkt23cf2lmnTR04yJbtu7tUkdSbfK8Mrr4M92cnpxbULg0q3yuDqy/Dffnw0ILapUHle2Vw9WW4b1yzkqGlSw5rG1q6hI1rVs7xDGkw+V4ZXH05oTo9EeQZANL8fK8Mrmiclt5do6OjOTY21u0yJKmvRMSOzJztGqT+HJaRJM3PcJekChnuklQhw12SKmS4S1KFDHdJqpDhLkkVMtwlqUKGuyRVyHCXpAr15WfLqG7eOUg6doa7eop3DpLaw2EZ9RTvHCS1h+GunuKdg6T2MNzVU7xzkNQehrt6incOktrDCVX1FO8cJLWH4a6es27ViGEuHSOHZSSpQh65q0peCKVBZ7irOl4IJTksowp5IZRkuKtCXgglOSyjCi0fHmJ8liDvxwuhnDvQ0fLIXdWp5UKo6bmD8ckpktfmDrbtHO92aeoDhruqs27VCNddfg4jw0MEMDI8xHWXn9N3R7zOHehYOCyjKtVwIZRzBzoWHrlLPcoPUdOxaOnIPSKeAn4EHARezczRiDgZuB1YATwFXJGZL0ZEAJ8BLgVeAT6QmQ+1v3Rp8fTCRObGNSsPO18fjm7uoBdei15vsX8vCzly/83MPDczR8vjTcB9mXk2cF95DHAJcHb52gDc1K5ipU7olYnMdswd9Mpr0eE68Xs5ljH3tcA7y/KtwLeAT5T22zIzgQciYjgiTs/MfcdSqNQp801kdvqI91jnDnrpteg1nfi9tHrknsA3ImJHRGwobac1BfZzwGlleQR4pum5e0vbYSJiQ0SMRcTYxMTEUZQuLY6aJjJrei016cTvpdVwf0dmrqYx5HJtRFzYvLIcpedCfnBm3pyZo5k5umzZsoU8VVpUNU1k1vRaatKJ30tL4Z6Z4+X7fuBO4Dzg+Yg4HaB831+6jwNnNj39jNIm9YVaLoKCul5LTTrxezliuEfECRHxpull4N3Ao8DdwPrSbT1wV1m+G7gqGs4HXnK8Xf2klougoK7XUpNO/F6iMaIyT4eIN9M4WofGBOy/ZOafR8QpwFeAXwCepnEq5A/LqZB/B1xM41TIqzNzbL6fMTo6mmNj83aRJM0QETuazmA8zBHPlsnMJ4G3ztL+AnDRLO0JXHsUdUqS2sSPH5Aq50VMg8lwlyrmXakGl58tI1XMT5YcXIa7VDEvYhpcDstUwnFVzaamu1JpYTxyr4AfDqW5eBHT4DLcK+C4qubiRUyDy2GZCjiuqvnUcFcqLZxH7hXww6EkzWS4V8Bx1fbbtnOcC66/n7M2fY0Lrr/f+Qv1HYdlKjD9J7dny7SHF/6oBoZ7JRxXbR/vXqQaOCwjzeAEtWpguEszOEGtGhju0gxOUKsGjrlLMzhBrRoY7tIsnKBWv3NYRpIqZLhLUoUMd0mqkOEuSRUy3CWpQoa7JFXIcJekChnuklQhw12SKmS4S1KFDHdJqpDhLkkVMtwlqUKGuyRVyHCXpAoZ7pJUIcNdkipkuEtShQx3SapQy/dQjYglwBgwnpnviYizgC8DpwA7gD/IzJ9ExBuA24C3AS8Av5+ZT7W9cqnHbds57k221TULOXL/CPBY0+MbgBsz8y3Ai8A1pf0a4MXSfmPpJw2UbTvH2bx1F+OTUyQwPjnF5q272LZzvNulaUC0FO4RcQbwO8DnyuMA3gXcUbrcCqwry2vLY8r6i0p/aWBs2b6bqQMHD2ubOnCQLdt3d6kiDZpWj9z/Bvg48NPy+BRgMjNfLY/3AtN/b44AzwCU9S+V/oeJiA0RMRYRYxMTE0dXvdSjnp2cWlC71G5HDPeIeA+wPzN3tPMHZ+bNmTmamaPLli1r5z8tdd3y4aEFtUvt1sqR+wXAZRHxFI0J1HcBnwGGI2J6QvYMYHowcRw4E6CsP5HGxKo0MDauWcnQ0iWHtQ0tXcLGNSu7VJEGzRHDPTM3Z+YZmbkCuBK4PzPfD3wTeG/pth64qyzfXR5T1t+fmdnWqqUet27VCNddfg4jw0MEMDI8xHWXn+PZMuqYlk+FnMUngC9HxJ8BO4FbSvstwBciYg/wQxr/IUgDZ92qEcNcXbOgcM/MbwHfKstPAufN0uf/gN9rQ22SpKPkFaqSVCHDXZIqZLhLUoUMd0mqUPTCWYoRMQE8vYg/4lTgB4v477eLdbZfv9Rqne3VL3XCsdX6i5k561WgPRHuiy0ixjJztNt1HIl1tl+/1Gqd7dUvdcLi1eqwjCRVyHCXpAoNSrjf3O0CWmSd7dcvtVpne/VLnbBItQ7EmLskDZpBOXKXpIFiuEtShaoI94g4MyK+GRHfjYjvRMRHZunzzoh4KSIeLl+f6katpZanImJXqWNslvUREZ+NiD0R8UhErO5CjSubttXDEfFyRHx0Rp+ubdOI+HxE7I+IR5vaTo6IeyPi8fL9pDmeu770eTwi1s/WZ5Hr3BIR3yu/2zsjYniO5867n3Sgzk9HxHjT7/fSOZ57cUTsLvvrpi7UeXtTjU9FxMNzPLeT23PWTOroPpqZff8FnA6sLstvAv4b+JUZfd4J/Fu3ay21PAWcOs/6S4GvAwGcDzzY5XqXAM/RuGCiJ7YpcCGwGni0qe0vgE1leRNwwyzPOxl4snw/qSyf1OE63w0cV5ZvmK3OVvaTDtT5aeBPW9g3ngDeDBwPfHvme2+x65yx/q+AT/XA9pw1kzq5j1Zx5J6Z+zLzobL8I+AxXrunaz9aC9yWDQ/QuOvV6V2s5yLgicxczKuIFyQz/5PG/QKaNd+cvfmm7c3WAPdm5g8z80XgXuDiTtaZmd/I1+4//ACNO5l11RzbsxXnAXsy88nM/AmNu7WtbWtxTearMyICuAL40mL9/FbNk0kd20erCPdmEbECWAU8OMvqt0fEtyPi6xHxq52t7DAJfCMidkTEhlnWH7rJeNF8A/JuuJK53zC9sk0BTsvMfWX5OeC0Wfr02rb9II2/0mZzpP2kEz5cho8+P8cQQi9tz98Ans/Mx+dY35XtOSOTOraPVhXuEfFzwL8CH83Ml2esfojGsMJbgb8FtnW4vGbvyMzVwCXAtRFxYRdrmVdEHA9cBnx1ltW9tE0Pk42/b3v6PN+I+CTwKvDFObp0ez+5Cfgl4FxgH40hj172PuY/au/49pwvkxZ7H60m3CNiKY2N+MXM3DpzfWa+nJn/W5bvAZZGxKkdLnO6lvHyfT9wJ6+/o9Whm4wXzTcg77RLgIcy8/mZK3ppmxbPTw9fle/7Z+nTE9s2Ij4AvAd4f3mTv04L+8miysznM/NgZv4U+Mc5fn6vbM/jgMuB2+fq0+ntOUcmdWwfrSLcy1jbLcBjmfnXc/T5+dKPiDiPxmt/oXNVHqrjhIh40/Qyjcm1R2d0uxu4KhrOB15q+lOu0+Y8GuqVbdqk+ebszTdtb7YdeHdEnFSGGd5d2jomIi4GPg5clpmvzNGnlf1kUc2Y5/ndOX7+fwFnR8RZ5a+8K2n8Hjrtt4DvZebe2VZ2envOk0md20c7MXO82F/AO2j8efMI8HD5uhT4EPCh0ufDwHdozOY/APx6l2p9c6nh26WeT5b25loD+HsaZyHsAka7VOsJNML6xKa2ntimNP7D2QccoDEmeQ1wCnAf8DjwH8DJpe8o8Lmm534Q2FO+ru5CnXtojKlO76v/UPouB+6Zbz/pcJ1fKPvfIzRC6fSZdZbHl9I4G+SJbtRZ2v95er9s6tvN7TlXJnVsH/XjBySpQlUMy0iSDme4S1KFDHdJqpDhLkkVMtwlqUKGuyRVyHCXpAr9PxFHKiG5QgWTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Question 3\n",
    "\n",
    "age = [2,2,2,6,9,10,13,10,14,14,12,7,11,11,14,20]\n",
    "freqHrz = [840, 670, 580, 470, 540, 660, 510, 520, 500, 480, 400, 650, 460, 500, 580, 500]\n",
    "\n",
    "plt.scatter(age, freqHrz)\n",
    "# this suggests a negative association, that frequency decreases with age. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "380.4375"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#B:\n",
    "globalmeanage  =(sum(age)/len(age))\n",
    "globalmeanfreq =(sum(freqHrz)/len(freqHrz))\n",
    "hyenas = pd.DataFrame({'age': age, 'freqhrz': freqHrz, 'meanage':globalmeanage, 'meanhrz': globalmeanfreq})\n",
    "hyenas\n",
    "sshyenaAge = sum((hyenas['age']-hyenas['meanage'])**2)\n",
    "sshyenaAge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "174175.0"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#C:\n",
    "sshyenafrz = sum((hyenas['freqhrz']-hyenas['meanhrz'])**2)\n",
    "sshyenafrz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4898.75"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#D:\n",
    "hyenas\n",
    "sumprod = sum((hyenas['age']-(hyenas['meanage']))*(hyenas['freqhrz'] - (hyenas['meanhrz'])))\n",
    "sumprod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.60179795461996"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#E:\n",
    "b = sumprod / sshyenaAge\n",
    "r = sumprod / ((sshyenaAge**(1/2))*(sshyenafrz**(1/2)))\n",
    "r\n",
    "# correlation coefficient is -0.60179795461996"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.6959612353906193"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#F\n",
    "# computed fischer's z transformation\n",
    "\n",
    "def ztransform(r):\n",
    "    z = (1/2)*(np.log(1+r)-np.log(1-r))\n",
    "    return z\n",
    "ztransform(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100+1.713913650100261j)"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# G:\n",
    "# computing the standard error\n",
    "((1-(z*z)/(16-2))**(1/2))\n",
    "# coputing the standard error of sampling distribution:\n",
    "sigma = ((1/16-3)**(1/2))\n",
    "sigma+100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# H:\n",
    "Zcrit = -1.645\n",
    "# at alpha = .05, the Z critical vlaue is -1.645. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-0.6959612353906195-3.3592707541965114j)\n",
      "(-0.695961235390619+3.3592707541965114j)\n",
      "(-0.6959612353906195-3.3592707541965114j)<Z<(-0.695961235390619+3.3592707541965114j)\n"
     ]
    }
   ],
   "source": [
    "# I:\n",
    "r_ = (1+r)/1.6017979546199599\n",
    "z2=0.5*np.log(r_)#\n",
    "z2\n",
    "\n",
    "low = z2-(1.96*(sigma))\n",
    "high = z2+(1.96*(sigma))\n",
    "\n",
    "print(low)\n",
    "print(high)\n",
    "\n",
    "print(str(low)+\"<Z<\"+str(high))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# J:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Year ['Sat', 'Sun', 'Mon', 'Tue', 'Wed', 'Thur', 'Fri']\n",
      "II1971 [0.83, 3.14, 4.2, 1.28, 1.16, 4.25, 2.08]\n",
      "I1972 [4.66, 4.15, 3.4, 1.74, 3.91, 5.15, 5.06]\n",
      "II1972 [3.03, 5.8, 2.29, 3.17, 3.5, 3.4, 3.04]\n",
      "I1973 [3.69, 3.72, 4.29, 2.06, 3.04, 2.3, 4.26]\n",
      "II1973 [2.35, 3.62, 3.56, 2.27, 4.46, 2.52, 3.36]\n",
      "I1974 [3.18, 3.28, 1.82, 3.75, 2.07, 3.54, 2.27]\n"
     ]
    }
   ],
   "source": [
    "# question #4:\n",
    "\n",
    "precipitationBoston = {\n",
    "\"Year\": [\"Sat\",\"Sun\", \"Mon\", \"Tue\", \"Wed\", \"Thur\", \"Fri\"],\n",
    "\"II1971\": [0.83, 3.14, 4.20, 1.28, 1.16, 4.25, 2.08],\n",
    "\"I1972\": [4.66, 4.15, 3.40, 1.74, 3.91, 5.15, 5.06],\n",
    "\"II1972\":  [3.03, 5.80, 2.29, 3.17, 3.50, 3.40, 3.04],\n",
    "\"I1973\":  [3.69, 3.72, 4.29, 2.06, 3.04, 2.30, 4.26],\n",
    "\"II1973\":  [2.35, 3.62, 3.56, 2.27, 4.46, 2.52, 3.36],\n",
    "\"I1974\":  [3.18, 3.28, 1.82, 3.75, 2.07, 3.54, 2.27],\n",
    "    \n",
    "    }\n",
    "#precipitationBoston.values[0]\n",
    "\n",
    "for key, value in precipitationBoston.items():\n",
    "    print(key, value)\n",
    "    \n",
    "df = pd.DataFrame(precipitationBoston, columns = ['Year',\"II1971\", \"I1972\",\"II1972\", \"I1973\", \"II1973\", \"I1974\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0    2.956667\n",
      "1    3.951667\n",
      "2    3.260000\n",
      "3    2.378333\n",
      "4    3.023333\n",
      "5    3.526667\n",
      "6    3.345000\n",
      "dtype: float64\n",
      "0    1.295001\n",
      "1    0.972552\n",
      "2    1.006837\n",
      "3    0.920335\n",
      "4    1.221354\n",
      "5    1.067027\n",
      "6    1.151308\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "df.set_index('Year')\n",
    "# compute the mean for each weekday\n",
    "print(df.mean(axis=1))\n",
    "# compute the std for each weekday\n",
    "print(df.std(axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ANOVA\n",
    "grandmean = df.mean()\n",
    "# Source:     \tSS  \tDF      \tMS      \tF\n",
    "# group: \tSST \tk−1 \tSST/(k−1) \tMST/MSE\n",
    "# Error:\tSSE \tN−k \tSSE/(N−k)\n",
    "# Total: SST + SSE, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "precipitationPittsburg = {\n",
    "\"Year\": [\"Sat\",\"Sun\", \"Mon\", \"Tue\", \"Wed\", \"Thur\", \"Fri\"],\n",
    "\"II1971\": [1.64,5.55,3.19,2.45,1.44,1.07,1.66],\n",
    "\"I1972\": [2.20,3.37,0.78,2.63,2.32,5.57,2.80],\n",
    "\"II1972\":  [2.75, 1.72,2.34,3.40,3.68,3.48,2.50],\n",
    "\"I1973\":  [2.23,4.31,2.02,1.83,4.35,4.07,2.66],\n",
    "\"II1973\":  [3.65,2.66,3.95,2.31,1.85,2.63,1.11],\n",
    "\"I1974\":  [4.96,3.00, 2.61,1.75,2.70,2.45,4.06], \n",
    "    }\n",
    "\n",
    "df2 = pd.DataFrame(precipitationPittsburg, columns = ['Year',\"II1971\", \"I1972\",\"II1972\", \"I1973\", \"II1973\", \"I1974\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0    2.905000\n",
      "1    3.435000\n",
      "2    2.481667\n",
      "3    2.395000\n",
      "4    2.723333\n",
      "5    3.211667\n",
      "6    2.465000\n",
      "dtype: float64\n",
      "0    1.212844\n",
      "1    1.339713\n",
      "2    1.077319\n",
      "3    0.601789\n",
      "4    1.107622\n",
      "5    1.541148\n",
      "6    1.017797\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "df2.set_index('Year')\n",
    "# compute the mean for each weekday\n",
    "print(df2.mean(axis=1))\n",
    "# compute the std for each weekday\n",
    "print(df2.std(axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ANOVA:\n",
    "grandmean = df2.mean()\n",
    "# Source:     \tSS  \tDF      \tMS      \tF\n",
    "# group: \tSST \tk−1 \tSST/(k−1) \tMST/MSE\n",
    "# Error:\tSSE \tN−k \tSSE/(N−k)\n",
    "# Total: SST + SSE, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f1638b1fb70>"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAASvElEQVR4nO3df6zddX3H8ed75eruAuM62016aa2J2IkwLLtBDMtGZsxFh9IxlsASUaJpYnRiZppRk0nG/gDXRDPFSRohwmLQRbtrdbg7MtjQJSK3vYVaapdGN+kpCQW94I+72db3/jineHu9t+d7es89557P9/lIbjjn+/30fN9fLn3xOZ/v5/v5RmYiSRp8v9LvAiRJ3WGgS1IhDHRJKoSBLkmFMNAlqRBn9evAq1evzg0bNvTr8JI0kHbv3v1sZq5ZaF/fAn3Dhg1MTU316/CSNJAi4n8W2+eQiyQVwkCXpEK0DfSI+NWI+FZEPB4R+yPirxdo89KI+EJEHIqIRyNiw7JUK0laVJUe+v8Bf5iZlwCvB66KiMvntXk38MPMfDXwceCjXa1SktRW20DPph+33g61fuYvAHMNcG/r9ReBN0VEdK1KSVJblWa5RMQqYDfwauBTmfnovCajwFMAmXk8Ip4HXg48O+9ztgBbANavX7+0yiVpwExMN9g+eZAjM7OsHRlm6/hGNm8a7drnV7oompknMvP1wPnAZRFx0ZkcLDN3ZOZYZo6tWbPgNEpJKtLEdINtO/fRmJklgcbMLNt27mNiutG1Y3Q0yyUzZ4CHgavm7WoA6wAi4izgXOC5LtQnSUXYPnmQ2WMnTtk2e+wE2ycPdu0YVWa5rImIkdbrYeDNwHfmNdsFvLP1+jrgoXShdUl60ZGZ2Y62n4kqPfTzgIcj4gngMeDBzPxqRNwWEW9vtbkbeHlEHAL+ArilaxVKUgHWjgx3tP1MtL0omplPAJsW2P6ROa//F/jTrlUlSYXZOr6RbTv3nTLsMjy0iq3jG7t2jL6t5SJJdXJyNstyznIx0CWpRzZvGu1qgM/nWi6SVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQZ/W7AKluJqYbbJ88yJGZWdaODLN1fCObN432uywVwECXemhiusG2nfuYPXYCgMbMLNt27gMw1LVkbYdcImJdRDwcEU9GxP6IuHmBNldGxPMRsbf185HlKVcabNsnD74Y5ifNHjvB9smDfapIJanSQz8OfCgz90TEOcDuiHgwM5+c1+7rmXl190uUynFkZraj7VIn2vbQM/PpzNzTev0j4ADgd0PpDKwdGe5ou9SJjma5RMQGYBPw6AK73xgRj0fE1yLidYv8+S0RMRURU0ePHu28WmnAbR3fyPDQqlO2DQ+tYuv4xj5VpJJUDvSIOBv4EvDBzHxh3u49wCsz8xLgk8DEQp+RmTsycywzx9asWXOGJUuDa/OmUW6/9mJGR4YJYHRkmNuvvdgLouqKyMz2jSKGgK8Ck5n5sQrt/xsYy8xnF2szNjaWU1NTHZQqSYqI3Zk5ttC+KrNcArgbOLBYmEfEK1rtiIjLWp/73JmXLEnqVJVZLlcA7wD2RcTe1rYPA+sBMvMu4DrgvRFxHJgFrs8qXX9JUte0DfTM/AYQbdrcCdzZraIkSZ1zLRdJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFaBvoEbEuIh6OiCcjYn9E3LxAm4iIT0TEoYh4IiIuXZ5yJUmLOatCm+PAhzJzT0ScA+yOiAcz88k5bd4CXND6eQPw6dY/VWMT0w22Tx7kyMwsa0eG2Tq+kc2bRvtdllSstj30zHw6M/e0Xv8IOADM/1t5DXBfNn0TGImI87perQbGxHSDbTv30ZiZJYHGzCzbdu5jYrrR79KkYnU0hh4RG4BNwKPzdo0CT815f5hfDn0iYktETEXE1NGjRzssVYNk++RBZo+dOGXb7LETbJ882KeKpPJVGXIBICLOBr4EfDAzXziTg2XmDmAHwNjYWJ7JZ2gwHJmZ7Wi7tBiH7qqr1EOPiCGaYf65zNy5QJMGsG7O+/Nb21RTa0eGO9ouLcShu85UmeUSwN3Agcz82CLNdgE3tma7XA48n5lPd7FODZit4xsZHlp1yrbhoVVsHd/Yp4o0iBy660yVIZcrgHcA+yJib2vbh4H1AJl5F/AA8FbgEPBT4KauV6qBcvIrsV+VtRQO3XWmbaBn5jeAaNMmgfd1qyiVYfOmUQNcS7J2ZJjGAuHt0N3CvFNU0orl0F1nKs9ykaRec+iuMwa6pBXNobvqHHKRpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCB1xIWhEmphs+mWiJDHRJfTcx3WDbzn3MHjsBQGNmlm079wEY6h0w0LvAnoW0NNsnD74Y5ifNHjvB9smD/l3qgIG+RPYspKU7MjPb0XYtzIuiS3S6noWkataODHe0XQsz0JfInoW0dFvHNzI8tOqUbcNDq9g6vrFPFQ0mA32J7FlIS7d50yi3X3sxoyPDBDA6Mszt117ssGWH2o6hR8Q9wNXAM5l50QL7rwS+DHyvtWlnZt7WxRpXtK3jG08ZQwd7FkvhBeb62rxp1N/1ElW5KPpZ4E7gvtO0+XpmXt2VigbMyf8ADaGl8wKztDRtAz0zH4mIDT2oZWDZs+jcQj1xp65JS9OtMfQ3RsTjEfG1iHjdYo0iYktETEXE1NGjR7t0aA2akz3xxswsyS964g0vMEtL0o1A3wO8MjMvAT4JTCzWMDN3ZOZYZo6tWbOmC4fWIFqsJ74qYsH2XmCWqlnyjUWZ+cKc1w9ExN9HxOrMfHapn12yOl/8W6zHfSKT4aFVXmCWztCSe+gR8YqIZtcqIi5rfeZzS/3cki025DAx3eh3aT2xWI/75FQ1p65JZ6bKtMX7gSuB1RFxGLgVGALIzLuA64D3RsRxYBa4PjNz2SouQN0v/p1uqqcXmKUzV2WWyw1t9t9Jc1qjKqr73aVO9ZSWh4tz9cHakeEFZ3TU6eKfPXGp+7z1vw9ct0LScrCH3gcOOUhaDgZ6nzjkIKnbHHKRpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhXDaojRHHVbBrMM51pWBLrXU4RF4dTjHOnPIRWo53SqYpajDOdaZgS61LPYIvMW2D6K6r/RZOgNdalnsEXiLbR9Ei63oWaeVPktmoEstJxZ5Lsti2weRK32WzUCXWkZP82i8UmzeNOpj/grmLBdVVvp0t9M9Gq8krvRZLgNdldRhupvr1GvQGeiqpC4Ptrb32pnSv7UNGgNdlTjdTfPV4VvboPGiqCpxupvm8yallcdAVyWDON1tYrrBFXc8xKtu+WeuuOMhJqYb/S6pKH5rW3kcclElg3bB0OGA5bd2ZHjBu2j91tY/BroqG6QLhnW5iNtPdZnmOUgMdBXJ4YDlN2jf2uqgbaBHxD3A1cAzmXnRAvsD+DvgrcBPgXdl5p5uFyp1wuGA3hikb211UOWi6GeBq06z/y3ABa2fLcCnl16WtDSDeBFXWqq2gZ6ZjwA/OE2Ta4D7sumbwEhEnNetAqUz4ZolqqNujKGPAk/NeX+4te3p+Q0jYgvNXjzr16/vwqGlxTkcoLrp6Tz0zNyRmWOZObZmzZpeHlqSiteNQG8A6+a8P7+1TZLUQ90I9F3AjdF0OfB8Zv7ScIskaXlVmbZ4P3AlsDoiDgO3AkMAmXkX8ADNKYuHaE5bvGm5ipUkLa5toGfmDW32J/C+rlUkSTojLs4lSYUw0CWpEK7losp8Oo20shnoqsTlaKWVzyEXVeLTaaSVzx66Kqn7crQON2kQ2ENXJXV+pujJ4abGzCzJL4abfKSdVhoDXZXUeTlah5s0KBxyUSV1fjpN3YebNDgMdFVW1+VoffqRBoVDLlIbdR5u0mCxhy61UefhJg0WA12qoK7DTRosDrlIUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVolKgR8RVEXEwIg5FxC0L7H9XRByNiL2tn/d0v1RJ0um0XT43IlYBnwLeDBwGHouIXZn55LymX8jM9y9DjZKkCqr00C8DDmXmdzPzZ8DngWuWtyxJUqeqBPoo8NSc94db2+b7k4h4IiK+GBHrulKdJKmybl0U/QqwITN/B3gQuHehRhGxJSKmImLq6NGjXTq0JAmqPYKuAcztcZ/f2vaizHxuztvPAH+70Adl5g5gB8DY2Fh2VCkwMd3wuY6StIgqPfTHgAsi4lUR8RLgemDX3AYRcd6ct28HDnSvxKaJ6Qbbdu6jMTNLAo2ZWbbt3MfEdKPtn5WkOmgb6Jl5HHg/MEkzqP8xM/dHxG0R8fZWsw9ExP6IeBz4APCubhe6ffIgs8dOnLJt9tgJtk8e7PahJGkgVRlyITMfAB6Yt+0jc15vA7Z1t7RTHZmZ7Wi7JNXNwNwpunZkuKPtklQ3AxPoW8c3Mjy06pRtw0Or2Dq+sU8VSdLKUmnIZSU4OZvFWS5lcMaS1H0DE+jQDHX/0g++kzOWTl7kPjljCfD3Ky3BwAy5qBzOWJKWh4GunnPGkrQ8DHT1nDOWpOVhoKvnnLEkLY+BuiiqMpzJjCVnxUjtGejqi05mLDkrRqrGIReteM6Kkaox0LXiOStGqsZA14rnrBipGgNdK56zYqRqvCiqFc91fKRqDHQNBNfxkdpzyEWSCmEPvaa8UUcqj4FeQ96oI5XJIZca8kYdqUwGeg15o45UJgO9hrxRRyqTgV5D3qgjlcmLojXkjTpSmQz0mvJGHak8DrlIUiEMdEkqRKVAj4irIuJgRByKiFsW2P/SiPhCa/+jEbGh65VKkk6rbaBHxCrgU8BbgAuBGyLiwnnN3g38MDNfDXwc+Gi3C5UknV6VHvplwKHM/G5m/gz4PHDNvDbXAPe2Xn8ReFNERPfKlCS1U2WWyyjw1Jz3h4E3LNYmM49HxPPAy4Fn5zaKiC3AltbbH0fEctxrvnr+cWukzucO9T5/z70+XrnYjp5OW8zMHcCO5TxGRExl5thyHmOlqvO5Q73P33Ov57nPV2XIpQGsm/P+/Na2BdtExFnAucBz3ShQklRNlUB/DLggIl4VES8Brgd2zWuzC3hn6/V1wEOZmd0rU5LUTtshl9aY+PuBSWAVcE9m7o+I24CpzNwF3A38Q0QcAn5AM/T7ZVmHdFa4Op871Pv8PXcRdqQlqQzeKSpJhTDQJakQAxnoEXFPRDwTEd8+TZsrI2JvROyPiP/oZX3Lqd25R8TW1nnvjYhvR8SJiPiNXte5XCqc/7kR8ZWIeLz1u7+p1zUulwrn/rKI+KeIeCIivhURF/W6xuUSEesi4uGIeLL1e715gTYREZ9oLUHyRERc2o9a+yozB+4H+H3gUuDbi+wfAZ4E1rfe/2a/a+7Vuc9r+zaaM476XncPf/cfBj7aer2G5kX6l/S77h6d+3bg1tbr3wb+rd81d/HczwMubb0+B/gv4MJ5bd4KfA0I4HLg0X7X3eufgeyhZ+YjNP+iLubPgJ2Z+f1W+2d6UlgPVDj3uW4A7l/GcnquwvkncE5r6YmzW22P96K25Vbh3C8EHmq1/Q6wISJ+qxe1LbfMfDoz97Re/wg4QPMO9bmuAe7Lpm8CIxFxXo9L7auBDPQKXgO8LCL+PSJ2R8SN/S6o1yLi14CrgC/1u5YeuxN4LXAE2AfcnJk/729JPfM4cC1ARFxG8xbx8/ta0TJorea6CXh03q6Flimp1VNcSg30s4DfBf4IGAf+KiJe09+Seu5twH9mZtXefCnGgb3AWuD1wJ0R8ev9LKiH7qDZK90L/DkwDZzoa0VdFhFn0+ykfDAzX+h3PStNqY+gOww8l5k/AX4SEY8Al9Acd6uL6ylsuKWim4A7sjmoeigivkdzPPlb/S1r+bUC7iZoXiAEvgd8t69FdVFEDNEM889l5s4FmlRZpqRopfbQvwz8XkSc1Rp6eAPNMbdaiIhzgT+g+e+hbr4PvAmgNX68kYJC7XQiYqS1PAfAe4BHSunFtv4HdTdwIDM/tkizXcCNrdkulwPPZ+bTPStyBRjIHnpE3A9cCayOiMPArcAQQGbelZkHIuJfgCeAnwOfycxFpzgOknbn3mr2x8C/tr6hFKXC+f8N8NmI2EdztsNfZmYRS6tWOPfXAvdGRAL7aT54phRXAO8A9rWGlKA5o2k9vHj+D9Cc6XII+Cmtbyt14q3/klSIUodcJKl2DHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUiP8HontRCUYYJp0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Question #5:\n",
    "\n",
    "facedim = [1.59, 1.67, 1.65, 1.72, 1.79, 1.77, 1.74, 1.74, 1.77, 1.78, 1.76, 1.81, 1.83, 1.83, 1.84, 1.87, 1.92, 1.95, 1.98, 1.99, 2.07]\n",
    "penaltyMin = [.44, 1.43, 1.57, .14, .27, .35, .85, 1.13, 1.47, 1.51, 1.99, 1.06, 1.20, 1.23, .8, 2.53, 1.23, 1.1, 1.61, 1.95, 2.95]\n",
    "#len(penaltyMin)\n",
    "plt.scatter(facedim, penaltyMin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#B:\n",
    "   The scater plot seems to yield some correlation between the face dimensions and the penalty minutes per game, the events seem positivly correlated where athlete with a bigger face seem to incur more penalty minutes per game.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.812857142857143\n",
      "1.2766666666666666\n"
     ]
    }
   ],
   "source": [
    "#C:\n",
    "facedim = np.array(facedim)\n",
    "penaltyMin = np.array(penaltyMin)\n",
    "\n",
    "# show varibale means\n",
    "print(facedim.mean())\n",
    "print(penaltyMin.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3.18898271]]\n",
      "[-4.50450342]\n"
     ]
    }
   ],
   "source": [
    "#D && E:\n",
    "facedim = facedim.reshape(-1, 1)\n",
    "penaltyMin = penaltyMin.reshape(-1, 1)\n",
    "\n",
    "# calculating the estimate of the regression slope and the estimate of the intercept \n",
    "lin = LinearRegression()\n",
    "lin.fit(facedim, penaltyMin)\n",
    "\n",
    "print(lin.coef_)      # A parameter\n",
    "print(lin.intercept_) # B parameter\n",
    "\n",
    "#yes an increase of 3 minute penalty time makes sense for the increase of 1 unit face size. \n",
    "A =lin.coef_\n",
    "B = lin.intercept_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "def yAxBline(A, B):\n",
    "    '''draws a line over the scatterplot that takes\n",
    "    the line's slope and intercept as arguments'''\n",
    "    getaxis = plt.gca()\n",
    "    x = np.array(getaxis.get_xlim())\n",
    "    y = B + A * x\n",
    "    plt.plot(x, y, '--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhsUlEQVR4nO3de3yV1Z3v8c8vV0ICSYAEciUgMVzlFhAdtShaqFq0lgo6raPtHNqe2tp2pC32dNo6Z6YXXtNOZ/SMh1rb2ulLgkqpenSorVbbWiGBgMhNKCg7CXJPuCWQ7KzzRzYY04TsJHvnyX729/168XLv51lm/5aRb1bWs571mHMOERGJfQleFyAiIpGhQBcR8QkFuoiITyjQRUR8QoEuIuITSV598IgRI1xJSYlXHy8iEpM2btx4xDmX09k5zwK9pKSEqqoqrz5eRCQmmdk7XZ3TlIuIiE8o0EVEfKLbQDezQWa2wcy2mNk2M/t2J21SzazCzPaY2XozK4lKtSIi0qVwRuhngeucc1OBacACM5vToc2ngOPOuXHAD4HvRbRKERHpVreB7tqcCr1NDv3puAHMLcDPQ6+fAuaZmUWsShER6VZYq1zMLBHYCIwDHnbOre/QpAAIADjnWsysARgOHOnwdZYCSwGKi4v7VrmIiEfWVteyYt0u6uobyc9KY9n8Mm6dXuB1WeFdFHXOBZ1z04BCYLaZTe7NhznnVjrnyp1z5Tk5nS6jFBEZ0NZW17J8zVZq6xtxQG19I8vXbGVtda3XpfVslYtzrh54GVjQ4VQtUARgZklAJnA0AvWJiAwoK9btorE5+L5jjc1BVqzb5VFF7wlnlUuOmWWFXqcBNwA7OzR7Bvi70OtFwEtOG62LiA/V1Tf26Hh/CmeEnge8bGZvAJXAi86558zsQTNbGGrzE2C4me0Bvgx8LTrlioh4Kz8rrUfH+1O3F0Wdc28A0zs5/o/tXjcBH4tsaSIiA8+y+WUsX7P1fdMuacmJLJtf5mFVbTzby0VEJBadX80yEFe5KNBFRHro1ukFAyLAO9JeLiIiPqFAFxHxCQW6iIhPKNBFRHxCgS4i4hMKdBERn1Cgi4j4hAJdRMQnFOgiIj6hQBcR8QkFuoiITyjQRUR8QoEuIuITCnQREZ9QoIuI+IQCXUTEJxToIiI+oUAXEfEJBbqIiE8o0EVEfEKBLiLiEwp0ERGfSPK6AJF4sba6lhXrdlFX30h+VhrL5pdx6/QCr8sSH1Ggi/SDtdW1LF+zlcbmIAC19Y0sX7MVQKEuEdPtlIuZFZnZy2a23cy2mdl9nbSZa2YNZrY59Ocfo1OuSGxasW7XhTA/r7E5yIp1uzyqSPwonBF6C/APzrlNZjYE2GhmLzrntndo9wfn3M2RL1Ek9tXVN/bouEhvdDtCd84dcM5tCr0+CewA9DuiSA/kZ6X16LhIb/RolYuZlQDTgfWdnL7CzLaY2QtmNqmLf3+pmVWZWdXhw4d7Xq1IjFo2v4y05MT3HUtLTmTZ/DKPKhI/CjvQzSwDeBr4onPuRIfTm4DRzrmpwH8Aazv7Gs65lc65cudceU5OTi9LFok9t04v4Du3TaEgKw0DCrLS+M5tU3RBVCLKnHPdNzJLBp4D1jnnfhBG+7eBcufcka7alJeXu6qqqh6UKiIiZrbROVfe2blwVrkY8BNgR1dhbmajQu0ws9mhr3u09yWLiEhPhbPK5W+ATwBbzWxz6NgDQDGAc+4RYBHwWTNrARqBJS6cob+IiERMt4HunPsjYN20eQh4KFJFiYhIz2kvFxERn1Cgi4j4hAJdRMQnFOgiIj6hQBcR8QkFuoiITyjQRUR8QoEuIuITCnQREZ9QoIuI+IQCXUTEJxToIiI+oUAXEfEJBbqIiE8o0EVEfEKBLiLiEwp0EZF+0hJs5eGX90Tt6yvQRUSipLXV8dqeI/zH73YDkJSYQLDV0doanSd0hvNMURER6YF3G5p4elMNFZUB9h87Q2ZaMnddUULm4GS+MK80ap+rQBcRiaDfbHuXz/zXRlodXDF2OF++4VIWTB7FoOTEqH+2Al1EpA/2HTlNRWWAKQWZ3HRZHrNKhvGZD1zC7eVFlIxI79daFOgiIj3UeC7IC28eYFVlgA37jpGYYHz6mrHcdFke2ekpfGXBeE/qUqCLiPTQJ39WyZ/3HqVk+GC+sqCMRTMKyR06yOuyFOgiIhfTcKaZX2+pZW11LT+9ZzaZacl8/rpx3Hd9KZePGYaZeV3iBQp0EZEOnHO8vvcYq6sCPL/1AGdbWpmUP5R3G5rITEvmynEjvC6xUwp0EZGQ1lZHQoKx98hp7vjx6wwZlMTt5UUsnlXE5IJMr8vrVreBbmZFwOPASMABK51zP+rQxoAfATcCZ4C7nXObIl+uiEhktQRbeXnXYSoq95ORmsS/LZnOJTkZPHZ3OVeMHUFaSvSXG0ZKOCP0FuAfnHObzGwIsNHMXnTObW/X5kNAaejP5cB/hv4pcsHa6lpWrNtFXX0j+VlpLJtfxq3TC7wuS+LUO0fblhs+tbGGQyfPkjMklTtmF184f934kR5W1zvdBrpz7gBwIPT6pJntAAqA9oF+C/C4c84Br5tZlpnlhf5dEdZW17J8zVYam4MA1NY3snzNVgCFuvSbpuYgyYkJJCYYFZUBHnnlL1w3PpfFs4qZW5ZDcmJs74bSozl0MysBpgPrO5wqAALt3teEjr0v0M1sKbAUoLi4GIkfK9btuhDm5zU2B1mxbpcCXaLuzdoGKioDrN1cy4+WTOO68SP55FVjuOuKEkZler/cMFLCDnQzywCeBr7onDvRmw9zzq0EVgKUl5dHZ3caGZDq6ht7dFziQzSn4c61tFJRuZ+KqgBv1p4gJSmBGyePYtTQNABGZKRG5HMGkrAC3cySaQvzXzrn1nTSpBYoave+MHRMBID8rDRqOwnv/Kw0D6qRgSAa03DOOWqON1I0bDBJCcYjr+xlaFoy3144iVunFZA5ODli9Q9E4axyMeAnwA7n3A+6aPYMcK+ZraLtYmiD5s+lvWXzy973lxcgLTmRZfPLPKxKvBTJabhDJ5t4emMtq6sCnGxq5rWvzSMlKYFf3/s3DE9PGVA3/0RTOCP0vwE+AWw1s82hYw8AxQDOuUeA52lbsriHtmWL90S8Uolp5/+CapWLnBeJabitNQ38+0u7eWnnIYKtjtklw7j32nEXzvtxWuViwlnl8kfgoj/eQqtbPheposSfbp1eoACXC3o7DffO0dMkJyaQn5XGqbMtVO+v5++vHsPt5UVckpMRrXJjgu4UFRFP9GQarqk5yH+/+S4VlQH+vPcod19ZwrcWTmLO2GH8efl1Mb/cMFIU6CLiiXCn4b733zv55evvcKKpheJhg7n/g5eyaGbbGgwzIzkxPubHw6FAFxHPdDYNd6KpmZd3HmLh1HzMjFNNLVw7PpfF5UXMGTuchAQFeFcU6CLiOecclW8fZ1Xlfp7feoCm5lZKc4cwMX8oD94yKW5WqfSVAl1EPLX74Ek+/YuN7D1ymozUJG6bUciSWUVMyBsCoDDvAQW6iPSrlmArr+4+TLAVbpg4ksLswRQNG8z/vHYcN04ZxeAUxVJv6b+ciPSL/UfPsLoqwJMbAxw8cZbZY4Zxw8SRpKUk8vNPzva6PF9QoItI1H33hZ088spfSDD4wKU5fHthMfMm5Hpdlu8o0EUk4rbXnWB1VYDPfOASRmUO4vKxw0hPSWRReSF5mdq/J1oU6CISESeamnl2Sx0VlQHeqGkgJTGBOWOHsSAzj2vLcrm2TCPyaFOgi0ifnTrbwpXfeYlTZ1sYP2oI3/zwRG6dVkB2eorXpcUVBbqI9Njhk2dZs6mGd46d4V8+MoWM1CSWzS9jWlEWlxVmaqmhRxToIhKWYKvj1d2HqdgQ4Lc7DtIS2t3wXEsrKUkJ/N2VJV6XGPcU6CISll/8+W2+9ex2hqen8Mmr2nY3HJcb37sbDjQKdBH5K03NQX6z/SAVlfv52Mwibp1ewIen5jMqcxDXjR9JSpJ2NxyIFOgicsHOd09QURngV9W11J9ppiArDUfb43+HZ6SyYHKexxXKxSjQReJcc7D1wn7iX1y1mb2HT/PBSSNZMquYKy/R7oaxRIEuEoecc2x85zirKgP8ftdhfr9sLhmpSfzr7VPJz0zTcsMYpUAXiSP1Z87xZFUNqyr385fDp0lPSWThtHzOnGshIzWJSfmZXpcofaBAF/G5YKvjVFMLmYOTOdDQxD8/v4MZxVl8/6OXcdNleaSnRj8G1lbX6gHh/UCBLuJTgWNneHJjDU9VBbh87HB+uHgaE/KG8vv751IyIr3f6lhbXfu+Z4fW1jeyfM1WAIV6hCnQ+5lGKhJtv9txkJ+99jZ/3HMEgGtKc7hxynurU/ozzKHtmaHtHwQN0NgcZMW6Xfp/P8IU6P1IIxWJlt0HT3JJTgYJCcb6fcfYe/g0980r5WPlRRRkebu7YV19Y4+OS+8p0PuRRioSSafOtvDcljpWVQbYHKjnF5+azdWlOdw3r5SvLhhP4gBZbpiflUZtJ+Gd7/EPGj9SoPcjjVQkEhrONPPPz2/nuTcOcOZckEtHZvCNmycyObRCpT8ucvbEsvll7/vNFCAtOZFl88s8rMqfBtZ33uc0UpHeOnrqLPuOnKa8ZBjpqYlUvXOchVPzuX1WEdOLsgb07obnf/vUtaPoM+fcxRuYPQbcDBxyzk3u5Pxc4NfAvtChNc65B7v74PLycldVVdXTemNaxzl0aBupfOe2KfqfuxPxfgE52Or4w+7DrK4K8OL2gwxLT+G1r80jMcFobXW6gzNOmdlG51x5Z+fCGaH/DHgIePwibf7gnLu5F7XFFY1UwhfvF5Cf33qA//3cduoamsgenMxdV5SweFbRhXlxhbl0pttAd869amYl/VBLXLh1ekFcBFJPdDYSj7cLyGdbgvx2+yEmFwxl9PB0hg5K5pLcDL5+00Sun5hLalKi1yVKDIjUHPoVZrYFqAPud85t66yRmS0FlgIUFxdH6KMllnU1Eu8Y5uf57QLyWwdPUlEZYM2mGo6faeYL80r58g2XclXpCK4qHeF1eRJjIhHom4DRzrlTZnYjsBYo7ayhc24lsBLa5tAj8NkS47oaiSeaEezk+o5fLiA757jzx+v5896jJCcaN0wcyeJZxVw1TiEuvdfnQHfOnWj3+nkz+z9mNsI5d6SvXzvexcNFwa5G3EHnSEtO9M1SN+ccm/bX88fdR7jv+lLMjFkl2cybkMtHphcwPCPV6xLFB/oc6GY2CjjonHNmNhtIAI72ubI4Fy8XBbtaylnQbi49ln+gHT11ll9V11JRGWD3oVMMTklkyewiRg4dxJc/GJs/nGTg6jbQzewJYC4wwsxqgG8CyQDOuUeARcBnzawFaASWuO7WQkq34uWi4MVuOon1C8h/2nOEu3+6geagY1pRFt+9bQo3T80nY4Dd+CP+Ec4qlzu6Of8QbcsaJYLi5a5SPy3lrK1v5MmqAMXDBnPbjEKmFWVx95UlLJpZRNmoIV6XJ3FAQ4UBKp7uKo3lkfj55YYVVQH+sPswzsEn5ozmthmFpKcm8fWbJnpdosQRBfoApf0vYsPnflnNb3ccJC9zEJ+/rpSPzSykaNhgr8uSOKVAH6D8NBXhF6fPtvD/3jjAUxtrePhvZ5AzJJW/v3oMH59TzNWlOQNmd0OJXwr0ASyWpyL8wjlHdaCe1ZUBnt1Sx+lzQS7JSaeuvpGcIanMGTvc6xJFLlCgi3TCOYeZ8e6JJj76n68xKCmRmy/LY8nsImYUZw/o3Q0lfinQRUJaWx1/+ssRVlUGCAYdj3xiJnmZaTx6VzmzxwxjyKBkr0sUuSgFusS9Aw2NPFlVQ0VlgNr6RrIGJ7NoRuGFUfq8CSO9LlEkLAp0iUvnWloxg+TEBNZsquUHL77FVeNG8LUPjeeGiSMZlKzdDSX2KNAlruw5dH53w1q+tXASH56az52zi1k4NV/LDSXmKdDF94KtjjWb2qZUqt45TlKCcf2EkRRmt92klZ2eQnZ6isdVivSdAl18yTlHXUMTBVlpJBg88spfcMDyD43nthmF5AzpfHdDP+1w6ae+SHgU6OIrx0+fu7C7YV19Ixu+fj1pKYk8sXQOORmpF11u6KcdLv3UFwlfgtcFiETCWwdP8vknqrn8X37Hg89tZ1ByAstvnMD5/M4dMqjbteMX2+Ey1vipLxI+jdAlZh1oaCTY6ijMHkxTc5BX3zrMnZcXs3hWERPyhvb463W2GdrFjg9k8bJbp7yfAl1iSnOwld/tOEhFZYBX3jrMopmFfH/RVKYUZLL+gXl9Wm7Y1WPvEmPwrtB42q1T3qNAl5jx0Eu7+dlrb3Pk1DlGDk3lc9eO42MziwAwsz6vHe8szC92fCDTbp3xSYEuA9aZcy38bschbpqSR0KC0dDYzIzibJbMLuKa0hySEiN7CajgIo/DizXarTM+KdClTyK9NM45x9baBlZVBnhmcx2nzrYwKnMQs0qG8cCNE6K6KZbfRrXarTP+KNCl1yK9NC5w7Az/4/Eqdr57kkHJCdw4JY/F5UWUj84GiPoOhxrVSqxToEuv9fVB1q2tjtf3HuVEUzMLJucxKnMQuUMH8fE5o1k4LZ+hHuxuqFHtX9MNSrFDgS691tulce82NPHUxgAVVQECxxqZkDeUBZPzSE5M4PFPzo5GqdJLukEptijQpdd6szTu4Zf38K+/2UWrgyvGDuf+D5Yxf9KoaJYpfdDX38KkfynQpdfCuYi49/ApKqoC3Dm7mNHD05lelMVn517C7eVFjB6eHvUaNV3QN7pBKbYo0KXXurqIOH/SKJ7e2La74Ya3j5GYYFyaO4TRw9O5ctwIrhw3ol/q03RB3+kGpdiiQJc+6XgR8WxLkCu/8xJHT5+jZPhgvrpgPB+dUUDu0EH9XpumC/rOb0s5/U6BLn3ScKaZtZtr2VbXwPcXTSU1KZEvXl/KpSOHMHvMME8fpqzpgr7TUs7Y0m2gm9ljwM3AIefc5E7OG/Aj4EbgDHC3c25TpAuVgaO11fH6vqNUVAZ44c13OdfSyuSCoZw+20J6ahKfuKLE6xIBTRdEipZyxo5w7p3+GbDgIuc/BJSG/iwF/rPvZclAtqa6ljt/vJ6Xdh5iyawinvv8VTz3+atJTx1Yv/Atm19GWof9XTRdIH7W7d9A59yrZlZykSa3AI875xzwupllmVmec+5ApIoU7zQHW3l55yEqKgNcP3Ekd8wu5oOTRvLDhKl8aHLegH6YsqYLJN5EYkhVAATava8JHfurQDezpbSN4ikuLo7AR0u07DtymorKAE9trOHIqbPkDknlugm5AAwdlMxHphd6XGF4NF0g8aRff0d2zq0EVgKUl5fH3p6kPtcSbL2wg+H9T25hc6Cea8tyWTKriLllkd/dUEQiKxKBXgsUtXtfGDomMeLN2obQBc4DvPilD5CdnsI/3TKZERkpniw3FJHeiUSgPwPca2argMuBBs2fD3wnm5pZW13LqsoA2+pOkJrUtrvhmeYg2cDE/J4/wk1EvBXOssUngLnACDOrAb4JJAM45x4BnqdtyeIe2pYt3hOtYqVvnHOcaGohMy2Zo6fO8Y1fb2NS/lD+6ZZJLJxWQGZa/+9uKCKRY86jx2uVl5e7qqoqTz473hw60cSTG2t4sipA2agh/N9PlANt+6yMzcnwuDoR6Qkz2+icK+/s3MBaOCwR9dqeIzz2p7d5edchgq2Oy8cM48YpeRfOK8xF/EWB7jNvHzlNQXYayYkJbHj7GFtq6ll6zVhuLy9izIjo724oIt5RoPtAU3OQF948wKoNAdbvO8YjH5/JgsmjWHrNWD537TiSo7jcUNvTigwcCvQYduZcC999YSe/qq7lZFMLo4cPZtn8MmYUZwEwOCW6315tTysysCjQY0xDYzNvHTzJrJJhpCUnsn7vMeaNz+X2WUXMGTOchIT+291Q29OKDCwK9BjgnGPDvmNUVAZ4/s0DDEpOZP0D80hNSuT5+64msR9DvL142Z5W00oSKxToA9zLuw7x4LPb2XfkNENSk1g0s5DF5cWkhObFvQpziI/taTWtJLFEgT7AtARbeeWtw4wens643AyGDkomJyOVe68dx41T8khLGTi7G8bD02w0rSSxRIE+QLxz9DSrq9p2Nzx44iyfumoM37h5IjNHZ7P6M1d4XV6n4mF72niZVhJ/UKAPAJ/+RRXrth0kwWBuWS4P3lLEdeNzvS4rLH7fnjYeppXEPxToHthed4LfbH+X++aVYmZMys9kSkEmi2YWMSpTuxsOJPEwrST+oUDvJyeamnlmcx0VlQG21jaQkpjAwqn5jM3J4AvzSr0uT7oQD9NK4h8K9H5Qvf84d/z4dZqaWxk/agjf+vBEbp1eQNbgFK9LkzD4fVpJ/EOBHgWHT57l6U01ZA9OZvGsYibmD2XJrGI+Mr2AywozMfNuqaGI+JcCPUJagq28uvswqzYEeGnnIVpaHR+ZXsDiWcWkJiXyrYWTvC5RRHxOgR4hX3nqDdZU1zIiI4VPXTWGj5UXMS5X29OKSP9RoPdCU3OQddveZXVVgO999DIKswfzt3OK+eCkUcybkBvV3Q1FRLqiQO+B7XUnWF0V4FfVtTQ0NlOYnUbN8UYKswczc/Qwr8sTkTinQO+Gcw4z4/jpcyx86I8kmDF/8iiWzCriirH9u7uhiMjFKNA74Zyj8u3jVFQGOH7mHI/dPYvs9BQe+fhMZo7OJjtdyw1FZOBRoLdz+ORZ1myqoaIywN4jp0lPSeSW6QUEWx2JCcb1E0d6XaKISJfiPtCDrY5gqyMlKYHn3qjjOy/spHx0Np+Zewk3TckjPTXu/xOJSIyI27QKHDvD6qoAT1bV8KUbSlk8q5jbZhRydWmOlhuKSEyKq0B3zvHsGweoqNzPn/YcJcHgmktzKB6WDkBmWjKZackeVyki0jtxEegHTzQxcuggzIxH/7CXo6fO8eUbLmXRzEJtgyoivuHbQD/Z1MyzWw5QURVg54ETbPj69WSmJfPoXeWMyEjVckMR8Z2wAt3MFgA/AhKBR51z3+1w/m5gBVAbOvSQc+7RCNYZtrePnObhl/fw3BsHaGwOUjZyCF9dMJ6kUIDnDtV+4yLiT90GupklAg8DNwA1QKWZPeOc296haYVz7t4o1NitI6fO0nguSNGwwbS0tvLCm+9y6/R8Fs8qZqp2NxSROBHOCH02sMc5txfAzFYBtwAdA71fBVsdr+4+zOrKAC9uP8j8yaN4+M4ZjMsdQtX/up5ByQPnYcoiIv0hnEAvAALt3tcAl3fS7qNmdg3wFvAl51ygkzYR8dM/7WPlq3s50NDEsPQU7r6yhMWzii6cV5iLSDyK1EXRZ4EnnHNnzezTwM+B6zo2MrOlwFKA4uLiXn/YsdPnKB05hG/cPJHrJ4wkJUm7G4qIhBPotUBRu/eFvHfxEwDn3NF2bx8Fvt/ZF3LOrQRWApSXl7seVdrOl2+4tEfz4mura/VMSBHxvXCGtpVAqZmNMbMUYAnwTPsGZpbX7u1CYEfkSvxrPQ3z5Wu2UlvfiANq6xtZvmYra6tru/13RURiSbeB7pxrAe4F1tEW1Kudc9vM7EEzWxhq9gUz22ZmW4AvAHdHq+CeWrFuF43Nwfcda2wOsmLdLo8qEhGJjrDm0J1zzwPPdzj2j+1eLweWR7a0yKirb+zRcRGRWOX7q4ld3dqvW/5FxG98H+jL5peR1mEZY1pyIsvml3lUkYhIdPh2L5fzzq9m0SqXgUurkEQiw/eBDm2hroAYmM6vQjp/4fr8KiRA3zORHvL9lIsMbFqFJBI5CnTxlFYhiUSOAl08pVVIIpGjQBdPaRWSSOTExUVRGbh6swpJq2JEOqdAF8/1ZBWSVsWIdE1TLhJTtCpGpGsKdIkpWhUj0jUFusQUrYoR6ZoCXWKKVsWIdE0XRSWmaG8eka4p0CXmaG8ekc5pykVExCc0Qpf30U07IrFLgS4X6KYdkdimKRe5QDftiMQ2BbpcoJt2RGKbAl0u0E07IrFNgS4X6KYdkdimi6JygW7aEYltCnR5H920IxK7NOUiIuITCnQREZ8IK9DNbIGZ7TKzPWb2tU7Op5pZRej8ejMriXilIiJyUd0GupklAg8DHwImAneY2cQOzT4FHHfOjQN+CHwv0oWKiMjFhTNCnw3scc7tdc6dA1YBt3Rocwvw89Drp4B5ZmaRK1NERLoTziqXAiDQ7n0NcHlXbZxzLWbWAAwHjrRvZGZLgaWht6fMbKDcUz6CDrX6VDz0Mx76COqnn/S0j6O7OtGvyxadcyuBlf35meEwsyrnXLnXdURbPPQzHvoI6qefRLKP4Uy51AJF7d4Xho512sbMkoBM4GgkChQRkfCEE+iVQKmZjTGzFGAJ8EyHNs8Afxd6vQh4yTnnIlemiIh0p9spl9Cc+L3AOiAReMw5t83MHgSqnHPPAD8BfmFme4BjtIV+LBlw00BREg/9jIc+gvrpJxHro2kgLSLiD7pTVETEJxToIiI+ETeBbmaPmdkhM3vzIm3mmtlmM9tmZq/0Z32R0l0/zWxZqI+bzexNMwua2bD+rrMvwuhjppk9a2ZbQt/Le/q7xkgIo5/ZZvYrM3vDzDaY2eT+rrGvzKzIzF42s+2h79V9nbQxM/v30NYib5jZDC9q7Ysw+znezP5sZmfN7P5efZBzLi7+ANcAM4A3uzifBWwHikPvc72uORr97ND2w7StSPK87gh/Lx8Avhd6nUPbhfoUr+uOQj9XAN8MvR4P/M7rmnvRxzxgRuj1EOAtYGKHNjcCLwAGzAHWe113lPqZC8wC/hm4vzefEzcjdOfcq7T9xe7KncAa59z+UPtD/VJYhIXRz/buAJ6IYjlREUYfHTAktP1ERqhtS3/UFklh9HMi8FKo7U6gxMxG9kdtkeKcO+Cc2xR6fRLYQdud5+3dAjzu2rwOZJlZXj+X2ifh9NM5d8g5Vwk09/Zz4ibQw3ApkG1mvzezjWZ2l9cFRZOZDQYWAE97XUsUPARMAOqArcB9zrlWb0uKii3AbQBmNpu2W8ILPa2oD0K7tE4H1nc41dn2IzH7FJaL9LPPFOjvSQJmAjcB84FvmNml3pYUVR8G/uScC3c0H0vmA5uBfGAa8JCZDfWyoCj5Lm2j1c3A54FqIOhpRb1kZhm0DS6+6Jw74XU90RLtfuoRdO+pAY46504Dp83sVWAqbXNdfrSEGJxuCdM9wHdd28TkHjPbR9sc8wZvy4qsUCDcA20XDoF9wF5Pi+oFM0umLeR+6Zxb00mTcLYfGfDC6GefaYT+nl8DV5lZUmg64nLa5rl8x8wygQ/Q1mc/2g/MAwjNKZcRg0HXHTPLCm3HAfD3wKuxNroN/SD6CbDDOfeDLpo9A9wVWu0yB2hwzh3otyIjIMx+9v1zQldXfc/MngDm0rZV5UHgm0AygHPukVCbZbSNeFqBR51z/+ZFrX0RZj/vBhY452Jtiwag+z6aWT7wM9pWFhhto/X/8qTYPgijn1fQ9hwCB2wDPuWcO+5Ntb1jZlcBf6DtWsf56xwPAMVwoZ9G23WRBcAZ4B7nXJUH5fZamP0cBVQBQ0NtTtG2EibsH9JxE+giIn6nKRcREZ9QoIuI+IQCXUTEJxToIiI+oUAXEfEJBbqIiE8o0EVEfOL/A2QEAxnSWdxIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(facedim, penaltyMin)\n",
    "yAxBline(3.18898271, -4.50450342)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from statsmodels.formula.api import ols"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}