5-Bayesian_Regression_with_PyMC3.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"5-Bayesian_Regression_with_PyMC3.ipynb","provenance":[],"authorship_tag":"ABX9TyO9GKG17WRQ3x1U4qWKuBW7"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"code","metadata":{"id":"1_8LroUeiXwB","executionInfo":{"status":"ok","timestamp":1616837694219,"user_tz":-270,"elapsed":4250,"user":{"displayName":"Erfan Salavati","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14Ggmz2YbEhjRBiZQITcg6T34VeDuIPFZ9DYQw8PaKg=s64","userId":"05241680325615283112"}}},"source":["import numpy as np\n","import matplotlib.pyplot as plt\n","import pymc3 as pm"],"execution_count":1,"outputs":[]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":282},"id":"60is82qEXeBX","executionInfo":{"status":"ok","timestamp":1616837694222,"user_tz":-270,"elapsed":2486,"user":{"displayName":"Erfan Salavati","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14Ggmz2YbEhjRBiZQITcg6T34VeDuIPFZ9DYQw8PaKg=s64","userId":"05241680325615283112"}},"outputId":"2dde429c-5d2e-4128-c436-6c2054f8e61f"},"source":["N = 10\n","x = np.random.uniform(-5,5,size=N)\n","epsilon = np.random.normal(0,.3,size=N)\n","y = -1 + .5 * x + epsilon\n","plt.plot(x,y,linestyle='None',marker='+',markersize=10)"],"execution_count":2,"outputs":[{"output_type":"execute_result","data":{"text/plain":["[<matplotlib.lines.Line2D at 0x7f9deee1f290>]"]},"metadata":{"tags":[]},"execution_count":2},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAMzElEQVR4nO3db4gc9R3H8c+nxmKxFh8kJWpMV+hFCGmalkUsfdDWM21sxKBFUPoHUUgLFa5BkEsDShEhEOhVqiBSxQcVpdCGSE7R0wp5oq0bG9NY9SKiNdKSE7G2CBXrtw92D+7Ou9xdZmZnvnvvFwTcvcvMl2jezs3+ZsYRIQBAXp+qewAAQDGEHACSI+QAkBwhB4DkCDkAJLeqjp2uXr06Wq1WHbsGgLQOHz78TkSsmft+LSFvtVrqdDp17BoA0rL95nzvc2oFAJIj5ACQHCEHgOQIOQD0wdjEZGXbJuQA0Ad3PX28sm0TcgBIjpADQHKEHACSq+WCIAAYZGMTk/OeE2+Njs96PTI8pF1bNxTen+t4sES73Q6u7ASwkrRGx/XG3u2FtmH7cES0577PqRUASI6QA0ByhBwAKlTlhUDTSgm57Qdsn7R9rIztAcCgmP7Qc2R4qLJ9lHVE/qCkbSVtCwAGThmrUxZSSsgj4pCkd8vYFgBgefp2jtz2Ttsd252pqal+7RYABl5p68httyQdjIhNi30v68gBDKKFLgSa63QvBFpoHTlXdgJASXZt3fCJQJdxIdBiWH4IAMmVtfzwYUnPSrrY9gnbN5WxXQAYNFWsKy9r1cr1EXFeRJwZEesi4v4ytgsAg6aKB0xwagUAKlTlhUDTCDkAVKjKC4GmEXIASI7lhwBQkX49YIIHSwBAHxVZV86DJQBgQBFyAEiOkANAcoQcAPqoinXlhBwA+qiKdeWEHACSI+QAkBwhB4DkCDkAJEfIASA5Qg4AyRFyAEiOkANATxWPYesHQg4APVU8hq0fCDkAJEfIASA5Qg4AyfGoNwArUr8ew9YPPOoNAHqKPIatH3jUGwAMKEIOAMkRcgBIjpADQE8Vj2HrB0IOAD1NX52yEEIOAMkRcgBpZb3JVdkIOYC0st7kqmyEHACSI+QAkBwhB4DkuGkWgBQG6SZXZeOmWQDSavpNrspW6U2zbG+z/art12yPlrFNAMDSFA657TMk3SPpCkkbJV1ve2PR7QIAlqaMI/JLJL0WEa9HxIeSHpG0o4TtAgCWoIyQXyDprRmvT/Tem8X2Ttsd252pqakSdgtgpct6k6uy9W35YUTcFxHtiGivWbOmX7sFMMBW2uqUhZQR8rclXTjj9breewCAPigj5M9LGrJ9ke1PS7pO0qMlbBcAsASFLwiKiI9s3yzpCUlnSHogIl4qPBkAYElKubIzIh6T9FgZ2wIALA/3WgGA5Ag5ACRHyAEgOUIOAMkRcgBIjpADQHKEHACSI+QAkBwhB4DkCDkAJEfIASA5Qg4AyRFyAEiOkANAcoQcAJIj5ACQHCEHgOQIOQAkR8gBIDlCDgDJEXIASI6QA0ByhBwAkiPkAJAcIQeA5Ag5ACRHyAEgOUIOAMkRcgBIjpADQHKEHACSI+QAkBwhB4DkCDkAJEfIASA5Qg4AyRUKue1rbb9k+2Pb7bKGAgAsXdEj8mOSrpF0qIRZgIEwNjFZ9whYYQqFPCJejohXyxoGGAR3PX287hGwwnCOHACSW7XYN9h+StLaeb60JyIOLHVHtndK2ilJ69evX/KAAIBTWzTkEXF5GTuKiPsk3SdJ7XY7ytgmAGAJIQewsLGJyXnPibdGx2e9Hhke0q6tG/o1FlYYR5z+wbHtqyX9WtIaSe9JOhIR31ns97Xb7eh0Oqe9X6DJWqPjemPv9rrHwACyfTgiPrHUu9AReUTsl7S/yDYAAMWwagUAkiPkAJAcIQdKNjI8VPcIWGEIOVAyVqeg3wg5ACRHyAEgOUIOAMkRcgBIjpADQHKEHACSI+QAkBwhB4DkCDkAJEfIASA5Qg4AyRFyAEiOkANAcoQcAJIj5ACQHCEHgOQIOQAkR8gBIDlCDgDJEXIASI6QA0ByhBwAkiPkAJAcIQeA5Ag5ACRHyAEgOUKewNjEZN0jAGgwQp7AXU8fr3sEAA1GyAEgOUIOAMkRcgBIblXdA2C2sYnJec+Jt0bHZ70eGR7Srq0b+jUWgAZzRPR9p+12OzqdTt/3m1VrdFxv7N1eybbHJib5HwKQhO3DEdGe+36hUyu299l+xfZR2/ttn1tke+g/VsQA+RU9Rz4haVNEbJY0KWl38ZEAAMtRKOQR8WREfNR7+ZykdcVHAgAsR5mrVm6U9PhCX7S903bHdmdqaqrE3Q6+keGhukcA0GCLfthp+ylJa+f50p6IOND7nj2S2pKuiSV8esqHnfVYaEXMXKyIAZppoQ87C69asX2DpB9LGo6ID5byewh5c1S5IgZAuRYKeaF15La3SbpV0jeWGnEAQLmKniO/W9I5kiZsH7F9bwkzAQCWodAReUR8saxBAACnh3utrHCsiAHyI+QrHKtTgPwIOQAkR8gBIDlCDgDJEXIASI6Ql4wn3gPoN0JeMu7vDaDfCDkAJEfIASA5Qg4AyRW618pKxxPvATRB4fuRn45Bvh859/cGUJWF7kfOqRUASI6QA0ByhBwAkiPkJeP+3gD6jZCXjNUpAPqNkANAcoQcAJIj5ACQHCEHgOQIOQAkR8gBIDlCDgDJEXIASI6QA0ByhBwAkiPkAJAcIQeA5Ag5ACRHyAEgOUIOAMkRcgBIjpADQHKEHACSI+QAkFyhkNu+w/ZR20dsP2n7/LIGAwAsTdEj8n0RsTkitkg6KOm2EmYCACxDoZBHxPszXp4tKYqNAwBYrlVFN2D7Tkk/kvQvSd86xfftlLRTktavX190twCAHkec+iDa9lOS1s7zpT0RcWDG9+2WdFZE3L7YTtvtdnQ6neXOCgArmu3DEdGe+/6ip1Yi4vKI2DTPrwNzvvUhSd8ra+CFjE1MVr0LAEil6KqVoRkvd0h6pdg4i7vr6eNV7wIAUil6jnyv7YslfSzpTUk/KT4SAGA5CoU8Iio/lQIAODWu7ASA5AovP6zS2MTkvOfEW6Pjs16PDA9p19YN/RoLABpl0eWHVSiy/LA1Oq439m4veSIAaL7TXn4IAGg2Qg4AyRFyAEguXchHhocW/yYAWEHShZzVKQAwW7qQAwBmI+QAkBwhB4DkarkgyPaUujfZOh2rJb1T4jhlYa6la+JMEnMtRxNnkgZ/ri9ExJq5b9YS8iJsd+a7sqluzLV0TZxJYq7laOJM0sqdi1MrAJAcIQeA5DKG/L66B1gAcy1dE2eSmGs5mjiTtELnSneOHAAwW8YjcgDADIQcAJJLHXLbt9gO26vrnkWSbN9h+6jtI7aftH1+3TNJku19tl/pzbbf9rkNmOla2y/Z/th27cvFbG+z/art12yP1j2PJNl+wPZJ28fqnmWa7QttP2P7b71/fyN1zyRJts+y/WfbL/bm+kXdM02zfYbtv9g+WNU+0obc9oWSvi3p73XPMsO+iNgcEVskHZR0W90D9UxI2hQRmyVNStpd8zySdEzSNZIO1T2I7TMk3SPpCkkbJV1ve2O9U0mSHpS0re4h5vhI0i0RsVHSpZJ+2pA/q/9Kuiwivixpi6Rtti+teaZpI5JernIHaUMuaUzSrZIa82ltRLw/4+XZashsEfFkRHzUe/mcpHV1ziNJEfFyRLxa9xw9l0h6LSJej4gPJT0iaUfNMykiDkl6t+45ZoqIf0TEC71//re6gbqg3qmk6PpP7+WZvV+1//2zvU7Sdkm/qXI/KUNue4ektyPixbpnmcv2nbbfkvR9NeeIfKYbJT1e9xANc4Gkt2a8PqEGxKnpbLckfUXSn+qdpKt3CuOIpJOSJiKiCXP9St0Dzo+r3MmqKjdehO2nJK2d50t7JP1c3dMqfXequSLiQETskbTH9m5JN0u6vQlz9b5nj7o/Gj/UlJmQk+3PSvq9pJ/N+Um0NhHxP0lbep8B7be9KSJq+3zB9pWSTkbEYdvfrHJfjQ15RFw+3/u2vyTpIkkv2pa6pwlesH1JRPyzrrnm8ZCkx9SnkC82l+0bJF0paTj6dPHAMv6s6va2pAtnvF7Xew/zsH2muhF/KCL+UPc8c0XEe7afUffzhTo/KP66pKtsf1fSWZI+Z/u3EfGDsneU7tRKRPw1Ij4fEa2IaKn7Y/BX+xHxxdie+Ry6HZJeqWuWmWxvU/fHu6si4oO652mg5yUN2b7I9qclXSfp0ZpnaiR3j57ul/RyRPyy7nmm2V4zvRrL9mckbVXNf/8iYndErOt16jpJf6wi4lLCkDfcXtvHbB9V99RPI5ZmSbpb0jmSJnpLI++teyDbV9s+IelrksZtP1HXLL0Pgm+W9IS6H979LiJeqmueabYflvSspIttn7B9U90zqXuU+UNJl/X+WzrSO+Ks23mSnun93Xte3XPklS33axou0QeA5DgiB4DkCDkAJEfIASA5Qg4AyRFyAEiOkANAcoQcAJL7P5b0jLfh0veEAAAAAElFTkSuQmCC\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"sdo-uYcFi15K"},"source":["# Probabilistic Model"]},{"cell_type":"code","metadata":{"id":"HOYHKRXc0Vp3","executionInfo":{"status":"ok","timestamp":1616837697804,"user_tz":-270,"elapsed":1056,"user":{"displayName":"Erfan Salavati","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14Ggmz2YbEhjRBiZQITcg6T34VeDuIPFZ9DYQw8PaKg=s64","userId":"05241680325615283112"}}},"source":["xs = np.linspace(-5,5,100)"],"execution_count":3,"outputs":[]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"sM236uTxi4m8","executionInfo":{"status":"ok","timestamp":1616837762773,"user_tz":-270,"elapsed":65539,"user":{"displayName":"Erfan Salavati","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14Ggmz2YbEhjRBiZQITcg6T34VeDuIPFZ9DYQw8PaKg=s64","userId":"05241680325615283112"}},"outputId":"d3b02093-f6b8-4d8f-ee9b-3587c1093c05"},"source":["with pm.Model() as linear_model:\n","    theta0 = pm.Normal(\"theta0\",mu=0,sigma=10)\n","    theta1 = pm.Normal(\"theta1\",mu=0,sigma=10)\n","    sigma = pm.Exponential(\"sigma\",.01)\n","    y_observed = pm.Normal(\"y_observed\",mu=theta0+theta1*x,sigma=sigma,observed=y)\n","    prior = pm.sample_prior_predictive()\n","    posterior = pm.sample()"],"execution_count":4,"outputs":[{"output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Sequential sampling (2 chains in 1 job)\n","NUTS: [sigma, theta1, theta0]\n","100%|██████████| 1000/1000 [00:01<00:00, 849.05it/s]\n","100%|██████████| 1000/1000 [00:01<00:00, 921.76it/s]\n","The acceptance probability does not match the target. It is 0.8890023634298677, but should be close to 0.8. Try to increase the number of tuning steps.\n"],"name":"stderr"}]},{"cell_type":"markdown","metadata":{"id":"FG8H5cvSt8iO"},"source":["# Plotting posterior distributions"]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":351},"id":"iF9nTPedmlBt","executionInfo":{"status":"ok","timestamp":1616832546834,"user_tz":-270,"elapsed":930,"user":{"displayName":"Erfan Salavati","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14Ggmz2YbEhjRBiZQITcg6T34VeDuIPFZ9DYQw8PaKg=s64","userId":"05241680325615283112"}},"outputId":"37a972e6-5c52-46b0-e4f8-d5b8f189562d"},"source":["#plt.hist(prior['theta0'])\n","#plt.hist(prior['theta1'])\n","plt.hist(prior['sigma'],bins=np.linspace(0,700,50))\n","plt.hist(posterior['theta0'])\n","plt.hist(posterior['theta1'])\n","plt.hist(posterior['sigma'])\n","plt.scatter(posterior['theta0'],posterior['theta1'])\n","plt.scatter(posterior['theta0'],posterior['sigma'])\n","plt.scatter(posterior['theta1'],posterior['sigma'])"],"execution_count":56,"outputs":[{"output_type":"execute_result","data":{"text/plain":["(array([  5.,  11.,  50., 125., 202., 248., 215.,  90.,  43.,  11.]),\n"," array([-1.07143548, -1.05221101, -1.03298654, -1.01376206, -0.99453759,\n","        -0.97531312, -0.95608864, -0.93686417, -0.9176397 , -0.89841522,\n","        -0.87919075]),\n"," <a list of 10 Patch objects>)"]},"metadata":{"tags":[]},"execution_count":56},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAQlklEQVR4nO3df6xfdX3H8efLFjFTjGBrrQUturoNE63mDlG3iZJMLMmKmWPwB6DR1GyQaGI2i2ZqlpBVM3Uxm2hVAmaKMn/RCf6oiFE3QQtBBDqkahl0BYpuqDNjKb73x/dcdyy33Hu/P+734uf5SL75nvM5n3POu+d+Lq97fny/pKqQJLXpUdMuQJI0PYaAJDXMEJCkhhkCktQwQ0CSGrZy2gUArFq1qtavXz/tMiTpEeX666+/r6pWj7KNZREC69evZ9euXdMuQ5IeUZLcMeo2vBwkSQ0zBCSpYYaAJDXMEJCkhs0bAkmOS3JNkluT3JLk9V3725PsS3Jj99rUW+eCJHuS3JbkZZP8B0iShreQp4MOAm+sqhuSHAVcn2Rnt+w9VfW3/c5JTgDOBJ4FPAX4cpJnVtWD4yxckjS6ec8Eqmp/Vd3QTf8U2A2se5hVNgMfr6oHquqHwB7gxHEUK0kar0XdE0iyHngucF3XdH6Sm5JcnOTorm0dcGdvtbuYIzSSbEmyK8muAwcOLLpwSdLoFhwCSR4HfAp4Q1X9BLgIeAawEdgPvGsxO66q7VU1U1Uzq1eP9IE3SdKQFvSJ4SRHMAiAj1bVpwGq6p7e8g8Cn+tm9wHH9VY/tmuTHpHWb71yKvvdu+20qexXbVnI00EBPgzsrqp399rX9rq9Ari5m94BnJnkyCTHAxuAb42vZEnSuCzkTOBFwNnAd5Pc2LW9GTgryUaggL3A6wCq6pYklwO3Mniy6DyfDJKk5WneEKiqbwCZY9FVD7POhcCFI9QlSVoCfmJYkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIatnHYBkua2fuuVU9nv3m2nTWW/mg7PBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWHzhkCS45Jck+TWJLckeX3XfkySnUlu796P7tqT5L1J9iS5KcnzJv2PkCQNZyFfG3EQeGNV3ZDkKOD6JDuBVwFXV9W2JFuBrcCbgJcDG7rX84GLundpaNP6CgXp1928ZwJVtb+qbuimfwrsBtYBm4FLu26XAqd305uBj9TAtcATkqwde+WSpJEt6p5AkvXAc4HrgDVVtb9bdDewppteB9zZW+2uru3QbW1JsivJrgMHDiyybEnSOCw4BJI8DvgU8Iaq+kl/WVUVUIvZcVVtr6qZqppZvXr1YlaVJI3JgkIgyREMAuCjVfXprvme2cs83fu9Xfs+4Lje6sd2bZKkZWYhTwcF+DCwu6re3Vu0Azi3mz4XuKLXfk73lNBJwP29y0aSpGVkIU8HvQg4G/hukhu7tjcD24DLk7wGuAM4o1t2FbAJ2AP8HHj1WCuWJI3NvCFQVd8AcpjFp8zRv4DzRqxLkrQE/MSwJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWHzhkCSi5Pcm+TmXtvbk+xLcmP32tRbdkGSPUluS/KySRUuSRrdQs4ELgFOnaP9PVW1sXtdBZDkBOBM4FndOu9LsmJcxUqSxmveEKiqrwE/XuD2NgMfr6oHquqHwB7gxBHqkyRN0Cj3BM5PclN3uejorm0dcGevz11dmyRpGRo2BC4CngFsBPYD71rsBpJsSbIrya4DBw4MWYYkaRRDhUBV3VNVD1bVL4AP8v+XfPYBx/W6Htu1zbWN7VU1U1Uzq1evHqYMSdKIhgqBJGt7s68AZp8c2gGcmeTIJMcDG4BvjVaiJGlSVs7XIcllwMnAqiR3AW8DTk6yEShgL/A6gKq6JcnlwK3AQeC8qnpwMqVLkkY1bwhU1VlzNH/4YfpfCFw4SlGSpKXhJ4YlqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlq2Lz/PwGpb/3WK6ddgqQx8kxAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGjZvCCS5OMm9SW7utR2TZGeS27v3o7v2JHlvkj1JbkryvEkWL0kazULOBC4BTj2kbStwdVVtAK7u5gFeDmzoXluAi8ZTpiRpEuYNgar6GvDjQ5o3A5d205cCp/faP1ID1wJPSLJ2XMVKksZr2HsCa6pqfzd9N7Cmm14H3Nnrd1fX9hBJtiTZlWTXgQMHhixDkjSKkW8MV1UBNcR626tqpqpmVq9ePWoZkqQhDBsC98xe5une7+3a9wHH9fod27VJkpahYUNgB3BuN30ucEWv/ZzuKaGTgPt7l40kScvMyvk6JLkMOBlYleQu4G3ANuDyJK8B7gDO6LpfBWwC9gA/B149gZolSWMybwhU1VmHWXTKHH0LOG/UoiRJS8NPDEtSwwwBSWqYISBJDZv3noCktqzfeuXU9r1322lT23erPBOQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWrYylFWTrIX+CnwIHCwqmaSHAN8AlgP7AXOqKr/HK1MSdIkjONM4CVVtbGqZrr5rcDVVbUBuLqblyQtQ5O4HLQZuLSbvhQ4fQL7kCSNwUiXg4ACvpSkgA9U1XZgTVXt75bfDayZa8UkW4AtAE996lNHLKMt67deOe0SJP2aGDUEfq+q9iV5ErAzyb/1F1ZVdQHxEF1gbAeYmZmZs48kabJGuhxUVfu693uBzwAnAvckWQvQvd87apGSpMkYOgSSPDbJUbPTwB8CNwM7gHO7bucCV4xapCRpMka5HLQG+EyS2e18rKq+kOTbwOVJXgPcAZwxepmSpEkYOgSq6gfAc+Zo/xFwyihFSZKWxqg3hiVpbKb15NvebadNZb/LgV8bIUkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWErp13AI9n6rVdOuwRJYzDN3+W9206b2r7BMwFJapohIEkNMwQkqWGGgCQ17BF/Y9ibs5I0vImdCSQ5NcltSfYk2Tqp/UiShjeREEiyAvgH4OXACcBZSU6YxL4kScOb1JnAicCeqvpBVf0v8HFg84T2JUka0qTuCawD7uzN3wU8v98hyRZgSzf7syS3LWC7q4D7xlLh+C3n2sD6RrGca4PlXd9yrg2WQX15x2EXLaS2p426/6ndGK6q7cD2xayTZFdVzUyopJEs59rA+kaxnGuD5V3fcq4Nlnd9S1XbpC4H7QOO680f27VJkpaRSYXAt4ENSY5P8mjgTGDHhPYlSRrSRC4HVdXBJOcDXwRWABdX1S1j2PSiLh8tseVcG1jfKJZzbbC861vOtcHyrm9JaktVLcV+JEnLkF8bIUkNMwQkqWFTD4Ekf5LkliS/SHLYx6EO9zUUSb6e5Mbu9R9JPtu1n5zk/t6yt06pvkuS/LBXx8auPUne2/W/KcnzplTfR7v2m5NcnOSIrn3k4zeG2o5Pcl3X/onuIQOSHNnN7+mWr19sbd12jkmyM8nt3fvRh+n3ju743JzkT3vtExt7Y6ht0uNu1PomNu7GVN/Ext4iantn9/uzu/uZJclRvWNzY5L7kvxd1/9VSQ70lr12QQVV1VRfwO8AvwV8FZg5TJ8VwPeBpwOPBr4DnDBHv08B53TTJwOfm3Z9wCXAK+dYZxPweSDAScB1U6pvU1dDgMuAPxvX8RtDbZcDZ3bT7+/V9ufA+7vpM4FPDFnfO4Gt3fRW4B1z9DkN2MngIYrHMnjy7fGTHnuj1rYE427U+iY27sZU38TG3gJreyHwL93vxwrgm8DJc/S7HviDbvpVwN8vtp6pnwlU1e6qmu/TwvN+DUWSxwMvBT67HOubw2bgIzVwLfCEJGuXur6quqqroYBvMfhMx1iMUluSMPh5frLrdylweje9uZunW35K13+x+tvpb7/vBOBrVXWwqv4buAk4td9hQmNvLLUdZrsjj7tR65vkuBu1viUYewuprYDHMPjD6EjgCOCefockzwSeBHx9kfv/FVMPgQWa62so1h3S53Tg6qr6Sa/tBUm+k+TzSZ41xfou7E6935PkyAWus5T10Z2Onw18ode8FMfvcLU9Efivqjo4R82/XKdbfn/Xf7HWVNX+bvpuYM0cfb7D4D8Mv5FkFfASfvWDkDCZsTeO2iY57sZy7CY47kapb9Jjb97aquqbwDXA/u71xarafUi32TOR/iOef9z9zD+Z5NBxOqcl+dqIJF8GnjzHordU1RVj2s1ZwId68zcAT6uqnyXZxOCvtA1TqO8CBj/oRzN47vdNwF8vZgNLdPzex+Cvotm/KhZ0/JaotqE9XH39maqqJA95XrqqvpTkd4F/BQ4wOC1/8JBuQ429Cdc20XE3hvpmDTXulrC+oYxaW5LfZHA5dfYMaWeS3+8dJxiEwNm9+X8GLquqB5K8jsFZxkvnLXax148m9eLhrxu/gEESzs5fAFzQm18F/Ah4zMNsfy+wahr19dpPprveCXwAOKu37DZg7ZSO39sY/LI9ahLHb5jaGFwrvg9YeWg/Bh9CfEE3vbLrlyHq+uUxB9YCty1gnY8BmyY99sZR2yTH3ZiO3cTG3Sj1TXrsLaQ24C+Av+rNvxX4y978c4DvPcw+VgD3L6SeR8rloPm+huKVDAb5/8w2JHny7LW6JCcyuPT1o6Wub/Z6a1fL6cDN3To7gHO6O/4nMfiB7X/opide32uBlzH4D8MvZldYwuM3Z201GMnXMPjZApwLzJ5Z7Ojm6ZZ/peu/WP3t9Lf/S0lWJHliN/1s4NnAl3pdJjX2RqptCcbdqPVNetwNXd8SjL15awP+HXhxkpXdJbMXA/3LQWcxuKHe//f07+380SH9D2+xCTvuF/AKBtfcHmBw42M2cZ8CXNXrtwn4HoMnSd5yyDa+Cpx6SNv5wC0MrvtdC7xwGvUBXwG+y+CX8B+Bx3XtYfA/3vl+t3zOv5SXoL6DXduN3eut4zp+Y6jt6QxuGu4B/gk4smt/TDe/p1v+9CGP3ROBq4HbgS8Dx3TtM8CHevu6tXtdC2xcirE3am1LMO5GrW9i425M9U1s7C2wthUMztp2d/W9+5Bt/AD47UPa/qZ37K45dPnhXn5thCQ17JFyOUiSNAGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWrY/wHkTBBZ7YP3KwAAAABJRU5ErkJggg==\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"VrkRM7dVt_jr"},"source":["# Confidence Interval"]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"OZ2j9qmyuA_8","executionInfo":{"status":"ok","timestamp":1616838356715,"user_tz":-270,"elapsed":2216,"user":{"displayName":"Erfan Salavati","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14Ggmz2YbEhjRBiZQITcg6T34VeDuIPFZ9DYQw8PaKg=s64","userId":"05241680325615283112"}},"outputId":"e301b80f-5e55-4ea1-d928-66a60ad17ead"},"source":["pm.hpd(posterior,.05)[0]['sigma']\n","pm.hpd(posterior,.05)[0]['theta0']\n","pm.hpd(posterior,.05)[0]['theta1']\n","sigma_mean = posterior['sigma'].mean()\n","theta0_mean = posterior['theta0'].mean()\n","theta1_mean = posterior['theta1'].mean()\n","map_estimate = pm.find_MAP(model=linear_model)\n","sigma_map = map_estimate['sigma']\n","theta0_map = map_estimate['theta0']\n","theta1_map = map_estimate['theta1']"],"execution_count":20,"outputs":[{"output_type":"execute_result","data":{"text/plain":["array(0.45120676)"]},"metadata":{"tags":[]},"execution_count":20}]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":282},"id":"tw5HJ_k4Ee7z","executionInfo":{"status":"ok","timestamp":1616838829745,"user_tz":-270,"elapsed":2911,"user":{"displayName":"Erfan Salavati","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14Ggmz2YbEhjRBiZQITcg6T34VeDuIPFZ9DYQw8PaKg=s64","userId":"05241680325615283112"}},"outputId":"1e4c5296-e0e4-4408-e7f2-d06d3fcd1c24"},"source":["xs = np.linspace(-5,5,100)\n","g = posterior['theta0']+posterior['theta1']*xs[:,np.newaxis]\n","g_map = theta0_map + theta1_map*xs\n","plt.plot(xs,g,color='gray',alpha=.5)\n","plt.plot(x,y,linestyle='None',marker='+',markersize=10)\n","plt.plot(xs,g_map)"],"execution_count":31,"outputs":[{"output_type":"execute_result","data":{"text/plain":["[<matplotlib.lines.Line2D at 0x7f9dd92c1590>]"]},"metadata":{"tags":[]},"execution_count":31},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeWyceZ7f9/dT98EqVhVvVvG+KeqgRFF9SLPrndnd2fUaRmIHSAIE8BpI2wgCJICBOM4GSP5MECBAYCeIB7bXiGFkvffuzM7u7PRke0bqntZJXRQp3mQVi3WwinXf9fzyh/p5RmxdPEpqUfq9gAFGbOnhUz09Hz79fb6/71cRQiBJkiSdXIZv+gYkSZKk45FBLkmSdMLJIJckSTrhZJBLkiSdcDLIJUmSTjjTN/FNW1tbRX9//zfxrSVJkk6s27dv7woh2r7+9W8kyPv7+7l169Y38a0lSZJOLEVRNp/3dVlakSRJOuFkkEuSJJ1wMsglSZJOOBnkkiRJJ9yxg1xRlB5FUf5GUZRHiqLMK4ry3zTixiRJkqSDaUTXSg34J0KIO4qiuIDbiqL8WAjxqAHXliRJkl7h2E/kQogdIcSdr/57FlgA/Me9riRJknQwDa2RK4rSD0wD1xt5XUmSpJOuWq2yvLxMrVZr+LUbdiBIUZQm4I+A/1YIkXnOX/8E+ASgt7e3Ud9WkiTprSaEIBaLsbKyQjqdxuPx0Nb2zOHMY2lIkCuKYuZJiP97IcQfP+/3CCG+B3wPYGZmRm6zkCTpnVcqlVhaWiIcDlMoFGhqauJ1LPM5dpAriqIA/xpYEEL878e/JUmSpJNNCMH29jbLy8skk0nMZjOdnZ0MDQ01/GkcGvNE/jHwXwAPFEW5+9XX/gchxA8bcG1JkqQTJZ/Ps7CwQDAYpFqt0tLSwtDQED09PUSjUYQQPHn+bZxjB7kQ4hrQ2LuSJEk6YVRVZXNzk4WFBVKpFG63m9HRUYaGhqjX63z66aesr6/zrW99i8nJyYZ+729k+qEkSdK7JJ1Oc+/ePba3tzGZTAwNDTE+Po7L5WJ+fp4bN26Qy+WwWq24XK6Gf38Z5JIkSUdUq9VYWlri0aNHlEolOjo6OH36NB0dHSQSCf7wD/+QSCQCgMfjwWazYbVaG34fMsglSZKOIBaLcfPmTXZ3d3G73Vy6dImBgQEUReFnP/sZDx48oFqt4nA48Hq9JJNJdnZ2GB4epr29vaH3IoNckiTpEMrlMnfu3GFlZQWDwcDk5CRTU1PY7XZWV1f57LPPSKfTWCwWuru7yefzbGxsUCwWAVhbW+PChQsNvScZ5JIkSQcghGB9fZ1bt25RLBbp7Ozk4sWL+Hw+MpkMP/zhDwkGgwgh8Pl8mEwmQqEQhUIBVVX161QqlYbfmwxySZKkV0in0/z85z8nEongcDi4fPkyg4OD1Ot1Pv/8c+bm5iiXyzidTlwuF/F4nHw+/9zj+I1uPQQZ5JIkSS9Ur9eZm5tjYWEBIQQTExOcP38ek8nE6uoqP/vZz0ilUphMJjo6OkilUmxsbLx0nko4HG74fcoglyRJ+hohBFtbW9y4cYN8Pk97ezsffvghHo+HRCLB1atX2dzcRFVV3G439XpdPwD0Oo7gv4oMckmSpKdkMhmuX7/O9vY2drudjz/+mOHhYXK5HNeuXePhw4cUi0UsFgsWi4VEIkGlUtlXB3+ZfD7f8HuWQS5JksSTl5Dz8/M8evSIer3OyMgIFy9eRFVV5ufnmZubY3d3FwC73U4ul2Nvb+/AAf4yxz22L4NckqT3mna0/u7du6TTaVpbW5mdncXn87G5ucnDhw8JBoNUKhVMJhPVapVEIkG9Xj/2927U3BUZ5JIkvZeEECSTSe7evUs4HMZisXDx4kVGR0eJRCJcvXqV9fV10uk0QgiEEGSzWVRVPfZT+NcD/LhhLoNckqT3jjahcG1tjWq1Sm9vL9PT05RKJW7fvs3W1haxWIxSqYQQgmq1Sr1eb0gZ5WmNakWUQS5J0nujWq2yvr7O0tISqVSK5uZmZmdncbvdLC8vE4vF2N7eJp1OU6vVUFVVD/BGdqO8dWNsJUmS3nZCCMLhMEtLS0SjUcxmMxMTE/T19bGzs8Pa2hqxWIxIJEKpVNKDuxFllDdBBrkkSe+0ZDLJ8vIykUiEcrlMW1sbQ0NDlEolHj58SD6fJxQKkUwm9dAWQjTkZeabIoNckqR3UqFQYHV1lXA4TCaTwel0MjAwgM1mIxgM8gcLeT5oSuoHebRyx+vYcv+6ySCXJOmdUqvV2NzcZGtri3Q6jaIodHV14fV6SafTJBIJSqUSf76m0t60hhACg8FwYsoozyODXJKkd4IQgp2dHdbX18lkMlSrVZqammhubqZWqxGPxzGZTPpLTZh+40/hZlGhqlgafl0Z5JIknXh7e3usrKyQyWQol8sYjUaam5sxmUzk83msViuxWIyVlRVqtRpGoxF4cwHuEAU+5DYXuce/E3+v4deXQS5J0olVLBZZXV0lHo/rbYJmsxmTyUS9XsdoNFKpVLhx4walUglFUfTTmW+CS2T5iFtc4AEmaswzRhFbw7+PDHJJkk6cWq3G1tYWwWCQer2Ooih6p4nJZMJoNGKxWLh16xZ7e3vcqXQxV+165jq/W5zZ9+tzpjDT5uOPmfWIFJe5yTnmAbjPBJ9zkYTiO/a1n0cGuSRJJ4YQgmg0ytraGuVyGYvFQr1eJ5vNYjAYaGpqwm63s7i4SCgUQlVVFEXhgjXCOdP2vmv9bnGG37bfauj9tYoEl7nBaRZRMXCH03zODGmluaHf5+tkkEuSdCKk02lWVlbIZrNYrVYsFgt7e3uUy2XcbjdNTU2Ew2FWVlaoVCooioKiKM/0g89Vuxvy1P20ThHjMteZZJkqJq5zni+4QE5pauj3eZGGBLmiKP8G+C0gJoSYasQ1JUmSAEqlkn7y0mKx0NzcTDQaJZVK0dTUpG/m0Q73aAOpXnSg526tcUEeEGGucJ1R1ilh5SqX+JLzFBV7Q65/UI16Iv+3wL8A/p8GXU+SpPdcvV7X6+AAra2tJJNJlpaWUBSF9vZ2SqUSDx480OeCa1MKX2s/uBAMEOQK1xkgSB47P+EyNzlLWbG+8o87nc6G31JDglwI8TNFUfobcS1Jkt5vQghisZheB29tbaVer7O4uEipVMLr9WIwGNjY2NC382jhfZjBVudMh3wqF4IR1rnCdXrYIYuTH/FL3OYMVcV84MvYbCe4a0VRlE+ATwB6e3vf1LeVJOkEyWQyej+4y+Wivb2dxcVFdnd3cTgceL1ednd32dvbo1Qq6RMKjzKZ8KDlFUWoTLDMFW7QSZw93PyA73CXSerK4SN0ZGTk0H/mVd5YkAshvgd8D2BmZubNbyeVJOmtVS6XWVtbIxqNYrFYGBgYYGtri/n5eRRFobm5mVwux+7uLvl8nnK5TL1ef2WAz1W7uVvrfubrB2k7NIg6p1nkMjdoZY9dvPwpv84DxlEV45E/q7YurpFk14okSd+Yer1OKBRic3MTgEAgQKlU4osvvqBYLGK326nX68TjcYrFIoVCgXK5fODrT5ufDein2w6f18FiFDXOMc/H3MRLhgit/AG/xQLDCMVw5M9qNBpxuVxcvHjxyNd4ERnkkiS9cUII4vE4a2trlEolWltbcbvd3L17l0QiAYDJZKJQKFCr1cjlchQKhYa/xHy6g8UsqlzgPh9xCxd5QnTyl/wKywzAMRZBmM1mnE4nFouFarX6Wl7ENqr98P8FfhloVRQlBPxPQoh/3YhrS5L0bslms6ysrJBOp3E6nYyOjrKyssLt27epVCoYDAYURaFWq1Eul0kmk691NrhVlJnlLpe4g5Mi6/TwJ/wG6/Q0JMC1kQDVahW73Y7ZfPAXowfVqK6V/6wR15Ek6d1VqVRYX19nZ2cHs9nMwMAA6XSan/70p2QyGYQQmEwmTKYnsRQOh6lUKq/tfuyiyD8x/T7/JX+JjTLL9PMzPiCkPFtTPyhFUbBYLNjtT/rItVVxXq+Xrq4uBgYG8Pv9jfoIOllakSTptVJVVa+Dq6pKV1cXRqORu3fvEo1GqVQqmM1mHA4Hdrtd35n5uvwt0wK/Jn7KBe5hMtZZZJirzBJROo58TYPBgMViwWw2oygKqqpisVjw+Xz4/X4GBgbo6urSf0g1mgxySZJeCyEEiUSC1dVVisUiPp8Pl8vF+vo6Gxsb5HI5DAYDbrcbr9dLNBrVw/512K7a+BbX+e+NP8UoVP5M/Yj/q/Z3WRX7n5APOjhLURR9OJc2DkBRFOx2O21tbfT29tLX10dbW1vDly1/nQxySZIaLpfLsbq6yt7eHna7nd7eXmKxGA8fPmR3dxdVVXG73XR2dpLNZrl///5rmw3eIpJc5gZnTAsIFO5yis+5yJ7Jw2rVf+jBWQaDAZPJpNfy4cmL2aamJrq6uujr66Onpwe32/3aA1wjg1ySpIapVCpsbGwQDocxmUx0d3dTKBS4d+8e29vbVCoVHA4H/f391Ot1Hj16RD6ffy330i7iXOE6p1iihokbnOMLZsgqriNdTxuPq1EUBZvNhsfjwe/309fXh9/v1+vjb5IMckmSjk1VVcLhMBsbG9Trddra2lBVlZWVFba2tigUCpjNZkZGRnC73dy/f59kMvla7sUvdrjCdcZYo4yFa8zyJecpKI5DX0tRFKxWK0ajkWq1ihACo9GIw+Ggra1ND/DXWf8+CBnkkiQdi1YHLxQKNDc3Y7FYCIfDBINBfflxT08Pvb293Lt3j7m5uSMdqX8pIegjxBWuM8QWBWz8DR9yg2lKyuFnmxgMBmw2GwaDgUqlQq1Ww2w209TURGdnp17/bm1tfWPlk5eRQS5J0pHk83lWV1dJJpNYrVZaW1vZ3d0lGAySTCZRVRWfz8fExATLy8v89V//deNfZArBMBtc4Tq9hMnh4K/5Frc4c6Alx08PztLWwNlsNoQQVCoVfXXc090nWv37bSKDXJKkQ6lWq3od3GAw4PV6yWQy+otMbXv9+Pg4u7u7fPrpp43fkSmEPsiqixgpXPyQX2GOU9QOMYlw2hxGURS9/bFWq1EqlRBCYDabaW9vp7e3l4GBAQKBwGuZXNgIMsglSToQIQThcJj19XVqtRput5tSqcTjx49JJBKUy2VsNhvDw8OoqsoXX3xBqVRq6D0oQmWKx1zmBu0kSODhz/g17jNx6EFW2stKh8NBpVLRl1KYzWba2toYHh5mYGCAzs7Ob7T+fRBv991JkvRWSCaTrK6uks/nsdvtmEwm1tbW9HngRqORnp4e7HY7Dx48oFAoNPT7G0WNszxiRszRpSSI0sIf8ps8YvTQg6wMBgN2ux273U65XNYPH1ksFjo7O5mYmHir6t8HIYNckqQXKhQKrK6ukkgk9PpxKBQikUhQq9UwGAy0tbVhs9lYW1sjl8s19EWmSVQ5zwM+4hbN5LgnBvmp8jGPGTr0HBSt28Rut1MsFkmlUsCTAPf7/Zw5c4be3t63rv59EDLIJUl6Rq1WY3Nzk1AoBIDVaiUWixGNRqnX65hMJjweDwaDgZ2dHbLZbEMD3CLKXOQeH3IbJ0U28fN9fo3/sfL3+G377UNdy2Qy4XA4cDgcFAoFfS2c3W6np6eH6elpenp63tr690HIIJckSSeEYGdnh/X1dX0GSiqVYmdnh1qths1mw+VyUa1W2dnZafhoWZso8gFzzDKHnTIr9HGVS2wpga9+x8GfwrXpgw6Hg3w+TyKR0AN8YGCAmZkZ/H7/vkM+J5UMckmSAEilUqysrJDNZjEYDBQKBSKRCOVyWV+zlsvlCAaDlMvlhga4U+T5kNtc5B4WqiwwxFUusaN0Huo62vTBpwNcGwlgt9sZGRlhZmaGzs7OE1P/PggZ5JL0nisWi6ytrRGPx1FVlUKhQCwWo1gs4nQ66e7uJpFIsLGxQalUauhscLfI8DG3mOYBRlTmGeUql4grrYda03besoPFYsHlcuF0OslkMsRiMYQQ2O12xsbGmJ2dpbW1tWH3/jaRQS5J76larcbW1hahUIhKpUKlUiEWi5HP57HZbAQCAfb29lhaWqJSqTQ0wL1ij8vc5CyPALjPBNeYJal49d/zqjVt8Isj9G53G01NTezt7REOP/kzDoeDiYkJLl26hMfjadi9v41kkEvSe0YIQTQaZW1tTd+Bubu7Sy6Xw2Kx0NXVRT6fZ2lpSV9y3ChtYlcfZKVi4DZn+JwZMsrhOkW0HvDm5macTid7e3tsb28jhKCpqYmJiQlmZ2dpbm5u2L2/zWSQS9J7JJ1O62vWCoUCyWSSbDaL0WikpaWFWq3G+vo6pVKpoWNlu0SEK1xnglXKmPk5F/g5F8grzkNfy+Fw4PP5sNlsJJNJ9vb2AHC5XIyNjTEzM/POP4F/nQxySXoPlEol1tbWiEajZLNZUqkU2WwWeBKABoOBSCRCoVBo6HH6HrHNt/iSYTYpYuX31F9lyzBMUTncqFeDwYDD4YAidHV1kUwm9SXNbreb0dFRpqen8fl879RLzIOSQS5J77B6vU4wGGRra4tUKkUqlSKXy1Gv13E4HJjNZjKZDLlcjkql0phecCEYZJMrXKefbfLY+fSrivi/rHx84EUOiqJgMBhoamqira0NIQQfFJPs7ESAJz+AhoeHOXv27BvZwvM2k0EuSe8gIQSxWIy1tTWSySTJZJJMJkOtVsNqteJwOCiVSqRSqcb1ggvBGKtc4Tp+omRo4q/4ZW5z+lCDrLQVai6Xi7a2Nmq1mj4KYNZhwOFwMTAwwNTUFF1dXRgMhzui/y6SQS5J75hMJsPKygrRaJTd3V0ymQzlchmz2Yzb7aZare57Mj8uRaicYonLXKeDBEma+T7f4R6T1JWDR4w2Rtbj8eDz+SgWi8Tjcer1ur7bs6enh8nJSQKBwDtxkKdRZJBL0juiXC6ztrZGMBgkGo2STqcpl8uYTCbcbjdCCHK5HNlstiF1cIOoc4YFLnODFlLE8fHHfJeHjB96kJXZbKalpYXm5mb9EI8W4C6Xi66uLsbHx+nr63vrJxF+Exryd0RRlO8C/wdgBP6VEOJ/acR1JUl6tXq9TigUYm1tjVAoRCqVolQq6QFuNBopFArk8/mGjJU1iSrTzPMRN/GQZYd2fp/fYoGRfYOsDnqg51u+HN8N1PVj9PCL+Sitra2MjY0xMDCAxfLqRRHvq2MHuaIoRuD/BH4VCAE3FUX5cyHEo+NeW5KkFxNCsLu7y+PHj9nc3GR3d5dKpaKXISwWC7lcjkKh0JAAN4sKM9znQ27jIs8W3fwF32aFgedOInzRgZ5/6LiN1Wqlo6NDv8d0+snTt9VqxWQy0dzczMjICMPDwyd6mNWb0ogn8llgRQixBqAoyu8BfxeQQS5Jr0k2m+Xx48esrKzo81C0F4RWq/WrcEw3ZCaKTZSYZY5LzOGgxCq9/BG/ySaBQ42S1bpKAoEAZrNZP4xkNpsxm80oioLL5WJwcJDR0VGczsP3mL+vGhHkfiD41K9DwKUGXFeSpK+pVCqsrq5y//59wuGwHuDa0uNsNks2m9UXBh+HQxT0QVZWKjxmkKtcYlvpOtR1FEXRD/Gw9GQ0gKqqOJ1ODAYDtVoNu91Of38/Y2NjJ3Ie+Dftjb01UBTlE+ATgN7e3jf1bSXpnaCqKpubm9y8eZPt7e1nAlx7iVmr1ahUKsf6Xi6R5SNucYEHmKgxzyjXuERUaTvUdbRDPF6vF0VR9B8s2qlLbbNQf38/ExMTeL3el11OeolGBPk20PPUrwNffW0fIcT3gO8BzMzMNG4CvSS9w7S5KJ9//jlbW1uUSiW9hGKxWMjn8+RyOWq1mv6ke1QekeIyNznHPAqCe0zyORdJKL5DXUdbpaYFuKqq2Gw2fD4f3y2r1Go1FEXRWwlP0kq1t1UjgvwmMKIoygBPAvw/Bf7zBlxXkt5ryWSSq1evsra2RqlU0k85Go1GKpUKxWKRWq1GtVo9VoC3igSXucFpFlExcIfTfM4MaeVwA6e0ANc2BwkhcDgceDwehBCkUimueFU6Op7sxezq6pIB3iDHDnIhRE1RlP8a+BFP2g//jRBi/th3Jknvqb29Pa5du8bKyooe4Ha7HaPRqJdOtAA/zoGeThHjMteZZJkqJq4zzRfMkFOaDnUdg8GAzWbD6/Xqh3RcLhetra2oqko8HqdardLa2sqpU6fw+/3yNGaDNaRGLoT4IfDDRlxLkt5HWivhzZs3WV5eplgs6rO2jUYjqqrqAa7956gCIswVrjPKOiUsXOUSX3L+SIOsbDYbHo8Hi8WCoii43W46Ojqo1+tEIhGKxSI+n4/JyUl6enrkYZ7XRP5dlaRvUL1eJxaLcffuXVZXV8nn8wghMJlM+lOrqqrU63V9NviRBlsJQT9BrnCdQYIUsPETPuYm5ygr1kNdSpsF7vV69b5vl8uF3++nVqsRDAYpFAo0Nzdz5swZBgYGMJsPPmtFOjwZ5JL0DahUKmxvb7O4uMjm5ibZbJZ6va7PGzEajXqnh/YkftQAH2GdK1ynhx2yOPkRv8RtzlA9xCAr+EWAa7PAzWYzzc3N+P1+6vU6GxsbZDIZ3G4358+fZ3h4GKv1cD8kpKORQS5Jb1AulyMUCrGxsUEoFCKdTuvtgiaTSX8Sr9Vq+mGeo9TBFaEywTJXuEEncVK4+QHf5i6nDjXICn6xTq21tRWHw4HFYsHj8RAIBFBVleXlZdLpNA6Hg7NnzzIyMvJkdrj0xsggl6TXTAhBMpkkFAoRiUSIRCKkUimKxSKqqmI0GvVaeKlU0rtQjhLgBlHnNItc5gat7LGLlz/l13nAOKpy+GmBNpuN9vZ2HA4HVqsVr9dLIBBACMHS0hKJRAKbzcbExARjY2O4XK5Dfw/p+GSQS9JrUq/XiUaj+pN3MpkkFouRy+X0AHc6nRiNRorFoj4XXFXVQ5dRjKLGOeb5mJt4yRChjT/gb7PAyKEnEQJYrVa6u7txOBx6Pbyn58lxkYWFBeLxOGazmeHhYSYmJt671WpvGxnkktRg5XKZcDhMOBymWCySy+XY2dkhnU5Tq9UwGo00NTXp80a0F5xaiB+GWVS5wH0+5BZu8oTo5C/5FZZfMMjqVaxWK36/H5fLpS837u3tRVEUFhcXiUQiKIpCb28vk5OTtLS0yF7wt4AMcklqkFwuRzAYJBaLUa1WqVarhMNh4vG4HuAulwuTyaQ/gWsBftgyilWUuchdPuAOToqs08Of8hus03PkAA8EAjQ3N+sthVq74MLCAqFQCFVV6erqYnJyko6ODhngbxEZ5JJ0DEIIEomEPgccngyFCoVCRKNRqtWqvhzh6SdwRVEQQhx6wYNdFPmAO8xyFxtllunnKpcIKv4j3b/FYqGvr08PcK/XS29vrx7gW1tbVKtV2traOHXqFN3d3fIwz1tIBrkkHYF24CUUClEsFjGZTAghCAaDRCIRfS64VkIpFot6gANUq9VD1cGbRI6PuM0F7mGhxiOGucolIkrHke7fbDbT39+P1+vVWwq1AF9aWtLHArS0tOiHeeRqtbeXDHJJOoRyucz29jbhcFgfvwqwvr5ONBqlXC4D4HQ6sVgsFAoFCoUC8KSN77AB3izSfMwtpnmIAZUHjHONWXaVliPdv8Viobe3l7a2NsxmM21tbfr6tOXlZZaXl8nn83g8Hqanp+nr65OHeU4AGeSSdADZbJZgMEg8HkcIgdvtJpPJsLS0RDwep1QqIYTAarVitVopFAoUi0XgaAHuE3tc5gZnWADgLqf4nIvsKUfrDjGbzQQCAbq6urBYLHR0dOg18NXVVR4/fkw2m8XtdnPx4kWGhobkarUTRAa5JL2ANv9Eax80Go14PB6SySTz8/Mkk0n9haXJZMJqtVIqlchms8DRArxdxLnCDU7xmBombnGWz5khqxytP9tsNtPV1YXf78dms9HV1aWXSTY2NlhYWCCVSuF0Ojl37hwjIyP6v2VIJ4cMckn6mlqtpte/S6USVqsVn89HPB7n/v37pNNpCoWC3olisVioVCp6DVw7Wn+YVkK/2OEK1xljjTIWPuciP+cCBeVoJyS1sklPTw9OpxO/308gEMBoNBIMBvUfRDabjVOnTjE2NkZT0+GmHkpvDxnkkvSVUqmk17/r9TpNTU14vV6i0Sirq6sUCgUymQyVSkU/tl6tVikUCvpLzHq9fvBWQiHoI8QVrjPEFgVs/A0fcoNpSsrRFg6bzWb99KXP5yMQCOhjY8PhMPPz8/phnpGRESYmJmhuPtzccentI4Nceu9lMhmCwSC7u7sAuN1uFEUhEomQSCQolUpkMhnK5bJeRlFVVX+xqSgK9Xr94E/gQjDEBt/iOr2EyeHgx1zhFmepKEerS2sDrLq7u+no6KC3t5euri4MBgPRaJT5+XkikQgGg4G+vj5OnTqlb/CRTj4Z5NJ7SQhBPB4nFAqRyWT0/ZdaW2EqlaJSqZDJZPSj89rWG20W+KFPYwrBOCtc4TrdxEjj4of8CnOconbISYQas9mM0+nUa9/9/f36YR2tlr+9vY2qqnR3dzM5OUl7e7sM8HeMDHLpvVKr1djZ2WF7e5tSqaQfgimVSkQiEbLZLNVqlXQ6rY+WhSehrfWGay8vD1pCUYTKFI+5zA3aSZDAw5/zq9xj8kiDrLRRt3a7ndbWVoaHhxkeHtZ3X6bTaR49esTm5ia1Wo329nampqbkarV3mAxy6b1QLBbZ3t5mZ2eHer2O0+nE4/GQy+X0TTblcplMJqPXwYUQ+lO39hLzME/gRlHjLI/4mJv4SBOjhT/iN5ln9EiDrLQA147Qj42N6dvnFUUhl8uxuLjI2toa5XKZlpYWpqamCAQC8jTmO04GufTOEkLsq38riqKPWc1ms5TLZf1lpRbgT2/h0Z68teP0Bw1wk6hynod8xE2aybFNB7/HL/GYoSPNQTEYDJhMJiwWCy6Xi/Hxcc6cOaNPHCwWiywtLbG8vEyhUMDr9XLhwgX9oI/07pP/K0vvHG3hbygUIpvNYjKZaG5uplqtkslk9KDO5XJ6CUXbSK89cRsMBj3AD1pCsYgyF7nHB9yhiQKb+Pk+v8YqfccKcLPZTFNTEyMjI1y4cEEP8MczRWoAACAASURBVHK5zOrqKouLi+RyOdxuN7OzswwODsrDPO8ZGeTSO6Narer173K5jNVqxe12UyqV9AM9BoNBnw2uDbDSttELITAYDBiNxkMFuE0UucQcl5jDTplV+vgDLrGlBI70ObR7sFgsOBwOBgcHuXTpEl6vF3hS519fX+fRo0ek02mcTifT09OMjIxgsx2tbVE62WSQSydesVjUt+9o9W+LxUI+n6dcLuNwODCZTEQiEZLJJKVSiVwuR6lU0sNaewIHDrzYwSnyfMhtZriHlSqLDHGVS4SVziN9jqcD3Gaz0dfXx+zsLB0dTwZj1et1/TBPIpHAbrdz+vRpxsbG5Gq195wMculEEkKQTqcJhUJ6/bupqUkvmSiKgs/nQ1VVtra2SCaTVCoVcrkc+Xx+39O29iLw6br4y7hFlo+4yXkeYERlnlGuMUtMaTvSZ9FeYmoB7vf7mZmZIRAI6C9Yd3Z2ePDgAfF4HIvFwvj4OBMTE3K1mgTIIJdOGFVVicVihEIhcrkcJpMJl8tFtVolm81iNpvx+/1Uq1VWV1dJJpP6C81sNqv3gAP6E/hBA9wrUlzmBmd5BMB9JrjGLEnFe6TPoigKRqMRm82G1Wqls7OTc+fO0d/fr4/FjUQiPHz4kGg0isFgYHBwkMnJSb3MIklwzCBXFOU/Af5nYAKYFULcasRNSdLXadt2tre3qVQqWK1WXC6XHtBOp5Ph4WGKxSKLi4skEglUVaVYLJJKpV4Y4AfRJna5zA2meIyKgduc5gsuklbc+37fXLWbaXP4lddTFAWDwYDNZsNut9PW1sbExIRe49aGdWmHeYQQ+P1+Tp06pfeKS9LTjvtE/hD4j4F/2YB7kaRnFAoFvf6tqioOhwOz2azXv30+H93d3WQyGebm5kgmk9TrdSqViv40rjlsgHeJKFe4zgQrVDDzJef5ORfIKc8fLnW39vIg1wLcarXidDppbW1lYGCA0dFRfd5JOp1mfn6era0tarUaHR0dTE1N0dnZKQNceqFjBbkQYgGQ/4BJDSWEIJVKEQwGSSaTKIqCw+FAVVUKhQIGgwG/309XVxeJRIKf//zn+hN4vV4nHo8/E+BPv8h8lR6xzRWuM8IGRaz8lEtc5zxF5WjjXbUAN5vNuN1ufD4ffr+f4eFh/Th9LpdjYWGB9fV1yuUyra2tTE1N6QOvJOll3liNXFGUT4BPAHp7e9/Ut5VOEK3+HQwGyefzGI1GnE6nPiLWarUyNDREe3s7kUiEzz77jN3dXf2wTiwWo1Kp7LvmgUNcCAbZ4grX6SdEHjufcpmbnKWiWI/0ebQAN5lMeDwevF4v7e3t9Pf309vbi9Fo1A/zLC0tUSwW8fl8XLx4Uf/rknQQrwxyRVE+BZ7XT/U7Qog/O+g3EkJ8D/gewMzMzMEn7UvvvEqlQjgcJhwOU6lU9P5pbc+l2+1mZGSE1tZWtre3+fGPf6xPKqzX688NcK2V75Ub6oVglDWucJ0AETI08Vf8Mrc5feRBVlqAG41GfD6f/p/Ozk6Ghoaw2WyUy2UWFxd5/Pixfpjn7NmzDAwMyNVq0qG9MsiFEN95EzcivX/y+by+bb5er2O327HZbBSLRX1zeyAQwOVyEQ6H+cu//Es9wFVVJRqNPjfALRYLtVrtpRvqFaEyyRJXuEEHu+zRzPf5DveYpK68+l9U56rd3K11P/P13y3O7Pv13/EZ+AeD7QwPD9Pc3EytVmNpaYlHjx6RyWRoamriwoULDA8PY7Ue7clfkmT7ofRGCSHY29sjFArp9W+bzUa9XqdUKmEymejr68Pv92OxWNjZ2eHatWvEYjHgyanGeDz+TIBrR9lrtRqlUumF398g6pxhgY+5SSt7xPHxx3yXh4wfapDVtDmsv9jUatj/On+ef9q3RWdnJw6HA5fLxeDgIB0dHaiqysbGBg8fPiSZTGK32zlz5gxjY2NytZp0bMdtP/yPgH8OtAF/oSjKXSHErzfkzqR3ilYCCQaDFAoFvX+6Wq1SLBZxOBz6LG2j0cjOzg5zc3N6gFerVWKx2L42Qngyj9tqtVIul/Vlx89jFDWmecjH3MRDlh3a+X1+iwVGjjQHBX4R4AaD4cn8kzx66aSnp4fe3l4URWF7e5sHDx6wu7uLxWJhYmKCiYkJuVpNapjjdq38CfAnDboX6R1UqVT09WnVahWz2YzNZtOPx2vryLRRrJFIhLm5OaLRKEIIqtUq8Xj8mQC3Wq2YzWbK5TK5XO6F398sKsxwnw+5jYs8Qbr4C77NCgPHCnCtDu5yuQgEAk9COVghEAgwODiI1WolGo3y8OFDIpEIRqORoaEhTp06JVerSQ0nSyvSa5HL5fT6txACq9WKxWKhXC5jNBrp7u7G7/fjdDoBiEajzM3N6f3ilUqF3d3dZ15UWq1WTCYTpVJJX7X2PFZR4hJ3ucQdHJRYo5c/4jfZJHDkANdOYmrtkD09PXg8HgwGA263G4gwMTFBIpHgyy+/JBwOI4Sgp6eHU6dO0dLSIlt1pddCBrnUMEIIkskkoVCIvb09FEXRXzxq0wiHhobo6urSOzNisRh37twhEolQq9UoFAqk0+lnWgUtFgtGo5FyufzSAHeIAh9wh4vcxUaFxwxylVm2lWdfTB6UFuDaaUy/3097ezuqqmKz2fQ6+D/eecAXX3zB1tYW9Xqdzs5OTp8+LVerSa+dDHLp2Or1OtFolFAopB/YsVqtVCoVyuUybrebQCBAa2urXleOx+PcuXNHbznM5/Pkcrl9py61YVImk+ml9W8Al8jyEbe5wH1M1HjEKFeZJaq0H+uzmUwmfSJhe3s7vb291Ot1fYlxb28vxWKRmzdvMpBfZ22vQmtrK2fOnNGXH0vS6yaDXDqycrn8TP1bK58IIWhra8Pv9++rCcdiMebm5tje3qZQKFAsFikUCvuu+/QTcKVSeWkboUek+ZibnGMeA6o+yCqh+I712bQAN5vN+Hw+BgcH9UmEHR0dDA4OIoTgwYMHLC8v64d5Ll26RE9PjzzMI71RMsilQ8tms4RCIWKxGKqqYrFYMJvNepj39vbi9/v3LTmIx+PMzc0RDAb1jTzPO4WpPcF+/eXm17WKBJe5wWkWUTFwl1N8zkVSyvFeJGptjAaDgebmZoaHh/UDPE1NTfqvV1ZW9M08Ho+Hc+fOMTAwIFerSd8I+U+ddCBCCBKJBKFQiFQqBaD3bVerVRwOBwMDA3R2du57Gt3d3WVubo61tTVyuRzlcplarfbM4CptO/2rtvJ0iBhXuM4ky1QxcZ1pvmDmhYOsDspkMmG1WvUXmYODg7S0tOjlnvHxcXw+HxsbG/pmHpfLxcWLFxkeHpar1aRvlAxy6aXq9TqRSIRQKESxWERRFMxms17y8Hq9BAIBfD7fvumCu7u73L59m9XV1WfWqWme/v2vGmYVEGGucJ1R1ilh4RqzfMl5CsrxNuMYjUZ9u47FYiEQCNDf308mk6FQKOiHk7a3t/nRj37E3t4edrud6elpRkdH5Wo16a0gg1x6rlKpxPb2Njs7O1SrVYxGI0ajkVqtRr1ep7u7m0AgoLcPwi8C/Msvv2RtbY1CoUC9Xn9mdZq21PiV42SFoJ8gV7jOIEEK2Pj/+JgbnKWsHC9AjUajfiDHYDDQ0dHBxMQEuVyOVCqlD7dKJpN8+umnJBIJLBYLk5OTTExM7PvckvRNk0Eu7ZPJZAiFQsTjcVRVxWQyYTAY9Fp4b28v3d3d+wY7CSGIRqNcu3aNra0tSqWSHt5PP2kfJsBHWOcK1+lhhyxOfsQvcZvTVJXjlTC02rd2316vlzNnzlCr1Ugmk7hcLiYnJymVSnz++edEo1GMRiPDw8OcOnXqq35xSXq7yCCX9CfpUChEOp1GCIHZbNYnB2qnF9va2va10wkh2Nra4tq1a3oboRbUzwvsVwW4IlQmWOEy1+kiTgo3P+Db3OXUgQZZvYx2jN5oNFKpVGhqauLMmTM4HA52d3exWq2MjY1hMpmYm5sjHH4yR6W3t5epqSl8vuN1wUjS6ySD/D1Wq9X0+rc2aEorn9RqNdrb2wkEArjd7n0HWur1OouLi9y4cUOff/KyAH8Vg6gzxWMuc4M2kuzi5U/5dR4wjqocr41PW8JsNpspFouYzWbOnz9PZ2cnkUiEUqlEf38/brebR48esbW1haqqdHV1cfr0adra2uRhHumtJ4P8PVQqlQiFQuzs7FCr1fTZIaqqoijKc9sH4cnatfv373P37l1SqZTeYaKVUJ5XB38Zo6hxjkd8zE28pInQyh/yt3nEyKEmEb6Iz+fDbreTz+dRFIXx8XFGR0eJRCKEw2Ha29vp6OhgZWWFL7/8Uh+dqx3mkQEunRQyyN8j6XT6mfo3PAliu91OIBB4pn1QK7vcu3dP75t+ukXweZ0or3oqN4kqF3jAR9zETZ4QnfwVv8wSg0eeg/I0j8eDx+MhnU5TKBTo7e3l3LlzJJNJtra2cLlcDA0NEQ6H+clPfkKpVKKlpYWpqSl6enrkaUzpxJFB/o4TQhCPxwmFQmQyGVRV1YNaVVV9+uDT7YPwi6f2Bw8e6KvXng7nowS4VZS5yD0+4DZOiqwT4E/5Luv0NiTA3W43ra2tpNNpUqkUbW1tzM7OUqlU2Nzc1Ge9pFIprl69Sj6fx+PxcP78efr7++VhHunEkv/kvqNqtRo7Ozv76t9aUCuKgt/vx+/375uJrQ29CgaDPH78mJ2dHUql0r5wft5hHu3PvohdFPmAO8xyFxtllunnKpcIKv6GfNampia6urpIp9MkEgmam5u5ePEidrud7e1tAHp6eiiXy9y6dYtMJoPL5eLSpUsMDQ3J1WrSiSeD/B1TLBb39X9rJya1SYTd3d10d3fvO4lYqVTY2dlhc3OTYDBIPB7X56XAk5B+UYC/TJPI8SG3meE+FqosMMxVLrGjdDTks2rloFwuRywWw+Fw8NFHH9Hd3c3W1ha7u7u0tbVhNBqZn59nb28Pp9PJ9PQ0Y2NjcrWa9M6QQf4OEELo9W9thrfRaNSDV5s+2N7ertd/hRCkUinC4TChUIhIJMLe3p4+IlY7Lv90CeUgLzABmkWaj7nFNA8xoPKQMa4xS1xpbcjn1frZK5UK8Xgci8XCmTNnmJycZHt7m5WVFZqammhpaWF1dVVvLzx9+jTj4+P6SU5JelfIID/BVFXV69/ZbFYfr6rp6OggEAjQ3Nysl1Wq1SqRSITt7W0SiQTxeJxUKrWvB1wLcNhfMnlViPvEHpe5wRkWALjHJNe4yJ7ibcjnNZlM9PT0oCgKiUQCg8HA0NAQFy5cYG9vj8XFRSwWC62trWxvb7OwsIDJZGJ0dJRTp07hcrkach+S9LaRQX4CVatVdnZ22N7e1k9RKoqiz+/u6urC7/frS321J/ZwOMz/fS3ItztKJBIJMpmMHuCqqurH77Un74O237WLOFe4wSkeU8PILc7yBTNklMYEp6IoBAIB7HY7yWRSHxFw8eJF6vU6y8vLAHi9XuLxOGtrayiKQn9/P1NTU3i9jflBIklvKxnkJ0ihUNDr39qYV60L5en2Qa37olarEY1GCYfDZDIZMpkMf7ZqYai4QaVS0U9ufr3+/XRt/GW6xQ5XuME4q5Qx8zkX+ZLz5JXGzSHp6OjA5/ORTCbJ5XL4fD5mZ2dxOp1sbm7qpzSz2Szz8/Ooqkp3dzenT5+mtbVV9oJL7wUZ5G85rZYdCoVIJBL6AR6tD1xrH3x6H2Q2myUcDhONRvXt8ru7u18tKR6iVCrppzfh4LXvr26IPkJc4QZDbFLEymd8yHXOUVLsDfvcXq+X7u5ukskkkUgEl8vF7OwsPT09bG1tEQ6HsVgsKIrC8vKyfhL1zJkzdHZ2ygCX3isyyN9SqqoSi8WeqX9r23O06YNa+6C2bi0cDpPNZvU1a1qAPx3cpVJJX+Lw9X7wFxKCITb4FtfpJUwOBz/mCrc4S+WYg6ye5nQ66evrI5PJsLOzg81mY3p6momJCaLRKPPz8/p7gK2tLSqVJ6vVTp8+jd/vl4d5pPeSDPK3TLVaJRwO76t/A/oeTL/fv699MJfL6U/ftVoNVVUpFArE43H9AFC1WtXr39q1XrXAQScE46xwhRt0EyWNix/yt/i3tV/jlDnesM9tsVgYGBigVCoRjUYxm82Mjo5y5swZCoUCjx490n/oxONxCoUCXq+X2dlZ+vr65Go16b12rCBXFOV/A/4OUAFWgd8WQqQacWPvm3w+TygUIhqN6ivQhBCYTCZ9+qDWPqgte9Bq30IIDAYDmUyGeDxOPp+nXq/zRcbLrdKzy4f/VW5636/PmcJMm8P7vqYIlVM85go3aCdBkmb+nF/lHpOoipEbtb6GBLmiKAwMDACQTCaBJ4d3pqenMRqNrK+v6z/Q9vb2KBQKuN1uPvjgA4aGhuRpTEni+E/kPwb+mRCipijK/wr8M+CfHv+23g9a/TsYDJJMJqlWq3q92mw209rauq99sFAoEA6HiUQi1Go1jEYjJpOJeDxOLBbT16hVKhVqtRrnTAWmXU8CWiur/G5xht+233rhPRlEnbM84jI38JEmRgt/xG8wz1hDBlk9rbu7G4fDoZeOtBq3z+djc3OTTCZDtVrVt/U4nU4uXLjA6OioXK0mSU85VpALIf76qV9+Cfz9493O+0FVVaLRKKFQiFwupwe4wWDAYrHQ2dmpt9tpveLhcFjflelwOBBCEIlE2N3dpVwuo6qqHuDwZByt1lJ4ECZR5TwP+YhbNJMlTAf/gW+xyHBD5qA8zefz4fP5yOVy7O3t4fV6mZycpLe3l0gkwsOHDymVSuRyOQqFAjabjTNnzjA+Pq63VEqS9AuN/PfSfwj8hwZe751TqVT0+rf29AxPlhjb7Xb8fj9dXV2YTCaKxSJra2v6UXur1YrH4yGVSrG6usre3p7eQvh0gGsdLV/fUP8iFlFhhnt8yG2aKLCJn+/zHVbpb3iAOxwOuru7yefzpFIpXC4XU1NTjIyMkMvlePjwoR7e2uzw8fFxpqam5Go1SXqJVwa5oiifAp3P+Uu/I4T4s69+z+8ANeDfv+Q6nwCfwJOtK++Tp+vfTx+BN5vNeDweAoEAra1Pjq8nEgnC4bBeL/Z6vSiKQiQSYWVlRS83aGH99AtM7cXmQdhEiUvMcYk72CmzSh9/wCW2lMAzv3eu2s3dWvczX//d4sy+Xz+v1g5PTmT6/X5qtRqpVAqHw8HQ0BBjY2MYDAaWl5dJp9PkcjlKpRJGo5HBwUGmpqbweDwH+jyS9D5TDjsI6ZkLKMo/AP4R8G0hROEgf2ZmZkbcuvXiOu27QAjB3t7evvo3oA+v0rbvuFwuyuUyOzs77OzsUC6XsVqt+Hw+KpUKGxsb7O7u6pvogX1P4PCkVPOqLfSataqLT0x/wQz3sFJlkSGucomw8ryf1S/2qlq7prOzE7PZTLlcxmw2093dzejoKB6PR/9sqVSKcrmM0WjE7/frdXLZCy5J+ymKclsIMfP1rx+3a+W7wH8H/NJBQ/xdp/Vzh0Ih8vm8XuIwmUzYbLZ90weTySQPHjwgkUgAT56+Ozo62Nvb4+HDhySTSb1jQ3varlar+05eHjTA3SLLR9zkvOkBJurMM8ZVZokpba/l70NzczMul0vvX+/u7mZoaIju7m52dna4c+eO/vkMBgNdXV2cPXuW9vZ2GeCSdEjHrZH/C8AK/Pir//N9KYT4x8e+qxOoUqmwvb1NOBzWT04KIbBarTQ1NREIBOjo6NDnhGuzvs1mM4FAALPZTDAYZGFhQR9iBU/CulKp6DNRtK8dNMC9IsVlbnCWRwDcZ4JrzJJs0CCrr7PZbLS0tOjH/1tbW+nt7WVgYIBMJsPc3JzeB240GvVOFb/fLwNcko7ouF0rw426kZMql8vtq39rEwStVuu+9sF0Os3CwgK7u7sIIfB4PPT29lIul1laWiIajeonMLUOlmKx+MrN9C/SJna5zA2meIyKgduc5gsuklbcr+Xvg9FopKWlRe+Wcbvd9PT06D3iCwsLhMNhcrkcBoOBtrY2pqam6Ovrk6cxJemY5GmKI3h6k442w1sLX7vdTmdnJ36/H7PZTCQSYWlpiWKxqL/08/l8xGIxbt68qS9x0IZfmUwmCoXCvs08B336BugUUa5wnUlWqGDmS87zcy6QU5pe/YcP4ZzpFy81PR6P3haobevp6+vD7XaztrbG1tYWe3t7GAwGfD4fp06dYnBwUB7mkaQGkf9POgTtROX29ja5XE7vQLFYLDgcDr18UiwW2dzc1JccNzc309vbi81mY2Njg9u3b7O3t6cPwDKbzQghyGaz+wL8ME/gPWKbK1xnhA1KWPkZl/iS8xQbOMjqadPmMA6HA6fTiclkwuFw6C9wu7q6CIVC3L17V6//ezweJicnGRkZkYd5JKnBZJAfQLlc1uvfxWJRn/lhs9nweDz4/X68Xi+xWIx79+6Rz+cxGo10dXXR2dlJoVBgaWmJ9fV1stks8KR33Ol0Uq/XSaVSr3wCn6t2P9vaJwSDbHGF6/QTIo+dn3CZm5ylrLy+NWYWi4WmpiZsNhs2m422tjba2toYGBggmUxy7do1IpGI/kNsfHyc8fFxbDbba7snSXqfySB/iWw2q9e/taBVFAWbzUZ7e7v+gi4cDrO4uIiqqrhcLsbGxvD5fEQiEa5evaq/ADUajTgcDkwmE6VSad9uTEVRXlpCuVt7KsiFYJQ1rnCdABEyOPkrfpk7nKaqvL5Fwtrcc7fbjc1mw+v10tLSwsDAAPV6nZs3bxIMBqlWq7hcLkZHR5mampKr1STpNZNB/jVCCBKJBMFgkFQqRbFYBH5x+rK7u5vOzk5SqRTLy8v6y7uOjg66u7sxm82sr69z7do1fX+m1hcOT16ORqNRvY1Qm61ykBKKIlQmWeYy1+lklz2a+T7f4R6T1JXX9z+l9sPL7XZjt9vxeDx4vV56e3tpamri4cOHrK2tUSwWaWpqYmJigtOnT+N2v54Xq5Ik7SeD/Cv1el1fn6adMNQCTGsfdDqdRKNRbt++Tb1ex+l0MjIyQkdHB4VCgYWFBRYXF0mn0xgMBlwuF06nUz/RqB2rBw4V4AZR5+8bf8p/xe/Tyh5xfPwJ3+UB4w0fZPV1ZrOZpqYm3G43zc3NeDwevWS0tLTEZ599Rj6fx2azMTk5yblz5+RqNUl6w977IC+VSvr6tEKhoLf/ORwOWlpa6Orq0nu/M5mM3jrX3d2Ny+UiHo/z2Wefsb6+TqFQwGKx0NXVhdVqpVAoEI1GSafT+hO4dpT+IAFuFDWmecjH3MJjzrBDG7/Pb7HI8GsPcK0MpD19Nzc34/V66e/vZ3t7mx/84AekUiksFgvDw8NMT0/T1tYme8El6Rvw3gZ5JpMhFAoRi8UoFAr67O+mpiY6Ozvx+XykUikeP35MrVbDbrczNDSkrxELBoP85Cc/0UfKak/nBoOBvb09fVOP9mJUC7iDtBIuVFuYEo/4xPQXtCsp7qjD/PPaP+Jv1HPAL4LyRbNNjkP5/9u7s9i4rvuO498zM/fOvnK4DMnhKnMRtVAVtVhLmzRuliatH/rUoi3SPhgFGjQB0gZN/NCHokCBFk0fWqAwCgQFGrhNXddu6zq1HEeunCiiZMsyJWrfuEjivs1CcpbTB2puJFm2uA1HQ/4/T+KQnDlXBn6++p///R+lrAFdlZWVBAIBfD4fLS0tJBIJjh07xujoKHa7ncbGRvbu3UssFpMAF6KEtlSQa60ZHx//WP3b6XTi9XqJxWKYpsno6Ch37txBKUU0GqWuro5gMMjCwgIXLlygr6+PqakplFLWmZkLCwuMjo4+dLADLA2zWm4Jxann2c+H/InjAzzMc4M4b/Ict2xxfpzfs6zZJmthGAaBQICamhqrHt7Y2Ihpmpw8eZKhoSG01sRiMXp6eojH4/IwjxBPgS0R5IXSyKP934WNu2g0aj1in8lkcLlctLS0UFNTg2maTE9P8+6773L16lWSySSmadLY2Eh9fT0TExPcvHmT8fFx0un0Q0ezLXeYlUenOMgH7ONDXCxymRZOsJ9h9fGJg8Vgt9vxer3U1tYSDoet8lA0GuXMmTNcv36dbDZLZWUlPT09tLS0yNFqQjxFNnWQz8/PMzQ0xN27d0kmk1b92+v1UllZidfrZXZ2lhs3bgAQjUatMAMYHBzk7NmzVsB7PB52795NTU0Nt27d4uzZs0xMTFgthPDzTczlBLhPJzjEGfbyEQZZ+nmG9zjAPfXx49mKobCZW11dTSwWwzAMIpEI8Xic/v5+jh8/zvz8POFw2DqZxzCK194ohFidTRfkWuvH1r8Ls7+j0ShKKSYmJhgbG8PpdNLU1GRtUGYyGfr6+h6aSlh4rNzn83Hp0iWOHz/O9PT0Q5MIV9KFEtIzHOY03VzARp6P6OQn7GNcVRT17+ZBpmkSjUaJx+M4nU48Hg/Nzc0MDg7y6quvkkgk8Pv9HD58mJ07d+J0Fu8BIyHE2myaIC+crl6of6dSS1N13W43fr8fv99PJpPh3r17wFI419bWUlFRgVKK2dlZent7uXz5MslkEsMwiMfjdHd3k8vl+OijjxgeHmZ2dvahDcyVBHiFnuQIveziInlsfEgXP2Ef0yr4xN99cLbJWtjtdkKhEE1NTXi9XgzDoLGxkUQiwZtvvsn09DQul4u9e/fS09MjR6sJUQbKPsgL9e+hoSHm5uasAVZer5dQKIRhGCQSCUZHRzEMg4aGBmKxGG63G601Q0NDnDt3jsHBQRYWFvB6vezcuZOdO3cyMTHBqVOnPvZofqFDY7lzUKr12P1BVlfI4OAUezjJXuaUf9nXudbuFJvNhtfrpaGhgYqKCmsGuGEYvPfee4yOjuJwOOjs7OTgwYPyMI8QZaRsgzydTlv178L4V4fDYd19K6WYm5uzRsa2trYSjUax2WwsLi5a5ZPJyUny+TzhXZiSigAAEJRJREFUcJiDBw/S2trKwMAAb7zxhvUIfaEDZaUBXqfvcJRe2rnBPCbvsZ+f8Quk1MY+su50Oqmvr6empgaHw0FFRQWhUIjTp08zNDSEUoqmpiaOHDlCRcXGlXeEEOujrIJca83MzAxDQ0NWm18+n8fpdBIIBKwa99zcnDUytra21jp1fnp6mgsXLnDlyhUSiQQOh4Oamhq6u7upra2lv7+fl19+2ZpMmMvlHuqPXlaAa00TgxzlFC0MksLFOxyil24W1MYOjXI4HFRVVVl1cJ/PR3V1Nf39/bz77rvkcjlqamo4cuQItbW10gsuRJkqqyC/evUqAwMDD9W/3W43drudxcVFkskkwWCQlpYWqqqqsNls5HI5bt++zfnz5xkaGmJhYQG32017ezu7d+8mFArR29vLD3/4Q6v+nc/nV3z3jdY8w02Ocoo4d5nDy1v8ImfYRUZt7NjWQn97YRaK0+kkFosxODjIG2+8weLiIuFwmEOHDtHa2ioBLkSZK6sgn5+fJ51O4/V6MU0TrbXVOVJbW0ssFsPnWzpAIZlMcv36dS5evMjExAS5XI5gMMiuXbvo7OzE4XBw4sQJLl26RCqV+ljP90oCvJOrHKWXGKNM4+cNPsdZuoo6yOqT+Hw+4vG41Q9eGPD1ox/9iGQyid/v59ChQ+zYsUN6wYXYJMoqyP1+P4FAgFwuRy6Xw+/3U1tbS1VVlXXE2MjIiDX7uzAbJRqN0tnZybZt20in0xw7doxbt25Zx6gtd/bJg5TOs5NLHKGXSiYZJ8xrfIE+OsirjQ9I0zSpr68nEolYc2Ky2SwnT55kenoat9vNgQMH6OnpkV5wITaZsgrywgM9hdq337/U9TE/P8/Nmze5fPky9+7dY35+HtM0aWpqoquri/r6eu7cucMPfvADazbKagPcrrPspp/DnCbCDPeI8gpfpp9nij7I6rHrsdupqakhGo1aUwoNw6Cvr4/x8XFM02THjh08++yzeL3eDV+fEKL4yirIGxsbaW5uxuFwWHPDBwYGuH79ulU+8Xg8dHZ20tnZSWVlJefPn+ett95iamrKah9cTYA7dIa99HGIMwRIMEw1L/MZrtACJagxK6UIBoPU1tbi9XoJBoN4PB6uXr3KyMgINpuN5uZmDh8+LJ0oQmxyZRXkpmmyuLjIwMAAt27dYnh4mJmZmYdaDNva2jBNk97eXl555RUSiYQV3IU2whV9pl5gP+c4yPt4SXOLel7nC9ygoSQBDkubvHV1dVavvM/nY3h4mOHhYasT5eDBg9TX18tGphBbQFkF+c2bN7l06RKjo6PMz89bJ/O0tbXR0tJCMpnkxIkTXL161TqabbUB7tZpDnCW/ZzFzQJXaeIEBxhUdUW4suUptEv6fD5CoRDBYJCJiQkuX77M4uIikUjEmokiUwmF2DrKKsivXbvGwMAATqeThoYGOjo6qK2tZWBggNdee43h4WHrBJ58Pr+qAPfqJIc4Qw8fYZLhIts4wQHuqur1vpxlKzxWHw6HCQQChEIhkskkFy5cIJVKEQgE2LdvH11dXTITRYgtaE1BrpT6c+B5IA+MAl/VWq/vSQcPKJyJ2dHRQTAYpK+vj3feeYfJyUlrA7PQB75SQT3LYU6zh/PYyHOedt5jP2MqWoQrWR6lFD6fj3A4jM/nIxKJkM1muXbtGqlUCrfbTXd3N93d3dbGrxBi61Er3fR76JeVCmitZ+//+Y+A7VrrP3jS7/X09OgzZ1Z+SEI+n2d8fJy+vj6uXLlincBTuPtezbVE9JQ1yArgHNt5j31MqdKeO+lyuQiFQni9XsLhMA6Hg5GRERKJhDUPfffu3XK8mhBbiFLqfa11z6Ovr+mOvBDi93mB1f9fYRmOHz9Of3+/9QBPLpcjm82u6r2q9BhH6KWLK+SwcYZd/JR9zK5gkFUxmKaJz+ezNjJdLhdTU1PMzc1ht9tpaGhgx44dNDQ0SB1cCAGsQ41cKfUXwO8CM8BnP+XnXgBeAGhoaFjVZz04IGs19W+AWn2Po5yig+ssYHCSvZxkL0n16T3WZzO1634+5oMcDgcejwePx4Pf78fr9ZJIJBgbG7M2dTs6OmhpaZE6uBDiIU8McqXU20DNY771otb6da31i8CLSqlvA18D/uxx76O1fgl4CZZKK6tZ7OTkpHVM20o16CGOcopt3CaNk+Mc5BR7mFfLm7f9YbY4QW4YBk6nE9M08Xq9+Hw+MpkMIyMjAFRUVNDc3ExbW5vUwYUQj/XEINdaP7fM9/o+8D98QpCvh0QisbJf0JpWbnOUUzQyTAIPb3OE03SzuMGDrB5lGAamaWK323G73Xi9Xmw2G3Nzc+TzeYLBIPX19bS1tVmnGgkhxOOstWvlGa311ftfPg9cWvuSPlkmk1neD2pNB9c4Si+1jDCDjzf5LB+wg6wq7ZwRh8OBaZo4HA4cDgdutxvDMKy2Sa/XS3V1Na2trdTX10sdXAjxRGutkf+lUqqdpfbD28ATO1bW4kldKUrn6eIyR+mligkmCfKf/Arn2F6SQVYPstlsuFwuDMNAa22VU5RSZLNZPB4P0WiUhoYGmpubMc3S/otBCFE+1tq18hvrtZC1sOkcu7jIEXqpYJpRKniVL3Ge9lUNsjqbqeXDbO3HXv9e+uGun27HnSfWzZVSeDweDMOwhn6ZpolhGNhsNtxuN5FIhKqqKrZt2yZ1cCHEipXVk52PcugMezjPYc4QZI47VPOv/BqX2LamOSh7jI8H9PfSPfyee2W97263G4/Hw+LiItlsFsMwrM1Nj8dDMBgkHA5bx9BJHVwIsRplGeSmXqSHczzL+/hIcZs6/ovnuE5TyQZZPcgwDGtuejqdxmazYRiG1ZVSOFe0sbFR6uBCiDUrqyB36XkOcJYDfICbBa7TyL9xgAFVX+qlAUt18EAggMPhIJVKobXG4XBYT2f6fD48Hg+xWEzq4EKIdVNWQf5FfsxuLnKJVk6wnzsqVuolAT+fieL3+5mbm7Puwj0eD5WVlQSDQZxOJ+FwmG3btlnH0QkhxHooqyD/Pw7yU3oYVZUb/tndjsdvarrdbioqKkgmk4yPj1uvVVVVEY1GrSc2pQ4uhCiWsgryyRIOsnp089MwDKLRKFprxsbGyOVymKZJVVWVVfd2OBxSBxdCFF1ZBfnTwGazEQwG8fl8TE5OMj8/j91up7KykpaWFkzTJJ/PSx1cCLFhJMiXyWazWZuWs7Oz3LlzB6UUfr+fjo4Oa0ZKIBCQOrgQYkNJkD+BUgqXy0UkEiGTyXD37l1yuRwul4v29nZqampIpVLY7XaZiyKEKAkJ8k+glMIwDCKRCA6Hg7GxMTKZDHa7nZaWFtrb20kmkywsLNDS0iJ1cCFEyUiQP0Iphd1ux+fz4fP5mJqaIp1Oo5SisrKSPXv2kMlkSCQSUgcXQjwVJMgfUGgV9Pl8pFIp7t69i9Yan89Hd3c3Ho+HdDpNKBSSOrgQ4qnx1AZ5YdLhRtSbC8OrfD4f2WyW0dFR8vk8pmnS1tZGPB4nkUigtaarq0vq4EKIp8pTGeQbFeI2mw2n04nP5yOfzzM9PW3VwePxOJ2dnaTTadLptNTBhRBPracyyIsd4IVRsoVTeRKJhHWwQyQSYfv27djtdtLptNTBhRBPvacyyIvJMAx8Ph92u51UKsXi4iJaa1wuF62trVRUVJDP56UfXAhRNrZMkBemEBqGQSqVIpPJkM/nMQyD6upq4vG4VWopBLrUwYUQ5WDTB3nhcGO3200qlSKdTpPP57HZbEQiEeLxOB6PB9M0ZS6KEKIsbdogL9xd+/1+0uk0MzMzaK2t8bKF6YROp1Pq4EKIsrYpg9zpdOL1esnlckxPT1sBbpomfr+fqqoqAoEAkUhE6uBCiLK3qYLcNE3cbjd2u51kMkk+n8fhcOBwODAMg3A4TGVlJaFQSOrgQohNY1MEeeFAY9M0SaVS5PN57HY7LpcLwCqlRCIRmpubqaurkzq4EGLTKKsgV0pZDwvB0kZmIcAXFxdJJBLWnJRcLofdbiccDhOLxYjH4zQ1NUkdXAix6axLkCulvgn8NVCptR5fj/d8HIfDQSaTAcDlcmGaJrlcjmQyid1ux+/3o5RicXERr9dLXV0ddXV1UgcXQmxqaw5ypVQc+DwwsPblfLpYLMbIyAiGYQCwsLBgHe7gdDpJpVLYbDZisRiNjY0888wzUgcXQmx663FH/l3gW8Dr6/BenyoSiTAxMUE2mwWwJhUmEgkSiQSBQIDm5mba2tqkDi6E2DLWFORKqeeBYa31uY24600mk2SzWdxuN6FQiPn5eaampjAMg+bmZrZv3y794EKILeeJQa6Uehuoecy3XgS+w1JZ5YmUUi8ALwA0NDSsYIk/F4vFyGQyaK2ZnJxEa01lZSVdXV20t7dLHVwIsSU9Mci11s897nWl1E6gGSjcjdcDHyil9mut7z3mfV4CXgLo6enRj35/OTweD3Nzc6TTaTweDx0dHezevVvq4EKILW3VpRWtdR9QVfhaKXUL6Clm18qNGzdYWFigsbGRAwcOWIOuhBBiKyurPvJnn32WmZkZOjs7pQ4uhBD3rVuQa62b1uu9PklTU9E/Qgghyo7UJYQQosxJkAshRJmTIBdCiDInQS6EEGVOglwIIcqcBLkQQpQ5CXIhhChzEuRCCFHm1IMn7mzYhyo1Btze8A9euyhQtBEET6Gtdr0g17xVlOs1N2qtKx99sSRBXq6UUme01j2lXsdG2WrXC3LNW8Vmu2YprQghRJmTIBdCiDInQb4yL5V6ARtsq10vyDVvFZvqmqVGLoQQZU7uyIUQosxJkAshRJmTIF8FpdQ3lVJaKRUt9VqKTSn1V0qpS0qpj5RS/6GUCpV6TcWilPqiUuqyUuqaUupPS72eYlNKxZVSP1ZK9SulLiilvl7qNW0EpZRdKXVWKfXfpV7LepEgXyGlVBz4PDBQ6rVskGPADq31LuAK8O0Sr6colFJ24O+BLwHbgd9USm0v7aqKLgt8U2u9HTgI/OEWuGaArwMXS72I9SRBvnLfBb4FbIldYq31W1rr7P0vfwbUl3I9RbQfuKa1vqG1XgT+BXi+xGsqKq31Xa31B/f/PMdSuNWVdlXFpZSqB74M/GOp17KeJMhXQCn1PDCstT5X6rWUyO8Db5Z6EUVSBww+8PUQmzzUHqSUagL2AKdKu5Ki+1uWbsTypV7Ielq3w5c3C6XU20DNY771IvAdlsoqm8qnXbPW+vX7P/MiS/8U//5Grk0Un1LKB/w78A2t9Wyp11MsSqmvAKNa6/eVUp8p9XrWkwT5I7TWzz3udaXUTqAZOKeUgqUSwwdKqf1a63sbuMR190nXXKCU+irwFeBzevM+eDAMxB/4uv7+a5uaUspgKcS/r7V+tdTrKbLDwK8rpX4VcAEBpdQ/a61/u8TrWjN5IGiVlFK3gB6tdTlOUFs2pdQXgb8BfklrPVbq9RSLUsrB0mbu51gK8NPAb2mtL5R0YUWklu5I/gmY1Fp/o9Tr2Uj378j/WGv9lVKvZT1IjVw8yd8BfuCYUupDpdQ/lHpBxXB/Q/drwP+ytOn3g80c4vcdBn4H+OX7/20/vH+3KsqM3JELIUSZkztyIYQocxLkQghR5iTIhRCizEmQCyFEmZMgF0KIMidBLoQQZU6CXAghytz/AzTOxZsMFvs6AAAAAElFTkSuQmCC\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]}]}