From ef33df6ff114125c9fbf7ef2c60c98d7443802c6 Mon Sep 17 00:00:00 2001 From: jeremyteitelbaum Date: Mon, 16 Apr 2018 12:21:24 -0400 Subject: [PATCH] Finished(?) 2.11.13 --- BDA 2.11.13.ipynb | 343 ++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 284 insertions(+), 59 deletions(-) diff --git a/BDA 2.11.13.ipynb b/BDA 2.11.13.ipynb index 46867cc..b198815 100644 --- a/BDA 2.11.13.ipynb +++ b/BDA 2.11.13.ipynb @@ -37,7 +37,7 @@ }, { "cell_type": "code", - "execution_count": 77, + "execution_count": 1, "metadata": {}, "outputs": [ { @@ -167,7 +167,7 @@ "9 1066 22 0.15 1985 7107.0" ] }, - "execution_count": 77, + "execution_count": 1, "metadata": {}, "output_type": "execute_result" } @@ -176,7 +176,7 @@ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", - "from scipy.stats import poisson\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", @@ -185,16 +185,23 @@ "airline_df" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Deaths vs miles (out of order -- this is part (3))" + ] + }, { "cell_type": "code", - "execution_count": 66, + "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_62649727c9442f9bc7cbfce75826d859 NOW.\n" + "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_d23d1d0b2e2b2dabfc71364b95c0d024 NOW.\n" ] } ], @@ -208,7 +215,8 @@ "}\n", "model {\n", "\n", - " // no prior here, what should we use?\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", @@ -218,31 +226,22 @@ }, { "cell_type": "code", - "execution_count": 67, - "metadata": {}, - "outputs": [], - "source": [ - "deaths=sm_simple.sampling(data=dict({'deaths':airline_df['Deaths']}))" - ] - }, - { - "cell_type": "code", - "execution_count": 68, + "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Inference for Stan model: anon_model_6ace3ad0d872dac6795ff2d2317d4760.\n", + "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.94 0.21 8.16 676.54 686.33 691.83 697.4 708.53 1553 1.0\n", - "lp__ 3.8e4 0.02 0.69 3.8e4 3.8e4 3.8e4 3.8e4 3.8e4 1551 1.0\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 Sun Apr 15 18:46:11 2018.\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" @@ -250,26 +249,26 @@ } ], "source": [ - "print(answer)" + "print(deaths)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "From the stan output reported above, the 95% interval for the poisson rate is (676.5,708.5)." + "From the stan output reported above, the 95% interval for the poisson rate is (675,707)" ] }, { "cell_type": "code", - "execution_count": 69, + "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_5173738091a79a17032e9cb1d2d99cd3 NOW.\n" + "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_9a1f6e45cdebbb8b75bf8564a74768dd NOW.\n" ] } ], @@ -283,10 +282,20 @@ " real theta ; \n", "}\n", "model {\n", - "\n", - " // no prior here, what should we use?\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" @@ -294,40 +303,50 @@ }, { "cell_type": "code", - "execution_count": 71, + "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "/home/jet08013/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", + "/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=sm_weights.sampling(data=dict({'deaths':airline_df['Deaths'],'miles':airline_df['Miles']}))" + "deaths_wtd=sm_weights.sampling(data=dict({'deaths':airline_df['Deaths'],'miles':airline_df['Miles']}))" ] }, { "cell_type": "code", - "execution_count": 72, + "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Inference for Stan model: anon_model_5173738091a79a17032e9cb1d2d99cd3.\n", + "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 4.4e-5 1.5e-3 0.12 0.12 0.12 0.12 0.12 1144 1.01\n", - "lp__ 3.8e4 0.02 0.71 3.8e4 3.8e4 3.8e4 3.8e4 3.8e4 1640 1.0\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 06:52:19 2018.\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" @@ -336,17 +355,18 @@ ], "source": [ "\n", - "print(deaths)" + "print(deaths_wtd)\n", + "predicted=deaths_wtd.extract()['predicted']\n" ] }, { "cell_type": "code", - "execution_count": 73, + "execution_count": 15, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHb1JREFUeJzt3Xl0lPd97/H3F4lFgEEsYtGCAQNiMXYEY8dL4niXwY5NbCBOcxviuOXeNDdN6lwa46Y3TdOeJJeem6Tn3pPUJ2nr9PQmGQHGDrFDHC/N0tjxCLHZWAZvzIxYxCJWAVp+94/nJxAgI6GR9Izm+bzO0ZmZ3zzS/OY58HxmnueZz5hzDhERiZ4BYU9ARETCoQAQEYkoBYCISEQpAEREIkoBICISUQoAEZGIUgCIiESUAkBEJKIUACIiEZUf9gQuZuzYsW7y5MlhT0NEpF+prq7e75wr6my5rA6AyZMnk0gkwp6GiEi/YmbvdWU57QISEYkoBYCISEQpAEREIkoBICISUQoAEZGIyuqzgEQkOtbVpFm1oZa6hkaKCwtYUVnOooqSsKeV0xQAIhK6dTVpVq7dSmNTCwDphkZWrt0KoBDoRdoFJCKhW7Wh9szGv01jUwurNtSGNKNoUACISOjqGhovaVx6hgJAREJXXFhwSePSMxQAIhK6FZXlFAzMO2esYGAeKyrLQ5pRNOggsIiEru1Ar84C6lsKABHJCosqSrTB72PaBSQiElEKABGRiFIAiIhElAJARCSiFAAiIhGlABARiSgFgIhIRHUaAGb2z2a2z8y2tRsbbWbPmdkOfznKj5uZ/aOZ7TSzLWY2r93vLPPL7zCzZb3zdEREpKu68g7gX4G7zht7FHjeOTcdeN7fBlgATPc/y4HvQRAYwFeBDwLXAl9tCw0REQlHpwHgnPs1cPC84fuAJ/z1J4BF7cZ/5AIvA4VmNhGoBJ5zzh10zh0CnuPCUBERkT7U3WMA451zuwH85Tg/XgIk2y2X8mPvNy4iIiHp6YPA1sGYu8j4hX/AbLmZJcwsUV9f36OTExGRs7obAHv9rh385T4/ngLK2i1XCtRdZPwCzrnHnXMx51ysqKiom9MTEZHOdDcAngbazuRZBjzVbvxT/myg64DDfhfRBuBOMxvlD/7e6cdERCQkndZBm9mPgZuBsWaWIjib55tA3MweBnYBS/zizwALgZ3ACeAhAOfcQTP7OvCqX+5vnXPnH1gWEZE+ZM51uCs+K8RiMZdIJMKehohIv2Jm1c65WGfL6ZPAIiIRpQAQEYkoBYCISETpO4FFRLLIm3uPUpVIUjZ6KJ+6fnKvPpYCQEQkZEdONvGzzXXEEyk2JxvIH2C9vvEHBYCISChaWx0vv3OAqkSKZ7ft5mRTKzPGD+crd8/iYxUljBk+uNfnoAAQEelD6YZG1lSnqKpOkjzYyGWD83lgXilLY2VcVToSs46ac3qHAkBEpJedbGrhudf3Ek8k+e3O/TgHN04bw5fuKKdyzgQKBuWFMi8FgIhIL3DO8VrdEeKJJE9tquNwYxMlhQX8+a3TWTy/lLLRQ8OeogJARKQnHTp+mnWb0sQTKbbvPsKg/AHcNWcCS2Nl3HDFGAYM6LtdPJ1RAIiIZKil1fGbHfVUJVI89/peTre0MrdkJF+/bw73Xl3CyKEDw55ihxQAIiLd9O7+46yuTrG6OsWeIycZNXQg/+W6y1kSK2XWxBFhT69TCgARkUtw4nQzz2zdQzyR5A/vHGSAwUdmFPHVj87m1lnjGJwfzgHd7lAAiIh0wjnHxl0NVCWSrN+ym2Onmpk8ZigrKst5YF4pE0YOCXuK3aIAEBF5H/uOnuTJjWniiSRv1R+nYGAed181kaWxMq6ZPKpPz9nvDQoAEZF2mlpaefGNfcQTKV6s3UdLqyN2+Si+9cBU7r6qmOGDc2ezmTvPREQkAzv2HiWeSPJkTZr9x05TdNlg/vTDU1kSK+WKouFhT69XKABEJLKOnGxi/ebdxBNJNvkStttmjWPJ/DJuLi8iPy+3G/MVACISKa2tjlfeOUhVIskz55WwLaooYWwflLBlCwWAiERC3ZkSthS7Dp4ItYQtWygARCRntZWwVVWn+M2OepyDG64YwyN3zAi1hC1bKACkX1lXk2bVhlrqGhopLixgRWU5iypKwp6WZJlt6cNUJZKs8yVsxSOH8Plbp7MkS0rYsoUCQPqNdTVpVq7dSmNTCxD0qq9cuxVAISAcOn6ap3wJ2+u+hK1yzgSWxkq54Yqx5GVRCVu2UABIv7FqQ+2ZjX+bxqYWVm2oVQBEVEur47c79xNPJHnutf5TwpYtFADSb9Q1NF7SuOSu9w6cLWHbfTgoYfvkdZNYMr+M2cXZX8KWLRQA0m8UFxaQ7mBjX1xYEMJspK+dON3Ms1v3UFWd5OW3gxK2m2YU8df3zOa2flbCli0UANJvrKgsP+cYAEDBwDxWVJaHOCvpTc45apJBCdvPNudOCVu2UABIv9G2n19nAeW++qOneLImRTyRYue+YzlXwpYtFADSryyqKNEGP0c1tbTyUm098USSF9/YR3OrY/7lo/jWA3NzroQtW2iNikiodu47SlUixZqNafYfO8XY4YN5+MNTWDK/jGnjcrOELVsoAESkzx092cT6LUEJW82usyVsS2NlfGRG7pewZQsFgIj0CeeCErZ4IskzW4MStunjolnCli0UACLSq3YfPlvC9t6BoITtfl/CdnVES9iyhQJARHrcqeYWfvX6PuKJJL/2JWzXTx3DF2+fzl1zJka+hC1bKABEpMe8VneYqkSKdZvSNJzwJWy3TGPx/DImjVEJW7ZRAIi0o7bRS9dw4jRPbaojnkjyWt0RBuUN4M4541kaK+PGaSphy2YKABFPbaNd11bCVpVI8ktfwnZlyQj+9r453Ht1MYVDB4U9RekCBYCIF2bbaH9557HrwAmqqpNnStgKhw7kjz44iSWxUuYUjwx7enKJMgoAM/sL4E8AB2wFHgImAj8BRgMbgT92zp02s8HAj4D5wAHg4865dzN5fJGeFFbbaLa/82g83cKz24Jz9l9++yBmcNP0Ir5y92xun60Stv6s2wFgZiXAnwOznXONZhYHHgQWAt92zv3EzL4PPAx8z18ecs5NM7MHgW8BH8/4GYj0kLDaRrPxew6cc2xKNhBPpFi/uY6jp5qZNHoo/+POGdw/r1QNrDki011A+UCBmTUBQ4HdwK3AH/n7nwD+hiAA7vPXAVYD/8fMzDnnMpyDSI8Iq200m77noP7oKdbVpIknkuzYd4whAwewcG5Qwnbt5NEM0AHdnNLtAHDOpc3sH4BdQCPwS6AaaHDONfvFUkDbS5gSIOl/t9nMDgNjgP3dnYNITwqrbTTs7znoqIStYlIh37h/LvdcNZHLhuhbtXJVJruARhG8qp8CNABVwIIOFm17hd/RS4cLXv2b2XJgOcCkSZO6Oz2RbgmjbTSsdx4dlbB95kNTWDK/lOnjL+vVx5bskMkuoNuBd5xz9QBmtha4ASg0s3z/LqAUqPPLp4AyIGVm+cBI4OD5f9Q59zjwOEAsFtPuIcl5ffnO4+jJJn7uS9g27mogb4Bx68yghO3m8iIGqoQtUjIJgF3AdWY2lGAX0G1AAngRWExwJtAy4Cm//NP+9u/9/S9o/79IoDffebSVsFUlUjyzdTeNTS1MGzecxxbOZFFFCeMu07dqRVUmxwBeMbPVBKd6NgM1BK/cfw78xMz+zo/90P/KD4F/M7OdBK/8H8xk4iJyceeXsA0fnM+iimKWxMqoKCtUCZtg2fwiPBaLuUQiEfY0RPqN9iVsv9lRT6uD66aOZmmsjLuunMDQQfrsZxSYWbVzLtbZcvrXINKJ/vAp3dfrjhBPJM+UsE0cOYTP3TKNxfNLuXzMsLCnJ1lKASByEdn8Kd22Eraq6iTb0kEJ2x2+hO1DKmGTLlAASJf0h1fBvSHbPqXb0ur43c79xNuVsM0pHsHX7g1K2EYNUwmbdJ0CQDqVza+Ce1u2fEp314ETrPYlbHWHTzKyIChhWzy/lCtLVMIm3aMAkE5l26vgvhTmp3QbT7fwi9d2E381xe/fPoAZfHh6EY/dPYvbZ41nyECVsElmFADSqWx5FRyGvv6UrnOOzanDxBNJfrbpbAnbl+6YwQPzVcImPUsBIJ0Ku6smTH31Kd39x07x5MY0VdVJ3tyrEjbpGzkZAFE9YNlbwuqqyRa99Snd5nYlbC/4ErYPlKmETfpOzgVAlA9Y9pawWjJz1c59x6iqTrJ2Y5r6o6cYO3yQStgkFDkXAFE+YNmbwmjJzCXHTjXz8y11xBMpqt87RN4A45bycSyNlXLLzHEqYZNQ5FwARPmApWQX5xyvvnuIeCLJz7cEJWxXFA1TCZtkjZwLgCgfsJTssOfwSdZsTFGVSPKuStgki+VcAET9gKWE41RzC89vD0rYfv1mUML2wSmj+fyt01kwVyVskp1y7l+lDlhKX9q+25ew1aQ5dKKJCSOG8Gc3ByVsk8eqhE2yW84FAOiApfSuwyeaeHpzmngixdb0YZWwSb+VkwEg0tNaWx2/e2s/8USKDa/t4XRzK7MmjuBvPjqb+z5QohI26ZcUACIXkTx4gqrqFGuqU6QbGhlZMJBPXFPGkliZStik31MAiJznZFMLv9i2h3giyX++FZSwfWjaWB5dMJM7ZquETXKHAkCE4Jz9Lb6E7enNdRw92UzZ6AIe8SVsJTqNWHKQAkAi7cCxUzxZk6YqkaJ271GGDBzAgisnsiRWynVTxqiETXKaAkAip7mllf94Myhhe3772RK2v//YlXz06mJGqIRNIkIBIJHxVv0xqhIp1mxMUX/0FGOGDeKhGyezJFbGDJWwSQQpACSnHTvVzDNbdhNPJEmcKWErYkmsjFtVwiYRpwCQnNNWwlaVSPLzrbs5cbqFqUXDWLlgJh+bpxI2kTYKAMkZbSVsq6tTvLP/OMMG5XHv1UEJ27xJKmETOZ8CQPq1082tPL99L/FEkv/wJWzXThnN526ZxkKVsIlclP53SL/0xp4jxF9NsW5TmoPHTzNhxBA+e/MVLJ5fxhSVsIl0iQJA+o22Eraq6hRbUocZmGfcOXsCS2KlfHh6kUrYRC6RAkCyWmur4z/fOkA8keQXvoRt5oTL+KovYRutEjaRblMASFZKHjzB6urggG66oZERQ/J58JoylsbKmFM8Qgd0RXqAAkCyxsmmFja8FpSw/W7n2RK2Ly+YyZ0qYRPpcQoACZVK2ETCowCQUJxfwjY4fwAL56qETaQvKQCkzzS3tPLrHfXEX03xq+17aW51XK0SNpHQKACk171df+zMt2rtUwmbSNZQAEivaCthq6pO8uq755aw3VI+jkH5KmETCZsCQHqMc47Ee4eIv3puCdujC2Zyf0UJ40aohE0kmygAJGN7jwQlbFWJsyVsH72qmKXXlDJv0iidsy+SpTIKADMrBH4AXAk44DNALfBTYDLwLrDUOXfIgq3Ad4GFwAng0865jZk8voTndHMrL7yxl3gixUu1+84pYVtw5QSGDdZrC5Fsl+n/0u8Cv3DOLTazQcBQ4DHgeefcN83sUeBR4MvAAmC6//kg8D1/Kf3IG3uOUJVI8WRNUMI2fsRglbCJ9FPdDgAzGwHcBHwawDl3GjhtZvcBN/vFngBeIgiA+4AfOecc8LKZFZrZROfc7m7PXvrE4cYmnt5cR1UieaaE7Y7Z41kSK+MmlbCJ9FuZvAOYCtQD/2JmVwPVwBeA8W0bdefcbjMb55cvAZLtfj/lxxQAWai11fH7t30J27Y9nPIlbP/zntksqlAJm0guyCQA8oF5wOedc6+Y2XcJdve8n45eJroLFjJbDiwHmDRpUgbTk+5IHQpK2KoSZ0vYPq4SNpGclEkApICUc+4Vf3s1QQDsbdu1Y2YTgX3tli9r9/ulQN35f9Q59zjwOEAsFrsgIKTntZWwVSVS/O6t/YBK2ESioNsB4JzbY2ZJMyt3ztUCtwGv+59lwDf95VP+V54G/ruZ/YTg4O9h7f8Pj3OObekjxBNJntqU5sjJZkpHFfDF22bwwPwSSkcNDXuKItLLMj0L6PPAv/szgN4GHgIGAHEzexjYBSzxyz5DcAroToLTQB/K8LGlGw4cO8W6TcEB3Tf2BCVsC66cwJJYGddPVQmbSJRkFADOuU1ArIO7butgWQd8LpPHk+5pbmnlNzv2E08k+dX2vTS1BCVsf7coKGEbWaASNpEo0qd1ctg7+49TlUiyZmOKvUeCErZl1wclbOUTVMImEnUKgBxz/FQzz2zdTVUixR/ePcgAg1vKx/G1e8u4daZK2ETkLAVADnDOUf3eIeKJJOu3+BK2scP48l0zuX9eCeNVwiYiHVAA9GP7jpxkzcY0VYkkb/sStnuumsjSWBnzL1cJm4hcnAKgnwlK2PZRlUjy0pv1tLQ6rp08ms/efAUL505UCZuIdJm2Fv1E7Z6jVCWSPFmT5sDx04y7bDD/9aapLJ5fytSi4WFPT0T6IQVAFjvc2MTPfAnbZl/Cdvus8SyNlfHh6WPJz9MBXRHpPgVAlmltdbzsS9iebVfC9tf3zGbRB4oZM3xw2FMUkRyhAMgSqUMnWFOdpqo6SepQI5cNyWdJrJSlsTLmlozUAV0R6XEKgBCdbGrhl6/vpSqR5Lc79+NcUMK2orKcyjkTVMImIr1KAdDHOiphKyks4Au3TeeBeaWUjVYJm4j0DQVAL1hXk2bVhlrqGhopLixgRWU5N80oYl1NmrgvYRvkS9iWqoRNREKiAOhh62rSrFy7lcamFgDSDY08Et+EmdHS6riqdCRfX3Ql915VzMihKmETkfAoAHrYqg21Zzb+bVodDBs0gDWfvYGZE0aENDMRkXMpAHpIWwlbuqGxw/tPnGrRxl9EsooCIAPOOTbuOkT81RTrt9Rx/HQLeQOCXT3nKy4sCGGGIiLvTwHQDfuOnGStP6D7dv1xhg7K4+65E1l6TRmpgyd47Mlt5+wGKhiYx4rK8hBnLCJyIQVAFzW1nC1he7E2KGG7ZvIo/ttHruDudiVs10wejZldcBbQooqSkJ+BiMi5FACdeHPvUeKvnlvCttyXsF3xPiVsiypKtMEXkaynAOjAkZNBCVs8kWJzsoH8Ab6E7ZpSbppepBI2EckJCgCvtdXx8jsHqEqkeHbbbk42tTJj/HC+cvcsPlZRohI2Eck5kQ+AdEMja6pTVFUnSR5s5LLB+TwwLyhhu6pUJWwikrsiHQBv7DnCgu/+Bufgxmlj+NIdQQlbwSCVsIlI7ot0AJSPv4zHFszirisnqIRNRCIn0gFgZvzpTVPDnoaISCh0OouISEQpAEREIkoBICISUQoAEZGIUgCIiESUAkBEJKIUACIiEaUAEBGJKAWAiEhEKQBERCJKASAiElEKABGRiFIAiIhEVMYBYGZ5ZlZjZuv97Slm9oqZ7TCzn5rZID8+2N/e6e+fnOlji4hI9/XEO4AvANvb3f4W8G3n3HTgEPCwH38YOOScmwZ82y8nIiIhySgAzKwUuBv4gb9twK3Aar/IE8Aif/0+fxt//22m71sUEQlNpu8AvgP8JdDqb48BGpxzzf52Cijx10uAJIC//7Bf/hxmttzMEmaWqK+vz3B6IiLyfrodAGZ2D7DPOVfdfriDRV0X7js74NzjzrmYcy5WVFTU3emJiEgnMvlKyBuBe81sITAEGEHwjqDQzPL9q/xSoM4vnwLKgJSZ5QMjgYMZPL6IiGSg2+8AnHMrnXOlzrnJwIPAC865TwIvAov9YsuAp/z1p/1t/P0vOOcueAcgIiJ9ozc+B/Bl4BEz20mwj/+HfvyHwBg//gjwaC88toiIdFEmu4DOcM69BLzkr78NXNvBMieBJT3xeCIikrkeCQARCayrSbNqQy11DY0UFxaworKcRRUlnf+iSAgUACI9ZF1NmpVrt9LY1AJAuqGRlWu3AigEJCupC0ikh6zaUHtm49+msamFVRtqQ5qRyMUpAER6SF1D4yWNi4RNASDSQ4oLCy5pXCRsCgCRHrKispyCgXnnjBUMzGNFZXlIMxK5OB0EFukhbQd6dRaQ9BcKAJEetKiiRBt86Te0C0hEJKIUACIiEaUAEBGJKAWAiEhEKQBERCJKASAiElEKABGRiFIAiIhElAJARCSiFAAiIhGlABARiSgFgIhIRCkAREQiSgEgIhJRCgARkYhSAIiIRJQCQEQkohQAIiIRpQAQEYkoBYCISEQpAEREIkoBICISUQoAEZGIUgCIiESUAkBEJKLyw56AiEhfWFeTZtWGWuoaGikuLGBFZTmLKkrCnlaoFAAikvPW1aRZuXYrjU0tAKQbGlm5ditApENAu4BEJOet2lB7ZuPfprGphVUbakOaUXZQAIhIzqtraLyk8ahQAIhIzisuLLik8ajodgCYWZmZvWhm283sNTP7gh8fbWbPmdkOfznKj5uZ/aOZ7TSzLWY2r6eehIjIxayoLKdgYN45YwUD81hRWR7SjLJDJu8AmoEvOedmAdcBnzOz2cCjwPPOuenA8/42wAJguv9ZDnwvg8cWEemyRRUlfOP+uZQUFmBASWEB37h/bqQPAEMGZwE553YDu/31o2a2HSgB7gNu9os9AbwEfNmP/8g554CXzazQzCb6vyMi0qsWVZREfoN/vh45BmBmk4EK4BVgfNtG3V+O84uVAMl2v5byY+f/reVmljCzRH19fU9MT0REOpBxAJjZcGAN8EXn3JGLLdrBmLtgwLnHnXMx51ysqKgo0+mJiMj7yCgAzGwgwcb/351za/3wXjOb6O+fCOzz4ymgrN2vlwJ1mTy+iIh0XyZnARnwQ2C7c+5/t7vraWCZv74MeKrd+Kf82UDXAYe1/19EJDyZVEHcCPwxsNXMNvmxx4BvAnEzexjYBSzx9z0DLAR2AieAhzJ4bBERyZAFJ+VkJzOrB97rg4caC+zvg8fpz7SOukbrqXNaR12TyXq63DnX6UHUrA6AvmJmCedcLOx5ZDOto67Reuqc1lHX9MV6UhWEiEhEKQBERCJKARB4POwJ9ANaR12j9dQ5raOu6fX1pGMAIiIRpXcAIiIRldMBYGZ5ZlZjZuv97Slm9oqvqv6pmQ3y44P97Z3+/snt/sZKP15rZpXhPJPeYWbvmtlWM9tkZgk/dsl13ma2zC+/w8yWvd/j9Ve+uHC1mb3h68+v13o6l5mV+39HbT9HzOyLWk/nMrO/8PX528zsx2Y2JNTtknMuZ3+AR4D/B6z3t+PAg/7694HP+ut/BnzfX38Q+Km/PhvYDAwGpgBvAXlhP68eXD/vAmPPG/tfwKP++qPAt/z1hcCzBJ1O1wGv+PHRwNv+cpS/Pirs59bD6+kJ4E/89UFAodbTRddXHrAHuFzr6Zz1UgK8AxT423Hg02Ful0JfKb24sksJvo/gVmC9/4e2H8j3918PbPDXNwDX++v5fjkDVgIr2/3NM8vlws/7BEAtMNFfnwjU+uv/BHzi/OWATwD/1G78nOX6+w8wwv+nNa2nLq+zO4HfaT1dsF7aGpFH++3MeqAyzO1SLu8C+g7wl0Crvz0GaHDONfvb7euoz1RV+/sP++W7VGHdjzngl2ZWbWbL/dil1nnn+jqaCtQD/+J3J/7AzIah9XQxDwI/9te1njznXBr4B4KKnN0E25lqQtwu5WQAmNk9wD7nXHX74Q4WdZ3c16UK637sRufcPIJva/ucmd10kWWjuo7ygXnA95xzFcBxzn7LXUeiup4A8Puv7wWqOlu0g7GcXk/++Md9BLttioFhBP/3ztdn26WcDACCorp7zexd4CcEu4G+AxSaWVsBXvs66jNV1f7+kcBBcrzC2jlX5y/3AU8C13Lpdd45vY4Inl/KOfeKv72aIBC0njq2ANjonNvrb2s9nXU78I5zrt451wSsBW4gxO1STgaAc26lc67UOTeZ4O3oC865TwIvAov9YudXVbedbbDYL+/8+IP+aPwUgu8z/kMfPY1eZWbDzOyytusE+223cel13huAO81slH+Fc6cfywnOuT1A0szavj38NuB1tJ7ezyc4u/sHtJ7a2wVcZ2ZDzcw4+28pvO1S2AdG+uDAy82cPQtoql9ROwneog7240P87Z3+/qntfv+vCI6y1wILwn4+PbhephKcSbAZeA34Kz8+huDg+Q5/OdqPG/B//brYCsTa/a3P+HW3E3go7OfWC+vqA0AC2AKsIzg7RevpwvU0FDgAjGw3pvV07jr6GvAGwYutfyM4kye07ZI+CSwiElE5uQtIREQ6pwAQEYkoBYCISEQpAEREIkoBICISUQoAEZGIUgCIiESUAkBEJKL+P42HafhjjkW6AAAAAElFTkSuQmCC\n", + "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": [ "
" ] @@ -358,53 +378,258 @@ "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 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": 90, + "execution_count": 111, "metadata": {}, "outputs": [ { "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VNeZ+P/PowJCCESHoaggAaYIgwsYbGwMBmLJRd6NEzts4hJ/SWynOtkkDvkl9m7YlO8W7/eb3STekmR3SRxnv4swBptm4x5s7BhkikEgUYQEAiGBUNec3x/nDBqEykia0WhmnvfrpZdm7rRz78w8z51z7n2OGGNQSikVu+LC3QCllFLhpYlAKaVinCYCpZSKcZoIlFIqxmkiUEqpGKeJQCmlYpwmgn5IRH4hIv9fJ7cbEcnuyzb1JyKySEQ+9rteIiK3hbNNrh2Xtaud2zPce5fQl+3qoC0PisibAd731yLyw1C3qYPXvmybichLIvJAD54nTURqRCQ++K2MfJoIgkREdojIOREZ2NvnMsZ80Rjz18FoV2fC+QVvj4g85b70X2mz/Gtu+VMAxpg3jDHTwtLITrRtV39JUNHEGHO7MeY3Xd2v7bY3xhwzxqQYY1pC28LIpIkgCEQkA1gEGOCusDYm8h0E2u7xfc4tVxFMLI05/ZC+KcHxOeCPwK9pE8REZJCI/J2IHBWRahF5U0QGudtuEpG3RaRKRI6LyINu+WV76iLylyJSJiInReThNs8/UET+VkSOicgp163ke/7FInJCRL4hIqfdczzkblsFrAS+5X4yb3DLvy0ipSJyQUQ+FpGlbVdWRG4QkXL/n9kico+I7HGX54nILhE579r0993Ylu8BySIy0z3XTGCQW+57rcUicqK9B4tInIh8R0QOi8hZEXleREa425JE5L/c8ioReU9ExnbVIBH5jYh8w12e4H6dPOauZ4tIpQtyl9olIv8JpAEb3Pb9lt9TrnTv1xkRWd3J6/5aRP7ZdYfUiMhbIjJORJ5xvz4PiMhcv/tPd79Mq0Rkr4jc5XfbSBF5wb0n7wJZbV7rKhHZ6tblYxH5VFfbxT3uQdeu/+s+3wf8PzOuPWtE5C2gFpgsIqki8m/u81gqIj/0fZZEJN59ns+IyBEgr83r7RCRR/yu/y8R2e8+r/tE5Jr2tr1c2cU03m2PShEpEpH/5fecT7nPzX+4590rItf53d7ldyTiGGP0r5d/QBHwGHAt0ASM9bvtn4AdwAQgHlgIDMR+UC8A9wOJwEhgjnvMr4EfusufAE4Bs4DBwG+xvzyy3e3PAC8AI4AhwAbgR+62xUAz8FfuNXKxX8bhbV/HXZ8GHAfGu+sZQFYH63wYWOZ3/Q/Ad9zld4DPusspwA0BbsengP8Cvgv8xC37KfCkW/6U33qd8HtcCXCbu/w1bFKe6LbzL4Hfudu+4LZPsnsvrgWGBtCuh4EN7vJn3Lr/3u+29V21y297GuBfsMntaqABmN7B6/4aOOPamQS8AhRjdzzigR8Cr7r7JmI/h98FBgBLsJ+vae7254Dn3WdoFlAKvOluG+ze94eABOAa97oz2/uctGnjg9jP2NddGz4NVAMj3O07gGPATPfciUCBe18GA2OAd4EvuPt/ETgATMJ+pl912yzB7/kecZfvdetxPSBANpDexbb3Pc9rwD+77ToHqACW+n0O67Hfl3jgR8Afu/sdiaS/sDcg0v+Am7DBf5S7fgD4urscB9QBV7fzuCeBdR0856UvHvDvwI/9bpvqPtDZ7sN/0f+DCCwAit3lxe71E/xuP40LzG2/4O45TwO3AYldrPcPgX93l4e4dqS7668DT/u2STe25VPYgJ/mgkei+z+JwBPBft8X2l33uPcnARu03wZmd7NdWUCVez9/gU0oJ9xtvwGe6Kpd7nqGe+8m+i17F7ivk8/Bv/hd/zKw3+96DlDlLi8CyoE4v9t/57ZpvNsGV/nd9je0JoJPA2+0ee1fAj9o73PS5n4PAicBabNOvh2BHcBf+d02Fpv8Bvktu5/WhPYK8EW/25bTcSLYDHy1g3Z1tO0T3OepBRjid/uPgF/7fQ63+d02A6jr7nckkv60a6j3HgC2GGPOuOu/pbV7aBR2j+NwO4+b1MHytsZj90B8jvpdHo3du33fdQdUAS+75T5njTHNftdrsXvpVzDGFGH3qJ8CTovIcyIyvoN2/Rb4M7GD438GfGCM8bXt89iEdcB1v9wRwHr6t+MYdu/2b4BDxpjjXTzEXzqwzm977Md+6ccC/4kNHs+J7Wb7qYgkBtCew0ANds9xEfAicFJEpgG3YPcuu6Pc73KH74dzyu9yXTvXfY8dDxw3xnj9bj+K/SU6GhsAO/ocpQPzfdvMbbeVwLgA1gWg1Lgo6ffc/p8b/9dNxyb4Mr/X+iX2l8Gl9eignW0F+h1qazxQaYy50OZ1Jvhdb/seJYlIQje/IxFDE0EviO2L/xRwi9g+83LsT+SrReRq7M/retr0xzrHO1jeVhn2A++T5nf5DDYYzDTGDHN/qcaYzgKLvytKzxpjfmuMuQn7hTXAT9p9oDH7sF+e27HdJb/1u+2QMeZ+7Jf7J8B/i8jgANvk8x/AN9z/7jgO3O63PYYZY5KMMaXGmCZjzNPGmBnYLro7sN0sgXgN+CQwwBhT6q5/DhgOfNjBY/qytO9JYJJcPhibhu06qcB233T0OToOvNZmm6UYYx4N8LUniIi0ee6Tftf9t8Nx7C+CUX6vNdQYM9Pd3tnnva3OvkOdbfuTwAgRGdLmdUo7eUzrEwf4HYkkmgh6Jx+7tzkDu7c4B5gOvAF8zu2d/Tvw925wKl5EFri96LXAbSLyKRFJcIN5c9p5jeeBB0VkhogkAz/w3eCe/1+AfxCRMXBpMHNFgO0/BUz2XRGRaSKyxLWvHptkOjvc7rfAV4CbsWMEvuf5CxEZ7dpX5RZ397C932O7BZ7v5uN+AawRkXTXltEicre7fKuI5LiByfPY7pIWd9tTIrKjk+d9DfgSttsLbBfFl7HdKx2t22XbN8R2YrvnviUiiSKyGLgTeM6173+Ap0QkWURmcPlBDS8CU0Xks+6xiSJyvYhMD/C1xwBfcY+7F/sd2NTeHY0xZcAW4O9EZKjYwf0sEbnF3eV591wTRWQ48J1OXvdfgW+KyLViZfvedzrZ9u4X5tvAj8QeQDAb+yt2bVcr2oPvSETQRNA7DwC/MvYY5XLfH/Az7JEhCcA3gULsUS+V2L2HONf9kYvd663E7lVe3fYFjDEvYQeEX8F2l7zS5i7fdsv/KCLngW3YAa1A/Bsww/1EL8AOrv4Y+0ujHPsF/24nj/8dtl/8Fb+uMbAD3HtFpAb4R2wfeD2AO4pjUVcNM8bUGWO2GWPqAlwXn3/EDp5vEZEL2IHj+e62ccB/Y5PAfmxw/y932yTgrU6e9zXsWIgvEbyJ7ZZ7vcNH2H7n77nt+81urke3GGMasYcu3459//4ZuzNywN3lS9hupHJsn/+v/B57AZt078PuLZdjP6eBnhOzE5jiXncN8EljzNlO7v857ID2PuAc9j3xuNv+Bdt9txv4AJvA2mWM+YN7vd9iB8YLsAPM0PW2vx87bnASWIcdD9na1YrS/e9IRJDLu/aUik0i8iF2kLmzAKbaEHvI8yOuq0RFqLCf6q5Uf2CMaa9bTqmYoF1DSikV47RrSCmlYpz+IlBKqRgXEWMEo0aNMhkZGeFuhlJKRZT333//jDFmdFf3i4hEkJGRwa5du8LdDKWUiigi0tmZ2Zdo15BSSsU4TQRKKRXjNBEopVSM00SglFIxThOBUkrFOE0ESikV4zQRKKVUjNNEoJRSMU4TgVJK9QfFa6EgA34bZ/8XdzlPTtBExJnFSikV1YrXwruroKXWXq89aq8DZK4M+cvrLwKllAq33atbk4BPS61d3gc0ESilVLjVHuve8iDTRKCUUuGWnNa95UGmiUAppcLt6jUQn3z5svhku7wPaCJQSqlwy1wJ856F5HRA7P95z/bJQDHoUUNKKdU/ZK7ss8Dflv4iUEqpGKeJQCml+pEzdWf6/DU1ESilYtrawrVkPJNB3NNxZDyTwdrCvjuj16euuY4Nhzfw+c2fZ+kflnL8wvE+fX0dI1BKxay1hWtZtWEVtU32ZK6j1UdZtcGe0bsyJ7T99cYYdlfspqCogM0lm6lpqmFiykQevfpRBicODulrt6WJQCkVs1ZvX30pCfjUNtWyevvqkCWC07Wn2XB4A+sPr6e4uphBCYNYlr6M/Ox8rh17LXHS9x01mgiUUjHrWHX7Z+52tLynGlsa2XF8BwVFBbx18i28xss1Y67hoYUPsTxjeZ//AmhLE4FSKmalpaZxtPpou8uD4UDlAQqKCth4ZCNVDVWMSR7Dw7Me5u6su8lIzQjKawSDJgKlVMxas3TNZWMEAMmJyaxZ2vMzes/Vn2NT8SYKigo4UHmAxLhElqQtIT87nwWeBcTHxQej6UGliUApFbN84wCrt6/mWPUx0lLTWLN0TbfHB5q9zbx98m0Kigp49firNHubmTFyBt+d/11yM3NJHZgaiuYHjRhjwt2GLl133XVm165d4W6GUkpdpri6mIKiAjYc3kBFXQXDBw4nb3Ie+dn5TBsxLdzNQ0TeN8Zc19X99BeBUkp1Q01jDS+XvExBUQG7K3YTL/EsmrCI/Ox8bp54M4nxieFuYrdpIlBKRZy1hWt73Z3THV7j5b3y9ygoKmDb0W3Ut9STlZrFN679Bndk3cGoQaNC9tp9QROBUiqi9OVJYKU1pawvWs8Lh1+gtKaUIYlDuCvrLvKz85k1ahYiEtTXCxcdI1BKRZSMZzLaPeQzPTWdkq+V9Pr565rr2HZ0G+uL1rOzfCeCMN8zn/zsfJamLSUpIanXr9FXdIxAKRWVQnESWEflHh6f8zh3Z92NJ8XT4+eOBJoIlFIRJZgngfXHcg/hoIlAKRVRensSWHvlHuaOmcvTC59mRcaKsJd7CAdNBEqpiNLTk8AipdxDOOhgsVIqakViuYdg0sFipVRMivRyD+GgiUApFRWOVB+hoKiAFw+/eKncw33T7us35R76s5AnAhGJB3YBpcaYO0QkE3gOGAF8AHzWGNMY6nYopaKPr9zDuqJ17KnYExXlHsKhL34RfBXYDwx1138C/IMx5jkR+QXweeDnfdAOpVQUiPZyD+EQ0kQgIhOBPGAN8ITY87GXAJ9xd/kN8BSaCJRSXYiVcg/hEOpfBM8A3wKGuOsjgSpjTLO7fgKY0N4DRWQVsAogLS04swUppSKLr9xDQVEB75a/e6ncw1fmfoUlaUsiqtxDfxayRCAidwCnjTHvi8hi3+J27tru8avGmGeBZ8EePhqSRiql+h3/cg8vl7zMxaaLMVXuIRxC+YvgRuAuEckFkrBjBM8Aw0Qkwf0qmAicDGEblFIRwlfuoaCogJLzJTFb7iEcQpYIjDFPAk8CuF8E3zTGrBSRPwCfxB459ACwPlRtUEr1b+2Ve7hmzDU8POthlmcsj8lyD+EQjvMIvg08JyI/BP4E/FsY2qCUCqMDlQdYd2gdG4s3Ut1QreUewqxPEoExZgeww10+Aszri9dVSvUfsV7uoT/TM4uVimJ9PaVjW1ruITJoIlAqSvXllI5ttVfu4f6r7ic/O5+pw6eG9LVV92kiUCpKrd6++rKa/QC1TbWs3r46JIlAyz1ELk0ESkWpUEzp2JaWe4gOmgiUilLBnNKxLS33EF00ESgVpXo7pWNbvnIP64vWs7N8p5Z7iCKaCJSKUj2d0tGff7mHzSWbqWmq0XIPUUinqlRKXcFX7mH94fUUVxdruYcIpVNVKqW6pamliR0ndrDu0LpL5R7mjpnL0wufZkXGCi33EMU0ESgV4w5UHqCgqICNRzZS1VCl5R5ikCYCpfqJQM8CDsbZwlX1VWws3qjlHhSgiUCpfiHQs4B7c7Zwe+Uepo+YzpPzniRvcp6We4hhOlisVD+Q8UxGu8f8p6emU/K1km7fz19xdTEFRQVsOLzhUrmHvMl55GfnM23EtGCtguqHdLBYqQgS6FnAgd7PV+6hoKiA3RW7tdyD6pQmAqXCyNffb9qfsfWKs4A7O1vYa7zsKt/FuqJ1Wu5BdYsmAqXCpG1/f1vtnQXc3tnCQ5OGcvuM28n9n9xL5R7uzLqT/Ox8ckblaLkH1SVNBEqFSXvVQX3SU9PbPRrI/2zh6uZqxg0bR0JiAm+VvcV8z3y+PPfLLE1bquUeVLdoIlAqTDrq7xek3YFfX7mHQxcO4RnjYUjTECakTCA/O5+7su5ifMr4ELdYRStNBEqFSaDVQStqK9hwZAMFRQVa7kGFhCYCpbohmFM/dlYd1FfuoaCogLdK36LFtGi5BxUymgiUClCwp35srzroVxZ8heN1x1nyhyW23MOgMTw06yEt96BCSk8oUypAPTmZKxDtlXu4ddKt3DPlHi33EEuK18Lu1VB7DJLT4Oo1kNm7KUX1hDKlgiyYUz/6l3vYcXwHTd6mS+UecjNzGZY0rLfNVZGkeC28uwpaXDdh7VF7HXqdDAKhiUCpAAVj6sf2yj18etqntdxDrNu9mtqmFl49fyNbz9/Ajyf+H5Kotb8QNBEo1X/0dOrH9so93DThJvKz87ll4i1a7iGG1TY288qB02zafz+vnr+OOpPEyPgqjjRMZMagYttN1Ac0ESgVoO5M/dheuYfJqZN54tonuDPrTi33EMMuNrjgX1jGqx+fpr7Jy6jE2fz5iO3kpr7J/MF7iRevvXNy4L82e0MTgVLdsDJnZadHCJXWlPJC0QusP7ye0ppSUhJTtNyD4mJDM9sPnGbTnjJ2HHTBP2Ug9147idwcD/N4mfhdv2kdIwCIT7YDxn1AE4FSvVTXXMe2o9tYX7SeneU7EUTLPShqGprZvv8UmwrL2PFxBQ3NXkYPGcinrrPB//qMEcTH+XYMVkIcQT9qKFCaCJTqAV+5h4KiAjaXbKamqYYJKRN4fM7jWu4hhvmC/8Y9Zbx20Ab/MUMGct/1Nvhfd1nwbyNzZZ8F/rY0ESjVDRW1Fbxw2Hb9aLkHBXChvont+0+zsdAG/8ZmL2OHDuT+eWk2+KcPJ66j4N9PaCJQqgta7kG1daG+iW37T7FxTzmvH2oN/p+Zl0bebA/XpvX/4O9PE4FSHThQeYCCogI2Htmo5R4U5+ub2LbP9vm/fvAMjS1exg1NYuX8NPJyPFwTYcHfnyYCpfy0V+5hSdoS8rPztdxDDKquaw3+bxyywX98ahKfXZBObo6HuZOGRWzw96eJQMU8/3IPrx5/lWZv86VyD3mT80gdmBruJqo+VF3XxNZLwb+CphbDhGGD+NyCdHJne5gzMTqCvz9NBCqiPfbztTx7ZDUtg48RfzGNVZPX8M+PBnbkxZHqIxQUFfDi4RcvlXu4b9p9nZZ7WLsWVq+GY8cgLQ3WrIGV4TnQQwVRdW0TW/aVs6mwjDeLzlwK/g8uzCA3x8OcScOi+hwQTQQqYj3287X8vHQVpNiTcFpSjtrrP6fDZNBeuYdFExaRn53PzRNv7rTcw9q1sGoV1Lpzfo4etddBk0EkqqptZIvb83/LL/g/dGMmuTkerp6YGtXB35+WoVYRK+EvM2hJubIIXHxNOs3/u+TSda/x8l75exQUFVxW7uGe7Hu4I+uOgMs9ZGTY4N9WejqUlFy5XPU/5y42snXfKTa64N/sNUwcPoi8HA+5OR5mR1nw1zLUKuq1DG6/IJdveWlNKeuL1vPC4RcorSllSOKQXpV7ONZB/a+Olqv+4dzFRjbvLWdjYRnvHD5Ls9cwacQgPr8ok7wcDzkToiv490TIEoGIJAGvAwPd6/y3MeYHIpIJPAeMAD4APmuMaQxVO1T0ir+YdsUvAhFhWPxkHtn8SNDLPaSltf+LIK1v6oKpbqh0wX9TYRlvHz5Li9eQNiKZRxZNJi/Hw6wJQ2M++PvrViIQkeHAJGPMngDu3gAsMcbUiEgi8KaIvAQ8AfyDMeY5EfkF8Hng591tuFKrJq+xYwKJtQwaOIjhKcNJTR5GfHwcpTWlPD7nce7OuhtPiicor7dmzeVjBADJyXa5Cr+zNQ1s3mv7/N85YoN/+shkVt1sg//M8Rr8O9JlIhCRHcBd7r4fAhUi8pox5onOHmfs4EONu5ro/gywBPiMW/4b4Ck0EageePqB5Rxf+zAfe7czMCker9cwzpvDj1c8EZJyD74BYT1qqP84W9PAy27P/49HKmnxGjJGJvOFmyeTq8E/YIH8Ikg1xpwXkUeAX7nunUB+ESAi8cD7QDbwT8BhoMoY0+zucgKY0MFjVwGrANL0t7dyrij3MKCFG8ZcR352fp+Ue1i5UgN/uJ2paeDlj3zB/yxeA5mjBvPFW2zwn+HR4N9dgSSCBBHxAJ8CVnfnyY0xLcAcERkGrAOmt3e3Dh77LPAs2KOGuvO6KvpouYfYVnHB7fnvKWNnsQ3+k0cN5rHF2eTmeJjuGaLBvxcCSQRPA5uBN40x74nIZOBQd17EGFPluphuAIaJSIL7VTARONnNNqsYoeUeYtvpC/Vs/sge7fNucSVeA1mjB/OlW7PJne1h2lgN/sESSCIoM8bM9l0xxhwRkb/v6kEiMhpocklgEHAb8BPgVeCT2COHHgDW96jlKippuYfYdvp8PS+54P9eSSXGQPaYFL60ZAp5OR6mjk3R4B8CgSSC/wtcE8CytjzAb9w4QRzwvDHmRRHZBzwnIj8E/gT8WzfbrKJQcXUxBUUFbDi8IeByDyo6nDpfz0uFZWwqLOe9ozb4TxmTwleWTCFvtoepY4eEu4lRr8NEICILgIXAaBHxP0JoKNDlb3J3iOncdpYfAeZ1v6kq2rRX7uGmCTdxT/Y9XZZ7UJGtvLqelz4qY1NhGbuOntPgH2ad/SIYAKS4+/i/K+exXTtKdZvXeNlVvot1ResuK/fwxLVPcGfWnQGXe1CRp7y6nk2FrcEfYOrYFL661Hb7TNHgHzYdJgJjzGvAayLya2NMO+dTKhW40ppSXiiyUzyW1pSSkpjSq3IPKjKUVdexqdAe6vm+C/5XjRvCE8umkpszjuwxGvz7g0DGCAaKyLNAhv/9jTFLQtUoFR3qmuvYdnQb64vWXyr3MM8zjy/N/RJL05YyKGFQuJuoQuBkVd2lPf8PjlUBNvh/Y9lUcmd7yBqdEuYWqrYCSQR/AH4B/CvQEtrmqEhnjGHPmT2sO7SOzSWbqWmqYULKBB6b8xh3Z93N+JTx4W6iCoHSqjpeKixjY2EZf3LBf7pnKN9cPpXcHA+TNfj3a4EkgmZjjJaAUJ2qqK1gw5ENFBQVUFxdzKCEQSxLX0Z+dn5Iyj2o8DtxrpaXCu2hnh8et8F/hmcof7liGrfPGqfBP4IEkgg2iMhj2DODG3wLjTGVIWuVighXlHswLcwdM5enFz7N8vTlpAzQQBBtjlfWXur22X2iGoCZ423wz8vxkDEqtCU+VGgEkggecP//0m+ZASYHvzkq3NYWrmX19tUcqz5GWmoaa5auYWXO5cV12iv38ODMB7k7+24yUzPD1HIVKscra9nogv8eF/xnTRjKtz4xjdxZGvyjQZeJwBij3+wYsbZwLas2rKK2ydZZPlp9lFUb7FyMeVPyrij3cOukW8nPzmfh+IVa7iHKHDvbGvwLS23wz5mQyrc/cRW5OeNIH6nBP5p0OVWliCRj5xBIM8asEpEpwDRjzIt90UDQqSr7SsYzGRytvvxI4ZRBKUwYPoHkpGSavE1MHzGd/Ox8cjNzGZY0LEwtVaFw9OzFS8H/o9LzAMyemEpujofcWR7SRiaHuYWqu4I5VeWvsKWkF7rrJ7BHEvVZIlB941i1nXNxQMIAhg8ZzrDBw0hMSKS5pZlPT/u0lnuIQiVnWoP/3pM2+F89MZUnb7+K3BwPk0Zo8I8FgSSCLGPMp0XkfgBjTJ3o2T9Rp6axhqzRWbTEtzA4aTDGGC7UXaCssozhicP59rxvh7uJUW/t2r6Z9Kb4zEU2FZaxcU8Z+8ps8J8zaRirc6dze844Jg7X4B9rAkkEja56qAEQkSz8jh5SkctX7qGgqICtR7eSNDiJxqZGyivLqbpYRXNLM8mJyaz5hM7FGGpr114+DebRo/Y6BCcZHKmoscG/sJz9LvjPTRvG9/Km84lZGvxjXSBjBMuA7wEzgC3AjcCDxpgdIW+do2MEwdVeuYfbM28nPzufPWV7+N4r3+v0qCEVfBkZNvi3lZ4OJSU9e87DFTVs2mNP8jpQfgGAa9KGkZvj4fYcDxOG6Znd0S7QMYIuE4F7spHYSWUE+KMx5kzvmxg4TQS9V99cz7Zj2ygoKmBnWWu5h/zsfC330A/ExUF7X0UR8HoDf56i0zWXjvP3Bf9r04fb4D9rHOM1+MeUXg8Wi0jb+QbK3P80EUkzxnzQmwaq0POVeygoKuDl4pe13EM/lpbW/i+CQKbrLjp9gY17bGG3j0/Z4H9d+nC+f8cMbs8ZhydVg7/qXGdjBH/n/icB1wG7sb8IZgM7gZtC2zTlL5ATvXx85R7WF63nSPWRbpV76KsBS3W5NWsuHyMASE62y9tz6NSFS0f7HDxVg4gN/j+4cwa3z/IwLjWpbxquokJnZahvBRCR54BVxphCd30W8M2+aZ6Czk/08iWDppYmXjvxGuuK1l0q9zBn9JxulXsI9YCl6phv+3aWhA+eusDGPTb4Hzptg//16SN46s4Z3J7jYexQDf6qZwIZLP7QGDOnq2WhFOtjBO2d6AWQnprO5s9tvlTu4VzDOcYMGsOdWXf2qNxDKAYsVc8ZYzh4qubSnn+RL/hnjCDP9fmP0eCvOhHME8r2i8i/Av+FPYT0L4D9vWyf6gbfiV4+8XHxpA5OJSE5gU9u+ORl5R4WjF9AQlwgb2s7r3Ose8tV8Blj+PjUhUtH+xyuuIgIzMsYwefunsknZmrwV8EXSMR4CHgU+Kq7/jqgZan7UFpqGkerj5IyKIXhKcMZkjyEOImjpbmFJ+c9GbRyD70ZsFQ9Z4zhQPkFd5x/GUcqLhInMD9zJA/emMmKmWMZM0SDvwqdQIrO1QP/4P5UHyupLmHZVct4s/RNEhISaG5ppvJCJTU19fx9Jl/rAAAfO0lEQVTijp/xmemfCdprdXfAUvWcMYb9ZRcuHep55IwN/jdMHsnDN2ayYuY4Rg8ZGO5mqhjR2eGjzxtjPiUihbiziv0ZY2aHtGUxrKaxhs0ltu//w4oPEeKpPz6NSj7iAsWY6kkkvvG3kLYScoL3uoEMWKqeM8awr+y8C/7lFLvgvyBrJJ9fZIP/qBQN/qrvdThYLCIeY0yZiKS3d3tfTmgfC4PFXuPl/VPvs+7QOrYd20Zdcx2ZqZnkZ+fz1/fdScne0Vc8Rgdx+z9jDHtPnr+0519ytpb4OGHB5JHk5nhYMXMsIzX4qxDp9WCxMcZ3AlkcUOa6iHB1h8YGpZWKkzUnWX94PeuLWss95E3OIz87n9mjZiMiPLKv/cfqIG7/5Av+vqN9jrrgvzBrJF+4JYsVM8cxYvCAcDdT9VdNdZDYtycBBjp5/UK/6y1u2fUhaVEM8C/38G7ZuxgM8z3z+dLcL7Vb7kEHcfs/YwwflbYG/2OVrcH/0VuyWK7BX3Wm8ggc3AKHNkPJW/DVD2Fo3535H0giSDDGNPquGGMaRUQ/0d3UUbmHR+c82mW5Bx3E7Z+MMRSWVl8K/scr60iIExZmj+LxW7NYPmMcwzX4q/Y0N8Kxt1uD/9kiu3xgHAyth5dvgOt/BJl9M0AXSCKoEJG7jDEvAIjI3UCfFp2LZGfqzrDh8AYKigo4Un2EpPgklqUv454p93RZ7sFHB3H7D2MMu09UX+rzP3HOBv8bs0fx5VunsHzmWIYla/BX7bhQDoe22sB/eAc0XoD4gZBxE2TNhbO/g3i3t9d8HN51p/X3QTII5MziLGAt4NtlPQF8zhhTFOK2XRJpg8W+cg8FRQW8WfrmpXIP+dn5rMhYEVC5B9V/GGP48HjVpaN9SqvqSIy3wT83x8PyGRr8VTu8Xjj5ARzcbIN/2W67fOgEmLLc/k2+BQYMhoIMqG2n/zc5HfJLetyEoJ1ZbIw5DNwgIinYxHGhx62Kch9XfnxFuYcHZz7Yo3IPKryMMfzpeBWb9pTx0ketwX/RlNF8fdlUlk0fS2pyYribqfqbuio4/Aoc2mL3/mvPgMTBxHmw9PswZQWMnWnri/ur7eDIj46WB1mXiUBE/gb4qTGmyl0fDnzDGPO9UDcuElTVV7GxeCPri9azv3I/iXGJLJ60mPzsfBaOX9jjcg+q73m9LvgXlvFSYRknq+sZEB/HoimjeGLZVG6bMZbUQcEN/lrttWf6zXYzBioOuL3+LXDsj2BaYNBwyL7NBv7spZA8ovPnSU7r4BdB3xwREkjX0J+MMXPbLPvAGNN2voKQ6W9dQy3eFt4++TYFRQW8evxVmrxNTB8xnbuz7yYvMy8o5R5U37DB/xwb95Tz0kdllLngf/NU2+1z24yxDE0KzZ5/22qvYA8CePZZTQadCft2a6qD4tdt4D+4BXy1wMbmwNTlNvhPvA7i4gN/zuK1dkygxW+l4pNh3rO9GiMI2gxlIrIHuN4Y0+CuDwJ2GWNm9rh13dRfEkFJdQkFRQVsOLyB03WnGTZwGHdMvoP87HymjZgW7uapAHm9hg+OnWNjYRkvFZZTft4X/EeTN3scS6eHLvj702qvPROW7VZ1rHWvv/h1aK6HxGSYvLi1vz91Qu9eo3gt7F5tu4OS0+DqNb0eKA5mIvgWcBfwK7foIeAFY8xPe9XCbghnImhb7iFe4rlpwk3kZ+dzy8RbSIzXfuJI4PUa3j92jo17ynj5Ixf8E+K4Zepo8nI8LJ0+hiF9EPz9BWt6yljTJ9utpQmO72zd669wBZeHZ8LUFTbwp98Iif27GGAwB4t/6n4V3IadoexloN2yE9HCa7zsKt9FQVHBZeUevn7t17lz8p2MTr6y3IPqf7xew66j52yf/0dlnDrfwICEOBZPHc2Ts69iyVV9H/z96YmCPROy7XbxTOvhnUWvQEM1xCVA+kKY+xc2AYzMvnKgNwoEOpJZDniBTwHFwP8LWYvCKJByD6p/a/EadpVUuuBfzukLDQxMiGPxtNHk5nhYOn0sKQP7xwB+T08U7DcDpWEStBMsvV4o3+1O6toCpe8DBgaPgel32v7+ybdC0tBgNr9f6qz66FTgPuB+4Czwe2xX0q191LY+Uddcx/Zj2wMu96D6nxav4d1iG/xf3ltOhQv+t04bQ+5sD0uuGtNvgr+/npwoqNOJ9vIEy4YLcPhVu9d/aCvUnAIEJlwDi5+EKcvAM8f2P8WQzqqPeoE3gM/7Th4TkSPGmMl92D4g+GMEHZV7uDv77i7LPaj+ocVr2Fl81gb/j05xpqaBpEQX/HNs8B/cD4N/b+kAczcZY8s3HNpiB3uPvg3eJhiYCtlLbF9/9jJIic7u3mCMEfw59hfBqyLyMvAcdowgYlXUVrDhyAbWF63vcbkHFT7NLV7eLa5kY2EZm/eWc6amkaTEOJZcZYP/rdOiM/j70+lEA9DcACVvtgb/c8V2+eir4IZHbV//pPmgB3pc0lkZ6nXAOhEZDOQDXwfGisjPgXXGmC191MZe8ZV7WFe0jrdK37pU7uGpBU9puYcI0NziZacv+H9UztmLjQxKjG8N/leNJnlAdAd/fzrA3IHzJ1uP8DmyA5ouQkISZN4MCx63e/7Do/oYl14J5Kihi9haQ2tFZARwL/AdoNNEICKTgP8AxmEHmp81xvyje47fAxlACfApY8y5XqxDh/7pw3/i9wd+r+UeIkxzi5c/HrHBf8tev+A/fQx5OR4WT4ut4O+vuwOlUTuw7G2BE7tsX//BLXCq0C5PnQRX32f3+jMWwYDkrp8rajdS4Lr1bTLGVAK/dH9dacaWovhARIYA74vIVuBBYLsx5sci8h1sUvl295odmJrGGq4fdz352fksGL9Ayz30Y80tXt45Yvv8N+89ReXFRpIH2D1/G/zHMGhAN87UjFLdGSiNuoHl2kpbx+fgZijaBnWVIPGQdgPc9rTd6x8zvXuHd0bdRuqZLk8oC9oLiawHfub+FrtpMD3ADmNMp6fl9nSw2Bijh3z2Y00tXt457Av+5ZyrbSJ5QDxLp48lL2cci6eNISlRg39PRfzAsjFwam/rXv+Jd8F4IXmkHeCduhyylti6Pj0V8Rupc0E7szhIjckAXgdmAceMMcP8bjtnjLninRSRVcAqgLS0tGuPtvdmqYjT1OLl7cNn2bSnjM37yqmqbWKwC/65rttHg39wROSZy40XbQmHg+7wzvMn7HLP1baGz9QVMH5u9+r4dCYiN1LggnZmcRAakoI9Ae1rxpjzge6hG2OeBZ4F+4sgdC1UodbU4uWtojNsKixjy75TVNU2kTIwgaXT7YDvLVM1+IdCxAwsVxa3ntFb/Aa0NMCAFFvHZ/G37d7/UE9oXjtiNlJohTQRiEgiNgmsNcb8j1t8SkQ8fl1Dp0PZBhUejc2XB//qOhv8b3PB/2YN/iHXb6c4bWmCY++0FnE7c9AuH5kN1z9iu3zSFkDCwI6fI1gDvP12I/WtkCUCsbv+/wbsN8b8vd9NLwAPAD92/9eHqg2qbzU2e3mzqIKNe8rZuq+c8/XNDBmYwG0zbLfPoimjNPj3oX41xWnN6dbj+o/sgIbzED/AFm677mE70DsyK7DnCuYAb7/aSOETsjECEbkJe2ZyIfbwUYDvAjuB54E04BhwrzsaqUP9pQy1ulJDcwtvHjrDxsIytu47xYX6ZoYkJbBsxlhyZ3lYNHUUAxM0+MccrxfK/tQ6OfvJP9nlQzw26E9dAZm3wMAenMcT5QO8wRT2MQJjzJt0fCby0lC9rgq9huYW3jhou3227m8N/stnjCNv9jhuzNbgH5Pqq93hnVugaCtcrAAEJl4PS75nB3vH5fS+eqeeXh10emC9Ckh9UwtvHLLBf9u+U1xoaGZoUgIrZo4jL8fDjdmjGJCgJTpiijG2f//SNI3vgLcZkoa5aRqX2/+DRwb3dXWAN+g0EagO1Te18PrBChv895+mpqGZ1EGJfGLWOHJne7gxS4N/zGmqs3V8fMG/ygXkMTNh4ZfdNI3XQ3wIQ4sO8AadJgJ1mfqmFl5zwX+7X/DPzRlHbo6HhRr8Y0/VcRv0D22BI69Bcx0kDLKHd974VbvnP2xS37VHB3iDThOBor6phR0f+4L/KS42tjAsOZG8HA+5sz0szBpJYrwG/5jR0mzP4vXt9Z/eZ5cPS4drPmv3+jNuCu80jStXauAPIk0EMcoG/9NsLCznFRf8hycncufV48nN8bBAg39suXjW1u85tBmKtkN9lZ2mMW0BLPtre5TPqKlROU2j0kQQU+oaW3j149NsKizjlQOnqW1sYcTgAdw1xwb/GyZr8I8ZxkD5ntbDO0/swk7TOBquyrPdPVm3QlJquFuq+oAmgihX29jMqwcq2PRRGa/sP01dkw3+d8+ZQF6OhxsmjyBBg39saLhgT+Y6tMWWdLhQZpePnwu3fNue0euZG3PTNCpNBFGptrGZVw7YPf9XD1RQ19TCyMEDuOcaG/znZ2rwjxlnD7u+fjdNY0sjDBxq9/anrLCHdw4ZG+5WqjDTRBAlLjb4Bf+PT1Pf5GVUygD+/NoJ5OZ4mJ85kvg47d+Nes0NNuD7yjlUHrbLR02DeatsX3/aAp2mUV1GE0EEu9jQzPYDp9m0p4wdB33BfyD3XjuJ3BwP8zJHaPCPBefL/A7v3AGNNRA/EDIXwfwvwpRlMEJn5VMd00QQYWoamtm+/xSbCsvY8XEFDc1eRg8ZyKeus8H/+gwN/hEr0Iqa3hYo/cBN2LLZDvoCDJ0AOfe6Oj43w4DB9jk/faseb98JnalSE0FE8AX/jXvKeO2gDf5jhgzkvutt8L9Og3/k66qiZt05e1jnoS32MM/asyBxMGk+LP2BDf5jZlx+eKdOw9gl3URWn01V2RuxWH30Qn0T2/efZmOhDf6NLvjn5nhs8E8fTpwG/+jRXkXN0XEwbzT82Vw4vhNMCwwaYbt6prhpGpNHdO85Qat0+on2TRT26qOq+y7UN7Ft/yk27inn9UM2+I8dOpDPzEsjb7aHa9M0+EetY8fstzEzAaa4v2FxQB00XoCbvm73+idcG/g0jVqls0u6iSxNBGF2vr6Jbftsn//rB8/Q2OJl3NAkVs5PIy/HwzUa/KPbuaO2u+fh4TC2CRIFGg0cbobXG6B+HOx7s2fPrVU6u6SbyNJEEAbVda3B/41DNvh7UpP4ixvSyZs9jrmTNPhHrZYm283jq+NTccAuH9QC77fAoSY42gIt2Iqaz/7oyucIdHRTq3R2STeRpYmgj1TXNbH1UvCvoKnFMD41ic8uSCc3x8PcScM0+Eermgo7UcvBzXD4VWiohrhEGDgZXmmBvXVQ6SbxEwGD7aRuL8B3Z3RTq3R2STeRpYPFIVRd28SWfeVsKizjzaIzNLUYJgwbdKmk85xJwxAt4hV9vF4o+9CWcTi02R7qiYGUcXagd+oKW8J5Wk73RyqjfXRTBZUOFodJdW0Tm13wf8sv+D+4MIO82eO5emKqBv9oVH8ejrzaOk1jzSnsNI3Xwa2rbQLwXH354Z09GanU0U0VApoIgqCqtpEte0+x0QX/Zq9h4vBBPHxjJrk5HmZr8I8+xsCZQ3aP/9AWOPoOeJtstc6spXavP/s2GDyq4+foyUiljm6qENBE0EPnLjayZV85GwvLedsF/0kjBvH5RZnk5XjImaDBP+o01cOv/hq2/QuMq4MRrnDfmBmw4HF7bP+k+YFP09iTkUod3VQhoImgGyovNrJlbzkbC8t4+/BZWryGtBHJPLJoMnk5HmZNGKrBP9pUl7q9/q1wcBuYRphioLgF3mmAE4nwt1+CZT0YXezJSKWObqoQ0MHiLpytaWDzXnu0zztHbPBPH5lMbo6HvBwPM8dr8I8qLc1Quqv18M5TH9nlqWnwx1PwwTkoaYZmv8foQK3qp3SwuBfO1jTw8l474PvHI5W0eA0ZI5P5ws2TydXgH31qK239noOb4fB2W9dH4t00jX9l6/aPngbx8XZsoC0dqFURThOBc6amgZc/8gX/s3gNZI4azBdvscF/hkeDf9Qwxu7p+/b6T7wHxgvJo2DqJ1rr+AwadvnjdKBWRamYTgQVF9ye/54ydhbb4D951GAeW5xNbo6H6Z4hGvyjRUMNFL/mgv9WuHDSLvfMgZv/0u71j+9imkYdqFVRKuYSwekL9Wz+yA74vltcaYP/6ME8fqsN/leN0+AfNSqPtE7OXvKmnaZxwBDIWgxTvmuP7R8yLvDn04FaFaViYrD49Pl6Xt5bzsY9ZbxbUokxkDV6MHk5HnJne5g2VoN/VGhuhGNvtwb/s0V2+cgp9rj+Kcttv3/CgPC2U6k+ooPFwO/fO8b/e7+U947a4D9lTApfWTKFvNkepo4dEu7mqWC4cMpN07gZDu+wJZvjB0DGIjtH75RlMGJyuFupVL8W1YngraKznKtt1OAfTbxeOPlB6+TsZR/a5ee9UJYEi78ID33fTtOolApIVCeCn35yNkmJAU7iofqvuio4/Irb898KtWfsNI0DM+H1FthXB6e8QA288DNIma399kp1Q1QnAk0CEcoYqPjYTc6+BY6946ZpHG7r90xZAdlLYcY1cPTi5Y+trbWDuZoIlApYVCcCFUGa6qD4jdYiblXuJK2xOXDT12zwn3jd5dM0aiVOpYJCE4EKn6rjrXv9xa9Dcx0kJtta/Tc9YY/ySZ3Q8eP1BC+lgkITgeo7Lc12mkZf8K/Yb5cPz4RrH7CBP/1GSEwK7Pn0BC+lgkITgQqti2cur+NTXw1xCZC+EOb+hT2+f2T25RO2BEpP8FIqKDQRqOAyBsp2tx7eWfo+dprGsTD9TrvXP/lWSBoanNdbuVIDv1K9pIlA9V7DBTspu+/wzppyQGDCNbD4SZi6HMZd3XkdH6VU2IQsEYjIvwN3AKeNMbPcshHA74EMoAT4lDHmXKjaoELoTJHr698MR9+20zQOTIXsJe7wztsgZXS4W6mUCkAofxH8GvgZ8B9+y74DbDfG/FhEvuOufzuEbVDdtXZt+33uzQ22cNuhLfav8oi9/+ir4IZHbV//pPkQnxje9iului1kicAY87qIZLRZfDew2F3+DbADTQT9x9q1lx+FU3kMfvYFOPkLaDwMTRchIQkyb4YbHrP9/cPTw9tmpVSv9fUYwVhjTBmAMaZMRMb08eurznzvuzCiAW4YCFMSYJw7eevUR7D4IbvXn7EIBiSHt51KqaDqt4PFIrIKWAWQpicIhU5tpa3jc3Az3FsJyYPBa+BYC2yrh4PNcMbA3/59uFuqlAqRvk4Ep0TE434NeIDTHd3RGPMs8CzY+Qj6qoH9Vkd9991lDJza23pS14l33TSNI+FkEnxYBYebod7vMena/aNUNOvrRPAC8ADwY/d/fR+/fmRq23d/9Ki9DoElg8aLtoSDb5rG8yfs8nGzYdE37FE+E66B3z0H61ZBfXPrY/VMXaWiXshmKBOR32EHhkcBp4AfAAXA80AacAy41xhT2dVz9XaGsoiXkdF+TZ30dCgpaf8xlcWtJ3WVvAktDTAgxdbxmboCspfBUM+VjwvWLw+lVNgFOkNZTExVGfHi4myXTlsidqIWsNM0Hv+j2+vfAmcO2uUjslqnaUxfCAkD+67dSqmw0qkqo0lHVTanTYQ//Zer4/Nq6zSN6TfCdQ/b4D8yq+/bq5SKKJoIIoF/lc3xcTAlEaYNAE81rH8chnhg1j22r3/yYhiYEu4WK6UiiCaC/q6+GuYkw1/dBKd3wiADBhiUCQtX2uA/Lqdn1TuVUgpNBP2PMbZ/39fXf+wd8DZD0jC4/s9td0/2bTB4ZLhbqpSKEpoI+oOmeih5ozX4V7nxgDEzYeGX3TSN10O8vl1KqeDTyBIuVcdbC7gdec1O05gwyPbx3/hVu+c/bFK4W6mUigGaCPpKS7M9i9d3UtfpvXb5sHS45rN2rz/jpsCnaVRKqSDRRBBKF8/aaRoPbYai7VBfZadpTFsAy39o9/pHTdWBXqVUWGkiCCZjoHyPO6N3C5x4DzAweDRclWcDf9atkJQa7pYqpdQlmgh6q6EGjuywe/2HtsKFMrt8/DWw+DswZRl45uo0jUqpfksTQU+cPdx6hM/Rt6ClEQYOtXv7U1bY4J+iUy0opSKDJoJANDfagO8r4lZ52C4fNQ3mrbK1fNIW6DSNSqmIpImgI+fL/A7v3AGNNRA/EDIXwfwv2r3+EZnhbqVSSvWaJgIfbwuUfuAmbNlsB30Bhk6AnHvtXn/mzTBgcHjbqZRSQRbbiaDunD2s89AWe5hn7VmQOJg0H5b+wAb/MTP08E6lVFSLrURgDJze3zpN4/GdYFpg0Ahbv2fqCshaAskjwt1SpZTqM9GfCBpr7TSNvsM7q4/b5eNy4Kav22P7J14HcfHhbadSSoVJdCeCDV+D3b+D5npIHGzr+Nz8TRv8h44Pd+uUUqpfiO5EMCwNrn3QBv6Mm3SaRqWUakd0J4JFT4S7BUop1e9p3QOllIpxmgiUUirGaSJQSqkYp4lAKaVinCYCpZSKcZoIlFIqxmkiUEqpGKeJQCmlYpwYY8Ldhi6JSAVwtIcPHwWcCWJz+rtYW1+IvXXW9Y1uwVzfdGPM6K7uFBGJoDdEZJcx5rpwt6OvxNr6Quyts65vdAvH+mrXkFJKxThNBEopFeNiIRE8G+4G9LFYW1+IvXXW9Y1ufb6+UT9GoJRSqnOx8ItAKaVUJzQRKKVUjIu4RCAiSSLyrojsFpG9IvK0W54pIjtF5JCI/F5EBrjlA931Ind7ht9zPemWfywiK8KzRoERkXgR+ZOIvOiuR/v6lohIoYh8KCK73LIRIrLVrfNWERnulouI/B+3bntE5Bq/53nA3f+QiDwQrvXpiogME5H/FpEDIrJfRBZE6/qKyDT3vvr+zovI16J1fX1E5OsuZn0kIr9zsax/fI+NMRH1BwiQ4i4nAjuBG4Dngfvc8l8Aj7rLjwG/cJfvA37vLs8AdgMDgUzgMBAf7vXrZL2fAH4LvOiuR/v6lgCj2iz7KfAdd/k7wE/c5VzgJffZuAHY6ZaPAI64/8Pd5eHhXrcO1vc3wCPu8gBgWDSvr996xwPlQHo0ry8wASgGBrnrzwMP9pfvcdg3UC83bjLwATAfeyZeglu+ANjsLm8GFrjLCe5+AjwJPOn3XJfu19/+gInAdmAJ8KJrf9Sur2tfCVcmgo8Bj7vsAT52l38J3N/2fsD9wC/9ll92v/7yBwx1QUJiYX3brONy4K1oX19sIjiOTVoJ7nu8or98jyOuawgudZN8CJwGtmKzYpUxptnd5QR2w0PrG4C7vRoY6b+8ncf0N88A3wK87vpIont9AQywRUTeF5FVbtlYY0wZgPs/xi3vaN0iZZ0nAxXAr1z337+KyGCid3393Qf8zl2O2vU1xpQCfwscA8qw38v36Sff44hMBMaYFmPMHOye8jxgent3c/+lg9s6Wt6viMgdwGljzPv+i9u5a1Ssr58bjTHXALcDj4vIzZ3cN9LXOQG4Bvi5MWYucBHbNdKRSF9fAFx/+F3AH7q6azvLImp93XjH3djunPHAYOxnu62wfI8jMhH4GGOqgB3YfsNhIpLgbpoInHSXTwCTANztqUCl//J2HtOf3AjcJSIlwHPY7qFniN71BcAYc9L9Pw2swyb8UyLiAXD/T7u7d7RukbLOJ4ATxpid7vp/YxNDtK6vz+3AB8aYU+56NK/vbUCxMabCGNME/A+wkH7yPY64RCAio0VkmLs8CLuB9wOvAp90d3sAWO8uv+Cu425/xdjOtReA+9zofCYwBXi3b9YicMaYJ40xE40xGdif0a8YY1YSpesLICKDRWSI7zK2H/kjLl+3tuv8OXd0yQ1Ateta2AwsF5Hhbo9suVvWrxhjyoHjIjLNLVoK7CNK19fP/bR2C0F0r+8x4AYRSRYRofU97h/f43APovRg0GU28CdgDzY4fN8tn+w2SBH2p+ZAtzzJXS9yt0/2e67V2PGFj4Hbw71uAaz7YlqPGora9XXrttv97QVWu+UjsYPmh9z/EW65AP/k1q0QuM7vuR5226IIeCjc69bJOs8BdrnPdQH2KJhoXt9k4CyQ6rcsatfXtfVp4ICLW/+JPfKnX3yPtcSEUkrFuIjrGlJKKRVcmgiUUirGaSJQSqkYp4lAKaVinCYCpZSKcZoIlPIjIkZE/tPveoKIVEhr1de7ROQ77vJTIvLNcLVVqWBJ6PouSsWUi8AsERlkjKkDlgGlvhuNMS9gT+pRKmroLwKlrvQSkOcuX3b2q4g8KCI/a/sAEckSkZddkbw3ROQqt/xeV39+t4i83ietV6qbNBEodaXnsKfxJ2HPZN/Zxf3BTjj+ZWPMtcA3gX92y78PrDDGXI0tsKZUv6NdQ0q1YYzZ42aEuh/Y1NX9RSQFW0DsD7aMDGDLBwC8BfxaRJ7HFhpTqt/RRKBU+17A1o9fjK2B05k4bF35OW1vMMZ8UUTmY7uaPhSROcaYs8FurFK9oV1DSrXv34G/MsYUdnVHY8x5oFhE7oVLc+xe7S5nGWN2GmO+j51lalInT6VUWGgiUKodxpgTxph/7MZDVgKfFxFfxdS73fL/LSKFIvIR8Dq2oqpS/YpWH1VKqRinvwiUUirGaSJQSqkYp4lAKaVinCYCpZSKcZoIlFIqxmkiUEqpGKeJQCmlYtz/D4R8lqFz/uPUAAAAAElFTkSuQmCC\n", "text/plain": [ - "{'copy_X': True, 'fit_intercept': True, 'n_jobs': 1, 'normalize': False}" + "
" ] }, - "execution_count": 90, "metadata": {}, - "output_type": "execute_result" + "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": [ - "lr=LinearRegression()\n", - "lr.fit(airline_df['Miles'].values.reshape(-1,1),airline_df['Deaths'].values.reshape(-1,1))\n", - "deaths_pred=lr.predict(np.linspace(4000,8000,10).reshape(-1,1))\n", - "ax.plot(np.linspace(4000,8000,10),deaths_pred)\n", + "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", "plt.show()\n", - "lr.get_params()" + "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": 84, + "execution_count": 100, "metadata": {}, "outputs": [ { - "ename": "AttributeError", - "evalue": "'LinearRegression' object has no attribute 'slope'", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mlr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mslope\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[0;31mAttributeError\u001b[0m: 'LinearRegression' object has no attribute 'slope'" - ] + "data": { + "text/plain": [ + "[]" + ] + }, + "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": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], - "source": [] + "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", @@ -430,7 +655,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.6.4" } }, "nbformat": 4,