Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 2.10.13\n",
"\n",
"Discrete data: Table 2.2 gives the number of fatal accidents and deaths on scheduled airline flights per year over a ten-year period. \n",
"We use these data as a numerical example for fitting discrete data models. \n",
"\n",
"1. Assume that the numbers of fatal accidents in each year are independent with a Poisson(theta) distribution. Set a prior distribution for theta and determine the posterior distribution based on the data from 1976 through 1985. Under this model, give a 95% predictive interval for the number of fatal accidents in 1986. You can use the normal approximation to the gamma and Poisson or compute using simulation.\n",
"2. Assume that the numbers of fatal accidents in each year follow independent Poisson distributions with a constant rate and an exposure in each year proportional to the number of passenger miles flown. Set a prior distribution for theta and determine the posterior distribution based on the data for 1976–1985. (Estimate the number of passenger miles flown in each year by dividing the appropriate columns of Table 2.2 and ignoring round-off errors.) Give a 95% predictive interval for the number of fatal iaccidents in 1986 under the assumption that 8 × 10 11 passenger miles are flown that year.\n",
"3. Repeat (1) above, replacing ‘fatal accidents’ with ‘passenger deaths.’\n",
"4. Repeat (2) above, replacing ‘fatal accidents’ with ‘passenger deaths.’\n",
"5. In which of the cases above does the Poisson model seem more or less reasonable? Why? Discuss based on general principles,without specific reference to the numbers in Table 2.2. Incidentally, in 1986, there were 22 fatal accidents, 546 passenger deaths, and a death rate of 0.06 per 100 million miles flown. We return to this example in Exercises 3.12, 6.2, 6.3, and 8.14.\n",
"\n",
"|Year |Fatal accidents |Passenger deaths |Death rate\n",
"|---|---|---|---| \n",
"|1976 | 24 | 734 | 0.19 \n",
"|1977 |25 |516 |0.12 \n",
"|1978 |31 |754 |0.15 \n",
"|1979 |31 |877 |0.16 \n",
"|1980 |22 |814 |0.14 \n",
"|1981 |21 |362 |0.06 \n",
"|1982 |26 |764 |0.13 \n",
"|1983 |20 |809 |0.13 \n",
"|1984 |16 |223 |0.03 \n",
"|1985 |22 |1066 |0.15 \n",
"\n",
"+ Table 2.2 Worldwide airline fatalities, 1976–1985.\n",
"+ Death rate is passenger deaths per 100 million passenger miles.\n",
"+ Source: Statistical Abstract of the United States.\n",
"\n",
"Gelman, Andrew; Carlin, John B.; Stern, Hal S.; Dunson, David B.; Vehtari, Aki; Rubin, Donald B.. Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science) (Page 60). CRC Press. Kindle Edition. \n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"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>Deaths</th>\n",
" <th>Fatal</th>\n",
" <th>Rate</th>\n",
" <th>year</th>\n",
" <th>Miles</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>734</td>\n",
" <td>24</td>\n",
" <td>0.19</td>\n",
" <td>1976</td>\n",
" <td>3863.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>516</td>\n",
" <td>25</td>\n",
" <td>0.12</td>\n",
" <td>1977</td>\n",
" <td>4300.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>754</td>\n",
" <td>31</td>\n",
" <td>0.15</td>\n",
" <td>1978</td>\n",
" <td>5027.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>877</td>\n",
" <td>31</td>\n",
" <td>0.16</td>\n",
" <td>1979</td>\n",
" <td>5481.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>814</td>\n",
" <td>22</td>\n",
" <td>0.14</td>\n",
" <td>1980</td>\n",
" <td>5814.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>362</td>\n",
" <td>21</td>\n",
" <td>0.06</td>\n",
" <td>1981</td>\n",
" <td>6033.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>764</td>\n",
" <td>26</td>\n",
" <td>0.13</td>\n",
" <td>1982</td>\n",
" <td>5877.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>809</td>\n",
" <td>20</td>\n",
" <td>0.13</td>\n",
" <td>1983</td>\n",
" <td>6223.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>223</td>\n",
" <td>16</td>\n",
" <td>0.03</td>\n",
" <td>1984</td>\n",
" <td>7433.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>1066</td>\n",
" <td>22</td>\n",
" <td>0.15</td>\n",
" <td>1985</td>\n",
" <td>7107.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Deaths Fatal Rate year Miles\n",
"0 734 24 0.19 1976 3863.0\n",
"1 516 25 0.12 1977 4300.0\n",
"2 754 31 0.15 1978 5027.0\n",
"3 877 31 0.16 1979 5481.0\n",
"4 814 22 0.14 1980 5814.0\n",
"5 362 21 0.06 1981 6033.0\n",
"6 764 26 0.13 1982 5877.0\n",
"7 809 20 0.13 1983 6223.0\n",
"8 223 16 0.03 1984 7433.0\n",
"9 1066 22 0.15 1985 7107.0"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy.stats import poisson, gamma, gengamma\n",
"from sklearn.linear_model import LinearRegression\n",
"import pystan\n",
"airline_df=pd.DataFrame(dict({'year':[x for x in range(1976,1986)],'Fatal':[24,25,31,31,22,21,26,20,16,22],'Deaths':[734,516,754,877,814,362,764,809,223,1066],'Rate':[.19,.12,.15,.16,.14,.06,.13,.13,.03,.15]}))\n",
"airline_df.set_index('year')\n",
"airline_df['Miles']=np.round(airline_df['Deaths']/airline_df['Rate'],0)\n",
"airline_df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Deaths vs miles (out of order -- this is part (3))"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_d23d1d0b2e2b2dabfc71364b95c0d024 NOW.\n"
]
}
],
"source": [
"stan_code='''\n",
"data {\n",
" int deaths[10];\n",
"}\n",
"parameters {\n",
" real<lower=0> theta ; \n",
"}\n",
"model {\n",
"\n",
" // this prior tends to say the number of deaths is between 600 and 800\n",
" theta~gamma(691,1);\n",
" deaths~poisson(theta);\n",
"}\n",
"\n",
"'''\n",
"sm_simple=pystan.StanModel(model_code=stan_code)\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Inference for Stan model: anon_model_d23d1d0b2e2b2dabfc71364b95c0d024.\n",
"4 chains, each with iter=2000; warmup=1000; thin=1; \n",
"post-warmup draws per chain=1000, total post-warmup draws=4000.\n",
"\n",
" mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\n",
"theta 691.9 0.25 8.18 675.85 686.54 691.81 697.26 708.22 1103 1.0\n",
"lp__ 4.2e4 0.02 0.76 4.2e4 4.2e4 4.2e4 4.2e4 4.2e4 1461 1.0\n",
"\n",
"Samples were drawn using NUTS at Mon Apr 16 10:55:25 2018.\n",
"For each parameter, n_eff is a crude measure of effective sample size,\n",
"and Rhat is the potential scale reduction factor on split chains (at \n",
"convergence, Rhat=1).\n"
]
}
],
"source": [
"print(deaths)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From the stan output reported above, the 95% interval for the poisson rate is (675,707)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_9a1f6e45cdebbb8b75bf8564a74768dd NOW.\n"
]
}
],
"source": [
"stan_code='''\n",
"data {\n",
" int deaths[10];\n",
" vector[10] miles;\n",
"}\n",
"parameters {\n",
" real<lower=0> theta ; \n",
"}\n",
"model {\n",
" // 691 is the average number of deaths\n",
" // 5715 is the averaage number of miles\n",
" // so gamma prior is about 691 deaths per 5715 miles\n",
" theta~gamma(691,5715);\n",
" deaths~poisson(miles*theta);\n",
"}\n",
"generated quantities {\n",
" int predicted[10] ; \n",
" int answer ; \n",
" \n",
" for(i in 1:10)\n",
" predicted[i]=poisson_rng(miles[i]*theta);\n",
" answer=poisson_rng(8000*theta) ; \n",
"}\n",
"\n",
"'''\n",
"sm_weights=pystan.StanModel(model_code=stan_code)\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jteitelbaum/anaconda3/lib/python3.6/site-packages/pystan/misc.py:399: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
" elif np.issubdtype(np.asarray(v).dtype, float):\n"
]
}
],
"source": [
"deaths_wtd=sm_weights.sampling(data=dict({'deaths':airline_df['Deaths'],'miles':airline_df['Miles']}))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Inference for Stan model: anon_model_37d152c75a3e2c95f06248eb3357d905.\n",
"4 chains, each with iter=2000; warmup=1000; thin=1; \n",
"post-warmup draws per chain=1000, total post-warmup draws=4000.\n",
"\n",
" mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\n",
"theta 0.12 3.7e-5 1.4e-3 0.12 0.12 0.12 0.12 0.12 1459 1.0\n",
"predicted[0] 467.72 0.36 22.23 425.0 453.0 468.0 483.0 512.0 3796 1.0\n",
"predicted[1] 520.0 0.41 23.67 474.0 503.0 520.0 536.0 566.0 3272 1.0\n",
"predicted[2] 608.23 0.43 25.42 559.0 591.0 608.0 625.0 658.0 3426 1.0\n",
"predicted[3] 663.18 0.45 26.81 611.0 645.0 663.0 681.0 716.0 3581 1.0\n",
"predicted[4] 704.22 0.47 27.73 651.0 686.0 704.0 722.0 761.0 3537 1.0\n",
"predicted[5] 729.83 0.47 28.63 672.0 710.0 730.0 749.0 787.0 3657 1.0\n",
"predicted[6] 711.57 0.48 28.39 657.4 692.0 711.0 730.0 769.0 3542 1.0\n",
"predicted[7] 753.78 0.46 28.29 698.4 734.0 754.0 773.0 810.59 3732 1.0\n",
"predicted[8] 898.95 0.56 31.12 840.0 877.0 898.0 920.0 962.0 3120 1.0\n",
"predicted[9] 860.63 0.51 30.58 802.0 840.0 860.0 881.0 922.0 3531 1.0\n",
"lp__ 3.6e4 0.02 0.72 3.6e4 3.6e4 3.6e4 3.6e4 3.6e4 1735 1.0\n",
"\n",
"Samples were drawn using NUTS at Mon Apr 16 10:56:28 2018.\n",
"For each parameter, n_eff is a crude measure of effective sample size,\n",
"and Rhat is the potential scale reduction factor on split chains (at \n",
"convergence, Rhat=1).\n"
]
}
],
"source": [
"\n",
"print(deaths_wtd)\n",
"predicted=deaths_wtd.extract()['predicted']\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4XNWZ+PHvO0Wj3i1bzdhyt2kGuRAwYJppDiaADSEBsiFZErYk2ZDApm3ab5OQ3WyySTabDUkgIUQCjIMpcegdbBljG+OKsVGXLKtLM5pyfn/cK1llxtJYo2a9n+fRo5lzyzn3ajTvvefcc44YY1BKKaWGyjHWBVBKKTWxaOBQSikVFQ0cSimloqKBQymlVFQ0cCillIqKBg6llFJR0cChEBEjIrPHuhzDISJPi8it9uvbROTVsS4T9C1XhOW/F5HvjWaZIhGRQyJyyRDWm2F/ZlyjUa4w+fecMxFZISJ7T3A/vxKRb8S2dJODBo5xxv7n7RSRVhFpEpHXReQOEYnJ30pEXhSR22Oxr9Fgf0HV9v6SEhGXiNSJSE8nJGPMFcaY+8emlJH1Ltd4CmgnC2PMK8aYeYOtF+7cG2PuMMZ8d+RKd/LSwDE+rTbGpACnAD8AvgrcN7ZFGlNNwBW93l8JNI5RWVQMjdVdixoeDRzjmDGm2RjzOLAOuFVETgUQEY+I/FhEPrSvxn8lIgn2sgwReUJE6kWk0X5dYC/7PrAC+LmItInIz3tld4mI7Le3+YWIiL3NbBF5SUSaReSIiJSEK6uI/FVE/qFf2nYR+ZhYfmLfJTSLyI7uYxmiPwC39Hp/C/BAv7wi3kmJyHwReUZEjorIXhFZ22vZlSLynn2HVykiXx6sMCIy074bdNjvfyMidb2W/1FEvtC7XCKyAPgVcI597pt67TJDRJ60y/CWiMyKkG93FdGnRKTc/lvdISJL7HPa1PtvKiIOEfm6iBy2z/0DIpLWa/kn7WUNIvK1fnk5RORuEXnfXl4qIpmDnRt720Mico99XhtF5HciEm8vu1BEKkTkqyJSA/zOTr9aRN7pdZd9eq/9LRaRt+3zUwLE91p2oYhU9HpfKCLr7c9/g4j8PNK5l37VhCLyGRE5YH9OHheRvF7LjH2uT/h/5KRijNGfcfQDHAIuCZP+IfA5+/V/AY8DmUAKsBH4d3tZFnAdkGgvexjY0Gs/LwK399u3AZ4A0oHpQD1wub3sIeBrWBcZ8cB5Ecp9C/Bar/cLse4UPMAqYKu9fwEWALlDPB8GOBWotbdPt1+fan18Bx4XcBvwqv06CSgHPgW4gLOAI8Aie3k1sMJ+nQGcNcRyfQicbb/eCxwEFvRatvh45eq1n98DR4GldvkeBP4cIc8Z9vn4lf23uAzwAhuAHCAfqAMusNf/O+AAUAQkA+uBP/T6+7QB59t/o/8EAtifPeALwJtAgb38f4GH+pXDdZzP8LtAIdZn9DXge/ayC+18fmjvN8H+m9QBywAncKu9Dw8QBxwGvgi4gesBf7/9VdivncB24Cf2373n83qcc9+9n4vsz8VZdr7/Dbwcy/+Rk+lH7zgmjiog077K+QzwRWPMUWNMK/D/gBsBjDENxphHjTEd9rLvAxcMYf8/MMY0GWM+BF4AzrTT/VhVZnnGGK8xJlId/WPAmSJyiv3+ZmC9McZn7yMFmA+IMWa3MaY6imP3YgXHdfZxPm6nDcXVwCFjzO+MMQFjzNvAo1hfQN3Ht1BEUo0xjfbyoXgJuEBEptnvH7HfzwRSsb7Ahmq9MWazMSaAFTjOHGT979p/i78B7Vhf6HXGmErgFWCxvd7NwH8aYw4aY9qAe4Abxaoeuh54whjzsv03+gYQ6pXH3wNfM8ZU2Mv/Dbhehl619HNjTLkx5ijWZ/CmXstCwLeMMT5jTCfW5/l/jTFvGWOCxmoT8gHL7R838F/GGL8x5hFgS4Q8lwJ5wF3GmPZBPq/93Qz81hjztn2892Ddoczotc5w/0dOGho4Jo58rCvTKVh3E1vt2/om4K92OiKSKCL/a1dBtAAvA+ki4hxk/zW9XndgXaECfAXrLmGziOwSkb8Lt7EdpJ7EDmD27wftZc8DPwd+AdSKyK9FJDWKYwerauoWwlRTDeIUYFn3ubLP181A9xf+dVhtJoft6oZzhrjfl7Cuds/HOscvYgXoC4BXjDGhiFsOFOncR1Lb63VnmPfd2+dhXa13O4x1VzPVXlbevcAY0w409Fr3FOCxXudsNxC0tx2K8l6vD9v5das3xvQO/KcA/9Lvb1Rob5MHVBr70r7X/sIpBA7bAThafc6VHWgbsP7vug3rf+RkooFjAhCRJVgf4Fexbqc7sapa0u2fNGNM94f4X4B5wDJjTCrWFxtYH2ywbrmHzBhTY4z5jDEmD+sq9JcS+dHdh4Cb7C/fBKyrsu79/MwYczawCJgL3BVNObCupHOxvriiuaIrB17qda7SjTHJxpjP2eXaYoy5BquqZwNQOsT9voTVXnSh/fpV4FyswPFShG1GeyjqKqwv5W7TsaqJarGq6Aq7F4hIIlY1Z7dy4Ip+5y3evqsZisJer6fbZenW/zyUA9/vl1eiMeYhu5z53e0JvfYXTjkwPcJd0WDnvs+5EpEkrPMx6PFG+T9yUtDAMY6JSKqIXA38GfijMWanfSX7f8BPRCTHXi9fRFbZm6VgBZYmuzHzW/12W4tV5z3UMtwgduM61pNMBuvKM5ynsP75vgOUdF912423y0TEjVW14j3OPsKyrzhXAx/td/U5mCeAuXZDsNv+WSIiC0QkTkRuFpE0Y4wfaOldLrtB9MII5dmPdZ4/gVUX3oJ1bq8jcuCoBQpEJC6K8g/HQ8AXxWrMT8aq0iyxr8gfAa4WkfPs8nyHvt8HvwK+3131KCJTROSaKPK+U0QK7M/gvwLHazD+P+AO+zMiIpIkIleJSArwBlaw+yexHsP+GFaVVDibsQLND+x9xIvIufaywc79n4BPiciZIuLBOldvGWMODXagUf6PnBQ0cIxPG0WkFesK6mtYDZef6rX8q1iNnm/a1VHPYt1lgNVwnoB1Z/ImVjVWbz/FqqtuFJGfDaEsS4C3RKQNq23hn40xH4Rb0a4bXg9cgvWP2C0V68uhEas6oAH4MYCI/KuIPD2EcmCM2WWM2TWUdXtt04rViHwj1lVlDccaZgE+CRyyz+MdWIEA+4ugDdh5nN2/BDTYdd7d7wXYFmH954FdQI2IHInmOE7Qb7GeSHsZ+AArYP8jWOcSuBPr71SN9bep6LXtT7H+3n+zP4tvYjVeD9WfgL9hPTRwEIjYydEYU4bVzvFzuxwHsBqzMcZ0AR+z3zditXOtj7CfINbFxWysBxQq7PVhkHNvjHkOq53nUazzMYtj1a6DGfL/yMlCort4U2pyEJFPYFUH3jPWZZloROQQ1pNkz451WdTI0M43SoVhjPnjWJdBqfFKq6qUUkpFRauqlFJKRUXvOJRSSkXlpGzjyM7ONjNmzBjrYiil1ISydevWI8aYKYOtd1IGjhkzZlBWVjbWxVBKqQlFRCL1yu9Dq6qUUkpFRQOHUkqpqGjgUEopFRUNHEoppaKigUMppVRUTsqnqpRSk8uGbZXcu2kvVU2d5KUncNeqeaxZnD/4huqEaOBQSk1oG7ZVcs/6nXT6rZHMK5s6uWe9NaixBo+RoVVVSqkJ7d5Ne3uCRrdOf5B7N+0doxKd/DRwKKUmtKqmzqjS1fBp4FBKTWh56QlRpavh08ChlJrQ7lo1jwS3s09agtvJXavmRdhCDZc2jiulJrTuBnB9qmr0aOBQSk14axbna6AYRVpVpZRSKioaOJRSSkVFA4dSSqmoaOBQSikVFQ0cSimloqKBQymlVFQ0cCillIqKBg6llFJRGbHAISK/FZE6EXm3V1qmiDwjIvvt3xl2uojIz0TkgIjsEJGzem1zq73+fhG5daTKq5RSamhG8o7j98Dl/dLuBp4zxswBnrPfA1wBzLF/Pgv8D1iBBvgWsAxYCnyrO9gopZQaGyMWOIwxLwNH+yVfA9xvv74fWNMr/QFjeRNIF5FcYBXwjDHmqDGmEXiGgcFIKaXUKBrtNo6pxphqAPt3jp2eD5T3Wq/CTouUPoCIfFZEykSkrL6+PuYFV0opZRkvjeMSJs0cJ31gojG/NsYUG2OKp0yZEtPCKaWUOma0A0etXQWF/bvOTq8ACnutVwBUHSddKaVGzaM1Ryl+fRe5L7xD8eu7eLSmfy385DLageNxoPvJqFuBv/RKv8V+umo50GxXZW0CLhORDLtR/DI7TSmlRsWjNUf58t5yKnx+DFDh8/PlveWTOniM5OO4DwFvAPNEpEJEPg38ALhURPYDl9rvAZ4CDgIHgP8DPg9gjDkKfBfYYv98x05TSqlR8e8Hq+kM9a0h7wwZ/v1g9RiVaOyN2EROxpibIiy6OMy6Brgzwn5+C/w2hkVTSqkhq/T5o0qfDMZL47hSSo1L03xhn8eJmD4ZaOBQSqnj+PweL/HBQJ+0+GCAz+/xjlGJxp4GDqWUOo7r217kR3vvpcBbg5gQBd4afrT3Xq5ve3Gsi9bHgbpWvv/ke/zhjUMjnteItXEopdRo+fqO57m/TvA70nGHmrg1x/C90y+Kyb7jXPeT1NFMZuUOulxOMgNBkjqaiHNtB+6JSR4nqs0X4IntVZSUlbPtwyZcDuGWc2aMeL4aOJRSE9rXdzzPffVJGKcHAL8zg/vqfbDj+ZgEj+edLfyuazrnvpJBktdJe3yQ++amEoorZ/Ww9x49YwxbDzdSsqWcJ3dW09EVZHZOMl+7cgHXnpVPdrJnxMuggUMpNaHdXyc9QaObcXi4v66D78Vg/w95T+Gc3Rk4glbNfrLXxTm7snhogWNUA0ddq5f1b1dSWlbOwfp2kuKcfPSMPNYuKWRxYToi4QbaGBkaOJRSE5rfkR5VerROez+1J2h0cwQdnPZ+akz2fzyBYIgX99ZTUlbO83vqCIYMxadkcMf1s7jqtFySPGPzFa6BQyk1obn8DQTissOmx0JcZ/ivyUjpsfDBkXZKy8p5dGsFda0+spPjuP28mdxQXMjsnOQRy3eoNHAopSa0pe+U8Hrx7eDoVV0V8rH0nRJYdcmw9+9O9pM0tZO8ZXW4kwP421xUvZVDe23CsPfdW0dXgKd21lC6pZzNh47iEFg5L4e1Swq5aH4Obuf4eQhWA4dSakJrSt5MylGhPX0tIWcWjmADSU2lNCVvjsn+C1ZUk5zbidNtdfiLSwlQeEE1bdXDDxzGGLZXNFOypZyN26to8wWYkZXIVy6fx3VnFTA1NX7YeYwEDRxKqQmtPilEfMcbxHe80S89Nvt3TzM9QaOb021wTzvxnuNH27t4bFslpVvK2VvbSrzbwZWn5bKuuJClMzNHtaH7RGjgUEpNaHGhBLqcnWHTYyHBPXDfx0uPJBgyvLK/nofLKvjbezX4g4YzCtL4/rWnsvqMPFLj3bEo7qjQwKGUmtCuTOvC7LqeU2rOpSvkIM4R4vC015BFj8Vk/0ckmykcCZs+FOVHO3i4rJxHtlZQ1ewlI9HNJ5afwrolhcyfNvJPZo0EDRxKqQntnMOrebp+Of+d1EWLw5AaEi6oX84VhwODbzwET/vWsC7u93jk2P58xsXTXWu4McI2Xn+QTbtqKC0r57UDDYjAijlT+NpVC7lkYQ4elzMmZRsrGjiUUhPa3w6dy9OJAQJ2s0CL0/B0YgDnoXO5Mgb7b3pnNyWz3VydFiTDaWgMCk80u2k6sBuu6LvuriqroXvDtkpavAHy0xP44iVzub64gPz02D6FNZY0cCilJrTnPSFmtu3jI41vkRJso9WZzOsZy3g+cV5M9l+bvoODHcLbHX2/+JPSdwDQ3OHnL9srKdlSzq6qFuJcDlYtmsa64kI+MisLh2N8N3SfCA0cSqkJbWrHXi5teBGnCQKQGmzj0oYXeQYgBvcc7QnBAWnGCM2hGfzzn7fx13dr8AVCLMhN5d9WL2TN4nzSE+OGne94poFDKTWhXdj0Wk/Q6OY0QS5seg344rD37+hKJuRpAyDkT8PffDb+prMx/iyer69jbXEh65YUcmp+2rDzmig0cCilJjRPIPyESpHSo9VWdxVmSohASzHSYBBAko/gcH7Alq99jnj3xG7oPhEaOJRSE5o4UiDUGj59GPbVtlKypZxA+2KkDcQjBIuSCeYngmcaznebJmXQAA0cahLZsK2Sezftpaqpk7z0BO5aNY81i/PHulhqmELJZ+JoeQPo/fiti1DymVHvq9Xr54kd1ZRsKeed8ibcToFMD12nJBPK9kCvHt2OOZOnaqo/DRxqUtiwrZJ71u+k02/VhVc2dXLP+p0AGjwmuKy8MpqCKzGdb1p3Ho4UJGE5WXmvDml7YwxbDlkTIz21s5pOf5A5Ocl8/aoFXLs4n0Wb3+sTMLp1JUzOuw3QwKEmiXs37e0JGt06/UHu3bRXA8cEN8PhovGsDSTMaCEuoYmuznQ6D9WQ0TDtuNvVtXp5dGslD5eVc/CINTHSmsV53FDcd2KktI4QzUkDg0RaR2hEjmci0MChJoWqpvDjCkVKVxNHelIrwbkVOFzWF7knsQn33BbSdwwc5TAQDPHC3npKtpTzwl5rYqQlMzL43IWzuOr0XBLjBn4lrtzZwBPF2QRcx4Y1dwVCrNzZAFeP3HGNZxo41KSQl55AZZggkXcS9eadrGrnHiLO1ffq3+EKUTv3UM/7g/VtlJZV8OjbFdS3+shO9nD7ipmsLS5k1pTjT4y0vP1RWvcu5a1ZC2jzJJDs62TZ+7tZ3r4ZuGwEjmj808ChJoW7Vs3r08YBkOB2cteq2PQuVmMnLqmTN6rO5rEDq2nwZpAV38i1szdyVs5OHi4rp7SsnC2HGnE6hJXzprC2uJCVUUyMVNuVhqOmHU9FDQHi8NCFw9VObZc2jit1Uutux9Cnqk4+b+0/mwc+vImukNVbu8GbyW/f/QT3SxC/2cHM7CS+evl8rjsrn5wTmBjpXXc+rwdmEMRq52jHw+uBGeA+8fk4JjoNHGrSWLM4XwPFSejhD9bQJX2H+AjhxBUKUnrHOSyZkTGsiZG2BQt7gka3IE62BQtPeJ8TnQYOpdSEFAwZXt5fT6OEn9OiS9wsnZk57HzaTByEiTtt5uQej+p4NHAopSaUDxs6eHirNTFSdXPkYUWS6IpJflPFQS0Dq6WmytDaSE5Gk/fIlTpZ7CiFn5wK/5Zu/d5ROtYlijmvP8hf3qnk4//3Juff+wI/f+EAc6em8Mubz2JBRh2m3zeZccApmU0xyftz+PD0S/PY6ZOV3nEoNZHtKOWF+57gvc7vYiQVqWph4aHHWPlp4PS1Y126YXu3spnSsmMTIxVkJPClS+dy/dkFPY9S311xCv6j8bj2tyLeICbeSWBOChWZ0TeEh3Od6/eYwOf4FSHqMOQg3IGD61y/B66NSR4TjQYOpSawF+5/lt3ej4MjzqqGlzTr/f3rWfkfEzNwNHf42fCONTHSe9XWxEhXnGpNjLS8aODESPWeDEyeg668vh3+6k1s+uhUu4QbzK+5PHgzQbJxcoQU54NUuYThDaM4cWngUGoC29d2GabfE0VG4tjXdhkrx6hMJyIUMrxxsIGSLeX8dVcNXYEQi/JS+c41i7jmjHzSEt0Rt030dtKeMLCXeKI3NqMCuIJ34jNvMCXuqzjlCEGTTXPgFlzBc2Ky/4lIA4dSMTTaI/AGJD2q9PGmqqmTh8sqeHhrORWNnaTGu7hpSSE3FA99YqQ5+wLsOTWA13ns6yw+GGDOvkAsJgDE3eWhU1bS2dU3FLuD2o9DKTVMYzECb5KjgfbQlLDp45UvEOTZ9+ooKSvnlf31GAPnzs7irlXzWLVoWtRzXHyvxrANP7+Y66A2XpjqNdy5z8/imth8sXcEWkhyDwxiHYGWmOx/ItLAoVSMxGIE3h07dvDcc8/R3NxMWloaF198MaeffnrE9Zcn/4GXWj5PgGMNwS68LE/+A7DuhI5jpOytsSZGemxbBY0dfnLT4vnHlbO5obiQwszEE95vDg6uqAlwRU2gT3ooRg+N7mh6mSVZl+NyHKsuC4T87Gh6mXmTdJTDMQkcIvJF4HbAADuBTwG5wJ+BTOBt4JPGmC4R8QAPAGcDDcA6Y8yhsSi3Uscz3BF4d+zYwcaNG/H7/QA0NzezceNGgIjBIzl+Nyv5JW+0fYK2UDbJjiOck/xHEuN3n8ARxF6r18/G7dWUlJWz3Z4Y6dKFU1lbXMiKOVNwOk68R3e32ngHud6Bdxe18Q6mD3vv8GHbe2AMp2dcQKIrlY5ACzsaX+LD9vFxjsfCqAcOEckH/glYaIzpFJFS4Eas2sifGGP+LCK/Aj4N/I/9u9EYM1tEbgR+yHi7lFKK4Y/A+9xzz1HYmcmSwCySiacNL1sC7/Pcc89FDBw/CqzlBwm/4dbEV3rSOkwcd/tv52cndhjDZoxh8wdHKSmzJkby+kPMm5rCN65eyLWL88lMim2P6zcyDNdWB5Few4IYgryR4WRJDPafkj2FD4/sHhAoUrIHVhFOFmPVAdAFJIiIC0gEqoGLgEfs5fcDa+zX19jvsZdfLMMZeEapEXLXqnkk9Kufj2YE3qyGeFYEFpBCAoKQQgIrAgvIaojcH+EpzuVu/+1UhLIJGaEilM3d/tt5inOHdSwnoq7Fyy9fPMBF//ES6379Jn/bVcu1iwvYcOe5/PULK/j0eTNjHjQALj3S3idoAAhOLj3SHpP9r7jxFlxxfbsAuuI8rLjxlpjsfyIa9TsOY0yliPwY+BDoBP4GbAWajDHdlZQVQHelcD5Qbm8bEJFmIAs4MqoFV2oQ4Ubg/XzRNJofPsQv/ncvyZkezrlmFnOXhZ+ZblloNocddZS5DtImXpJNPMWBIpaFZkfMMzQnlY3eDJ5NLsQRSiTkyKKzLQOJDz9+U6z5gyFe2FNHaVk5L+ytJxgyLJ2ZyZ0rZ3PladPCTowUa8n+8O0jkdKjtWCF9TTVK39+gNaGI6RkZbPixlt60iejsaiqysC6i5gJNAEPA1eEWbW70jLc3cWACk0R+SzwWYDp02NRs6lU9HqPwLvvrRpeeHAPgS5rkqG2oz5eeHAPQNjgUUUje/KeYkHR23g87fh8Sew+eBam6koi3bMEs3fy7T3PclPNh8RxlC7aeSj5Wb5VkAScPxKHCMD79W2Ubinn0bcrOdLmIyfFw2fPL+KGswsoGmRipFgLxDfg9maHTY+VBStWTupA0d9YNI5fAnxgjKkHEJH1wEeAdBFx2XcdBUCVvX4FUAhU2FVbacDR/js1xvwa+DVAcXHx5H3AWo0bb/zlfZ7MDLHD+DC+EOJxcLp4iP/L+2EDx8G8v3Ik1ctvNt/VMyHRmqInOchfWUn4XuDf3vs0iYlpnHvqL6n05JDvq+MrB+/j23ufhtV3xvR42n0BntxZTemWcsoOWxMjXTQ/h3XFhVw4bwquIU6MFGt1sx8l973bcISOVSeFHD7qZj/KzEk6JMhIG4vA8SGwXEQSsaqqLgbKgBeA67GerLoV+Iu9/uP2+zfs5c8bYzQwqHHv0TgfO32dYOzbZl+IndIJHusD3d+HSX4e2rOuz4REf9izjpvmrI+Yhyshjbvm/iM+l9UOUhE/jbvm/gvf3/ffMTkGYwzbypso3VLOxu1VtHcFKcpO4u4r5vOxs/LJSYnNeFDDUZmQDAt/R86B63F5swjEN1A3+xEqEybrgCAjbyzaON4SkUewHrkNANuw7hSeBP4sIt+z0+6zN7kP+IOIHMC607hxtMus1FC1b6ujZdMhgk0+dju99L/EMQZ2+8MPBb7x8GWEkneRNGUT4m7C+NPx1a9i4+HL+HaE/H4049aeoNHN54rnRzNu5RPDOI4jbT42bLPGi9pf10aC28nVp+eybkkhZ58yvImRYm1V0XfYdOAbNK/4Cg4JETIO2o5exKqi74x10U5aY9KPwxjzLeBb/ZIPAkvDrOsFbhiNcik1HO3b6mh8eDeErCobf8CEbaHzB8LfMLfEHeJTnj+xPLML0oGmDt5seZDf8fGIedYlhH8kNFL68QRDhpf31VOypZxnd9cSCBkWT0/nBx87javPyCPZMz77CyctzmEV36Vl06cJNvlwpntIXTWDpMU5Y120k9b4/CQoNQE1b3gXQsd6FwthnuIg/NMeALfFlbJ8URd0P7GaAcuTujDvlgLfCLtNalsTLSkZYdOH6sOGDkrLrImRalq8ZCbFcdtHZrB2SSFzp06M6p6kxTkaKEaRBg6lYiTk6/vvZDCECxMmbDiBc2Z2HAsa3eLs9AiW7NzCS8suItBrgD9XMMCSnVvgo5GfAvL6g/z13RpKtpTzxsEGHAIXzJ3Ct1Yv5OIFU4lz6RxvKjINHGpUjPaosWPBST1Bjl31ZgWDNLgG/oslSCj8DtLB7F5NQcUVJJkE2qWTioKnkfkbI+Y5t9MLe7fxVtEi2jwJJPs6WXZwl5UexruVzZRsKWfDO5W0egNMz0zky5fN5bqzC8hNi838Ferkp4FDjbixGDV2LOxJ+AvUL2XP0dfpCLawsmgGG0KXEZBjvZpdJsjyGTVhtzfbVzOrbg1uuxd0sklkVvka3vdhPXsYbhtniDn1lcypr+yXfux1U0cXG7ZVUlpWwXvVLXjsiZHWLilk+cyBEyMpNRgNHGrExWLU2Ing13KIsxqO4rBvKOIdbVxQ9xKb05fQ6komJdDG0qYt5Lubsfuq9pFff2VP0Ojmxkl+feRJJRxBDyHXwLmvJeDh1f1HKCkrZ5M9MdKp+al895pFfPTMfNISIk+MpNRgNHCoETfcUWMnivnvJ7Kn6HReWXYpLcnpJLe1cMFbm7jtwB/7rNf2wcDZ6gCSTfg+EZHSAaZ25RGa9gLT7d7m1c15bNpzNTs7F/C7+94iLcHNx5dO54biAhblDW1iJKUGo4FDjbjhjho7URwuOItn5l2OvN2Bx1tFV7yT/XlnkpGUSlKXlw63B29nO6fsfyfs9q3iIzVMkGiVgXcU3c7MP0TjnC280zCfVyvPYVfDPAwOTvO08JnrLuCyhVOjnhhJqcFo4IhgMjTmjpa7Vs3r08YB0Y0aO1G8VHQlyewGAAAgAElEQVQhrt0tGGO1GRR11XJu3WFcdmN4kt9HXFwc++cVh93+7WA7nQ74DX7qMOQg3I6bhGA7C8Osv6emhfskwOuvfpM2fzKZ8Ue5umgT5+W/xTSEC864aWQOVE16GjjCmCyNuaMl3KixJ2UgPtTVEzQAznZV9gSNbu5QiFRP+Kqn1zsTeCbJi9+e7aAWw4/wcmlnQk8v8Bavn43bqyjdUs72imZcUszinB2cl/8mC7P24hDrUd8IfQyVigkNHGFMlsbc0dR71NiTVdAPy6Zt4WNzniArvhGfL4mNDet4Nntln0dlZ/d7AqrbS0kd+Pt15PDj4MWkDt482EDplnKeeteaGGn+tBS+efVCsqtuICm7ccC+HC2eAWlKxYoGjjAmS2PuZNO8cSN1P/kvAtXVuHJzyfniF0hbvTpm+z9/2pvcvOhRXM4uALbGL+ap/IvpEusOoy0+kZfnn4nDhG9zaCf8k04duLnx12+S4nHxsbMKWFdcyOkFaYgIWz/RQvMnHJi4Y3c20uUg9dEAOjCsGikaOMKYLI25k0nzxo1Uf+ObGK/VMS5QVUX1N74JELPgcdOcJ3uCBkApN/cEjW5+h5uy2eHbdjzuZnz+9AHpIkF+fP3ZXHlaLglxfYNO3K4M0v7YSOs1EMwE51FI+YsQt2vgMCRKxYqOKxDGcKcAVeNP3U/+qydodDNeL3U/+a+Y5REX39Ln/REGTi4E0BQ3cGa6A3VtFCVX0H90K7d08fG5f+K6swsGBA2AA6flkrTNwdRvxJF3ZxxTvxFH0jYHB07LPfEDUWoQescRxqRpzJ1E/NXVYQcX9FdXxy4TbzokHBtcMJsjHGHgwHtZ1AP2xEg7qikpK2fr4UacsoBTUso56s2g1Z9MVnwj187eyPLcrRGzXDpjB8mODup2pBDocOJKDJJzeisp03fE7riU6kcDRwSToTF3MulISCAt41Q8i65FEjIxnUfx7XqM5sZ3Y5bHwQ/O4JR5r+J0Wg9WrOVBfmM+16e6Ki7UycojL/HVR6bxxA57YqQpSdxzxXzS69aSnTVwVNuu1sj/prnSwOMFH+HeaeuoIos8GrjLVcJH5fWYHZdS/WngUJNC9TmXk5O0EofDetpIErPwnP1JqttfiFkeFUcK6ZLlzJj5Dh5PO0v8mwm5HTzCTRzxZZBafYTEiiY2dnyExLiqnomRzppuTYx039cyyTy/GYf7WHVVyC9UbZ4O14TP8/7AZfwoeBOdWMdVyRTuCXyGRpPMp2J2ZEr1pYFDTQozM5fi8Pd9RNXh8DAzc8DcYSesw+mlvr6I+voiAEIG2lI6mWW20Xl0AV3GSZ7Dy1n1W/nZL38wYGIk3+FLOLjlJQoWl+OJ78DnTaRiWyHBwxdEzPPewDo6pe9xdeLh3sA6DRxqxAwpcIjIj4DvYc0R/lfgDOALxpg/HndDpcaJNH8Gf4z7Gxv8s6k1mUyVo6xxH+Dmrktjlke1M4GiYIi2UAIHgtkcCGbT4YsjPuTl9NadLGzdTaa/CZ87GHY2vYL0vbzdtZTqLb2CWQjOSt8bMc8OwncmjJSuVCwM9Y7jMmPMV0TkWqACayrXFwANHGpC+L3naR70LcFrV+nUmGx+35WC3/M03+D8Ye/f6w+S/EEiz6XPpdLpRjAUSDPntWxh9pF3cWL1swhJiAPZqWH3sT8rY+C8Tw47PYKUQBut7oGz9KUE2k74WJQazFADR3fPpCuBh4wxR8fTZPVKDebJrvk9QaObFw9Pds2PMCnr4Iwx7LQnRnp8exWtWYtIDworOp0s6nKRYhIJ+JoIyvuYUBvx/gBFtUd5P68u7P5aSA775FcLyRHLcGb7dt5IXUbAcazzoCvk58z27cCNJ3hkSh3fUAPHRhHZg1VV9XkRmQKEn2JMqTC2/uw/eevlZ+l0CAkhw7LzL+Hsf/rSqOVfazKjSj+exvYuNrxTScmWcvbUtOJxObjytFzcr9QxPeBEen39uzwLcMXN56KX/qEn7eOvht9vu/GQLF1h0yNpK57DBW++wua0XnN+NG/h6PJFUR+XUkM1pMBhjLlbRH4ItBhjgiLSTsTnPJTqa+vP/pNXXnmWoNPqb9rpFF555VmAUQseOcGj1DoHdsjLCR4d0vahkOHVA9bESM/sqqUrGOL0gjS+t+ZUPnpmHqnxbn7w5gu8m+fmhdMTaE50kNYRYuWOTor39e0rktUSPo99/nROjzvSZ2DEgHGwL0xv8m5L5+7gueCpXLt5I6ltzbQkp/H60pVcPFf7caiRE81TVQuAGSLSe5sHYlwedRJ66+VjQaNb0OHgrZefHbXAsXbP8/zPgo8RcBz7+LpCAdbueR64NeJ2FY0dPFxWwSNbK6hs6iQ90c3Hl01nbXEhC/P6tlVsPaORN/Nn4HdZx9qc5OTJJYk4Ahn4PD/H4zvKrIOP4/DuDpvXpxwb+R//9ZzuqiVJumg3cewITOVzjkeAz4fdZmngCfxzGimdcytHyCabI6zlQZYGXgNi1yteqd6G+lTVH4BZwDtA97CxBg0cagg6I8xpHSl9JLjSXfhPzcAc7EC8QUy8E39RCq7qgQMLev1BnnmvltKycl49cASA82Znc/cV87n0OBMjHZiaTLCuk7j9rT15BOak8NxZGSyqasYXn8WeeR/HmfMs4R6wvci1AwnCvV19O/OtdEa+e4j35HKu71XO5dV+6XlDPzlKRWmodxzFwEJjjI7yr6KWEDL44ufjSliBOFIwoVYCna/giXDlPRJ+c+Mn8SWmQmFqv/RP8M/26/eqWigtK2fDO5U0dfjJT0/gny6aww3FBRRkDBxfqr+6hkTc7zUjIevfRLxB3Luaae01C1PI6cHdujzs9s0kc6VzM2tcx3p9dxkXNWQTadLXollfZs+erxEKHRuU0+FIoGjWlwctr1InaqiB411gGhDDgX3UZJE4Yy2h1lxErKt7cabiTrqMxGmnjVoZjiSE/+o94krhD28epnRLOTsrm4lzOrhs0VTWLSnkI7OycUZxV+Te32L1+utFQoa4/c3Qq1+Fryt8g3xCqIPfT7uGXxWto9ozhVxfPXccLOHammci5pk7zWpqPPj+j/H6qon35FI068s96UqNhOMGDhHZiFUllQK8JyKbgZ4JkI0xHx3Z4qmTgb9jOv2f3hZx4++YPmplyKL+2ICDxuBo7MJZ0Y6ztpNvhGqYPy2Fb61eyJoz88lIijv+ziLxhRjYEaM7/ZhkR33YzX9/2lf5ZWYxnU4ryFTFT+Xf532W5ilncddxss2ddo0GCjWqBrvj+PGolEKd1LpChnBfqF2h0av5vNb/KPcHbyFYFcJZ2Y6jM4hxCXOnHeLHH7uZ0/KtiZGGI9XfSot7YOe+lEA7YFV1ufByZsJDhOtjUVKwkk6fv09apzOekoKVxw0cSo224wYOY8xLACLyQ2PMV3svsx/PfWkEy6ZOEq1iKHcFeSU+QIvDkBoSVnhdFAZGfjqYrkCI5/fU8uobxTi9jTgQghlxeGbBDVMeYdqBLk4vuDMmed2w/xnun//RAZ3xLmhuBGcWyY4jLE16kMyW8BOCVfYLGoOlKzVWhtrGcSnw1X5pV4RJU2qAvyR5qXcKAfuCvsVp2JToZ0rQcPcI5XmgrpWSLeWsf7uShvYuEiSbZSm7WTX/KQrTP8TnS+LQgTOpaJgRszxXV+5AgNI5l9LiSiE10Mra/c+wtuI15qyuwd/hpH5HGq9mX8V1YbbPFSdVJhg2XanxZLA2js9hPUBeJCK9nwlMAXTAfzUkda4QQfp++QXESo+lNl+AJ7ZXUVJWzrYPm3A5hIsX5LBuSSGv/+KrBNyFHN55AYe7NwgFSWk8fLxdRuXgrMtZs2sD1+17oyct6HSya9H1BEvW96Rl1YWfmOnO/T6+O9OJ13msyiw+aLjzAx9cGLNiKjVsg91x/Al4Gvh36HNx2GqMGVqXWzXpBSPMUBwpPRrGGLYebqRkSzlP7qymoyvIrClJ/OuV87l2cQFTUqzhOrbV1dPlDdE1JR/jjkP8XcTVVxLX0jDsMvS4/H32BG6k6P2niPc14vVkcHDWlTgu3wW9LrtyOgdO1gSw6v1OQu0ufjHXQ228MNVruHOfj1U1gdiVUakYGKyNoxloBm4CEJEcrOcKk0Uk2Rjz4cgXUU10Hk8nPt/AfhAeT2eYtYemvtXH+rcrKCkr52B9O4lxzgETI/UmrkTiWo4S13J0QHqspKWW0bwatuz8AoGOLFyJDUw57TEykzcDx57U8qeED5jOdA9X1Pi4ol+gcKZHHqtKqbEw1J7jq4H/BPKAOuAUYDegI6mpQbXPSsG5J9TTMQ7AOIT2WQOHAz+eQDDEi3vrKS0r5/k9dQRChrNPyeBH183iytNzw85x0c3hvoBg4Fmg95eyC4c78iRJ0Xp1+7mc95FXSDtl87FEH6Q8eKyaLuhws7/gDM4Is33qqhk0rd+P8R+rwhO3g9RVM2JWRqViYaiN498DlgPPGmMWi8hK7LsQpQYTKJhKyNmBq99QHKHcoV3tf3CkndKych7dWkFdq4/s5Dj+7ryZrC0uYHbO0IKPK34BIkLA+yqEWsGRgiv+PJye+cM5tD7+6r4Z8zqsOOM1yAghjQ6CL81HdjZgsKuuij5KdXb4PJMWW/1MWjYdItjkw5nuIXXVjJ50pcaLoQYOvzGmQUQcIuIwxrxgP46r1KDiQw1487Lpykvqmx48EnGbzq4gT+2spqSsnM0fHMUhsHJeDjcUF3LxghzczujaRwIOL27PAlyeBX3S/Y4Try7r7++Nh3/338Dvy24A4MtN8dYQ6+f0Xc9B5P4rSYtzNFCocW+ogaNJRJKBV4AHRaSOvvf8SkV0c6CEPzg+TZccG3Yjzni5OVACXNKTZoxhR0UzJWXlbHynilZfgBlZidy1ah7Xn13A1NQTnw51ifePbHfdhunVx0JCfpZ0/RG46oT329slnUKIeP4XH3UY2sWQbAZ2KvQk6iRoamIbauC4BmsSpy8ANwNpwHdGqlDq5LLS/Tweuig1N/cZ+vsjbmtE16PtXTy2rZKHy6yJkeLdDq48NZe1SwpZNjNz2D26Adrr69k/p4qczgJSjINWCVGXXMXs/eGH/zgRznQPlzXBZXZDeHl8kG2dAUyvp8ccTjh/3cJIu1BqQhjqRE7tInIKMMcYc7+IJAIn3CtJRNKB3wCnYo2F9XfAXqAEmAEcAtYaYxrF+tb4Kda0tR3AbcaYt080bzX6Aq0uzk3tO/R3yAjbyxfx2INv88x74SdGiqXSBat4wZFFMO7YDHtOsmhdsIrLYpRH/8btQo8TcQl7jIP2Nj/JmR7OuWYWc5dNi1GOSo2NoT5V9Rngs0Am1rwc+cCvgItPMN+fAn81xlwvInFYA/n8K/CcMeYHInI3Vr+Rr2L1UJ9j/ywD/sf+rSaIys1TmH5BDQ634UhnJq9WLuO1yuUc9WWQnniEm5dbEyMtyB04zlOsvOGYNaATYhAnbzhmxSyPcI3bp62awXJts1AnmaFWVd0JLAXeAjDG7Lf7dERNRFKB84Hb7H11AV0icg3H+sfeD7yIFTiuAR6w5wJ5U0TSRSTXGKNDvE8QDR9MYYd7MXsTZrG3yfqinuEsZ2ldGT/99X/gcY38kBrthB/xNlL6idLGbTUZDDVw+IwxXd11zfb0sSc6tGkRUA/8TkTOALYC/wxM7Q4GxpjqXoEpHyjvtX2FnaaBY5zbVdXMw2UV/HH6bQR8LpJ8Ps5wVTPbeYRk08nhtNxRCRoAyfhoY2DjevKxWQKUUkM01MDxkoj8K5AgIpdijV+1cRh5ngX8ozHmLRH5KRx3rLtwLaMDgpaIfBarOo3p00dvngfVV3Onn8ffqaSkrJx3K1uIczrIdbaywFlHrqPl2Lwc4iQrZeCAfiPls64n+WlgTZ/qKidBPut6EsIOOaiUimSogeNu4NPATuDvgaewGrdPRAVQYYx5y37/iL3/2u4qKBHJxeqh3r1+Ya/tC4Cq/js1xvwa+DVAcXGxTnE7ikIhw5sfNFC6pZyn363BFwixIDeVf1u9kDWL8/nJD78f9smoJBm9q/1PyF/BBf8XuJJW4knBy2dcT1npSqmoDPWpqpCIbAA2GGOG9fyiMaZGRMpFZJ4xZi9WA/t79s+twA/s33+xN3kc+AcR+TNWo3iztm+MDzXNXh7ZWk5pWQUfHu0gJd7FDcUFrCuezqn5qT3BIpUArQx8Sip1FLsCbWtfwB0pG/gn16M9aV3GySuti0/4CQ+lJqvBhlUX4FvAP2BVGYmIBIH/NsYMpx/HP2J1JIwDDgKfAhxAqYh8GvgQuMFe9ymsR3EPYD2O+6lh5DshbNhWyb2b9lLV1EleegJ3rZrHmsX5Y10swJoY6bndtZSUlfPyvnpCBs4pyuKLl87hilNziXcPbLNY4p/JS+4qgnJsDCancbDEP3PUyp137dd57k9f4yMZB0h1+2jxe3i9cTYzPv71USuDUieLwe44vgCcCywxxnwAICJFwP+IyBeNMT85kUyNMe8AxWEWDbj4s5+mis0UbRPAhm2V3LN+J51+q/6/sqmTe9bvBBjT4LG/1poY6bFt1sRI01Lj+fyFs7mhuIBTspKOu21RaC4hfyplroO0iZdkE09xoIii0NRRKj0sWLES+D4lf36A1oYjpGRls+Ljt9jpSqloDBY4bgEuNcb0DCpkjDkoIp8A/gacUOBQkd27aW9P0OjW6Q9y76a9ox442nwBNm6vorTXxEiXLJjKuiWFnD93Ck7H0Hp0h0wbs0O5zO7K7ZfeOhLFjmjBipUaKJSKgcECh7t30OhmjKkXkdh27VUAVDZ2QJiG5MrGjlHJ3xhDWffESDuq6fQHmZ2TzNeuXMC1Z+WTnRz93BD1u0rIWXQrvT8yxvip31XCdK6MZfGVUqNgsMDRdYLL1AnK8jfREJcRNn0k1bV6Wf92JaX2xEhJcU6uOTOPtUsKWVyYPqzxouIPbMbbCZ5F1yIJmZjOo/h2PUZ85ebBN1ZKjTuDBY4zRKQlTLpAmN5UatjObXmT+kUZrJn7NFnxjTR4M9iw7wqm7GoEPhHTvLonRiqxJ0YKhgzFp2Rwx/WzuOq0XJKOMzFSNI6mQHblZgL9AsXR6OZxUkqNE4NNHTs63XpVjzNy3mH6ohocLqsrSnZCI3+36CE+rI/dwHgH69soLavg0bcrqLcnRrr9vJncUFzI7JzkmOXTbd+5QsqzBk+vp299Lit9RcxzU0qNtNhcUqqYKVxe2xM0ujlchsLltcPab0dXgKd21lC6pZzNh47idAgr501hbXEhK+dHPzFSNGas9rLVkcDcV0Jktlp3GvtWOJhxVewmUVJKjR4NHOOMIykUVfrxGGN4p7yJ0rIKNm6vos0XYGZ2El+5fB7XnTW8iZGi0RlKZsbV7XRdDTV22gygMxj7uxul1MjTwBHJjlJ47jvQXAFpBXDxN+H0tSOebdArbG44mw0Hr6bBm0FWfCNrip5gadbWIe/jaHsX69+uoLSsnH21bdbESKflsq64kKUxmhgpGr91fpp1+x6gdksW/jY37mQ/U5c0UDL3Fq4e1ZIopWJBA0c4O0qpfv2LHJznxuvJJN7XRtHrXyQXRjx4vF67hIcOrKUrZA333eDN5A97biQw28Gq42wXDBle2V9PaVk5z7xXiz9oOKMwne9feyqrz4j9xEjRSH0tiYpd+WCsuyZ/WxwVL+aTWp8EF41ZsZRSJ0gDRxjVW7/JnlkeQk7rytwb72TPLAds/Sa5Ixw4Nn54ZU/Q6NYVimPjh1fy7TDrlx/t4OGych7ZWkFVs5eMRDefXD6DtUsKmD9t5CZGisapu5/vCRo9TMhK5/YxKZNS6sRp4Ajj4FQvIWffB8pCTuHgVC+5EbaJlaPegX04+qd7/UE27aqhZEs5r7/fgAismDOFr121kEsW5ozaHBdDJcFwT3RHTldKjW8aOMLwesJ/8UZKj6VUaaPFDOzgkCptvFvZTGlZORu2VdLiDVCQkcAXL5nL9cUF5KcnjHjZTpTDlUooMDBIOFzj445IKRUdDRxhuHwZBOIbw6aPtGlpAZpboNdAshgBnyOBq//7VeJcDi5fNI11Swo5pygLxxDHixpLZ1y6lm1P/x76DKPu4oxLR/5hA6VU7GngCCN73w1syNzO+oNX9DzZ9LGip1lz9AxGemilg2cU4W8I4drbivit6CEGfAkevn3RfK45M4/0xNjOkz3SLrptDQDbnyklFGjB4UrljEvX9qQrpSYWDRxhPFOzmN/WLCSI1SmuwZvJb9+7iSQ8LBrBfKuaOumsCOCq6sDhD2FcQjA3kWBBIibFza0fmTGCuY+si25bo4FCqZOEBo4wfmG8BKVve0YQB78wXr4Q47y6AiGe3V1LyZZyXtlfj9tAMNND1+xUQlMTwH6yK7NlYNWZUkqNBQ0cYfgJP/xGpPQTsa/XxEhH27vITYvnzpWz+du2/2X3mTcSch+rjorzdzF7fwmgc0kopcaeBo5wIrU3D7MdutXrZ+P2akrKytle3oTbKVy6cCo3FBdy/hxrYqTfVT/DFbuyeWvuZdQmOJnaGWTZvmd5bcprw8tcKaViRANHGMbtYHnWZj4254meoc3X77+aNxuWRr8vY9hyyJoY6amd1sRIc6cm8/WrFnDt4nyy+k2MdGHrhXyhbinxtccGAPTKUtzu6mEfl1JKxYIGjjCKz9zFbal/Is5pTeGandDIbYv+RKAliaE+VlXX6uXRrZU8XFbOwSPtJHtcrFmcx9riQs48zsRId9SvIt70DSbxxsMd9ccbcEQppUaPBo4wbkp5sCdodItzBrkp5UHgrojbBYIhXthbT8mWcl7Ya02MtGRGBp+7cBZXnZ5LYtzgpzvNF75TXKR0pZQabRo4wkh3tkaV3n9ipCkpHj6zooi1xQUUTYlu6PBQsAGnMztsulJKjQcaOMJwejN4rbGIxw6s7ukAeO3sjZybcXDAui/vq+eW3262J0bKYd2SQi6cN+WEJ0ZytpYQSr0dh+NYdVUo5MPZWgJce6KHpJRSMaOBI4xd22/ngbapfYY2f+C9m0hPrmXlVX3XXTozk3+9cj5rzswnJwYTIzWk7iC39Vd4425G4rMw3gbiux6kOvU9Coa9d6WUGj4NHGE80FxIV7+2665QHA80F/IP/daNdzv57PmzYpZ34xIfzsAbzDr0HPG+EN50B+/PSKTRNTqz9Sml1GA0cIRRJ4ZwnTas9JEVzIRaiad2ar9AMfJZK6XUkMSuK/TJxBPhtERKjyFne1JU6UopNdo0cITRNScN02+4cuMQuuakjXjeUw7ejPj73giK38WUgzePeN5KKTUUWlUVhmOqB7+k4drfiniDmHgngTkpOHI8g288TJnelbAbjsx5lEB8Ay5vFtn7r7PSlVJqHNDAEcayD97lraJT6co7Vj3kCgZYdvBdoHhE8w4tbyDlb2eTVvORY2kOH6HLtB+HUmp80MARxmlVHxDv9/NW0SLaPAkk+zpZdnAXc+orRjzvQ+7/h3vhdKYcuB6XN4tAfAP1sx/B7/6QfK4Y8fyVUmowGjjCSDXtzKmvZE595YD0keb1VePNq6I1782+C3zjf4pYpdTkoI3jYVzCy7jx90lz4+cSXh7xvOM9uVGlK6XUaNPAEcY0+YDVPEMaLYAhjRZW8wzT5IMRz7to1pdxOBL6pDkcCRTN+vKI562UUkOhVVVh/JCP8l3zOKfL3p60DhPHN+Sj/McI55077RoADr7/Y7y+auI9uRTN+nJPulJKjTUNHGE8470WvyOHr7hKyZMGqkwWPwqs5cXQuaOSf+60azRQKKXGLQ0cYbQiPB46j8e7zuuTrs3TSimlbRxhZUUIEZHSlVJqMhmzwCEiThHZJiJP2O9nishbIrJfREpEJM5O99jvD9jLZ4x02dbkpNK/j7jHTldKqcluLO84/hnY3ev9D4GfGGPmAI3Ap+30TwONxpjZwE/s9UbU1790Hp/MSWMKggBTED6Zk8bXv3TeoNsqpdTJTowZ/fG6RaQAuB/4PvAlYDVQD0wzxgRE5Bzg34wxq0Rkk/36DRFxATXAFHOcghcXF5uysrKRPxCllDqJiMhWY8yg4yqN1R3HfwFfAUL2+yygyRgTsN9XAPn263ygHMBe3myv34eIfFZEykSkrL6+fiTLrpRSk9qoBw4RuRqoM8Zs7Z0cZlUzhGXHEoz5tTGm2BhTPGXKlBiUVCmlVDhj8TjuucBHReRKIB5IxboDSRcRl31XUQBU2etXAIVAhV1VlQYcHf1iK6WUgjG44zDG3GOMKTDGzABuBJ43xtwMvABcb692K/AX+/Xj9nvs5c8fr31DKaXUyBpP/Ti+CnxJRA5gtWHcZ6ffB2TZ6V8C7h6j8imllGKMe44bY14EXrRfHwSWhlnHC9wwqgVTSikVkQ45otQ4smFbJfdu2ktVUyd56QnctWoeaxbnD76hUqNIA4dS48SGbZXcs34nnf4gAJVNndyzfieABg81roynNg6lJrV7N+3tCRrdOv1B7t20N8IWSo0NDRxKjRNVTZ1RpSs1VjRwKDVO5KUnRJWu1FjRwKHUOHHXqnkkuJ190hLcTu5aNW+MSqRUeNo4rtQ40d0Ark9VqfFOA4dS48iaxfkaKNS4p1VVSimloqKBQymlVFQ0cCillIqKBg6llFJR0cChlFIqKho4lFJKRUUDh1JKqaho4FBKKRUVDRxKKaWiooFDKaVUVDRwKKWUiooGDqWUUlHRwKGUUioqGjiUUkpFRQOHUkqpqGjgUEopFRUNHEoppaKigUMppVRUNHAopZSKigYOpZRSUdHAoZRSKioaOJRSSkVFA4dSSqmoaOBQSikVFQ0cSimlouIa6wIopdR4t2FbJfdu2ktVUyd56QnctWoeaxbnj3WxxowGDqWUOo4N2yq5Z/1OOv1BACqbOrln/U6ASRs8tLdY2EIAAAr6SURBVKpKKaWO495Ne3uCRrdOf5B7N+0doxKNPQ0cSil1HFVNnf+/vXMPtrqq4vjnOxdEwOKCoINXB0SRRFReEoSaUwZpDZKDCZMFmqmpM6mlgjqUmhHSlNrDx5hppjw0H0QYOgg6PriGylNFUUlR1KuW5CtRV3/sdbjnXs65hyNwfvexPjO/Ofu39v7t397rnvtbZz9+a5UlbwuE4QiCIGiCPao7liVvC1TccEjaS9IiSU9LWi3pRy7vJuk+Sc/5Z1eXS9JVktZKWiFpcKXbHARB2+Xc0f3o2L6qgaxj+yrOHd0voxZlTxYjjo+BH5vZ/sBw4AxJ/YHJwEIz6wss9HOAo4C+fpwCXF35JgdB0FYZO6iGacceSE11RwTUVHdk2rEHttmFcchgV5WZbQA2ePq/kp4GaoBjgCO82E3AYuB8l//ZzAxYIqlaUk+vJwiCYIczdlBNmzYUjcl0jUNSb2AQUAvsnjMG/rmbF6sBXs67bL3LGtd1iqSlkpbW1dXtyGYHQRC0aTIzHJJ2Af4KnGVmG5sqWkBmWwjMrjOzoWY2tEePHturmUEQBEEjMjEcktqTjMYtZnaHi1+X1NPzewJvuHw9sFfe5XsCr1aqrUEQBEFDsthVJeCPwNNm9uu8rLnARE9PBO7Ok3/Pd1cNB96J9Y0gCILsyMLlyEjgu8BKSctcdgHwS2COpO8DLwHHed584GhgLfA+cGJlmxsEQRDko7RZqXUhqQ74V4Vv2x14s8L3bEmEfkoTOipN6Kg026KjXmZWcpG4VRqOLJC01MyGZt2O5kropzSho9KEjkpTCR2Fy5EgCIKgLMJwBEEQBGURhmP7cV3WDWjmhH5KEzoqTeioNDtcR7HGEQRBEJRFjDiCIAiCsgjDEQRBEJRFGI4mkFQl6UlJ8/x8b0m1HjNktqSdXN7Bz9d6fu+8Oqa4fI2k0dn0ZMchaZ2klZKWSVrqsrJjq0ia6OWfkzSx2P1aGu7N+XZJz3gMmhGhn3ok9fPvTu7YKOms0FFDJJ3t8YtWSZopaedMn0dmFkeRAzgHuBWY5+dzgPGevgb4oadPB67x9Hhgtqf7A8uBDsDewPNAVdb92s46Wgd0byS7HJjs6cnAdE8fDdxDclw5HKh1eTfgBf/s6umuWfdtO+nnJuBkT+8EVId+iuqqCngN6BU6aqCXGuBFoKOfzwEmZfk8ylwpzfUgOVNcCHwFmOdf1DeBdp4/Aljg6QXACE+383ICpgBT8urcXK61HEUMxxqgp6d7Ams8fS0woXE5YAJwbZ68QbmWegCf9394hX62Sl+jgIdDR1voJRdaops/X+YBo7N8HsVUVXGuAM4DPvXzXYH/mNnHfp4fF2RzzBDPf8fLb1UskRaOAfdKelzSKS4rN7ZKa9VTH6AO+JNPeV4vqTOhn2KMB2Z6OnTkmNkrwK9IPvw2kJ4vj5Ph8ygMRwEkfRN4w8wezxcXKGol8rYqlkgLZ6SZDSaF+D1D0uFNlG1remoHDAauNrNBwHvUh0QuRFvTz2Z8fn4McFupogVkrVpHvr5zDGl6aQ+gM+n/rTEVex6F4SjMSGCMpHXALNJ01RVAtaScR+H8uCCbY4Z4fhfgbdpALBEze9U/3wDuBIZRfmyV1qqn9cB6M6v189tJhiT0syVHAU+Y2et+Hjqq50jgRTOrM7NNwB3Al8jweRSGowBmNsXM9jSz3qTh8/1m9h1gETDOizWOGZLbxTHOy5vLx/suh72BvsBjFerGDkdSZ0mfy6VJc9SrKD+2ygJglKSu/utqlMtaNGb2GvCypH4u+irwFKGfQkygfpoKQkf5vAQMl9RJkqj/HmX3PMp64ae5H8AR1O+q6uOKXksaUndw+c5+vtbz++RdfyFp98Ia4Kis+7OdddOHtEtjObAauNDlu5I2Fjznn91cLuD3ro+VwNC8uk5y/a0FTsy6b9tRRwOBpcAK4C7Sjp/QT0MddQLeArrkyUJHDXV0MfAM6YfZzaSdUZk9j8LlSBAEQVAWMVUVBEEQlEUYjiAIgqAswnAEQRAEZRGGIwiCICiLMBxBEARBWYThCCqGJJN0c955O0l1qvc+PEbSZE//TNJPPH2jpHGevl5S/wzaPsO9k85oJJ/kfVgm6SlJP6h023YU/jcwSfvmyc522VA/ny+p2tPvZtXWoLK0K10kCLYb7wEDJHU0sw+ArwGv5DLNbC7pJaWimNnJO7aJRTkV6GFm/yuQN9vMzpS0G7Ba0lyrfwO6xSCpysw+aSReSXoJ9ud+Po708hkAZnZ0hZoXNCNixBFUmnuAb3i6wdvC/uv9d01dLGlx3q/dCUqxQFZJmp5X5l1Jl0laLmmJpN1dfpyXXS7pwQJ1y0cWq7ze410+l+QfqDYnK4QltyvPA70kDZP0iDs3fCT39rikAyQ95iOUFZL6+hv4f/d2rcq77xBJDyg5kFyQ54JjsaTpXs+zkg5zeSdJc7ze2UqxGHK6GiXpUUlPSLpN0i4uXydpqqSHgOMKdOsukp8kJPUhOcyry9PZOkndC+jyXEn/9LZc7LKC/QxaHmE4gkozi+T2YGfgIKC2RPmCSNoDmE7yIzYQOETSWM/uDCwxs4OBB4Hc9NFUYLTLxxSo9liv62CSf6AZknqa2RjgAzMbaGazm2hTH9LbvGtJb/kebsm54VTgF17sNOBKMxsIDCX5D/o68KqZHWxmA4B/SGoP/BYYZ2ZDgBuAy/Ju187MhgFnAT912enAv83sIOBSYIi3qztwEXCkJYeUS0mxZnJ8aGaHmtmsAt3aSHKbMoBk6Iv2P08Po0juLIaR9DlEyfnlFv0sVVfQPAnDEVQUM1sB9CY9hOZvQ1WHAIstOX77GLgFyHnm/YgUswCS++nenn4YuNHXIaoK1HkoMNPMPvGppgf8PqU4XtIy0ujpVDN7m+RY7jZJq4DfAAd42UeBCySdD/TyKbuVwJE+ijjMzN4B+gEDgPu87otITuly3FGgf4eSDDNmtork5gRSwKP+wMNe10RSsKQcpYzBLNJ01ViSI8tSjPLjSeAJ4AskQ1Kon0ELJNY4giyYS4ovcATJJ9FnoZCL6BybrN6Xzif499zMTpP0RdJU2TJJA83sra2ssylmm9mZjWSXAovM7FtKoTsXextulVTrbVgg6WQzu1/SEFJ0u2mS7iU9oFeb2Ygi98yttWzuXxPtF3CfmU0okv9ek72DvwEzgKVmtlEqqSYB08zs2i0yGvXTzC4pVVnQ/IgRR5AFNwCXmNnKbaijFviypO6SqkgjmAeaukDSPmZWa2ZTSVHR9mpU5EHS6KFKUg/SCOazejPuQv3C/6S8NvQBXjCzq0gG9CCfdnvfzP5CMqiDSU7oekga4de1l3QATfMQ8G0v3x840OVLgJHy3VG+FrLf1nbER0Xn03CqrCkWACflraPUSNqtSD+DFkiMOIKKY2brgSu3sY4NkqaQXEsLmG9md5e4bIakvl5+Icmrbz53kkJwLicFuDnPkmv0z8LlwE2SzgHuz5MfD5wgaRMpvvYlpOmwGZI+BTaRYkd/pLQF+SpJXUj/q1eQvBAX4w9+zxWkaaIVJLfjdZImATMldfCyFwHPbm1niqx/FCt7r6T9gUd9dPIucAKwb+N+bm2dQfMivOMGQSvBR17tzexDSfuQjON+ZvZRxk0LWhkx4giC1kMnYJHvyBI+csm4TUErJEYcQRAEQVnE4ngQBEFQFmE4giAIgrIIwxEEQRCURRiOIAiCoCzCcARBEARl8X9EDQOve6MPygAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig,ax=plt.subplots(1)\n",
"ax.scatter(airline_df['Miles'],airline_df['Deaths'] )\n",
"x=np.linspace(4000,8000,10)\n",
"ax.plot(np.linspace(4000,8000,10),(.12)*np.linspace(4000,8000,10))\n",
"for x in range(200):\n",
" ax.scatter(airline_df['Miles'],predicted[x])\n",
"ax.set_title('Deaths vs. Miles, with model predictions')\n",
"ax.set_xlabel(\"Millions of Passenger Miles\")\n",
"ax.set_ylabel(\"Deaths\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The plot above shows that there is far more variation in the deaths per passenger mile than can be accounted for by a poisson model. For example, just look at rows 5 and 7, both had comparable number of miles but the number of deaths is way outside the 95% interval. \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Now let's look at accidents vs miles."
]
},
{
"cell_type": "code",
"execution_count": 101,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_9e55af8f55ddc462276e2095df7f3378 NOW.\n"
]
}
],
"source": [
"stan_code='''\n",
"data {\n",
" int accidents[10];\n",
" vector[10] miles;\n",
"}\n",
"parameters {\n",
" real<lower=0> theta ; \n",
"}\n",
"model {\n",
" \n",
" theta~gamma(24,5000);\n",
" accidents~poisson(miles*theta);\n",
"}\n",
"generated quantities {\n",
" int predicted[10] ; \n",
" int answer ; \n",
" \n",
" for(i in 1:10)\n",
" predicted[i]=poisson_rng(miles[i]*theta);\n",
" answer=poisson_rng(8000*theta) ; \n",
"}\n",
"\n",
"'''\n",
"sm_weights=pystan.StanModel(model_code=stan_code)"
]
},
{
"cell_type": "code",
"execution_count": 102,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jteitelbaum/anaconda3/lib/python3.6/site-packages/pystan/misc.py:399: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
" elif np.issubdtype(np.asarray(v).dtype, float):\n"
]
}
],
"source": [
"accidents=sm_weights.sampling(data=dict({'accidents':airline_df['Fatal'],'miles':airline_df['Miles']}))"
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Inference for Stan model: anon_model_9e55af8f55ddc462276e2095df7f3378.\n",
"4 chains, each with iter=2000; warmup=1000; thin=1; \n",
"post-warmup draws per chain=1000, total post-warmup draws=4000.\n",
"\n",
" mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\n",
"theta 4.2e-3 6.9e-6 2.6e-4 3.7e-3 4.0e-3 4.2e-3 4.4e-3 4.8e-3 1440 nan\n",
"predicted[0] 16.24 0.07 4.14 9.0 13.0 16.0 19.0 25.0 3615 1.0\n",
"predicted[1] 18.09 0.07 4.38 10.0 15.0 18.0 21.0 27.0 3688 1.0\n",
"predicted[2] 21.16 0.08 4.8 12.0 18.0 21.0 24.0 31.0 3765 1.0\n",
"predicted[3] 23.08 0.08 5.07 14.0 20.0 23.0 26.0 34.0 3623 1.0\n",
"predicted[4] 24.49 0.09 5.17 15.0 21.0 24.0 28.0 35.0 3676 1.0\n",
"predicted[5] 25.42 0.09 5.37 15.0 22.0 25.0 29.0 36.0 3223 1.0\n",
"predicted[6] 24.76 0.08 5.2 16.0 21.0 24.0 28.0 36.0 3796 1.0\n",
"predicted[7] 26.24 0.09 5.39 16.0 23.0 26.0 30.0 38.0 3455 1.0\n",
"predicted[8] 31.27 0.1 5.96 20.0 27.0 31.0 35.0 44.0 3601 1.0\n",
"predicted[9] 29.86 0.1 5.87 19.0 26.0 30.0 34.0 42.0 3581 1.0\n",
"answer 33.7 0.1 6.04 23.0 30.0 33.0 38.0 46.0 3748 1.0\n",
"lp__ 354.66 0.02 0.72 352.61 354.51 354.94 355.11 355.16 1713 1.0\n",
"\n",
"Samples were drawn using NUTS at Mon Apr 16 12:12:57 2018.\n",
"For each parameter, n_eff is a crude measure of effective sample size,\n",
"and Rhat is the potential scale reduction factor on split chains (at \n",
"convergence, Rhat=1).\n"
]
}
],
"source": [
"print(accidents)"
]
},
{
"cell_type": "code",
"execution_count": 112,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Wl4VNeV6P3/0gBCCMQMxSAJJMAMwuABDDY2BgOx5EHujh07dOIhviSOMzrpxA55E7s7dIbbg/vedCdxpztJd0gcp28jjMFmsvEcbOzYyAwGgYRASCAQEgjNqv1+2LtQITSUpCqVqmr9nkePqk5N+5yqWuvU3uesLcYYlFJKxa64cDdAKaVUeGkiUEqpGKeJQCmlYpwmAqWUinGaCJRSKsZpIlBKqRiniaAfEpGfi8j/18ntRkSy+rJN/YmILBaRj/2uF4vILeFsk2vHJe1q5/YM994l9GW7OmjLAyLyRoD3/bWI/CDUbergtS/ZZiLyoojc34PnSRORGhGJD34rI58mgiARkZ0iclZEBvb2uYwxXzDG/G0w2tWZcH7B2yMiT7ov/VfaLP+aW/4kgDHmdWPM9LA0shNt29VfElQ0Mcbcaoz5TVf3a7vtjTElxpgUY0xLaFsYmTQRBIGIZACLAQPcEdbGRL6DQNs9vs+65SqCiaUxpx/SNyU4Pgv8Cfg1bYKYiAwSkX8QkaMiUi0ib4jIIHfbDSLylohUicgxEXnALb9kT11E/lpEykTkhIg81Ob5B4rI34tIiYicdN1KvudfIiLHReQbInLKPceD7rbVwCrgW+4n80a3/NsiUioi50XkYxFZ1nZlReQ6ESn3/5ktIneJyB53eb6I7BaRc65N/9iNbfkukCwis9xzzQIGueW+11oiIsfbe7CIxInI4yJyWETOiMhzIjLC3ZYkIr91y6tE5F0RGdtVg0TkNyLyDXd5gvt18kV3PUtEKl2Qu9guEfkvIA3Y6Lbvt/yecpV7v06LyJpOXvfXIvKvrjukRkTeFJFxIvK0+/V5QETm+d1/hvtlWiUie0XkDr/bRorI8+49eQfIbPNaV4jINrcuH4vIPV1tF/e4B1y7/q/7fB/w/8y49qwVkTeBWmCKiKSKyL+7z2OpiPzA91kSkXj3eT4tIkeA3Davt1NEHva7/r9EZL/7vO4Tkava2/ZyeRfTeLc9KkWkUET+l99zPuk+N//pnneviFzjd3uX35GIY4zRv17+AYXAF4GrgSZgrN9t/wLsBCYA8cAiYCD2g3oeuA9IBEYCc91jfg38wF3+BHASmA0MBn6H/eWR5W5/GngeGAEMATYCP3S3LQGagb9xr5GD/TIOb/s67vp04Bgw3l3PADI7WOfDwHK/638EHneX3wY+4y6nANcFuB2fBH4LfAf4sVv2E+AJt/xJv/U67ve4YuAWd/lr2KQ80W3nXwC/d7d93m2fZPdeXA0MDaBdDwEb3eVPu3X/g99tG7pql9/2NMC/YZPblUADMKOD1/01cNq1Mwl4GSjC7njEAz8AXnH3TcR+Dr8DDACWYj9f093tzwLPuc/QbKAUeMPdNti97w8CCcBV7nVntfc5adPGB7Cfsa+7NnwKqAZGuNt3AiXALPfciUC+e18GA2OAd4DPu/t/ATgATMJ+pl9x2yzB7/kedpfvdutxLSBAFpDexbb3Pc+rwL+67ToXqACW+X0O67Hfl3jgh8CfuvsdiaS/sDcg0v+AG7DBf5S7fgD4urscB9QBV7bzuCeA9R0858UvHvAfwI/8bpvmPtBZ7sN/wf+DCCwEitzlJe71E/xuP4ULzG2/4O45TwG3AIldrPcPgP9wl4e4dqS7668BT/m2STe25ZPYgJ/mgkei+z+JwBPBft8X2l33uPcnARu03wLmdLNdmUCVez9/jk0ox91tvwEe66pd7nqGe+8m+i17B7i3k8/Bv/ld/zKw3+96NlDlLi8GyoE4v9t/77ZpvNsGV/jd9ne0JoJPAa+3ee1fAN9v73PS5n4PACcAabNOvh2BncDf+N02Fpv8Bvktu4/WhPYy8AW/21bQcSLYAny1g3Z1tO0T3OepBRjid/sPgV/7fQ63+902E6jr7nckkv60a6j37ge2GmNOu+u/o7V7aBR2j+NwO4+b1MHytsZj90B8jvpdHo3du33PdQdUAS+55T5njDHNftdrsXvplzHGFGL3qJ8ETonIsyIyvoN2/Q74C7GD438BvG+M8bXtc9iEdcB1v9wWwHr6t6MEu3f7d8AhY8yxLh7iLx1Y77c99mO/9GOB/8IGj2fFdrP9REQSA2jPYaAGu+e4GHgBOCEi04GbsHuX3VHud7nD98M56Xe5rp3rvseOB44ZY7x+tx/F/hIdjQ2AHX2O0oEFvm3mttsqYFwA6wJQalyU9Htu/8+N/+umYxN8md9r/QL7y+DienTQzrYC/Q61NR6oNMacb/M6E/yut32PkkQkoZvfkYihiaAXxPbF3wPcJLbPvBz7E/lKEbkS+/O6njb9sc6xDpa3VYb9wPuk+V0+jQ0Gs4wxw9xfqjGms8Di77LSs8aY3xljbsB+YQ3w43YfaMw+7JfnVmx3ye/8bjtkjLkP++X+MfDfIjI4wDb5/CfwDfe/O44Bt/ptj2HGmCRjTKkxpskY85QxZia2i+42bDdLIF4FPgkMMMaUuuufBYYDH3TwmL4s7XsCmCSXDsamYbtOKrDdNx19jo4Br7bZZinGmEcCfO0JIiJtnvuE33X/7XAM+4tglN9rDTXGzHK3d/Z5b6uz71Bn2/4EMEJEhrR5ndJOHtP6xAF+RyKJJoLeycPubc7E7i3OBWYArwOfdXtn/wH8oxuciheRhW4veh1wi4jcIyIJbjBvbjuv8RzwgIjMFJFk4Pu+G9zz/xvwTyIyBi4OZq4MsP0ngSm+KyIyXUSWuvbVY5NMZ4fb/Q74CnAjdozA9zx/JSKjXfuq3OLuHrb3B2y3wHPdfNzPgbUiku7aMlpE7nSXbxaRbDcweQ7bXdLibntSRHZ28ryvAl/CdnuB7aL4MrZ7paN1u2T7htgubPfct0QkUUSWALcDz7r2/Q/wpIgki8hMLj2o4QVgmoh8xj02UUSuFZEZAb72GOAr7nF3Y78Dm9u7ozGmDNgK/IOIDBU7uJ8pIje5uzznnmuiiAwHHu/kdX8JfFNErhYry/e+08m2d78w3wJ+KPYAgjnYX7HrulrRHnxHIoImgt65H/iVsccol/v+gJ9ijwxJAL4JFGCPeqnE7j3Eue6PHOxebyV2r/LKti9gjHkROyD8Mra75OU2d/m2W/4nETkHbMcOaAXi34GZ7id6PnZw9UfYXxrl2C/4dzp5/O+x/eIv+3WNgR3g3isiNcA/Y/vA6wHcURyLu2qYMabOGLPdGFMX4Lr4/DN28HyriJzHDhwvcLeNA/4bmwT2Y4P7b91tk4A3O3neV7FjIb5E8Aa2W+61Dh9h+52/67bvN7u5Ht1ijGnEHrp8K/b9+1fszsgBd5cvYbuRyrF9/r/ye+x5bNK9F7u3XI79nAZ6TswuYKp73bXAJ40xZzq5/2exA9r7gLPY98Tjbvs3bPfdh8D72ATWLmPMH93r/Q47MJ6PHWCGrrf9fdhxgxPAeux4yLauVpTuf0ciglzatadUbBKRD7CDzJ0FMNWG2EOeH3ZdJSpChf1Ud6X6A2NMe91ySsUE7RpSSqkYp11DSikV4/QXgVJKxbiIGCMYNWqUycjICHczlFIqorz33nunjTGju7pfRCSCjIwMdu/eHe5mKKVURBGRzs7Mvki7hpRSKsZpIlBKqRiniUAppWKcJgKllIpxmgiUUirGaSJQSqkYp4lAKaVinCYCpZSKcZoIlFKqPyhaB/kZ8Ls4+7+oy3lygiYizixWSqmoVrQO3lkNLbX2eu1Rex1g8qqQv7z+IlBKqXD7cE1rEvBpqbXL+4AmAqWUCrfaku4tDzJNBEopFW7Jad1bHmSaCJRSKtyuXAvxyZcui0+2y/uAJgKllAq3yatg/jOQnA6I/T//mT4ZKAY9akgppfqHyav6LPC3pb8IlFIqxmkiUEqpfuR03ek+f01NBEqpmLauYB0ZT2cQ91QcGU9nsK6g787o9alrrmPj4Y18bsvnWPbHZRw7f6xPX1/HCJRSMWtdwTpWb1xNbZM9meto9VFWb7Rn9K7KDm1/vTGGDys+JL8wny3FW6hpqmFiykQeufIRBicODulrt6WJQCkVs9bsWHMxCfjUNtWyZseakCWCU7Wn2Hh4IxsOb6CouohBCYNYnr6cvKw8rh57NXHS9x01mgiUUjGrpLr9M3c7Wt5TjS2N7Dy2k/zCfN488SZe4+WqMVfx4KIHWZGxos9/AbSliUApFbPSUtM4Wn203eXBcKDyAPmF+Ww6somqhirGJI/hodkPcWfmnWSkZgTlNYJBE4FSKmatXbb2kjECgOTEZNYu6/kZvWfrz7K5aDP5hfkcqDxAYlwiS9OWkpeVx0LPQuLj4oPR9KDSRKCUilm+cYA1O9ZQUl1CWmoaa5et7fb4QLO3mbdOvEV+YT6vHHuFZm8zM0fO5DsLvkPO5BxSB6aGovlBI8aYcLehS9dcc43ZvXt3uJuhlFKXKKouIr8wn42HN1JRV8HwgcPJnZJLXlYe00dMD3fzEJH3jDHXdHU//UWglFLdUNNYw0vFL5FfmM+HFR8SL/EsnrCYvKw8bpx4I4nxieFuYrdpIlBKRZx1Bet63Z3THV7j5d3yd8kvzGf70e3Ut9STmZrJN67+Brdl3saoQaNC9tp9QROBUiqi9OVJYKU1pWwo3MDzh5+ntKaUIYlDuCPzDvKy8pg9ajYiEtTXCxcdI1BKRZSMpzPaPeQzPTWd4q8V9/r565rr2H50OxsKN7CrfBeCsMCzgLysPJalLSMpIanXr9FXdIxAKRWVQnESWEflHh6d+yh3Zt6JJ8XT4+eOBJoIlFIRJZgngfXHcg/hoIlAKRVRensSWHvlHuaNmcdTi55iZcbKsJd7CAdNBEqpiNLTk8AipdxDOOhgsVIqakViuYdg0sFipVRMivRyD+GgiUApFRWOVB8hvzCfFw6/cLHcw73T7+035R76s5AnAhGJB3YDpcaY20RkMvAsMAJ4H/iMMaYx1O1QSkUfX7mH9YXr2VOxJyrKPYRDX/wi+CqwHxjqrv8Y+CdjzLMi8nPgc8DP+qAdSqkoEO3lHsIhpIlARCYCucBa4DGx52MvBT7t7vIb4Ek0ESiluhAr5R7CIdS/CJ4GvgUMcddHAlXGmGZ3/Tgwob0HishqYDVAWlpwZgtSSkUWX7mH/MJ83il/52K5h6/M+wpL05ZGVLmH/ixkiUBEbgNOGWPeE5ElvsXt3LXd41eNMc8Az4A9fDQkjVRK9Tv+5R5eKn6JC00XYqrcQziE8hfB9cAdIpIDJGHHCJ4GholIgvtVMBE4EcI2KKUihK/cQ35hPsXnimO23EM4hCwRGGOeAJ4AcL8IvmmMWSUifwQ+iT1y6H5gQ6jaoJTq39or93DVmKt4aPZDrMhYEZPlHsIhHOcRfBt4VkR+APwZ+PcwtEEpFUYHKg+w/tB6NhVtorqhWss9hFmfJAJjzE5gp7t8BJjfF6+rlOo/Yr3cQ3+mZxYrFcX6ekrHtrTcQ2TQRKBUlOrLKR3baq/cw31X3EdeVh7Thk8L6Wur7tNEoFSUWrNjzSU1+wFqm2pZs2NNSBKBlnuIXJoIlIpSoZjSsS0t9xAdNBEoFaWCOaVjW1ruIbpoIlAqSvV2Sse2fOUeNhRuYFf5Li33EEU0ESgVpXo6paM//3IPW4q3UNNUo+UeopBOVamUuoyv3MOGwxsoqi7Scg8RSqeqVEp1S1NLEzuP72T9ofUXyz3MGzOPpxY9xcqMlVruIYppIlAqxh2oPEB+YT6bjmyiqqFKyz3EIE0ESvUTgZ4FHIyzhavqq9hUtEnLPShAE4FS/UKgZwH35mzh9so9zBgxgyfmP0HulFwt9xDDdLBYqX4g4+mMdo/5T09Np/hrxd2+n7+i6iLyC/PZeHjjxXIPuVNyycvKY/qI6cFaBdUP6WCxUhEk0LOAA72fr9xDfmE+H1Z8qOUeVKc0ESgVRr7+ftP+jK2XnQXc2dnCXuNld/lu1heu13IPqls0ESgVJm37+9tq7yzg9s4WHpo0lFtn3krO/+RcLPdwe+bt5GXlkT0qW8s9qC5pIlAqTNqrDuqTnpre7tFA/mcLVzdXM27YOBISE3iz7E0WeBbw5XlfZlnaMi33oLpFE4FSYdJRf78g7Q78+so9HDp/CM8YD0OahjAhZQJ5WXnckXkH41PGh7jFKlppIlAqTAKtDlpRW8HGIxvJL8zXcg8qJDQRKNUNwZz6sbPqoL5yD/mF+bxZ+iYtpkXLPaiQ0USgVICCPfVje9VBv7LwKxyrO8bSPy615R4GjeHB2Q9quQcVUnpCmVIB6snJXIFor9zDzZNu5q6pd2m5h1hStA4+XAO1JZCcBleuhcm9m1JUTyhTKsiCOfWjf7mHncd20uRtuljuIWdyDsOShvW2uSqSFK2Dd1ZDi+smrD1qr0Ovk0EgNBEoFaBgTP3YXrmHT03/lJZ7iHUfrqG2qYVXzl3PtnPX8aOJ/4ckau0vBE0ESvUfPZ36sb1yDzdMuIG8rDxumniTlnuIYbWNzbx84BSb99/HK+euoc4kMTK+iiMNE5k5qMh2E/UBTQRKBag7Uz+2V+5hSuoUHrv6MW7PvF3LPcSwCw0u+BeU8crHp6hv8jIqcQ5/OWIHOalvsGDwXuLFa++cHPivzd7QRKBUN6zKXtXpEUKlNaU8X/g8Gw5voLSmlJTEFC33oLjQ0MyOA6fYvKeMnQdd8E8ZyN1XTyIn28N8XiJ+929axwgA4pPtgHEf0ESgVC/VNdex/eh2NhRuYFf5LgTRcg+KmoZmduw/yeaCMnZ+XEFDs5fRQwZyzzU2+F+bMYL4ON+OwSqII+hHDQVKE4FSPeAr95BfmM+W4i3UNNUwIWUCj859VMs9xDBf8N+0p4xXD9rgP2bIQO691gb/ay4J/m1MXtVngb8tTQRKdUNFbQXPH7ZdP1ruQQGcr29ix/5TbCqwwb+x2cvYoQO5b36aDf7pw4nrKPj3E5oIlOqClntQbZ2vb2L7/pNs2lPOa4dag/+n56eRO8fD1Wn9P/j700SgVAcOVB4gvzCfTUc2abkHxbn6Jrbvs33+rx08TWOLl3FDk1i1II3cbA9XRVjw96eJQCk/7ZV7WJq2lLysPC33EIOq61qD/+uHbPAfn5rEZxamk5PtYd6kYREb/P1pIlAxz7/cwyvHXqHZ23yx3EPulFxSB6aGu4mqD1XXNbHtYvCvoKnFMGHYID67MJ2cOR7mToyO4O9PE4GKaF/82TqeObKGlsElxF9IY/WUtfzrI4EdeXGk+gj5hfm8cPiFi+Ue7p1+b6flHtatgzVroKQE0tJg7VpYFZ4DPVQQVdc2sXVfOZsLynij8PTF4P/Aogxysj3MnTQsqs8B0USgItYXf7aOn5WuhhR7Ek5LylF7/Wd0mAzaK/eweMJi8rLyuHHijZ2We1i3Dlavhlp3zs/Ro/Y6aDKIRFW1jWx1e/5v+gX/B6+fTE62hysnpkZ18PenZahVxEr46wxaUi4vAhdfk07z/y6+eN1rvLxb/i75hfmXlHu4K+subsu8LeByDxkZNvi3lZ4OxcWXL1f9z9kLjWzbd5JNLvg3ew0Thw8iN9tDTraHOVEW/LUMtYp6LYPbL8jlW15aU8qGwg08f/h5SmtKGZI4pFflHko6qP/V0XLVP5y90MiWveVsKijj7cNnaPYaJo0YxOcWTyY320P2hOgK/j0RskQgIknAa8BA9zr/bYz5vohMBp4FRgDvA58xxjSGqh0qesVfSLvsF4GIMCx+Cg9veTjo5R7S0tr/RZDWN3XBVDdUuuC/uaCMtw6focVrSBuRzMOLp5Cb7WH2hKExH/z9dSsRiMhwYJIxZk8Ad28AlhpjakQkEXhDRF4EHgP+yRjzrIj8HPgc8LPuNlyp1VPW2jGBxFoGDRzE8JThpCYPIz4+jtKaUh6d+yh3Zt6JJ8UTlNdbu/bSMQKA5GS7XIXfmZoGtuy1ff5vH7HBP31kMqtvtMF/1ngN/h3pMhGIyE7gDnffD4AKEXnVGPNYZ48zdvChxl1NdH8GWAp82i3/DfAkmghUDzx1/wqOrXuIj707GJgUj9drGOfN5kcrHwtJuQffgLAeNdR/nKlp4CW35/+nI5W0eA0ZI5P5/I1TyNHgH7BAfhGkGmPOicjDwK9c904gvwgQkXjgPSAL+BfgMFBljGl2dzkOTOjgsauB1QBp+ttbOZeVexjQwnVjriEvK69Pyj2sWqWBP9xO1zTw0ke+4H8Gr4HJowbzhZts8J/p0eDfXYEkggQR8QD3AGu68+TGmBZgrogMA9YDM9q7WwePfQZ4BuxRQ915XRV9tNxDbKs47/b895Sxq8gG/ymjBvPFJVnkZHuY4Rmiwb8XAkkETwFbgDeMMe+KyBTgUHdexBhT5bqYrgOGiUiC+1UwETjRzTarGKHlHmLbqfP1bPnIHu3zTlElXgOZowfzpZuzyJnjYfpYDf7BEkgiKDPGzPFdMcYcEZF/7OpBIjIaaHJJYBBwC/Bj4BXgk9gjh+4HNvSo5SoqabmH2HbqXD0vuuD/bnElxkDWmBS+tHQqudkepo1N0eAfAoEkgv8LXBXAsrY8wG/cOEEc8Jwx5gUR2Qc8KyI/AP4M/Hs326yiUFF1EfmF+Ww8vDHgcg8qOpw8V8+LBWVsLijn3aM2+E8dk8JXlk4ld46HaWOHhLuJUa/DRCAiC4FFwGgR8T9CaCjQ5W9yd4jpvHaWHwHmd7+pKtq0V+7hhgk3cFfWXV2We1CRrby6nhc/KmNzQRm7j57V4B9mnf0iGACkuPv4vyvnsF07SnWb13jZXb6b9YXrLyn38NjVj3F75u0Bl3tQkae8up7NBa3BH2Da2BS+usx2+0zV4B82HSYCY8yrwKsi8mtjTDvnUyoVuNKaUp4vtFM8ltaUkpKY0qtyDyoylFXXsbnAHur5ngv+V4wbwmPLp5GTPY6sMRr8+4NAxggGisgzQIb//Y0xS0PVKBUd6prr2H50OxsKN1ws9zDfM58vzfsSy9KWMShhULibqELgRFXdxT3/90uqABv8v7F8GjlzPGSOTglzC1VbgSSCPwI/B34JtIS2OSrSGWPYc3oP6w+tZ0vxFmqaapiQMoEvzv0id2beyfiU8eFuogqB0qo6XiwoY1NBGX92wX+GZyjfXDGNnGwPUzT492uBJIJmY4yWgFCdqqitYOORjeQX5lNUXcSghEEsT19OXlZeSMo9qPA7fraWFwvsoZ4fHLPBf6ZnKH+9cjq3zh6nwT+CBJIINorIF7FnBjf4FhpjKkPWKhURLiv3YFqYN2YeTy16ihXpK0gZoIEg2hyrrL3Y7fPh8WoAZo23wT8320PGqNCW+FChEUgiuN/9/2u/ZQaYEvzmqHBbV7CONTvWUFJdQlpqGmuXrWVV9qXFddor9/DArAe4M+tOJqdODlPLVagcq6xlkwv+e1zwnz1hKN/6xHRyZmvwjwZdJgJjjH6zY8S6gnWs3ria2iZbZ/lo9VFWb7RzMeZOzb2s3MPNk24mLyuPReMXabmHKFNypjX4F5Ta4J89IZVvf+IKcrLHkT5Sg3806XKqShFJxs4hkGaMWS0iU4HpxpgX+qKBoFNV9pWMpzM4Wn3pkcIpg1KYMHwCyUnJNHmbmDFiBnlZeeRMzmFY0rAwtVSFwtEzFy4G/49KzwEwZ2IqOdkecmZ7SBuZHOYWqu4K5lSVv8KWkl7krh/HHknUZ4lA9Y2Sajvn4oCEAQwfMpxhg4eRmJBIc0szn5r+KS33EIWKT7cG/70nbPC/cmIqT9x6BTnZHiaN0OAfCwJJBJnGmE+JyH0Axpg60bN/ok5NYw2ZozNpiW9hcNJgjDGcrztPWWUZwxOH8+353w53E6PeunV9M+lN0ekLbC4oY9OeMvaV2eA/d9Iw1uTM4NbscUwcrsE/1gSSCBpd9VADICKZ+B09pCKXr9xDfmE+245uI2lwEo1NjZRXllN1oYrmlmaSE5NZ+wmdizHU1q27dBrMo0ftdQhOMjhSUWODf0E5+13wn5c2jO/mzuATszX4x7pAxgiWA98FZgJbgeuBB4wxO0PeOkfHCIKrvXIPt06+lbysPPaU7eG7L3+306OGVPBlZNjg31Z6OhQX9+w5D1fUsHmPPcnrQPl5AK5KG0ZOtodbsz1MGKZndke7QMcIukwE7slGYieVEeBPxpjTvW9i4DQR9F59cz3bS7aTX5jPrrLWcg95WXla7qEfiIuD9r6KIuD1Bv48hadqLh7n7wv+V6cPt8F/9jjGa/CPKb0eLBaRtvMNlLn/aSKSZox5vzcNVKHnK/eQX5jPS0UvabmHfiwtrf1fBIFM11146jyb9tjCbh+ftMH/mvThfO+2mdyaPQ5PqgZ/1bnOxgj+wf1PAq4BPsT+IpgD7AJuCG3TlL9ATvTy8ZV72FC4gSPVR7pV7qGvBizVpdauvXSMACA52S5vz6GT5y8e7XPwZA0iNvh///aZ3Drbw7jUpL5puIoKnZWhvhlARJ4FVhtjCtz12cA3+6Z5Cjo/0cuXDJpamnj1+KusL1x/sdzD3NFzu1XuIdQDlqpjvu3bWRI+ePI8m/bY4H/olA3+16aP4MnbZ3JrtoexQzX4q54JZLD4A2PM3K6WhVKsjxG0d6IXQHpqOls+u+ViuYezDWcZM2gMt2fe3qNyD6EYsFQ9Z4zh4Mmai3v+hb7gnzGCXNfnP0aDv+pEME8o2y8ivwR+iz2E9K+A/b1sn+oG34lePvFx8aQOTiUhOYFPbvzkJeUeFo5fSEJcIG9rO69T0r3lKviMMXx88vzFo30OV1xABOZnjOCzd87iE7M0+KvgCyRiPAg8AnzVXX8N0LLUfSgtNY2j1UdJGZTC8JThDEkeQpzE0dLcwhPznwhauYfeDFiqnjPGcKD8vDvOv4wjFReIE1gweSQPXD+ZlbPGMmaIBn8VOoEUnasH/sn9qT5WXF3M8iuW80bpGyQkJNDc0kzl+UpqaupFGeDlAAAfjElEQVT5+W0/5dMzPh201+rugKXqOWMM+8vOXzzU88hpG/yvmzKSh66fzMpZ4xg9ZGC4m6liRGeHjz5njLlHRApwZxX7M8bMCWnLYlhNYw1bim3f/wcVHyDEU39sOpV8xHmKMNWTSHz97yFtFWQH73UDGbBUPWeMYV/ZORf8yylywX9h5kg+t9gG/1EpGvxV3+twsFhEPMaYMhFJb+/2vpzQPhYGi73Gy3sn32P9ofVsL9lOXXMdk1Mnk5eVx9/eezvFe0df9hgdxO3/jDHsPXHu4p5/8Zla4uOEhVNGkpPtYeWssYzU4K9CpNeDxcYY3wlkcUCZ6yLC1R0aG5RWKk7UnGDD4Q1sKGwt95A7JZe8rDzmjJqDiPDwvvYfq4O4/ZMv+PuO9jnqgv+izJF8/qZMVs4ax4jBA8LdTNVfNdVBYt+eBBjo5PWL/K63uGXXhqRFMcC/3MM7Ze9gMCzwLOBL877UbrkHHcTt/4wxfFTaGvxLKluD/yM3ZbJCg7/qTOUROLgVDm2B4jfhqx/A0L478z+QRJBgjGn0XTHGNIqIfqK7qaNyD4/MfaTLcg86iNs/GWMoKK2+GPyPVdaRECcsyhrFozdnsmLmOIZr8FftaW6Ekrdag/+ZQrt8YBwMrYeXroNrfwiT+2aALpBEUCEidxhjngcQkTuBPi06F8lO151m4+GN5Bfmc6T6CEnxSSxPX85dU+/qstyDjw7i9h/GGD48Xn2xz//4WRv8r88axZdvnsqKWWMZlqzBX7XjfDkc2mYD/+Gd0Hge4gdCxg2QOQ/O/B7i3d5e8zF4x53W3wfJIJAzizOBdYBvl/U48FljTGGI23ZRpA0W+8o95Bfm80bpGxfLPeRl5bEyY2VA5R5U/2GM4YNjVReP9imtqiMx3gb/nGwPK2Zq8Fft8HrhxPtwcIsN/mUf2uVDJ8DUFfZvyk0wYDDkZ0BtO/2/yemQV9zjJgTtzGJjzGHgOhFJwSaO8z1uVZT7uPLjy8o9PDDrgR6Ve1DhZYzhz8eq2LynjBc/ag3+i6eO5uvLp7F8xlhSkxPD3UzV39RVweGX4dBWu/dfexokDibOh2Xfg6krYewsW1/cX20HR350tDzIukwEIvJ3wE+MMVXu+nDgG8aY74a6cZGgqr6KTUWb2FC4gf2V+0mMS2TJpCXkZeWxaPyiHpd7UH3P63XBv6CMFwvKOFFdz4D4OBZPHcVjy6dxy8yxpA4KbvDXaq8902+2mzFQccDt9W+Fkj+BaYFBwyHrFhv4s5ZB8ojOnyc5jXXbjrLmOSg5DWmjYO09sGp53xwREkiUutUY8x3fFWPMWRHJwc5aFpNavC28deIt8gvzeeXYKzR5m5gxYgaPz3+c3Mm5QSn3oPqGDf5n2bSnnBc/KqPMBf8bp43imyunc8vMsQxNCs2ev1Z77Zmwb7emOih6zQb+g1vBVwtsbDbc8DUb/CdeA3HxAT/luuIcVv/yZ9S6w3KOnobVvwSm5tAXqxTIGMEe4FpjTIO7PgjYbYyZ1QftA/rPGEFxdTH5hflsPLyRU3WnGDZwGLdNuY28rDymj5ge7uapAHm9hvdLzrKpoIwXC8opP+cL/qPJnTOOZTNCF/z9abXXngnLdqsqad3rL3oNmushMRmmLGnt70+d0OOnz8jI4Gg7K5Wenk5xL1YqmNVHfwvsEJFfuesPAr/pccsiTNtyD/ESzw0TbuCJrCe4aeJNJMZrP3Ek8HoN75WcZdOeMl76yAX/hDhumjaax7OvYNmMMQzpg+DvT6u99kyfbLeWJji2q3Wvv8IVXB4+Ga5+wAb+9OshMTjFAEs6aHxHy4MtkMHin7hfBbdgZyh7CWi37ES08Bovu8t3k1+Yf0m5h69f/XVun3I7o5MvL/eg+h+v17D76Fnb5/9RGSfPNTAgIY4l00bzxJwrWHpF3wd/f3qiYM+EbLtdON16eGfhy9BQDXEJkL4I5v0VTFsJI7MuH+gNgrS0tHZ/EaT10Ych0JHMcsAL3AMUAf8vZC0Ko0DKPaj+rcVr2F1c6YJ/OafONzAwIY4l00eTk+1h2YyxpAzsHwP4PT1RsN8MlIZJ0E6w9Hqh/EN3UtdWKH0PMDB4DMy4HaatgCk3Q9LQYDa/XWvXrmX16tXU+q1UcnIya/vorNHOqo9OA+4F7gPOAH/Ajinc3Cct6yN1zXXsKNkRcLkH1f+0eA3vFNng/9Lecipc8L95+hhy5nhYesWYfhP8/fXkRMGwD5T2A706wbLhPBx+xe71H9oGNScBgQlXwZInYOpy8MyFuK5P9AymVa7xa9asoaSkhLS0NNauXXtxeah1Vn3UC7wOfM538piIHDHGTOmTlvkJ9mBxR+Ue7sy6s8tyD6p/aPEadhWdscH/o5OcrmkgKdEF/2wb/Af3w+DfWzrA3E3G2PINh7bawd6jb4G3CQamQtZS29eftRxSorO7NxiDxX+J/UXwioi8BDyLHSOIWBW1FWw8spENhRt6XO5BhU9zi5d3iirZVFDGlr3lnK5pJCkxjqVX2OB/8/ToDP7+dIA5AM0NUPxGa/A/W2SXj74CrnvE9vVPWgB6oMdFnZWhXg+sF5HBQB7wdWCsiPwMWG+M2dpHbewVX7mH9YXrebP0zYvlHp5c+KSWe4gAzS1edvmC/0flnLnQyKDE+Nbgf8VokgdEd/D3pwPMHTh3ovUInyM7oekCJCTB5Bth4aN2z394VB/j0iuBHDV0AVtraJ2IjADuBh4HOk0EIjIJ+E9gHHag+RljzD+75/gDkAEUA/cYY872Yh069C8f/At/OPAHLfcQYZpbvPzpiA3+W/f6Bf8ZY8jN9rBkemwFf3/dHSiN2oFlbwsc3237+g9uhZMFdnnqJLjyXrvXn7EYBiR3/VxRu5EC161vkzGmEviF++tKM7YUxfsiMgR4T0S2AQ8AO4wxPxKRx7FJ5dvda3ZgahpruHbcteRl5bFw/EIt99CPNbd4efuI7fPfsvcklRcaSR5g9/xt8B/DoAGBn6kZrbozUBp1A8u1lbaOz8EtULgd6ipB4iHtOrjlKbvXP2ZG9w7vjLqN1DNdnlkctBcS2QD81P0tcdNgeoCdxphOT8vt6WCxMUYP+ezHmlq8vH3YF/zLOVvbRPKAeJbNGEtu9jiWTB9DUqIG/56K+IFlY+Dk3ta9/uPvgPFC8kg7wDttBWQutXV9eiriN1LnAh0s7pNEICIZwGvAbKDEGDPM77azxpjL3kkRWQ2sBkhLS7u6vZMtVORpavHy1uEzbN5TxpZ95VTVNjHYBf8c1+2jwT844uJsLG1LxB5C3y81XrAlHA66wzvPHbfLPVfaGj7TVsL4ed2q49OpiNxIgQtmiYneNiQFewLa14wx5wLdQzfGPAM8A/YXQehaqEKtqcXLm4Wn2VxQxtZ9J6mqbSJlYALLZtgB35umafAPhYgZWK4saj2jt+h1aGmAASm2js+Sb9u9/6Ge0Lx2xGyk0AppIhCRRGwSWGeM+R+3+KSIePy6hk6Fsg0qPBqbLw3+1XU2+N/igv+NGvxDrt9OcdrSBCVvtxZxO33QLh+ZBdc+bLt80hZCwsCOnyNYA7z9diP1rZAlArG7/v8O7DfG/KPfTc8D9wM/cv83hKoNqm81Nnt5o7CCTXvK2bavnHP1zQwZmMAtM223z+KpozT496F+NcVpzanW4/qP7ISGcxA/wBZuu+YhO9A7MjOw5wrmAG+/2kjhE7IxAhG5AXtmcgH28FGA7wC7gOeANKAEuNsdjdSh/lKGWl2uobmFNw6dZlNBGdv2neR8fTNDkhJYPnMsObM9LJ42ioEJGvxjjtcLZX9unZz9xJ/t8iEeG/SnrYTJN8HAHpzHE+UDvMEU9jECY8wbdHwm8rJQva4KvYbmFl4/aLt9tu1vDf4rZo4jd844rs/S4B+T6qvd4Z1boXAbXKgABCZeC0u/awd7x2X3vnqnnl4ddHpgvQpIfVMLrx+ywX/7vpOcb2hmaFICK2eNIzfbw/VZoxiQoCU6Yooxtn//4jSNb4O3GZKGuWkaV9j/g0cG93V1gDfoNBGoDtU3tfDawQob/PefoqahmdRBiXxi9jhy5ni4PlODf8xpqrN1fHzBv8oF5DGzYNGX3TSN10J8CEOLDvAGnSYCdYn6phZedcF/h1/wz8keR062h0Ua/GNP1TEb9A9thSOvQnMdJAyyh3de/1W75z9sUt+1Rwd4g04TgaK+qYWdH/uC/0kuNLYwLDmR3GwPOXM8LMocSWK8Bv+Y0dJsz+L17fWf2meXD0uHqz5j9/ozbgjaNI09smqVBv4g0kQQo2zwP8WmgnJedsF/eHIit185npxsDws1+MeWC2ds/Z5DW6BwB9RX2Wka0xbC8r+1R/mMmhaSaRpV+GkiiCF1jS288vEpNheU8fKBU9Q2tjBi8ADumGuD/3VTNPjHDGOgfE/r4Z3Hd2OnaRwNV+Ta7p7MmyEpNdwtVX1AE0GUq21s5pUDFWz+qIyX95+irskG/zvnTiA328N1U0aQoME/NjSctydzHdpqSzqcL7PLx8+Dm75tz+j1zOvzaRpV+GkiiEK1jc28fMDu+b9yoIK6phZGDh7AXVfZ4L9gsgb/mHHmsOvrd9M0tjTCwKF2b3/qSnt455Cx4W6lCjNNBFHiQoNf8P/4FPVNXkalDOAvr55ATraHBZNHEh+n/btRr7nBBnxfOYfKw3b5qOkwf7Xt609bqNM0qktoIohgFxqa2XHgFJv3lLHzoC/4D+TuqyeRk+1h/uQRGvxjwbkyv8M7d0JjDcQPhMmLYcEXYOpyGKGz8qmOaSKIMDUNzezYf5LNBWXs/LiChmYvo4cM5J5rbPC/NkODf8QKtKKmtwVK33cTtmyxg74AQydA9t2ujs+NMGCwfc5P3azH23dCZ6rURBARfMF/054yXj1og/+YIQO591ob/K/R4B/5uqqoWXfWHtZ5aKs9zLP2DEgcTFoAy75vg/+YmZce3qnTMHZJN5HVZ1NV9kYsVh89X9/Ejv2n2FRgg3+jC/452R4b/NOHE6fBP3q0V1FzdBzMHw1/MQ+O7QLTAoNG2K6eqW6axuQR3XtO0CqdfqJ9E4W9+qjqvvP1TWzff5JNe8p57ZAN/mOHDuTT89PInePh6jQN/lGrpMR+GycnwFT3NywOqIPG83DD1+1e/4SrA5+mUat0dkk3kaWJIMzO1TexfZ/t83/t4GkaW7yMG5rEqgVp5GZ7uEqDf3Q7e9R29zw0HMY2QaJAo4HDzfBaA9SPg31v9Oy5tUpnl3QTWZoIwqC6rjX4v37IBn9PahJ/dV06uXPGMW+SBv+o1dJku3l8dXwqDtjlg1rgvRY41ARHW6AFW1HzmR9e/hyBjm5qlc4u6SayNBH0keq6JrZdDP4VNLUYxqcm8ZmF6eRke5g3aZgG/2hVU2Enajm4BQ6/Ag3VEJcIA6fAyy2wtw4q3SR+ImCwndTtBfjujG5qlc4u6SaydLA4hKprm9i6r5zNBWW8UXiaphbDhGGDLpZ0njtpGKJFvKKP1wtlH9gyDoe22EM9MZAyzg70TltpSzhPz+7+SGW0j26qoNLB4jCprm1iiwv+b/oF/wcWZZA7ZzxXTkzV4B+N6s/BkVdap2msOYmdpvEauHmNTQCeKy89vLMnI5U6uqlCQBNBEFTVNrJ170k2ueDf7DVMHD6Ih66fTE62hzka/KOPMXD6kN3jP7QVjr4N3iZbrTNzmd3rz7oFBo/q+Dl6MlKpo5sqBDQR9NDZC41s3VfOpoJy3nLBf9KIQXxu8WRysz1kT9DgH3Wa6uFXfwvb/w3G1cEIV7hvzExY+Kg9tn/SgsCnaezJSKWObqoQ0ETQDZUXGtm6t5xNBWW8dfgMLV5D2ohkHl48hdxsD7MnDNXgH22qS91e/zY4uB1MI0w1UNQCbzfA8UT4+y/B8h6MLvZkpFJHN1UI6GBxF87UNLBlrz3a5+0jNvinj0wmJ9tDbraHWeM1+EeVlmYo3d16eOfJj+zy1DT400l4/ywUN0Oz32N0oFb1UzpY3Atnahp4aa8d8P3TkUpavIaMkcl8/sYp5Gjwjz61lbZ+z8EtcHiHresj8W6axr+xdftHT4f4eDs20JYO1KoIp4nAOV3TwEsf+YL/GbwGJo8azBdussF/pkeDf9Qwxu7p+/b6j78LxgvJo2DaJ1rr+AwadunjdKBWRamYTgQV592e/54ydhXZ4D9l1GC+uCSLnGwPMzxDNPhHi4YaKHrVBf9tcP6EXe6ZCzf+td3rH9/FNI06UKuiVMwlglPn69nykR3wfaeo0gb/0YN59GYb/K8Yp8E/alQeaZ2cvfgNO03jgCGQuQSmfsce2z9kXODPpwO1KkrFxGDxqXP1vLS3nE17yninuBJjIHP0YHKzPeTM8TB9rAb/qNDcCCVvtQb/M4V2+cip9rj+qStsv3/CgPC2U6k+ooPFwB/eLeH/vVfKu0dt8J86JoWvLJ1K7hwP08YOCXfzVDCcP+mmadwCh3faks3xAyBjsZ2jd+pyGDEl3K1Uql+L6kTwZuEZztY2avCPJl4vnHi/dXL2sg/s8nNeKEuCJV+AB79np2lUSgUkqhPBTz45h6TEACfxUP1XXRUcftnt+W+D2tN2msaBk+G1FthXBye9QA08/1NImaP99kp1Q1QnAk0CEcoYqPjYTc6+FUredtM0Drf1e6auhKxlMPMqOHrh0sfW1trBXE0ESgUsqhOBiiBNdVD0emsRtyp3ktbYbLjhazb4T7zm0mkatRKnUkGhiUCFT9Wx1r3+oteguQ4Sk22t/hses0f5pE7o+PF6gpdSQaGJQPWdlmY7TaMv+Ffst8uHT4ar77eBP/16SEwK7Pn0BC+lgkITgQqtC6cvreNTXw1xCZC+COb9lT2+f2TWpRO2BEpP8FIqKDQRqOAyBso+bD28s/Q97DSNY2HG7Xavf8rNkDQ0OK+3apUGfqV6SROB6r2G83ZSdt/hnTXlgMCEq2DJEzBtBYy7svM6PkqpsAlZIhCR/wBuA04ZY2a7ZSOAPwAZQDFwjzHmbKjaoELodKHr698CR9+y0zQOTIWspe7wzlsgZXS4W6mUCkAofxH8Gvgp8J9+yx4HdhhjfiQij7vr3w5hG1R3rVvXfp97c4Mt3HZoq/2rPGLvP/oKuO4R29c/aQHEJ4a3/UqpbgtZIjDGvCYiGW0W3wkscZd/A+xEE0H/sW7dpUfhVJbATz8PJ34OjYeh6QIkJMHkG+G6L9r+/uHp4W2zUqrX+nqMYKwxpgzAGFMmImP6+PVVZ777HRjRANcNhKkJMM6dvHXyI1jyoN3rz1gMA5LD206lVFD128FiEVkNrAZI0xOEQqe20tbxObgF7q6E5MHgNVDSAtvr4WAznDbw9/8Y7pYqpUKkrxPBSRHxuF8DHuBUR3c0xjwDPAN2PoK+amC/1VHffXcZAyf3tp7UdfwdN03jSDiRBB9UweFmqPd7TLp2/ygVzfo6ETwP3A/8yP3f0MevH5na9t0fPWqvQ2DJoPGCLeHgm6bx3HG7fNwcWPwNe5TPhKvg98/C+tVQ39z6WD1TV6moF7IZykTk99iB4VHASeD7QD7wHJAGlAB3G2Mqu3qu3s5QFvEyMtqvqZOeDsXF7T+msqj1pK7iN6ClAQak2Do+01ZC1nIY6rn8ccH65aGUCrtAZyiLiakqI15cnO3SaUvETtQCdprGY39ye/1b4fRBu3xEZus0jemLIGFg37VbKRVWOlVlNOmoyub0ifDn37o6Pq+0TtOYfj1c85AN/iMz+769SqmIookgEvhX2RwfB1MTYfoA8FTDhkdhiAdm32X7+qcsgYEp4W6xUiqCaCLo7+qrYW4y/M0NcGoXDDJggEGTYdEqG/zHZfeseqdSSqGJoP8xxvbv+/r6S94GbzMkDYNr/9J292TdAoNHhrulSqkooYmgP2iqh+LXW4N/lRsPGDMLFn3ZTdN4LcTr26WUCj6NLOFSday1gNuRV+00jQmDbB//9V+1e/7DJoW7lUqpGKCJoK+0NNuzeH0ndZ3aa5cPS4erPmP3+jNuCHyaRqWUChJNBKF04YydpvHQFijcAfVVdprGtIWw4gd2r3/UNB3oVUqFlSaCYDIGyve4M3q3wvF3AQODR8MVuTbwZ94MSanhbqlSSl2kiaC3GmrgyE67139oG5wvs8vHXwVLHoepy8EzT6dpVEr1W5oIeuLM4dYjfI6+CS2NMHCo3dufutIG/xSdakEpFRk0EQSiudEGfF8Rt8rDdvmo6TB/ta3lk7ZQp2lUSkUkTQQdOVfmd3jnTmisgfiBMHkxLPiC3esfMTncrVRKqV7TRODjbYHS992ELVvsoC/A0AmQfbfd6598IwwYHN52KqVUkMV2Iqg7aw/rPLTVHuZZewYkDiYtgGXft8F/zEw9vFMpFdViKxEYA6f2t07TeGwXmBYYNMLW75m2EjKXQvKIcLdUKaX6TPQngsZaO02j7/DO6mN2+bhsuOHr9tj+iddAXHx426mUUmES3Ylg49fgw99Dcz0kDrZ1fG78pg3+Q8eHu3VKKdUvRHciGJYGVz9gA3/GDTpNo1JKtSO6E8Hix8LdAqWU6ve07oFSSsU4TQRKKRXjNBEopVSM00SglFIxThOBUkrFOE0ESikV4zQRKKVUjNNEoJRSMU6MMeFuQ5dEpAI42sOHjwJOB7E5/V2srS/E3jrr+ka3YK5vujFmdFd3iohE0BsistsYc02429FXYm19IfbWWdc3uoVjfbVrSCmlYpwmAqWUinGxkAieCXcD+lisrS/E3jrr+ka3Pl/fqB8jUEop1blY+EWglFKqE5oIlFIqxkVcIhCRJBF5R0Q+FJG9IvKUWz5ZRHaJyCER+YOIDHDLB7rrhe72DL/nesIt/1hEVoZnjQIjIvEi8mcRecFdj/b1LRaRAhH5QER2u2UjRGSbW+dtIjLcLRcR+T9u3faIyFV+z3O/u/8hEbk/XOvTFREZJiL/LSIHRGS/iCyM1vUVkenuffX9nRORr0Xr+vqIyNddzPpIRH7vYln/+B4bYyLqDxAgxV1OBHYB1wHPAfe65T8HHnGXvwj83F2+F/iDuzwT+BAYCEwGDgPx4V6/Ttb7MeB3wAvuerSvbzEwqs2ynwCPu8uPAz92l3OAF91n4zpgl1s+Ajji/g93l4eHe906WN/fAA+7ywOAYdG8vn7rHQ+UA+nRvL7ABKAIGOSuPwc80F++x2HfQL3cuMnA+8AC7Jl4CW75QmCLu7wFWOguJ7j7CfAE8ITfc128X3/7AyYCO4ClwAuu/VG7vq59xVyeCD4GPO6yB/jYXf4FcF/b+wH3Ab/wW37J/frLHzDUBQmJhfVts44rgDejfX2xieAYNmkluO/xyv7yPY64riG42E3yAXAK2IbNilXGmGZ3l+PYDQ+tbwDu9mpgpP/ydh7T3zwNfAvwuusjie71BTDAVhF5T0RWu2VjjTFlAO7/GLe8o3WLlHWeAlQAv3Ldf78UkcFE7/r6uxf4vbsctetrjCkF/h4oAcqw38v36Cff44hMBMaYFmPMXOye8nxgRnt3c/+lg9s6Wt6viMhtwCljzHv+i9u5a1Ssr5/rjTFXAbcCj4rIjZ3cN9LXOQG4CviZMWYecAHbNdKRSF9fAFx/+B3AH7u6azvLImp93XjHndjunPHAYOxnu62wfI8jMhH4GGOqgJ3YfsNhIpLgbpoInHCXjwOTANztqUCl//J2HtOfXA/cISLFwLPY7qGnid71BcAYc8L9PwWsxyb8kyLiAXD/T7m7d7RukbLOx4Hjxphd7vp/YxNDtK6vz63A+8aYk+56NK/vLUCRMabCGNME/A+wiH7yPY64RCAio0VkmLs8CLuB9wOvAJ90d7sf2OAuP++u425/2djOteeBe93o/GRgKvBO36xF4IwxTxhjJhpjMrA/o182xqwiStcXQEQGi8gQ32VsP/JHXLpubdf5s+7okuuAate1sAVYISLD3R7ZCresXzHGlAPHRGS6W7QM2EeUrq+f+2jtFoLoXt8S4DoRSRYRofU97h/f43APovRg0GUO8GdgDzY4fM8tn+I2SCH2p+ZAtzzJXS90t0/xe6412PGFj4Fbw71uAaz7ElqPGora9XXr9qH72wuscctHYgfND7n/I9xyAf7FrVsBcI3fcz3ktkUh8GC4162TdZ4L7Haf63zsUTDRvL7JwBkg1W9Z1K6va+tTwAEXt/4Le+RPv/gea4kJpZSKcRHXNaSUUiq4NBEopVSM00SglFIxThOBUkrFOE0ESikV4zQRKOVHRIyI/Jff9QQRqZDWqq93iMjj7vKTIvLNcLVVqWBJ6PouSsWUC8BsERlkjKkDlgOlvhuNMc9jT+pRKmroLwKlLvcikOsuX3L2q4g8ICI/bfsAEckUkZdckbzXReQKt/xuV3/+QxF5rU9ar1Q3aSJQ6nLPYk/jT8Keyb6ri/uDnXD8y8aYq4FvAv/qln8PWGmMuRJbYE2pfke7hpRqwxizx80IdR+wuav7i0gKtoDYH20ZGcCWDwB4E/i1iDyHLTSmVL+jiUCp9j2PrR+/BFsDpzNx2Lryc9veYIz5gogswHY1fSAic40xZ4LdWKV6Q7uGlGrffwB/Y4wp6OqOxphzQJGI3A0X59i90l3ONMbsMsZ8DzvL1KROnkqpsNBEoFQ7jDHHjTH/3I2HrAI+JyK+iql3uuX/W0QKROQj4DVsRVWl+hWtPqqUUjFOfxEopVSM00SglFIxThOBUkrFOE0ESikV4zQRKKVUjNNEoJRSMU4TgVJKxbj/Hz5boXPd/mdcAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"predicted range for 8000 is from 23.0 to 46.0\n",
"Gelman says the true number was 22\n"
]
}
],
"source": [
"predicted=accidents.extract()['predicted']\n",
"answer=accidents.extract()['answer']\n",
"fig,ax=plt.subplots(1)\n",
"ax.scatter(airline_df['Miles'],airline_df['Fatal'],color='blue')\n",
"x=np.linspace(3000,8000,10)\n",
"ax.plot(np.linspace(3000,8000,10),(.0042)*np.linspace(3000,8000,10))\n",
"ax.plot(x,.0037*x-6)\n",
"ax.plot(x,.0048*x+6)\n",
"l=np.percentile(predicted,2.5,axis=1)\n",
"ax.scatter(airline_df['Miles'],np.percentile(predicted,2.5,axis=0),color='red')\n",
"ax.scatter(airline_df['Miles'],np.percentile(predicted,97.5,axis=0),color='green')\n",
"ax.scatter([8000],answer.mean(),color='orange')\n",
"ax.scatter([8000],np.percentile(answer,2.5),color='orange')\n",
"ax.scatter([8000],np.percentile(answer,97.5),color='orange')\n",
"ax.set_title('Accidents vs. Miles, with model predictions')\n",
"ax.set_xlabel(\"Miles\")\n",
"ax.set_ylabel(\"Accidents\")\n",
"ax.scatter([8000],[22],color='black')\n",
"plt.show()\n",
"print('predicted range for 8000 is from ',np.percentile(answer,2.5),' to ',np.percentile(answer,97.5))\n",
"print('Gelman says the true number was 22')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is more promising, because the blue dots lie inside the sampled region.\n"
]
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"For completeness, here is the Gamma prior we are using. "
]
},
{
"cell_type": "code",
"execution_count": 100,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x1a1ac97748>]"
]
},
"execution_count": 100,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl0XGed5vHvr7RauyzJsizbkrwkzooDBgIxzIGwBEiTzGnoDs2BMB1O6GmaYek5Q5imZ3pm+vQAc5rtNEMTCNOBQ0No6B7SrBOykDhssbM4CbZjy0ssydps7Xup3vmj7nXKimRJdlXdpZ7POTqqunXr3p+vyo9evfe99zXnHCIiEl+JoAsQEZHcUtCLiMScgl5EJOYU9CIiMaegFxGJOQW9iEjMKehFRGJOQS8iEnMKehGRmCsOugCAxsZG197eHnQZIiKRsm/fvkHnXNNy64Ui6Nvb29m7d2/QZYiIRIqZnVjJeuq6ERGJOQW9iEjMKehFRGJOQS8iEnMKehGRmFPQi4jEnIJeRCTmFPQxtefwIHuPnwm6DBEJgVBcMCXZk5xP8b9+doivPHyUuooS9nz89VSV6ccsUsjUoo+R0+MzvPfrv+UrDx/lTZc3Mzw5x92/PB50WSISMAV9jHzgm/vYe2KIz7zjau587y6u37GOOx8+yuj0XNCliUiAVhz0ZlZkZk+Y2Q+95x1m9hszO2xm95hZqbe8zHt+xHu9PTelS6bB8Rn2nhjiQ6/bxh/s2gTAR994CSNTc/zDo8eDLU5EArWaFv2HgQMZzz8NfM45tx0YAm7zlt8GDDnntgGf89aTHHv0yCAAr73khRvZXdlayxsvb+arjxxlZEqtepFCtaKgN7ONwNuAr3nPDXg98D1vlbuBm73HN3nP8V6/3ltfcuiRw4PUrinhytbac5Z/5A3bGZtOcteeYwFVJiJBW2mL/vPAfwJS3vMGYNg5l/SedwGt3uNW4CSA9/qIt77kiHOOPYcH2b2tkaLEub9Tr9hQy/U71vHdx07inAuoQhEJ0rJBb2Y3Av3OuX2ZixdZ1a3gtczt3m5me81s78DAwIqKlcV1DozTOzrN7u2Ni77+by5tond0mu7hqTxXJiJhsJIW/XXA283sOPAd0l02nwfqzMwfoL0R6PEedwGbALzXa4EXXbnjnLvTObfLOberqWnZCVLkPB5+Lt0/v3vb4kH/srZ6APYeH8pbTSISHssGvXPuE865jc65duAW4AHn3LuBB4F3eKvdCvzAe3yv9xzv9Qec+gxyas+RQToaK9m0tmLR13esr6GqrJjHdKWsSEG6mHH0Hwc+ZmZHSPfB3+Utvwto8JZ/DLjj4kqU85lNpvj10dNLtuYBihLGS9vq2XdCLXqRQrSqa+Odcw8BD3mPjwKvWGSdaeCdWahNVuDx54eYnJ3nNUv0z/te3lbPZ3/+HCOTc9RWlOSpOhEJA10ZG3F7Dg9SlDCu3Xr+gU272tfiXPoXg4gUFgV9xD1yZJCdm+qoKT9/K33npjqKE6Z+epECpKCPsOR8igM9o+zyRtWcz5rSIq5orWWv+ulFCo6CPsJODk0xO59i67qqFa3/8rZ6njo5zExyPseViUiYKOgjrLN/HIBtKwz6Xe1rmUmmeKZ7NJdliUjIKOgjrHMgHfRbm1Ya9P6FU+qnFykkCvoIO9I/TlN1GbVrVjZcsrGqjI7GSvXTixQYBX2EdQ6Ms7WpclXveVlbPU9oiKVIQVHQR5Rzjs6BiRV32/h2rK9mcHyW0+MzOapMRMJGQR9Rg+OzjEzNrfhErO+S5moAnusbz0VZIhJCCvqIWu2JWJ8f9If7x7Jek4iEk4I+os4G/Spb9M01ZVSXF/Ncn4JepFAo6CPqSP84FaVFtNSUr+p9ZsYlzdXquhEpIAr6iOocmGBLUyWJxOqn472kuYrn+sY0taBIgVDQR1Rn//iq++d929dVMzw5x4BG3ogUBAV9BE3OJukenmLbBQb92ROy6r4RKQgK+gg6OjABrP5ErO+S9en36YSsSGFQ0EeQP+JmtWPofU1VZdRVlOiErEiBUNBHUGf/OAmDtobFJwNfjplxybpqDqtFL1IQFPQR1Dkwwea1FZQVF13wNrZr5I1IwVDQR1DnwPgFd9v4LmmuZnQ6Sf+YRt6IxJ2CPmJSKcfRwQm2XOCIG9/2Zp2QFSkUCvqIGRyfYTaZYlP9movajj/E8lCvgl4k7hT0EdM9PAXAhrqLC/rGqjLWVpZqLL1IAVDQR4wf9K0X2aIH2L6uiud0F0uR2FPQR0z3kBf0F9mih3T3zeG+cY28EYk5BX3E9AxPUV1eTHX5yuaJPZ8tTZWMzyR1zxuRmFPQR0z38FRWWvMA7Y3p+WaPD05mZXsiEk4K+ojpGppiYxb65wG2nA36iaxsT0TCSUEfMd3DUxc94sbXWreG4oRx7LSCXiTOFPQRMjo9x9h0MmtdN8VFCTavrVCLXiTmFPQR0pPFoZW+9sZKjinoRWJNQR8h/tDKbHXdALQ3VHLi9KSGWIrEmII+QvwW/cYsBn1HYwVTc/P0jWqIpUhcKegjpGt4itKiBI1VZVnbpj/EUt03IvGloI+Q7qEpWurKSSQsa9tsb1DQi8Sdgj5CerJ4sZRvQ90aSosTHNcQS5HYUtBHSDbH0PuKEkbb2gq16EViTEEfEbPJFP1jM1lv0UO6n15j6UXia9mgN7NyM/utmT1lZs+a2X/zlneY2W/M7LCZ3WNmpd7yMu/5Ee/19tz+EwpD78g0zmV3DL2vo7GSE2cmSaU0xFIkjlbSop8BXu+cewmwE7jBzK4FPg18zjm3HRgCbvPWvw0Ycs5tAz7nrScXqWs4feOxnLToGyqZTaboGZnK+rZFJHjLBr1L86chKvG+HPB64Hve8ruBm73HN3nP8V6/3syyN0ykQGXzPvQLtTdWALqLpUhcraiP3syKzOxJoB+4D+gEhp1zSW+VLqDVe9wKnATwXh8BGrJZdCHyZ5ZqqSvP+rY7/LH0GnkjEksrCnrn3LxzbiewEXgFcNliq3nfF2u9v6jz18xuN7O9ZrZ3YGBgpfUWrJ7hKZqqyygrLsr6tpuryykvSeiErEhMrWrUjXNuGHgIuBaoM7Ni76WNQI/3uAvYBOC9XgucWWRbdzrndjnndjU1NV1Y9QUkmxOOLJRIGO0NGnkjElcrGXXTZGZ13uM1wBuAA8CDwDu81W4FfuA9vtd7jvf6A053zLpo3UNTORlx4+torFTXjUhMraRF3wI8aGb7gceA+5xzPwQ+DnzMzI6Q7oO/y1v/LqDBW/4x4I7sl11YUilHz8h0zlr0kB5Lf/LMJMn5VM72ISLBKF5uBefcfuCaRZYfJd1fv3D5NPDOrFQnAJyemGU2mWJDbfZPxPraGyqYm3ecGplm09qKnO1HRPJPV8ZGQN/oNADrcxj0bd7NzXTPG5H4UdBHQO+IH/S567ppa0i34k+c1lh6kbhR0EdAr9+ir8ldi765upyy4gQn1KIXiR0FfQT0jkyTMGisKs3ZPhIJo62hQi16kRhS0EdA7+g0TdVlFBfl9se1eW2lgl4khhT0EdA3Op3T/nlfe0MFJ85MaKJwkZhR0EdA78g062uyN0/sUtoaKpieS9/3XkTiQ0EfAb2j0zk9Ees7O8RSt0IQiRUFfchNzCQZm07SnMMx9L6zQyzPqJ9eJE4U9CHnD61syUPQt9atoThhGmIpEjMK+pDr8y6Was5D101xUYLW+jUaeSMSMwr6kMvHxVKZ2ho0xFIkbhT0IXdqJPf3ucnUtraC46c1xFIkThT0Idc3Ok1NeTEVpcveaDQr2hoqGJtOMjw5l5f9iUjuKehDrndkOm+teXhhiKVG3ojEh4I+5PpGp/NyItbXfvYulhp5IxIXCvqQOzWSn4ulfP6kIzohKxIfCvoQS86nGByfycsYel95SREtteWagEQkRhT0ITYwPkPKkZerYjO1NVTwvFr0IrGhoA+xszNL5bHrBqBtbSXHFfQisaGgD7HePF4Vm6mtsYLB8RkmZpJ53a+I5IaCPsTyeZ+bTG1rNVG4SJwo6EOsd3Sa0qIEaytzN4XgYtobNfJGJE4U9CHWNzLNupoyzCyv+233Lpo6pvvSi8SCgj7E8j2G3ldZVkxTdZkmIBGJCQV9iKXnis1/0AN06C6WIrGhoA8p51zephBcTHtjBcd0MlYkFhT0ITU6lWR6LhVYi76toZKBsRnGNcRSJPIU9CHlD63M9xh6X0ejdxdLtepFIk9BH1KnRqaA/I+h9/kjb44Pqp9eJOoU9CHVF3CLvs27XbEumhKJPgV9SPWOzADBBX1lWTHrqss0ll4kBhT0IdU7OkVDZSmlxcH9iNobK9VHLxIDCvqQyvcUgotpb6jgmProRSJPQR9SvaMzgY2h97U3VjI4PsPYtCYKF4kyBX1I9Y1O533CkYU6/InCdYWsSKQp6ENoem6eMxOztATcom9r0O2KReJAQR9C/aPeiJug++i92xXr5mYi0aagDyH/qtig++grSotprinTtIIiEbds0JvZJjN70MwOmNmzZvZhb/laM7vPzA573+u95WZmXzSzI2a238xemut/RNz4V8UGPeoG0t03atGLRNtKWvRJ4M+dc5cB1wIfNLPLgTuA+51z24H7vecAbwG2e1+3A1/OetUx518VG4ag72ioVB+9SMQtG/TOuVPOuce9x2PAAaAVuAm421vtbuBm7/FNwDdc2q+BOjNryXrlMdY7MkNFaRHVZcVBl+INsZzVEEuRCFtVH72ZtQPXAL8Bmp1zpyD9ywBY563WCpzMeFuXt2zhtm43s71mtndgYGD1lcdYn3cf+nxPIbiYdv+eN7pwSiSyVhz0ZlYFfB/4iHNu9HyrLrLMvWiBc3c653Y553Y1NTWttIyCcGpkKhTdNgBb11UBcHRwPOBKRORCrSjozayEdMh/yzn3z97iPr9Lxvve7y3vAjZlvH0j0JOdcgtDXwiuivW1NVSQMOjsV9CLRNVKRt0YcBdwwDn32YyX7gVu9R7fCvwgY/l7vdE31wIjfhePLC+VcqG4KtZXVlzE5rUVdA7ohKxIVK3kbN91wHuAp83sSW/ZfwY+BXzXzG4Dngfe6b32Y+CtwBFgEvh3Wa045gYnZkimXGha9ABbm6o4oha9SGQtG/TOuT0s3u8OcP0i6zvggxdZV8Hq8+5DH5Y+eoBt66p45PAg8ylHUSL4E8Qisjq6MjZkwnJVbKatTVXMzqfoGtLIG5EoUtCHTG+ILpbybV2XvrlZ54C6b0SiSEEfMr0jUxQljMaqsqBLOWtLY3qIZWe/TsiKRJGCPmR6R2ZYV10Wqr7w+spSGipL1aIXiSgFfcj0jU4HNiH4+WxtqlLQi0SUgj5kTo1MhepErG/rukoNsRSJKAV9yPSNzoTqRKxva1MVQ5NznJmYDboUEVklBX2IjM8kGZ9JhjPovXveqPtGJHoU9CHSOxK+MfS+bU3+yBsFvUjUKOhDpGc4PbNUSwhb9Bvq1lBWnFCLXiSCFPQh4gd9a/2agCt5saKE0dFYqZubiUSQgj5EuoenSBihHF4J6X56tehFokdBHyLdw+mhlSVF4fyxbGuq4uSZSabn5oMuRURWIZyJUqC6h6bYUBe+bhvf1nVVpByaLFwkYhT0IdIzMhXK/nmfP/LmuT5134hEiYI+JOZTjlPD0yFv0VdSnDAO9Z5vymARCRsFfUgMjKVnlgpz0JcVF7G1qYqDp8aCLkVEVkFBHxLd3tDKjSEOeoAdLdUc7FXQi0SJgj4k/KAPc4seYMf6GrqHpxiZmgu6FBFZIQV9SPScDfpwjqH37WipBuCQWvUikaGgD4nuoSlqyoupLi8JupTzumx9DQAHdUJWJDIU9CHRMzxFa31F0GUsq7mmjLqKEg7ohKxIZCjoQ6J7eIrWkHfbAJgZO9ZXq0UvEiEK+pBIB324T8T6dqyv4VDvGKmUC7oUEVkBBX0IjE7PMTadDP2IG99lLdVMzs7z/JnJoEsRkRVQ0IdAmG9PvJgdOiErEikK+hDoHorGGHrfJc3VmKETsiIRoaAPgZ6IXBXrW1NaREdDpVr0IhGhoA+BruEpSosSNFaVBV3KiulWCCLRoaAPgZ7haVrqykkkLOhSVmzH+hpOnJ5kYiYZdCkisgwFfQj0DE+xoTYa3Ta+Heu9WyH0qVUvEnYK+hDoHgr3hCOLuawlPfLm2e6RgCsRkeUo6AM2N5+ibyzcE44sZmP9GhoqS3nypIJeJOwU9AHrHZnGOSJx+4NMZsbOTXU8eXIo6FJEZBkK+oB1eWPoW+vCf0Ozha7ZXEfnwITuTS8Scgr6gJ04PQFAW0P0gn7npnoA9ncNB1yJiJyPgj5gxwYnKC1KRK6PHuDqTbWYwZPPK+hFwkxBH7BjgxO0NVRQFKEx9L6a8hK2NlXx5EkFvUiYKegDdmxwgo7GyqDLuGDpE7LDOKdbFouE1bJBb2ZfN7N+M3smY9laM7vPzA573+u95WZmXzSzI2a238xemsvio24+5ThxejLyQX96YvbsSWURCZ+VtOj/AbhhwbI7gPudc9uB+73nAG8BtntftwNfzk6Z8dQzPMXsfCryQQ/whLpvREJr2aB3zj0MnFmw+Cbgbu/x3cDNGcu/4dJ+DdSZWUu2io2bY4PpETdRDvod66spL0nohKxIiF1oH32zc+4UgPd9nbe8FTiZsV6Xt+xFzOx2M9trZnsHBgYusIxoOxv0TdEN+uKiBFe11urCKZEQy/bJ2MWGjix6ls45d6dzbpdzbldTU1OWy4iGY4MTVJYW0RSh2xMvZuemOp7pGWU2mQq6FBFZxIUGfZ/fJeN97/eWdwGbMtbbCPRceHnxdmxwgo6mSsyiN7Qy085N9cwmU5qIRCSkLjTo7wVu9R7fCvwgY/l7vdE31wIjfhePvFh6aGVV0GVctJ2b0ydkHz+h7huRMFrJ8MpvA78CLjWzLjO7DfgU8EYzOwy80XsO8GPgKHAE+CrwpzmpOgZmkvN0DUV7aKWvtW4Nm9dW8MvO00GXIiKLKF5uBefcu5Z46fpF1nXABy+2qEJw8swkKQcdjdG7x81irtvWyA/395CcT1FcpOvwRMJE/yMDcmxwEiAWXTcAu7c1Mjad5GlNRCISOgr6gBwbHAegoyH6XTcAr9ragBnsOTwYdCkisoCCPiDHBidoqCyltqIk6FKyYm1lKVdsqGHPEQW9SNgo6ANydGCC9hiciM103bZGHn9+iMnZZNCliEgGBX1Ajp+O9l0rF7N7WyNz847fHlt4xwwRCZKCPgATM0n6RmdiF/Qvb19LaXGCR9V9IxIqCvoA+Pe42RKzoC8vKWJXWz17jmg8vUiYKOgD0DngjbiJ8M3MlnLdtkYOnBplcHwm6FJExKOgD8Az3SOUFifY2hSPMfSZdm9rBFD3jUiIKOgDsL9rhMtbaiiJ4RWkV7bWsraylJ8f6F9+ZRHJi/glTcilUo5nuke4emNt0KXkRFHCuOHK9dx/oI+p2fmgyxERFPR5d3RwnInZea7eWBd0KTlz49UtTM7O88BBtepFwkBBn2f7u9L3golrix7glR0NNFWX8cP9mopAJAwU9Hm2v2uENSVFsTwR6ytKGG+9cj0PHOxnfEZXyYoETUGfZ093j3Blaw1FiWjPKrWcG1+ygZlkivsP9AVdikjBU9DnUXI+xbM9I7Hun/e9bHM962vK+eF+TTAmEjQFfR4d7h9nei4V6/55XyJhvO3qFn5xaIDR6bmgyxEpaAr6PHraOxF7VWv8gx7So29m51Pc96y6b0SCpKDPo6e6hqkuL6Y9JpONLGfnpjraGir49m+fD7oUkYKmoM+jp7tHuKq1lkTMT8T6zIz3vbqdvSeGeOL5oaDLESlYCvo8mUnOc+DUKFcVQP98pnfu2kR1eTFf23Ms6FJECpaCPk8O9Y4xN+94SQGMuMlUVVbMH71yMz95+hQnz0wGXY5IQVLQ58njJ9JdF4VyIjbT+17dTsKMu395POhSRAqSgj5Pfn6gny1NlWysXxN0KXnXUruGt13dwnceO6mhliIBUNDnwfDkLL86epo3X7Ees8I4EbvQ+3dvYXwmyXc0Akck7xT0eXD/gX7mU44brlgfdCmBuWpjLa/Z3siXHuzktGafEskrBX0e/PTZXlpqywviitjz+S83Xs7ETJLP/PRQ0KWIFBQFfY5NziZ5+LmBgu628W1vrua23R3cs/ckj2tcvUjeKOhz7BeHBphJpnjTFc1BlxIKH7p+O+tryvnL//sM8ykXdDkiBUFBn2M/e7aX+ooSXtG+NuhSQqGqrJhP3ngZz/aM8s1fHQ+6HJGCoKDPodlkivsP9vOGy5opjuFE4BfqbVe18NpLmvibnxxUF45IHih9cuhXR08zNp3khisLd7TNYsyML/zhTtbXlHP7N/bRMzwVdEkisaagzxHnHF975Ci1a0q4bltj0OWETn1lKXfduouZuXnef/deJmc15aBIrijoc+ShQwM8cniQD1+/nfKSoqDLCaXtzdV88Y+u4WDvKB/81uMKe5EcUdDnwNx8ir/+0e/Y0ljJe17VFnQ5ofa6S9fxP26+koeeG+APv/Jr+kangy5JJHYU9Dnwj795ns6BCT7x1sso0UnYZb37lW189T276BwY5+YvPcrvekaDLkkkVpRCWTYyOcfnfv4cr97awBsuWxd0OZHxhsub+e4HXoVzcPOXHuUzPz2orhyRLFHQZ9HU7Dx//k9PMTI1xyffdnnBXwm7Wle21vKvH9rNjS9p4X8/1Mn1f/sL/uWJLmaTqaBLE4m0nAS9md1gZofM7IiZ3ZGLfYRN/9g0t9z5K+4/2Md/vfFyLt9QE3RJkdRUXcZn/2An3/uTV1FfUcpH73mKV3/qfj7904McH5wIujyRSDLnsnsZupkVAc8BbwS6gMeAdznnfrfUe3bt2uX27t2b1TryxTnHb4+d4WPffYozE7N84ZadvKmA71KZTfMpx8OHB/jWr5/ngYN9pBx0NFaye1sjr97awOUbathUX1Ewc/CKLGRm+5xzu5ZbrzgH+34FcMQ5d9Qr5DvATcCSQR8lM8l5ekem6R6a4tHOQX7wZA9dQ1M0VZdxzweu5eoCmyowl4oSxusuXcfrLl3HqZEpfvJ0L3uODPL9x7v45q9PAFBRWsT2dVVsrK9gQ10562vX0FBZSm1FCfUVpVSVFVFRWkxlaTFlJQlKixL6xSAFJxdB3wqczHjeBbwyB/vhu4+d5KuPHM3KtjL/rnHOpZ87SDnH3LxjJpliJjnP2PQLJwiLEsbubY189A2X8OYr11NVlovDKZCepeqPd3fwx7s7mE2meLZnhEO9YxzsHeNI/zgHTo1y/8E+pueW788vThglRQmKE0ZxkVGUMBLmf6Wv3E0kwDDMwP+1YGZnH5PxuyLz18Zqz8voV478h+u383sv2ZDTfeQimRb77L6of8jMbgduB9i8efMF7aiuooTtzVUX9N7F2IL/vQnvP3ZJUYKykgRlxQnq1pTSWr+GDXXlXNpcTUNVWdb2LytTWpzgms31XLO5/pzlzjmGJ+c4MznL8OQsw5NzTMzOMzmTZHJ2/uwv65lkivmUY24+RXLeMe8czjnmU46UA+f9gj/7C5/0Mpexn7P7PKeA1f073GrfILFUu6Yk5/vIRdB3AZsynm8Eehau5Jy7E7gT0n30F7KjN12xXv3hcpaZUV9ZSn1ladCliIRKLkbdPAZsN7MOMysFbgHuzcF+RERkBbLeonfOJc3sz4CfAUXA151zz2Z7PyIisjI5OXvonPsx8ONcbFtERFZHV8aKiMScgl5EJOYU9CIiMaegFxGJOQW9iEjMZf2mZhdUhNkAcOIC394IDGaxnGxRXaujulYvrLWprtW5mLranHNNy60UiqC/GGa2dyV3b8s31bU6qmv1wlqb6lqdfNSlrhsRkZhT0IuIxFwcgv7OoAtYgupaHdW1emGtTXWtTs7rinwfvYiInF8cWvQiInI+zp9gIcAv4AbgEHAEuGOR18uAe7zXfwO0Z7z2CW/5IeDNy20T6PC2cdjbZmm+6iJ9n/4HgQPAs8CHM9b/K6AbeNL7emuej9dx4Glv33szlq8F7vOO131AfR6P16UZx+NJYBT4SL6OF9Dg/bzGgb9b8J6XecfrCPBFXvjrOOfHa6m6gArgR8BB7/P1qYzX3gcMZByv9+f5eD3kbdPf/7rlPhN5OF7VCz5fg8Dn83i83gjs8z5H+4DXZ/PzdU4NK1kpl1+kb2XcCWwBSoGngMsXrPOnwN97j28B7vEeX+6tX0Y6wDu97S25TeC7wC3e478H/n0e62oBXprxIXsuo66/Av5jEMfLe+040LjI/j7jf3iBO4BP57OuBdvvJT1uOF/HqxLYDfwJLw6u3wKvIj2j2k+At+TxeC1aF+mgf533uBR4JKOu9y38N+T5eD0E7Fpkf4tuK191LXj/PuC1eTxe1wAbvMdXAt3Z+nwt/ApD183ZycSdc7OAP5l4ppuAu73H3wOut/TknDcB33HOzTjnjpH+7feKpbbpvef13jbwtnlzvupyzp1yzj0O4JwbI92yb13hccpZXcvsL3NbeT1eC957PdDpnFvthXUXXJdzbsI5tweYzlzZzFqAGufcr1z6f9w3eOG45Px4LVWXc27SOfeg93gWeJz0DG+rkfW6lrHUZyKvdZnZdmAd6V+Oq3ExdT3hnPNn33sWKDezsix9vs4RhqBfbDLxheF3dh3nXBIYIf3n2FLvXWp5AzDsbWOpfeWyrrPMrJ30b/TfZCz+MzPbb2ZfN7NzJ0TNfV0O+H9mts+bz9fX7Jw75W3rFOn/DPmsy3cL8O0Fy3J9vJbS6m1nsW3m43gty8zqgN8D7s9Y/Pve8fqemW1a4q25rOv/mNmTZvaXGWG+0m3l9HgB7yLd0s4cnZLP4/X7wBPOuRmy8/k6RxiCfiWTiS+1TraW56uu9JvMqoDvk+5vHvUWfxnYCuwETgF/m+e6rnPOvRR4C/BBM3vtEvtfSi6PVynwduCfMl7Px/FaymrXv9BtXNB+zKyY9C/FLzrnjnqL/5V03/DVwM95oVWYr7re7Zy7CniN9/WeVW4rZ8fLs7CbBIr/AAACbklEQVQhkbfjZWZXAJ8GPrCKba5KGIJ+JZOJn13H+xDXAmfO896llg8Cdd42ltpXLuvCzEpIh/y3nHP/7K/gnOtzzs0751LAV1m6SyUndfl/Qjrn+oF/ydh/n/enpN9l0Z/PujxvAR53zvX5C/J0vJbSxbldIpnbzMfxWs6dwGHn3Of9Bc65015rEdLH62X5rMs51+19HwP+kRd+XivdVs6Ol5m9BCh2zu3LqDcvx8vMNpL+//Ze51xnxvoX+/k6RxiCfiWTid8L3Oo9fgfwgPcn1r3ALV6/VgewnfRJjEW36b3nQW8beNv8Qb7q8v5cvQs44Jz7bOaG/B+e598Cz+Sxrkozq/bqqATelLH/zG3l9XhlvO9dLOi2ydPxWpT3J/OYmV3r/UzfywvHJR/Ha0lm9tekg+QjC5ZnHq+3kz4/lJe6zKzYzBq9xyXAjSz++TrftnJyvDzLfb5ycry87rUfAZ9wzj3qr5ylz9e5ljpLm88v4K2kR6B0An/hLfvvwNu9x+Wk/2w/QjoAtmS89y+89x3COzO91Da95Vu8bRzxtlmWr7pIn/l3wH4WDAsEvkl6ONV+74fZkse6tpAeLfAU6ZNCmcergXQ/72Hv+9o8/xwrgNNA7YJ95et4HSfd+hon3dLyR0ntIh1WncDf8cLwt3wdrxfVRbrl50iH0jnDAoH/6f1snyLd2NmRx7oqSY9o2e/V8AVeGO215Lby8XP0Xju68Hjk43gBnwQmOHeIpz/s9KI/X5lfujJWRCTmwtB1IyIiOaSgFxGJOQW9iEjMKehFRGJOQS8iEnMKehGRmFPQi4jEnIJeRCTm/j+xwD5PR/yRAwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rv=gamma(a=24,scale=1/5000)\n",
"plt.plot(np.linspace(0,0.02,100),rv.pdf(np.linspace(0,0.02,100)))"
]
},
{
"cell_type": "code",
"execution_count": 99,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"23.8"
]
},
"execution_count": 99,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.mean(airline_df['Fatal'])"
]
},
{
"cell_type": "code",
"execution_count": 110,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"14.8"
]
},
"execution_count": 110,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
".0037*4000\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}