Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
BDA/BDA 3.10.8 - Stan.ipynb
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
299 lines (299 sloc)
66.1 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "code", | |
| "execution_count": 44, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stderr", | |
| "output_type": "stream", | |
| "text": [ | |
| "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_b2ef2a6e4a8e372508b34c24c518abeb NOW.\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "\n", | |
| "# coding: utf-8\n", | |
| "\n", | |
| "# #### Problem 3.10.8\n", | |
| "# \n", | |
| "# Analysis of proportions: a survey was done of bicycle and \n", | |
| "# other vehicular traffic in the neighborhood of the campus of the \n", | |
| "# University of California, Berkeley, in the spring of 1993. \n", | |
| "# Sixty city blocks were selected at random; each block was observed \n", | |
| "# for one hour, and the numbers of bicycles and other vehicles traveling \n", | |
| "# along that block were recorded. The sampling was stratified into six \n", | |
| "# types of city blocks: busy, fairly busy, and residential streets, with \n", | |
| "# and without bike routes, with ten blocks measured in each stratum. \n", | |
| "# Table 3.3 displays the number of bicycles and other vehicles \n", | |
| "# recorded in the study. For this problem, restrict your attention \n", | |
| "# to the first four rows of the table: the data on residential streets. \n", | |
| "# \n", | |
| "# (a) Let $y_1$ , . . . , $y_{10}$ and $z_1$ , . . . , $z_8$ be the \n", | |
| "# observed proportion of traffic that was on bicycles in the residential \n", | |
| "# streets with bike lanes and with no bike lanes, respectively \n", | |
| "# (so $y_1 = 16/(16 + 58)$ and $z_1 = 12/(12 + 113)$, for example). \n", | |
| "# Set up a model so that the $y_i$ ’s are independent and identically \n", | |
| "# distributed given parameters $\\theta_y$ and the $z_i$ ’s are \n", | |
| "# independent and identically distributed given parameters $\\theta_z$ . \n", | |
| "# \n", | |
| "# (b) Set up a prior distribution that is independent in \n", | |
| "# $\\theta_y$ and $\\theta_z$ . \n", | |
| "# \n", | |
| "# (c) Determine the posterior distribution for the parameters \n", | |
| "# in your model and draw 1000 simulations from the posterior distribution. \n", | |
| "# (Hint: $\\theta_y$ and $\\theta_z$ are independent in the posterior \n", | |
| "# distribution, so they can be simulated independently.) \n", | |
| "# \n", | |
| "# (d) Let $\\mu_y = E(y_i |\\theta_y )$ be the mean of the distribution \n", | |
| "# of the $y_i$ ’s; $\\mu_y$ will be a function of $\\theta_y$. \n", | |
| "# Similarly, define $\\mu_z$ . Using your posterior simulations from (c), \n", | |
| "# plot a histogram of the posterior simulations of $\\mu_y-\\mu_z$, the \n", | |
| "# expected difference in proportions in bicycle traffic on residential \n", | |
| "# streets with and without bike lanes. We return to this example in \n", | |
| "# Exercise 5.13.\n", | |
| "# \n", | |
| "# Gelman, Andrew; Carlin, John B.; Stern, Hal S.; Dunson, David B.; \n", | |
| "# Vehtari, Aki; Rubin, Donald B.. Bayesian Data Analysis, \n", | |
| "# Third Edition (Chapman & Hall/CRC Texts in Statistical Science) (Page 81). \n", | |
| "# CRC Press. Kindle Edition. \n", | |
| "# \n", | |
| "# #### Data\n", | |
| "# |Type |Bike lane? |Counts of Bikes/others|\n", | |
| "# |--- |----------|----|\n", | |
| "# |Residential |yes |16/58, 9/90, 10/48, 13/57, 19/103, 20/57, 18/86, 17/112, 35/273, 55/64 |\n", | |
| "# |Residential |no |12/113, 1/18, 2/14, 4/44, 9/208, 7/67, 9/29, 8/154|\n", | |
| "# \n", | |
| "# Gelman, Andrew; Carlin, John B.; Stern, Hal S.; \n", | |
| "# Dunson, David B.; Vehtari, Aki; Rubin, Donald B.. \n", | |
| "# Bayesian Data Analysis, Third Edition (\n", | |
| "# Chapman & Hall/CRC Texts in Statistical Science) \n", | |
| "# (Page 81). CRC Press. Kindle Edition. \n", | |
| "\n", | |
| "# #### Probably best to first do 3.10.6\n", | |
| "# For that problem see the reference \n", | |
| "# [Raftery, 1988](https://www.stat.washington.edu/raftery/Research/PDF/bka1988.pdf)\n", | |
| "import pystan\n", | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "from scipy.stats import beta\n", | |
| "\n", | |
| "stan_code=\"\"\"\n", | |
| "data {\n", | |
| " int<lower=1> N;\n", | |
| "\n", | |
| " int bikes[N];\n", | |
| " int others[N];\n", | |
| "}\n", | |
| "\n", | |
| "parameters {\n", | |
| " real<lower=0> theta_b;\n", | |
| " real<lower=0> theta_v;\n", | |
| "\n", | |
| "}\n", | |
| "\n", | |
| "model {\n", | |
| " //theta_b~uniform(0,100);\n", | |
| " //theta_v~uniform(50,300);\n", | |
| " \n", | |
| " bikes~poisson(theta_b);\n", | |
| " others~poisson(theta_v);\n", | |
| "}\n", | |
| "\n", | |
| "generated quantities {\n", | |
| " real b_ppc;\n", | |
| " real o_ppc;\n", | |
| " real p ; \n", | |
| "\n", | |
| " o_ppc=poisson_rng(theta_v);\n", | |
| " b_ppc=poisson_rng(theta_b);\n", | |
| "\n", | |
| " p=o_ppc/(o_ppc+b_ppc);\n", | |
| "}\n", | |
| "\"\"\"\n", | |
| "sm=pystan.StanModel(model_code=stan_code)\n", | |
| "fit=sm.sampling(data=dict({'N':10,'bikes':[16,9,10,13,19,20,18,17,35,55],'others':[58, 90, 48, 57, 103, 57, 86, 112, 273, 64] }))\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 69, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "fit=sm.sampling(data=dict({'N':10,'bikes':[16,9,10,13,19,20,18,17,35,55],'others':[58, 90, 48, 57, 103, 57, 86, 112, 273, 64] }))\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 70, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "b_ppc=fit.extract()['b_ppc']\n", | |
| "o_ppc=fit.extract()['o_ppc']\n", | |
| "theta_b=fit.extract()['theta_b']\n", | |
| "theta_v=fit.extract()['theta_v']\n", | |
| "p_lane=fit.extract()['p']\n", | |
| "%matplotlib inline\n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 74, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAELCAYAAAA1AlaNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xm4HVWZ7/HvzyCTyJz2AUJMFOxLVEQ4RJwiLYJBbxNUkGB7JS1907bNvc7d8aqI0e4OoNL0I90SIc3gAIjKTUskIpN9BTSHKSHESIgxHMIjkSBKI8Mh7/1jrRM3m73PqTPs4Zz1+zxPntSuWlX11tmr3r127VWrFBGYmVkZntfpAMzMrH2c9M3MCuKkb2ZWECd9M7OCOOmbmRXESd/MrCBO+mZmBXHSNzMriJO+mVlBtut0APX23nvvmDZtWqfDsAnstttu+01ETG73fl23rZWq1uuuS/rTpk2jt7e302HYBCbpV53Yr+u2tVLVeu3LO2ZmBXHSNzMriJO+mVlBnPTNzAripG9mVhAnfSuapNmS1kpaJ2lBg+WzJN0uqV/SCQ2W7yrpAUlfaU/EZqPjpG/FkjQJOA84FpgBnCxpRl2xjcA84JtNNvN54KZWxWg21pz0rWQzgXURsT4ingIuA+bUFoiIDRGxEthav7Kkw4AXAT9sR7BmY8FJ30q2H3B/zeu+PG9Ikp4HfAn4RAviMmuZrrsjt0TTFly9bXrDord3MJLiqMG8qLjuB4FlEXG/1GgzeQfSfGA+wNSpU4cd4EQyUM9dxzvLSd9K1gfsX/N6CrCp4rqvBd4o6YPALsD2kh6LiGf9GBwRi4HFAD09PVU/UMxaxknfSrYCOFDSdOABYC7wniorRsRfDExLmgf01Cd8s27ka/pWrIjoB04DlgNrgCsiYrWkhZKOA5B0uKQ+4ETgfEmrOxex2ei5pd8htdfxrXMiYhmwrG7e6TXTK0iXfQbbxkXARS0Iz2zMuaVvZlYQJ30zs4I46ZuZFcRJ38ysIP4h18xayp0Wuotb+mZmBXFLv8t4SAYzayW39M3MCuKkb2ZWEF/e6WK+1GMTmet3Z7ilb2ZWECd9M7OCOOmbmRXESd/MrCCVfsiVNBs4F5gEXBARi+qWfxT4K6Af2Ay8PyJ+lZc9A6zKRTdGxHFjFLuZjUO+Q7ezhkz6kiYB5wFHkx4vt0LS0oi4p6bYHaQnBz0u6W+As4CT8rI/RMQhYxy3mZmNQJXLOzOBdRGxPiKeAi4D5tQWiIgbIuLx/PJWhnjohJmZdUaVpL8fcH/N6748r5lTgR/UvN5RUq+kWyUdP4IYzcxsjFS5pq8G86JhQem9QA/wpprZUyNik6SXANdLWhUR99WtNx+YDzB16tRKgY9HvpZpZp1WpaXfB+xf83oKsKm+kKS3AJ8CjouIJwfmR8Sm/P964Ebg1fXrRsTiiOiJiJ7JkycP6wDMzKy6Kkl/BXCgpOmStgfmAktrC0h6NXA+KeE/VDN/D0k75Om9gdcDtT8Am5lZGw2Z9COiHzgNWA6sAa6IiNWSFkoa6H55NrAL8G1Jd0oa+FA4COiVdBdwA7CortePWcdImi1praR1khY0WD5L0u2S+iWdUDP/EEm3SFotaaWkk+rXNetWlfrpR8QyYFndvNNrpt/SZL2bgVeOJkCzVqjYFXkjMA/4eN3qjwPvi4h7Je0L3CZpeUT8tg2hm42KR9m0Um3rigwgaaAr8rakHxEb8rKttStGxC9qpjdJegiYDDjpW9fzMAxWquF2RW5I0kxge+C+ocqadQMnfStV5a7ITTcg7QNcCvxlRGxtUmZ+vk+ld/PmzSMI02xsOelbqSp1RW5G0q7A1cCnI+LWZuXcHdm6jZO+lWrIrsjN5PLfAy6JiG+3MEazMeekb0Wq0hVZ0uGS+oATgfMlrc6rvxuYBczLXZTvlORBBW1ccO8dK1aFrsgraDB4YER8Hfh6ywM0awG39M3MCuKW/jhRO1jbhkVv72AkZjaeuaVvZlYQJ30zs4I46ZuZFcRJ38ysIE76ZmYFcdI3MyuIk76Zjdq0BVf7GdDjhPvpj0Pus29mI+WWvpl1nL8ptI+TvplZQXx5x8zGnFvt3cstfTOzgjjpm5kVxEnfzKwgTvpmZgVx0jczK4h775jZmHGvne7nlr6ZWUGc9K1okmZLWitpnaQFDZbPknS7pH5JJ9QtO0XSvfnfKe2L2mzkKiX9CifGRyXdI2mlpOskvbhmmU8M60qSJgHnAccCM4CTJc2oK7YRmAd8s27dPYHPAq8BZgKflbRHq2M2G60hk37FE+MOoCciDgauBM7K6/rEsG42E1gXEesj4ingMmBObYGI2BARK4Gtdeu+Fbg2IrZExCPAtcDsdgRtNhpVWvpVTowbIuLx/PJWYEqe9olh3Ww/4P6a1315XqvXNeuYKkl/uJX7VOAHw1lX0nxJvZJ6N2/eXCEkszGhBvNiLNd13bZuU6XLZuUTQ9J7gR7gTcNZNyIWA4sBenp6qp5044K7sHW1PmD/mtdTgE3DWPfIunVvrC80keu2jU9VWvqVTgxJbwE+BRwXEU8OZ12zDlkBHChpuqTtgbnA0orrLgeOkbRH/p3qmDzPrKtVSfpDnhiSXg2cT0r4D9Us8olhXSsi+oHTSHVyDXBFRKyWtFDScQCSDpfUB5wInC9pdV53C/B50vmxAliY55l1tSEv70REv6SBE2MSsGTgxAB6I2IpcDawC/BtSQAbI+K4iNgiaeDEgEJODF/SGT8iYhmwrG7e6TXTK/hjx4T6dZcAS1oaoNkYqzQMQ4UT4y2DrOsTw8ysS/iOXDOzgjjpm5kVxEnfzKwgTvpm1jWmLbjaHSFazEnfzKwgTvpmZgXxk7PMbMR8KWb8cUvfzKwgTvpmZgVx0jczK4iv6ZvZsPg6/vjmlr6ZWUGc9M3MCuKkb2ZWECd9M7OC+Ifcca72R7UNi97ewUjMbDxwS9/MrCBO+mZmBXHSNzMriJO+mVlBnPStaJJmS1oraZ2kBQ2W7yDp8rz8p5Km5fnPl3SxpFWS1kj6ZLtjNxsJJ30rlqRJwHnAscAM4GRJM+qKnQo8EhEHAOcAZ+b5JwI7RMQrgcOAvx74QDDrZk76VrKZwLqIWB8RTwGXAXPqyswBLs7TVwJHSRIQwAskbQfsBDwF/K49YZuNnPvpjxEPQjUu7QfcX/O6D3hNszIR0S/pUWAv0gfAHOBBYGfgIxGxpeURm42SW/pWMjWYFxXLzASeAfYFpgMfk/SS5+xAmi+pV1Lv5s2bRxuv2ag56VvJ+oD9a15PATY1K5Mv5ewGbAHeA1wTEU9HxEPAT4Ce+h1ExOKI6ImInsmTJ7fgEMyGx0nfSrYCOFDSdEnbA3OBpXVllgKn5OkTgOsjIoCNwJuVvAA4Avh5m+I2GzEnfStWRPQDpwHLgTXAFRGxWtJCScflYhcCe0laB3wUGOjWeR6wC3A36cPj3yNiZVsPwGwEKv2QK2k2cC4wCbggIhbVLZ8F/DNwMDA3Iq6sWfYMsCq/3BgRx2HWJSJiGbCsbt7pNdNPkLpn1q/3WKP5NrYGOkh4MMGxM2TSr+nLfDTp+uYKSUsj4p6aYhuBecDHG2ziDxFxyBjEamaFcG+41qnS0t/WlxlA0kBf5m1JPyI25GVbWxCjmZmNkSrX9Bv1Zd5vGPvYMXdZu1XS8Y0KuFubmVl7VEn6VfoyD2ZqRPSQurj9s6SXPmdj7tZmZtYWVS7vVOnL3FREbMr/r5d0I/Bq4L5hxGgV+SlaZjaUKi39Kn2ZG5K0h6Qd8vTewOup+S3AzMzaa8ikX6Uvs6TDJfWRurCdL2l1Xv0goFfSXcANwKK6Xj9mZtZGlfrpV+jLvIJ02ad+vZuBV44yRjMzGyO+I9fMrCBO+mZmBXHSNzMriJO+mVlBnPTNzAripG9mVhAnfTOzgjjpm5kVxEnfzKwgle7INTPzg00mBrf0zcwK4qRvZlYQX96ZoDy2fjWSZgPnApOACyJiUd3yHYBLgMOAh4GTah4PejBwPrArsBU4PD9I3axruaVvxZI0CTgPOBaYAZwsaUZdsVOBRyLiAOAc4My87nbA14EPRMTLgSOBp9sUutmIOelbyWYC6yJifUQ8BVwGzKkrMwe4OE9fCRwlScAxwMqIuAsgIh6OiGfaFLfZiDnpW8n2A+6ved2X5zUskx8o9CiwF/AyICQtl3S7pL9rQ7xmo+Zr+lYyNZgXFctsB7wBOBx4HLhO0m0Rcd2zVpbmA/MBpk6dOuqAzUbLLX0rWR+wf83rKcCmZmXydfzdgC15/k0R8ZuIeJz0ZLlD63cQEYsjoicieiZPntyCQzAbHid9K9kK4EBJ0yVtD8wFltaVWQqckqdPAK6PiCA9M/pgSTvnD4M3ARPy+c/TFlztG7MmEF/esWJFRL+k00gJfBKwJCJWS1oI9EbEUuBC4FJJ60gt/Ll53UckfZn0wRHAsohwZmwRd0EeO076VrSIWEa6NFM77/Sa6SeAE5us+3VSt02zccNJfxT8ldfMxhtf0zczK4iTvplZQZz0zWxccW+i0XHSNzMriJO+mVlBKiV9SbMlrZW0TtKCBstn5fFH+iWdULfsFEn35n+n1K9rZmbtM2TSrzj87EZgHvDNunX3BD4LvIY0ouFnJe0x+rDNzGwkqrT0hxx+NiI2RMRK0oMkar0VuDYitkTEI8C1wOwxiNvMzEagStKvMvxsK9Y1M7MxViXpVxl+dlTrSpovqVdS7+bNmytu2sxK5q6bI1Ml6VcZfnZU63r4WTOz9qiS9KsMP9vMcuAYSXvkH3CPyfPMzKwDhkz6+RFxA8PPrgGuGBh+VtJxAJIOl9RHGo3wfEmr87pbgM+TPjhWAAvzPDMz64BKo2xWGH52BenSTaN1lwBLRhGjmZmNEQ+tbGYN+UfSicnDMJiZFcRJvwDu2mZmA5z0zcwK4qRvZlYQJ30zs4K4984w+dr4xCJpNnAuMAm4ICIW1S3fAbgEOAx4GDgpIjbULJ8K3AOcERFfbFfcZiPllr4Vq+Kw4acCj0TEAcA5wJl1y88BftDqWM3GipO+lWzIYcPz64vz9JXAUZIEIOl4YD2wuk3xmo2ak76VrMrQ39vK5CFJHgX2kvQC4O+Bz7UhTrMx46RvJasy9HezMp8DzomIxwbdgYcNty7jH3KtZFWG/h4o0ydpO2A3YAvpEaAnSDoL2B3YKumJiPhK7coRsRhYDNDT01P1ORRmLeOkbyXbNmw48ABp2PD31JVZCpwC3AKcAFwfEQG8caCApDOAx+oTvlk3ctK3YkVEv6SBYcMnAUsGhg0HeiNiKXAhcKmkdaQW/tzORWw2ek76VrQKw4Y/QXpOxGDbOKMlwZm1gH/INTMriJO+mVlBnPTNzAripG9mVhAnfTOzgjjpm5kVxF02K5gowynXHseGRW/vYCRmY2egXrtOV+Okb2bPMlEaOdaYL++YmRXESd/MrCBO+mZmBXHSNzMriJO+mVlBKiV9SbMlrZW0TtKCBst3kHR5Xv5TSdPy/GmS/iDpzvzvq2MbvpmZDceQXTYlTQLOA44mPUVohaSlEXFPTbFTgUci4gBJc4EzgZPysvsi4pAxjtvM7Fl8H0o1VVr6M4F1EbE+Ip4CLgPm1JWZA1ycp68EjpLU6NmiZmbWQVWS/n7A/TWv+/K8hmUioh94FNgrL5su6Q5JN0l6I2Zm1jFV7sht1GKvf8BzszIPAlMj4mFJhwFXSXp5RPzuWStL84H5AFOnTq0QkpmNJd+FW44qLf0+YP+a11OATc3KSNoO2A3YEhFPRsTDABFxG3Af8LL6HUTE4ojoiYieyZMnD/8ozMyskipJfwVwoKTpkrYnPRh6aV2ZpcApefoE4PqICEmT8w/BSHoJcCCwfmxCNzOz4Rry8k5E9Es6DVgOTAKWRMRqSQuB3ohYClwIXCppHbCF9MEAMAtYKKkfeAb4QERsacWBmJnZ0CqNshkRy4BldfNOr5l+AjixwXrfAb4zyhitBdy9LZE0GziX1KC5ICIW1S3fAbgEOAx4GDgpIjZIOhpYBGwPPAV8IiKub2vwZiPgO3KtWDX3oBwLzABOljSjrti2e1CAc0j3oAD8BvjziHgl6dLmpe2J2mx0nPStZCO+ByUi7oiIgQ4Nq4Ed87cCs67mpG8lG+09KAPeBdwREU/W70DSfEm9kno3b948ZoGbjZSTvpVsNPegpIXSy0mXfP660Q7cHdm6jZO+lWzE96Dk11OA7wHvi4j7Wh6tVTZtwdW+4awJPyPXSu7Js+0eFOABUlfj99SVGbgH5RaefQ/K7sDVwCcj4idtjNmGwQ9Nfy4n/TpuHZRjlPegnAYcAHxG0mfyvGMi4qH2HoXZ8DjpW9FGcQ/KF4AvtDxAszHma/pmZgVxS9+sYL6cWR639M3MCuKkb2ZWECd9M7OCOOmbmRXESd/MrCDuvWPPUvDduUVxr51yuaVvZlYQJ30zs4I46ZuZFcRJ38ysIE76ZjbheXz9P3LSNzMriLts4u5rZqXwQ1Wc9G0Q7rM/8biBY768Y2ZWELf0zSY4t+6fq+RvscUmfZ8Iw1PySTJeuY4PTynX+4tK+j4JrASu5zaYStf0Jc2WtFbSOkkLGizfQdLleflPJU2rWfbJPH+tpLeOXejWKQN9nidCcnHdtnoTpW43M2RLX9Ik4DzgaKAPWCFpaUTcU1PsVOCRiDhA0lzgTOAkSTOAucDLgX2BH0l6WUQ8M9YHYjZcrtsG5X0zqnJ5ZyawLiLWA0i6DJgD1J4Yc4Az8vSVwFckKc+/LCKeBH4paV3e3i1jE7512ji/1j+h6nZpyavVGv09x2Edf44qSX8/4P6a133Aa5qViYh+SY8Ce+X5t9atu9+Io63TrJJPhDdmPBqrpNPG969r67ZZq1RJ+mowLyqWqbIukuYD8/PLJyXdXSGupnTmaNYeU3sDv+l0EGOkbcfShvfvTwd21WBZV9ftNhtP9bctsY5h3WxFvC+uUqhK0u8D9q95PQXY1KRMn6TtgN2ALRXXJSIWA4sBJPVGRE+V4Ludj6U7SerNk67bgxhP8Y6nWKGz8VbpvbMCOFDSdEnbk368WlpXZilwSp4+Abg+IiLPn5t7QEwHDgR+Njahm42a67YVZ8iWfr6OeRqwHJgELImI1ZIWAr0RsRS4ELg0/5i1hXTykMtdQfphrB/4W/dusG7hum0lUmq0dA9J8/NX4nHPx9KdOnUs4+1vOJ7iHU+xQmfj7bqkb2ZmreNRNs3MCtLRpC9piaSHaruxSTpD0gOS7sz/3tbJGKuStL+kGyStkbRa0ofy/D0lXSvp3vz/Hp2OdTCDHMe4e18k7SjpZ5LuysfyuTx/eh5S4d48xML2Ldj3R/I+75b0rRxLy/c7UpI+lGNdLenDeV7X1N0muaJhfEr+RWmIjJWSDu2CWE/Mf9utknrqyrd3OI+I6Ng/YBZwKHB3zbwzgI93Mq4RHss+wKF5+oXAL4AZwFnAgjx/AXBmp2Md4XGMu/eF1Jd+lzz9fOCnwBHAFcDcPP+rwN+M8X73A34J7JRfXwHMa/V+RxHvK4C7gZ1JnTt+ROqN1DV1t0muaBgf8DbgB/n9PwL4aRfEehDp/pAbgZ6a+TOAu4AdgOnAfcCkVsbX0ZZ+RPyY1CNi3IuIByPi9jz9e2AN6eSfA1yci10MHN+ZCKsZ5DjGnUgeyy+fn/8F8GbSkArQuvdkO2Cn3Ld/Z+DBNu13JA4Cbo2IxyOiH7gJeAddVHeb5Ipm8c0BLsnv/63A7pL2aU+kjWONiDURsbZB8W3DeUTEL4GB4Txapluv6Z+Wv5Yt6fbLIY0ojcT4alLL8kUR8SCkhAr8SeciG56644Bx+L5ImiTpTuAh4FpSS+q3OblBC4ZPiIgHgC8CG0nJ/lHgtlbvdxTuBmZJ2kvSzqSW8v50f91tFl+j4TW65W9dr+2xdmPS/zfgpcAhpBPmS50NZ3gk7QJ8B/hwRPyu0/GMVIPjGJfvS0Q8ExGHkO6YnUlq1T6n2FjuM38gziF9Xd8XeAFwbKv3O1IRsYY0eui1wDWkyw39g67U3SoNkdEl2h5r1yX9iPh1PlG3Al+jxV91xpKk55MS5Tci4rt59q8Hvlrm/x/qVHxVNTqO8fy+AETEb0nXU48gfd0fuDGx4fAJo/QW4JcRsTkinga+C7yuDfsdsYi4MCIOjYhZpEsT99L9dbdZfJWGyOgSbY+165J+3bW3d5C+enY9SSLdvbkmIr5cs6j2Nv5TgP/b7tiGo9lxjMf3RdJkSbvn6Z1IyXgNcANpSAVozXuyEThC0s7573kU6c7dVu93xCT9Sf5/KvBO4Ft0f91tFt9S4H25F88RwKMDl4G6UPuH82jnr9oNfuX+FulSwdOkT7xTgUuBVcDK/AfZp5MxDuNY3kD6WrYSuDP/extpGN7rSC2n64A9Ox3rCI9j3L0vwMHAHTnmu4HT8/yXkE6sdcC3gR1asO/PAT/P+72U1Duj5fsdRbz/Sfpgugs4Ks/rmrrbJFc0jI90yeQ80u83q6jpLdPBWN+Rp58Efg0sryn/qRzrWuDYVsfnO3LNzArSdZd3zMysdZz0zcwK4qRvZlYQJ30zs4I46ZuZFcRJ38ysIBMu6UvaXdIHa14fKen7w9zGPEn7jnD/N9YPndpk+18Zyfa7UbPhmPOypkPKVlz/VZJukbRK0n9I2rUdx9SNOl23K25/yPo/nkjaQ9L38phTP5P0irrlkyTd0ex9kPRiSdfl9W+UNKVm2VRJP8z1/p481lXLTbikD+wOfHDIUoObRxozxarpBz4WEQeRhjn4W0kz8rK7SXd4/niE619AGj73lcD3gE+04gDGCdft9vs/wJ0RcTDwPuDcuuUfIt3l3cwXSSN+HgwsBP6pZtklwNm53s+kTcNcTMSkvwh4qdKDPs7O83aRdKWkn0v6Rr41HkmHSbpJ0m2SlkvaR9IJQA/wjbyNnSSdLmmF0kMmFg+sP4j3Sro5l282Rs2+kq5RegDEWQMzJf2bpF7VPPQjz1+UWwMrJX2xfmNKDzm5OLccNkh6p6Szcgv5GqXxdBoec57/P/Mx3iXpO0qjLSLpIqUHUtwsaX3++zxLDDIcczQfUrbS+qQxyAc+MK4F3jXYtia4jtVtSQdJ+lnN62mSVjaJ88TcKv6FpDfWlP9PSbfnf6/L8/eR9OMcz90D5ev2vUHSPyp94+uVdGg+pvskfaCm3CfysaysO3euyn+H1ZLm18x/TNI/5Dp/q6QXNTiWGaS7fYmInwPTBsoptdrfTmqYNLNtfdIwHHPyujOA7SLi2rztxyLi8UG2M3Y6dVt1C2+BnsazH15wJGlo2ymkD7lbSEMNPB+4GZicy50ELMnTN/LsBx3sWTN9KfDng+z/RuBreXpWbSw1ZeYB64HdgB2BXwH71+4LmJS3dTCwJ+kW7YE7qHdvsM0zgP+Xj+tVwOPkW7pJLeTjhzjmvWq29QXgf+Xpi0hDBjyPVIHXVfj7bwR2bfB3GfJ2+Pr1c7xz8vRHgd93uo4VXLfvBF6Sp/8e+HST+v+lPP024Ed5emdgxzx9INCbpz8GfKqmzr+wwTY3kB84A5xDGlbjhcBk4KE8/xhgMWkIhucB3wdm1Z1TO5G+ee6VX8fA8ZIeyNLoeP4R+HKenkn6VnpYfn0lcFh+H77f5G/2TeBDefqdeZ97kc7H75MG47sDOJsWPzxl4N/AiH8T3c8iog9AaWz1acBvSU8MujY3biaRxsto5M8k/R2p4u4JrAb+Y5D9fQvSwxQk7Spp90ijPNa6LiIezTHdA7yYNK72u3NrZDvSU6xmkMZEeQK4QNLVpMrSyA8i4mlJq/LxXJPnr8rH/KeDHPMrJH2BdAlhF2B5zXavijS65j1NWkPk4xjVsNJN1n8/8C+STieN+fPUcLc7wbWzbl8BvJv0jeOk/K+RgRFmb8vxQPog+oqkQ4BngJfl+SuAJfmb6FURcWeTbS7N/68iPQ3t98DvJT2hNKjeMfnfHbncLqQPlx8D/1vSO/L8/fP8h0l1aeBcug04usF+FwHn5r/tqrz9fkn/nfSBc5ukI5vEDPDxfNzzciwPkD44tgPeSHpexUbgclJj8MJBtjUmSkn6T9ZMP0M6bgGrI+K1g60oaUfgX0mto/slnUFqnQ+mfkCjRgMcPScmpVH2Pg4cHhGPSLqI1DrqV7pMdBQwFziN9BSmhtuMiK2Sno7cvAC2MvQxXwQcHxF35Qp6ZJNYm339bzSsdGXN1o/0lfqYXOZlpK/T9kftrNuXA9+W9F3Sg8nuHSKmgXgAPkIaaOxVpJb4E7CtYTSL9L5eKunsiLhkkG1u5dnHXFu3/ykizq87xiNJo6u+NiIel3RjzTHWniO1sW6TGx9/mbcl0mMwf0k6D49Telb0jsCukr4eEe+tW38TqYU/0Kh5V0Q8KqkPuCMi1udlV5F+z2p50p+I1/R/T/rqN5S1wGRJr4WUdCS9vME2BirIb/Kb9pxr2g2clLf5BtKwro9WjH1X4L+AR3OL+ti8nV2A3SJiGfBh0oNMRmKwY34h8GBOvn8xnI3mk6HRsNKjXl9/HPL3ecCnSc+WLVVH63ZE3EdKjp8hfQAMx27Ag/kb4/8gfftA0otJLeavkerASB9ivhx4fz4OJO2X685uwCM54f83UmKtTKnH1MAD7P8K+HFE/C4iPhkRUyJiGukD4Pr6hJ/X3zvXXYBPAkvy9ApgD0mT8+s3k77Rt9yES/oR8TDwk/yj0NmDlHuKVMnPlHQX6Xrl6/Lii4Cv5q90T5IeGrIKuIr0Zg3lEUk3kxLUqcOI/S7S18fVpMrxk7zohcD38w9nN5FaTcM2xDF/hvRYxGtJQwIPx+tJJ/Kb8w9yd+YWEJLekVs1rwWulrQ8z99X0rKh1gdOlvSLHNMm4N+Hf+QTQ5fU7cuB95Iu9QzHvwKnSLqVdGnnv/L8I4E7Jd1B+pG+vndMJRHxQ9L181vy5c0rSefNNaRv0SuBzwO3DnPTBwGrJf2c1Aj70BDlkbRQ0nH55ZHA2lyHXwT8Q473GdK3+utyvCK9Fy3noZXNzAoy4Vr6ZmbWXCk/5I45SeeRLkvUOjciir38YBPNmHrjAAAALElEQVSD6/bE5ss7ZmYF8eUdM7OCOOmbmRXESd/MrCBO+mZmBXHSNzMryP8HfInw3aYlIMUAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 2 Axes>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XucXHV9//HXZ3eTEJLNzWzInQ1yEUSqNCD+sEoFK0UE688LPgTFH8rD+mvV1mqD8quotdrWWlurtaki1iKi2AqCqFRFULktyCUh3HPPXrObbHaXXHbn8/vjfHeZLDs7szPnMjP7fj4e+8iZOWfO+ZyZyXu/+z3nfI+5OyIiUvsasi5ARETioUAXEakTCnQRkTqhQBcRqRMKdBGROqFAFxGpEwr0acDMNprZWSlv08zsG2bWZ2b3TjD/UjP71SSvv9XM3pVQbWeZ2Y4k1l3Ctt9hZj+dZH6WtbmZHVtg3mF1T7bsFLd5u5m9p9L1SKQp6wJqkZltAY4CRoBB4EfAn7r7QJZ1AZjZNcAOd79y9Dl3f3EGpbwSeC2w0t0Hp/pid//D+EvKnrtfC1w7+tjMHDjO3Z/Krqrixtct1Ukt9PK9wd3nAqcCpwFXjl8gtFJTe4/NrDGtbZXgaGBLOWEuIuVRoFfI3XcCtwInw9ifkJ8xs18DQ8AxZrbczG4ys14ze8rM3jv6ejO7ysxuMLPrzWyfmT1gZr+TN//EsM49oevkgrx515jZv5rZj8xsELgMeAfwUTMbMLMfhuW2mNk5YXqWmX3RzHaFny+a2aww7ywz22FmHzazLjNrN7N3F9r3QvtlZpcBXwNeEer4ZOFV2JfMbK+ZPWZmZ+fNOOxPcTN7r5ltCu/Ro2Z2qpl9xMy+P26FXzKzL4bpRaHbZ1fo+vnBJPvxfTPrNrPNZvaBvHmnm1mbmfWbWaeZfaHAOn5pZv87TL8ydEmcFx6fY2YPhumxriYzuyO8/KHwPr0tb32lfga3m9mnzezX4b35qZktzpt/Qfje7AnLnlhoXcF5ZvaMmfWY2d+PNkhski6ysL/bzez3w+MXmdlt4XvxuJm9tcg2R9fzQjP7uZntDtu/1swW5M3fYmZ/YWYPh+/M9WZ2RN78883swbCvvzGzU0rZbl1xd/1M8QfYApwTplcBG4FPh8e3A9uAFxN1ac0Afgl8BTgCeCnQDZwdlr8KOAS8OSz7F8DmMD0DeAr4GDATeA2wDzghvPYaYC9wJtEv5yPCc389Sb2fAu4GlgAtwG/yaj8LGA7LzADOI/qltLDA+zDZfl0K/GqS9/DSsK0/C9t6W9iXRXnv43vC9FuAnUR/CRlwLNFfAMuIurwWhOWagC7gd8PjW4DrgYVhG6/O288dYboBuB/4q/AeHwM8A7wuzL8LuCRMzwXOKLA/nwK+FKY/BjwN/G3evH+a6H0BHDg27/FUP4Pbw7aOB2aHx58L844P789rw7o+SvR9mllgXQ78AlgErAaeyPsMJqwbeB2wHTg9PD8nPH53+DxOBXqAF09S/+g2jg21ziL6bt4BfHHc9/heYHmocRPwvjDv1PDZvxxoBN4Vlp+VdV6kmk1ZF1CLP+GLMgDsAbYShdrsMO924FN5y64i6mtvznvus8A1Yfoq4O68eQ1AO/B74acDaMibfx1wVZi+BviPcbVdw+SB/jRwXt681xF1jUAUJs8CTXnzu5ggxErYr8MCYILXXwrsAizvuXt5Ljzz/6P/BPhggfXcCrw3TJ8PPBqmlwE5JghCDg/0lwPbxs2/AvhGmL4D+CSwuMh34mzg4TD9Y+A9o58r0S++N030vjBxoJf0GeS9T1fmPX4/8OMw/f+A7477bu0EziqwLgfOHbeun01S9xVE3/+X5D3/NuDOcev9N+ATk9T/ngLz3gj8dtz3+OK8x38HfDVM/yuhYZI3/3HCL/Hp8qMul/K90d0XuPvR7v5+d382b972vOnlQK+778t7biuwYqLl3T0H7AivWw5sD88VfW2Jlod15K9ved7j3e4+nPd4iKhlOtF6iu1XMTs9/M8rUMuoVUS/iCbyTeDiMH0x8K281/S6e1+RGo4Gloc/0/eY2R6iFvZRYf5lRC3dx8zsPjM7v8B67gKON7OjiP5a+Q9gVej+OJ3oF0OpSv0MRnUUWPawzzp8j7Yz+WeU/30q9HmM+hDRL4xH8p47Gnj5uPfzHcDSSdYDgJktMbPvmNlOM+sH/hNYPG6xQvt6NPDhcdtdVaT+uqNAT0Z+SO0CFplZc95zq4laSqNWjU6EPsuV4XW7iEKhYZLXjh8us9jwmbuIvvz569tV5DWF1lNsv4pZYWZWQi3bgRcWWMcPgFPM7GSiFvq1ea9ZlN8HW8B2YHP45Tz60+zu5wG4+5Pu/naiLqq/BW4wsznjV+LuQ0RdNx8ENrj7QaLurD8Hnnb3niJ1JOGwzzq816uY/DNalTdd7LvxFuCNZvahvOe2A78c937Odfc/LqHezxJ9f09x93lEv6Bt8pcctt3PjNvuke5+XYmvrwsK9IS5+3ai/9ifNbMjwoGayzj8FLDfNbM3mVkTUavnAFE/9z1EfaAfNbMZFp1L/gbgO5NsspOoH7iQ64ArzawltB7/iqgllMR+FbME+EDYt7cAJxKdAjre14C/MLPftcixZnZ0qGM/cAPwbeBed98Wnm8n6o75ipktDNt41QTrvhfoN7O/NLPZZtZoZieb2WkAZnaxmbWE1u2e8JqRAvvzS+BPwr8QdSfkP55Isc+rEt8FXm9mZ5vZDODDRN+t30zymo+E92sV0S+n6ydZdhdRV9MHzOz94bmbif5SuSS85zPM7LQSDsYCNBO6Ms1sBfCREl4z6t+B95nZy8N3ZI6ZvX5cg6PuKdDT8Xagleg/wH8T9Sfeljf/RqK+xz7gEqL+1kOhlXcB8IdEB5a+ArzT3R+bZFtfB04Kf3ZOdFbHXwNtwMPAI8AD4bkk9quYe4DjiPbtM8Cb3X33+IXc/Xth/reJDgr/gOig2KhvAi/hue6WUZcQHXB+jKgf+kPj5uPuI0S/JF9KdDC6h+gXyPywyLnARjMbAP4JuCj8EpnIL4lC6Y4CjydyFfDN8HmVdDZIqdz9caJW7peI9usNRKfbHpzkZTcS/aXxINFB5a8X2cY2olD/SzN7T+iC+wPgIqLvRQfRXzazSij5k0QHN/eGbf9XCa8ZraMNeC/wL0T/j54i6vefVuzwLkxJm5ldRXRQ7OJiy8rEzGw1UWgvdff+rOsRyYpa6FLTwvGFPwe+ozCX6U6X/kvNCgcnO4nOxjg343JEMqcuFxGROqEuFxGROpFql8vixYu9tbU1zU2KiNS8+++/v8fdW4otVzTQzexqogs2utx9dACqvyc6Beog0RV873b3PYXXEmltbaWtra3YYiIiksfMthZfqrQul2t4/gGn24CT3f0UogF8rphSdSIiEruige7udwC94577ad5YE3cTXaouIiIZiuOg6P8husRaREQyVFGgm9nHicZuLjh+h5ldbtENAtq6u7sr2ZyIiEyi7EC36Aa+5wPv8ElOZnf39e6+1t3XtrQUPUgrIiJlKuu0RTM7F/hLosHjh+ItSUREylG0hW5m1xEN3n+CRfebvIxoRLNm4LZwD7+vJlyniIgUUbSFHgb3H2/SITVFRCR9uvRfRKROKNAz0LruFlrX3ZJ1GSJSZxToIiJ1QoEuIlInFOgiInVCgS4iUicU6CIidUKBnrXcSNYViEidUKBn5Ej2w7f+CP5mOTz83azLEZE6kOot6OQ5/7fpB/D0z2HhGrjxTzjj2/vp4AVs+dzrsy5NRGqUWugZmMVBLmm8DV78R/DOGyF3iHc1/TTrskSkxinQM3BWw4PMs2fh1HfCwqPhuD/gDY13AQVHIRYRKUqBnoFXNzxEvx8Jra+KnjjxDay0Hl5c2n1gRUQmpEDPwOkNj3Nf7gRoDIcwjj0HgDMaNmZYlYjUOgV62ga6ObZhF/fmXvTcc81L2ZI7itMansiuLhGpeQr0tG27CyBqoedp8xNY2/A4FL6bn4jIpBToaWt/iGFvYKO3HvZ0W+54Fls/7H46m7pEpOYp0NPW9SjP+DIOMPOwMdEfzh0TTXQ8lFFhIlLrFOhp69zI477qeU8/5Ss45I3QqQOjIlIeBXqaDuyDPVt5LLf6ebMOMoOnfTl0bMigMBGpBwr0NHU9BjBhCx1gk6+GTgW6iJRHgZ6mnui0xCd9xYSzH8uthv6dMNSbZlUiUicU6Gnq2wzWyE5fPOHssaDf/VSKRYlIvVCgp6l3M8xfyXCBQS43+7JoQqcuikgZFOhp6tsMi9YUnL3dWxj2BrXQRaQsCvQ09W6Oxj8vYJgmtvkSBbqIlEWBnpZn98CzvZO20CF0u6jLRUTKUPSORWZ2NXA+0OXuJ4fnFgHXA63AFuCt7t6XXJl1oG9z9O+4Fnr+1aIAW3wp9P4yGtPFLK3qRKQOlNJCvwY4d9xz64CfuftxwM/CY5lMbwj0oi30pXBoCPa1p1CUiNSTooHu7ncA40+MvhD4Zpj+JvDGmOuqP2Mt9NZJF3tm9EyXnieTrUdE6k65fehHuXs7QPh3SaEFzexyM2szs7bu7u4yN1cH9myH2YtgVvOki23NHRWW192LRGRqEj8o6u7r3X2tu69taWlJenPVq38XzJ/4CtF8HSwCa4Q921IoSkTqSbmB3mlmywDCv13xlVSn+nfCvJVFFxuhEeatUKCLyJSVG+g3Ae8K0+8CboynnDrWv7OkFjoAC1ZFXTQiIlNQNNDN7DrgLuAEM9thZpcBnwNea2ZPAq8Nj6WQg0PwbB9/95t9pS2/YLVa6CIyZUXPQ3f3txeYdXbMtdSv/p0A7PIXlLb8/FWwbxeMHILGGQkWJiL1RFeKpiEEegeLSlt+wWrw3NjrRERKoUBPw94pttAXhDsaqdtFRKZAgZ6G/l0AdPrC0pZfEO5opEAXkSlQoKehfwc9Po8DzCxt+XkrAeOL3/9ZomWJSH1RoKehfxcdXmL/OUDTTGhexkrrSa4mEak7CvQ07Gsvvbtl1ILVrLRpPFSCiEyZAj0NA110+YKpvWbBKlagFrqIlE6BnrTcCAx208UUA33ecpZYH+RyydQlInVHgZ60wR7wHN1TbaHPW8EsG4ah3cnUJSJ1R4GetIFOALp9/tRe1xzGRdfFRSJSIgV6wi790g8BymqhA7pzkYiUTIGesBbbA1BWHzqgFrqIlEyBnrAW9gLQM9Uul7lLGPaGsatMRUSKUaAnrMX20O+z2c+sqb2woTFq1fery0VESqNAT1iL7Zl6/3nQ4YvU5SIiJVOgJ6zF9tI91f7zoN0XqctFREqmQE/YEvqmfspi0Dka6O4xVyUi9UiBnrAW21t2l0u7L4JDg3CgP+aqRKQeKdCTdGCAubZ/6uO4BGMjNKrbRURKoEBP0mAXUMZFRYECXUSmQoGepH3hsn/K60NvR4EuIqVToCdpMBrPPP+iotZ1t5T88q7RMdQV6CJSAgV6kkKg7/Z5Zb38IDPgyMWwT4EuIsUp0JMUhr7to7n8dcxbrha6iJREgZ6kwR76fXbU0i5X8zKNuCgiJVGgJ2moh94yu1vGzFum8VxEpCQVBbqZ/ZmZbTSzDWZ2nZkdEVdhdWGwh95KulsgaqEP9cDwwXhqEpG6VXagm9kK4APAWnc/GWgELoqrsLowtLvsA6Jjwp2Lzrzy21M6Q0ZEpp9Ku1yagNlm1gQcCejoXb7BnooD/dLvbwfgKOuLoyIRqWNlB7q77wQ+D2wD2oG97v7TuAqree4wtJteKgv0znC1qAJdRIqppMtlIXAhsAZYDswxs4snWO5yM2szs7bu7u7yK601+/dC7hC7vbI+9I5wcZECXUSKqaTL5Rxgs7t3u/sh4L+A/zV+IXdf7+5r3X1tS0tLBZurMeEc9ErPcumjmQPexFIFuogUUUmgbwPOMLMjzcyAs4FN8ZRVB8JVopV2uYDR5QtZokAXkSIq6UO/B7gBeAB4JKxrfUx11b7BHoCKu1wAOlnIUnorXo+I1LemSl7s7p8APhFTLfVlKAr0ii8sIupHP9G2VbweEalvulI0KaMt9Iq7XKJRF3VQVESKUaAnZWg3zJjDAWZWvKoOX8hc289chmIoTETqlQI9KYM9MOcFsayqU6cuikgJKupDl0kM9cCc8k7THH+Jfye6uEhEilMLPSmDPdHNKWIw1kJHgS4ihSnQkzK0G+bEG+i6uEhEJqNAT4J7aKHH04c+xBH0+2xdXCQik1KgJ+HAPhg5EFsLHaJBupaaLi4SkcIU6EkIFxXF1YcO0OkLdFBURCalQE/CYDQwV6wtdBYp0EVkUgr0JCTSQl/IEvZALhfbOkWkvijQkxAu+4/rwiKIrhadYSPP/bIQERlHgZ6EhFroAOxrj22dIlJfFOhJGOyBptm0fuL22FY5eis6+hXoIjIxBXoSxi4qsthW2aEWuogUoUBPQowXFY3qYT45NwW6iBSkQE/CUE+spywCDNNED/MV6CJSkAI9CTEOzJWv0xeoD11EClKgJ2Ew/hY6QIcvUgtdRApSoMft4CAMPxt7HzpEt6JToItIIQr0uI1dVJREC31hdAbN8IHY1y0itU+BHrcELioa1cnoqYsdsa9bRGqfAj1uYwNzlXf7ucmMXVykbhcRmYACPW5D8Y/jMkqX/4vIZBTocRtMrstl7GpRnbooIhNQoMdtqAcaZ8Ks5thXvYe50DhLLXQRmZACPW6Du2kfnkvrFT9KYOUGzUsV6CIyoYoC3cwWmNkNZvaYmW0ys1fEVVjNGuym1+NvnY9pXqazXERkQpW20P8J+LG7vwj4HWBT5SXVuKEedvu85NY/bxn070pu/SJSs8oOdDObB7wK+DqAux909z1xFVazBnvoJYUWunty2xCRmlRJC/0YoBv4hpn91sy+ZmZzxi9kZpebWZuZtXV3d1ewuRoxtJveBFvon7lzDxwahAP9iW1DRGpTJYHeBJwK/Ku7vwwYBNaNX8jd17v7Wndf29IS/8U2VeXQfjg4QE+Cgf7cuegdtK67JbHtiEjtqSTQdwA73P2e8PgGooCfvsJFRb0kGeijt6JTP7qIHK7sQHf3DmC7mZ0QnjobeDSWqmpVuKgoybNcOtDVoiIysaYKX/+nwLVmNhN4Bnh35SXVsBDoSZ7lcvjl//MT246I1J6KAt3dHwTWxlRL7Uuhy2U/s+CI+eHy/xclth0RqT26UjROYy30BE9bhHDqorpcRORwCvQ4DfVAQxP9PO/szXgp0EVkAgr0OI3dHNqS3Y4u/xeRCSjQ4zS0O5Fbzz3PvCjQG8glvy0RqRkK9DgNdidyc+jnaV4GPsIL0NWiIvIcBXqcBnvSaaE3LwPgKOtNflsiUjMU6HEa2p3IvUSfZyzQ+5LflojUDAV6XIYPwIF+Pv+r3clva14U6EsV6CKSR4Eel6EoyJO8qGjMnCVgDSxRoItIHgV6XAajoYETv6gIoLEJ5ixhKQp0EXmOAj0uYwNzpdBCB2heqj50ETmMAj0uY10uKbTQAeYt11kuInIYBXpcQgu9x1MaAVEtdBEZR4Eel6EesEb6OTKd7TUvZ5ENRHdJEhFBgR6fcJWop/WWNi+N/tUgXSISKNDjMpjSOC6jwrnoGqRLREYp0OMy1JPOOC6jmpdH/+7TvUVFJFLpLehk1GAPLDsllU21rruF+Qzw0BGohS4iY9RCj8tQTzrjuAR7mcN+nwH9aqGLSESBHofhg7B/b7i5RVosumG0WugiEijQ4xAuKrrytnTPOOlgkc5yEZExCvQ4DI3eHDqly/6DLl+gQBeRMQr0OKQ9jkvQ4Yugvx3cU92uiFQnBXocQpfL7rTGcQk6fSEMPxv134vItKdAj8NAFwDdviDVzXb6wmhC3S4iggI9HgOdHPRG9jIn1c2OBbpOXRQRFOjxGOymh/mApbrZnR5Ok9y7I9Xtikh1qjjQzazRzH5rZjfHUVBNGuhMvbsFwmmL1gh7tqW+bRGpPnG00D8IbIphPbVroDO9cdDzjNAI81co0EUEqDDQzWwl8Hrga/GUU6MGuunOINABWHC0Al1EgMpb6F8EPgrkCi1gZpebWZuZtXV3d1e4uSqUG4HBbrpJv8sFgAWrFegiAlQQ6GZ2PtDl7vdPtpy7r3f3te6+tqUlvcGrUjPUCz6SWQv9H+/bT66/HYYPZLJ9EakelbTQzwQuMLMtwHeA15jZf8ZSVS0Z6ARSvJfoODu8hQZznekiIuUHurtf4e4r3b0VuAj4ubtfHFtltWIwm4uKRm338FfPnq2ZbF9EqofOQ6/U6FWiZNdCB9SPLiLx3LHI3W8Hbo9jXTUn4y6XThZyyBuZoUAXmfbUQq/UQBc0zWaA2ZlsfoRG2n2RWugionuKVmygC+YugYF0L/vPt8Nb6HroId583y0AbPnc6zOrRUSyoxZ6pQY6Ye5RmZaww1tYaXV4jr+ITIkCvVKD3VELPUPbvYWl1sdMDmVah4hkS4FeqYHOzAN99EwXtdJFpjcFeiVGDkV3K8q4y2WrR9tvtY5M6xCRbCnQKxHuJfrx2zozLWOzLwVgjQJdZFpToFdiIArQzEZaDPpoZq8fqRa6yDSnQK9Ef3Qvzw5flHEhxmZfqkAXmeYU6JXYF93LM/tAhy2+lDUNCnSR6UyBXon+doa9IdxPNFtbfCnL2c0sDmZdiohkRIFeif5ddLGAXBW8jZtzS2kwZ5V1ZV2KiGQk+ySqZft2VUV3C0QtdNCZLiLTmQK9Ev3tVRPoo6cu6sCoyPSlQK9E/y46fWHWVQDQz1x6fa5a6CLTmAK9XPv74eC+aOjaKrFFpy6KTGsK9HLtq5Zz0J+zWacuikxrCvRy9UfnoHdWUaA/k1vOMuuFA/uyLkVEMqBAL1doobdTPYH+pK8A4I1XXU3rulsyrkZE0qZAL1f/ToCqOSgK8ISvBOC4hh0ZVyIiWVCgl6u/nV6fywFmZl3JmG1+FPt9BsebAl1kOlKgl2tfe1X1nwPkaOApX6FAF5mmdJPoKRrtm7515gZ2+OKMq3m+J3wlr2h4NOsyRCQDaqGXxVlh3WO3fqsmT+ZWssx6mcdg1qWISMoU6GWYxyDz7NkqbaFHZ7ocazszrkRE0qZAL8Mqi249V40t9NEzXY7XmS4i007ZgW5mq8zsF2a2ycw2mtkH4yysmq20bqA6A32Ht7DPZ3Oibc26FBFJWSUHRYeBD7v7A2bWDNxvZre5e90fkVsZxhyvxkB3GnjUj+YlDZuzLkVEUlZ2C93d2939gTC9D9gErIirsGq20nrY57PZy5ysS5nQhtwaTrRtMDKcdSkikqJY+tDNrBV4GXDPBPMuN7M2M2vr7u6OY3OZW2nd4YCoZV3KhB7JrWG2HYSeJ7IuRURSVHGgm9lc4PvAh9y9f/x8d1/v7mvdfW1LS/V1UZQjCvQlWZdR0AZvjSbaH8q0DhFJV0WBbmYziML8Wnf/r3hKqnbOqrEWenV6xpcz5LOg/cGsSxGRFFVylosBXwc2ufsX4iupurWwh7m2f+yWb9UoFw6MqoUuMr1U0kI/E7gEeI2ZPRh+zouprqp1TLgj0DO+PONKJvdIbg2DWx/ghetuyroUEUlJ2actuvuvqNajggk6piG6scXmXPW20AF+mzuWdzf9hBfZ9qxLEZGU6ErRKVpjHRzwGeziBVmXMqm23AkArG14PONKRCQtCvQpWmPtbPaleJW/dbtYzE5/Aacp0EWmjepOpSp0TAj0WtCWOyFqobtnXYqIpECBPhUjw6y2Lp7xZVlXUpL7ciew1Pqgb0vWpYhIChToU9G3hRk2wuYaCnQAtt2dbSEikgoF+lR0bQTg8dyqjAspzRO+kj6fC5vvyLoUEUmBAr2I1nW3jN12jo4NjLiNjTle7ZwGfpU7GZ76H8jlsi5HRBKmQJ+Kzg1s9mUcYGbWlZTs9pGXwmAXdD6SdSkikjAF+lR0bGCTr866iim5I3dKNPHkbdkWIiKJU6CXav9e2LuNTbmjs65kSrpZwIZca9TtIiJ1TYFeqs7ogOijNdZCB7g99zsMb72Hl677TtaliEiCFOilCiMX1loLHeDWkdNpshznNt6XdSkikiAFeqm23wvzVtLJoqwrmbKN3sozuaWc33BX1qWISIIU6CXaseEObu6rjdMVn8+4OXcGr2h4FPZ1Zl2MiCREgV6Co+hlpfXwQO64rEsp240jZ9JoDg+rH12kXinQSzA6YuH9NRzoT/sK7sm9iC0/+bIuMhKpUwr0Eryy4RH2+pFs8DVZl1KRa4fPprWhk0uu/Nvnrn4VkbqhQC/KeWXjBn6TezEjNGZdTEV+nDudbp/P5Y03Z12KiCRAgV7EMdbOSuvhV7mXZF1KxQ4yg/XDr+f3GjfwMnsy63JEJGYK9CLObbgXgJ+NvCzjSuJx7cg57PZmPtz0Xd34QqTOKNCLOL/xHtpyx9NR5fcQLdUQR/DPw2/ilY0bYdMPsy5HRGKkQJ9M1yZOatjKj0ZennUlsfrPkXPYlFsFP74Cnt2TdTkiEhMF+iS+9c9XcsBn8N8jZ2ZdSqxGaGTdoffCQAfc/CF1vYjUCQV6IYO7eVPjnfww9wr6mJd1NbF7yI+F3/8YbPxvuPMfsi5HRGIw7QP9sDsShccA3Pl5juAgXx0+P6PKUnDmn8FL3go//zTc9ZWsqxGRCk37QJ/Qzgfg3vV8d+QsnqqR282VpaEBLvwyt46cBj+5An70ETi0P+uqRKRMTVkXUG2W0AffuwLmLuWzXW/PupxEjf410sAH+Zhfy3vuXQ9P/wJe9zdw3GvBLOMKRWQqKmqhm9m5Zva4mT1lZuviKiozO9r43sxPMtDXyYXdl9PP3KwrSkWOBv56+BIuObgORg7Ct98CX/09+M2/QN8WHTQVqRFlt9DNrBH4MvBaYAdwn5nd5O6PxlVcWdzDTw6YYJrw2B0ODnCqPcFJDVvhW1+Dp3/ODFvEJQeviA4aTjN35k7huI6TePIte6Dtavjpx6Of5uWw4lRYtAYWroHmpXDEAl73bw/T73O468o/hMYZ0U9D+Fete5HUmZfZ+jKzVwBXufvrwuMrANz9s4Ves3btWm9ra5v6xn78MbggOJtBAAAEZUlEQVT/msLBPDpdgW25Fr438mq+MXIuAxxZ0brqxTG2izMbNnBaw+NccNRu6NsKIwdKe7E15oV6+LfsxyJ14G3fghe+pqyXmtn97r622HKV9KGvALbnPd4BPO8KHDO7HLg8PBwws8fL3N5ioKfM15agH3gauDq5TUxdwvs8ua3AL9LfbKb7nBHt83Tw8bMr2eeS7n1ZSaBP1Hx6XnPf3dcD6yvYTrQxs7ZSfkPVE+3z9KB9nh7S2OdKDoruAFblPV4J7KqsHBERKVclgX4fcJyZrTGzmcBFwE3xlCUiIlNVdpeLuw+b2Z8APwEagavdfWNslT1fxd02NUj7PD1on6eHxPe57LNcRESkuujSfxGROqFAFxGpE1UX6MWGEzCzWWZ2fZh/j5m1pl9lvErY5z83s0fN7GEz+5mZlXROajUrddgIM3uzmbmZ1fQpbqXsr5m9NXzOG83s22nXGLcSvterzewXZvbb8N0+L4s642RmV5tZl5ltKDDfzOyfw3vysJmdGmsB7l41P0QHV58GjgFmAg8BJ41b5v3AV8P0RcD1Wdedwj7/PnBkmP7j6bDPYblm4A7gbmBt1nUn/BkfB/wWWBgeL8m67hT2eT3wx2H6JGBL1nXHsN+vAk4FNhSYfx5wK9F1PGcA98S5/WproZ8OPOXuz7j7QeA7wIXjlrkQ+GaYvgE426ymrxEvus/u/gt3HwoP7yY657+WlfI5A3wa+Dug1sf0LWV/3wt82d37ANy9K+Ua41bKPjuM3T1mPnVwHYu73wH0TrLIhcB/eORuYIGZLYtr+9UW6BMNJ7Ci0DLuPgzshZq+g3Mp+5zvMqLf8LWs6D6b2cuAVe5+c5qFJaSUz/h44Hgz+7WZ3W1m56ZWXTJK2eergIvNbAfwI+BP0yktU1P9/z4l1TYeeinDCZQ05EANKXl/zOxiYC3w6kQrSt6k+2xmDcA/ApemVVDCSvmMm4i6Xc4i+gvsTjM72d1r9S7epezz24Fr3P0fwmB/3wr7XNlIe9Ut0fyqthZ6KcMJjC1jZk1Ef6pN9idOtStpCAUzOwf4OHCBu5c45GHVKrbPzcDJwO1mtoWor/GmGj4wWur3+kZ3P+Tum4HHiQK+VpWyz5cB3wVw97uAI4gG7apniQ6ZUm2BXspwAjcB7wrTbwZ+7uFoQ40qus+h++HfiMK81vtWocg+u/ted1/s7q3u3kp03OACdy9j7OWqUMr3+gdEB78xs8VEXTDPpFplvErZ523A2QBmdiJRoHenWmX6bgLeGc52OQPY6+7tsa0966PCBY4CP0F0hPzj4blPEf2HhuhD/x7wFHAvcEzWNaewz/8DdAIPhp+bsq456X0et+zt1PBZLiV+xgZ8AXgUeAS4KOuaU9jnk4BfE50B8yDwB1nXHMM+Xwe0A4eIWuOXAe8D3pf3OX85vCePxP291qX/IiJ1otq6XEREpEwKdBGROqFAFxGpEwp0EZE6oUAXEakTCnQRkTqhQBcRqRP/H8soW8eH02OOAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "fig,ax=plt.subplots(1,2)\n", | |
| "ax[0].set_xlim(15,30)\n", | |
| "ax[1].set_xlim(80,115)\n", | |
| "j=ax[0].hist(theta_b,density=True,bins=40)\n", | |
| "j=ax[1].hist(theta_v,density=True,bins=40)\n", | |
| "ax[1].set_xlabel('theta_v has mean '+str(np.round(np.mean(theta_v),2)))\n", | |
| "ax[0].set_xlabel('theta_b has mean '+str(np.round(np.mean(theta_b),2)))\n", | |
| "plt.show()\n", | |
| "fig,ax=plt.subplots(1)\n", | |
| "s=ax.hist(1-p_lane,bins=50,density=True)\n", | |
| "ax.plot(np.linspace(0,1,1000),beta.pdf(np.linspace(0,1,1000),b=94.91,a=21.3))\n", | |
| "ax.set_title('Proportion of bicycles with no bike lane')\n", | |
| "plt.show()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 60, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "bikes=[12,1,2,4,9,7,9,8]\n", | |
| "other=[113,18,14,44,208,67,29,154]\n", | |
| "fit=sm.sampling(data=dict({'N':8,'bikes':bikes,'others':other }))\n", | |
| "b_ppc=fit.extract()['b_ppc']\n", | |
| "o_ppc=fit.extract()['o_ppc']\n", | |
| "theta_b=fit.extract()['theta_b']\n", | |
| "theta_v=fit.extract()['theta_v']\n", | |
| "p_nolane=fit.extract()['p']\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 76, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAELCAYAAAA1AlaNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAH8JJREFUeJzt3XuYHVWZ7/Hvj0RARO4ZD7eQKDhDVERooniJHBEm6AzxEiQ4HMnInIzjcMa74lERo84EUBk8MiMRMtxURFSmRwMRuTlHQdPcAiEiIcbQhEciwSiDXELe+WOths1md3d19751r9/nefpJ7apVtd/aWfXutatWrVJEYGZmZdiq0wGYmVn7OOmbmRXESd/MrCBO+mZmBXHSNzMriJO+mVlBnPTNzAripG9mVhAnfTOzgkzudAD1dtttt5g2bVqnw7AJ7KabbvptRExp9/u6blsrVa3XXZf0p02bRl9fX6fDsAlM0q878b6u29ZKVeu1T++YmRXESd/MrCBO+mZmBXHSNzMriJO+mVlBnPStaJJmS7pL0mpJJzdYPkvSzZI2S5rbYPkOku6T9JX2RGw2Nk76VixJk4CzgaOAGcBxkmbUFVsHzAe+MchmPgtc36oYzZrNSd9KNhNYHRFrIuJx4BJgTm2BiFgbESuALfUrSzoYeAHww3YEa9YMTvpWsj2Be2te9+d5w5K0FfBF4CMtiMusZbrujlxLpp38g2e8XrvozR2KZEJTg3lRcd33Aksj4l6p0WbyG0gLgAUAU6dOHXGAE53refs56VvJ+oG9a17vBayvuO6hwOskvRfYHtha0sMR8YyLwRGxGFgM0NPTU/ULxaxlnPStZMuB/SRNB+4D5gHvrLJiRPzVwLSk+UBPfcI360Y+p2/FiojNwEnAMmAVcGlErJS0UNLRAJIOkdQPHAOcI2ll5yI2Gzu39LtI/flNa72IWAosrZt3Ss30ctJpn6G2cT5wfgvCM2s6t/TNzAripG9mVhAnfTOzgjjpm5kVxBdyzayt3GGhs9zSNzMriFv640Rt68i3qpvZaLmlb2ZWECd9M7OC+PTOOOSRCa0Eruet4Za+mVlBnPTNzAripG9mVhAnfTOzglS6kCtpNnAWMAk4NyIW1S3/IPA3wGZgA/DuiPh1XvYkcHsuui4ijm5S7GY2wfhu3dYbNulLmgScDRxBerzcckm9EXFnTbFbSE8OekTS3wGnA8fmZX+MiAObHLeZmY1CldM7M4HVEbEmIh4HLgHm1BaIiGsj4pH88kaGeeiEmZl1RpWkvydwb83r/jxvMCcCV9S83lZSn6QbJb1lFDGamVmTVDmnrwbzomFB6XigB3h9zeypEbFe0guBayTdHhH31K23AFgAMHXq1EqBTwQ+f2lm7Valpd8P7F3zei9gfX0hSW8EPgEcHRGPDcyPiPX53zXAdcAr6teNiMUR0RMRPVOmTBnRDpiZWXVVkv5yYD9J0yVtDcwDemsLSHoFcA4p4T9QM39nSdvk6d2A1wC1F4DNzKyNhk36EbEZOAlYBqwCLo2IlZIWShrofnkGsD3wbUm3Shr4Utgf6JN0G3AtsKiu149Zx0iaLekuSaslndxg+SxJN0vaLGluzfwDJd0gaaWkFZKOrV/XrFtV6qcfEUuBpXXzTqmZfuMg6/0UeNlYAjRrhYpdkdcB84EP163+CPCuiLhb0h7ATZKWRcTv2hC62Zh4lE0r1VNdkQEkDXRFfirpR8TavGxL7YoR8cua6fWSHgCmAE761vU8DIOVaqRdkRuSNBPYGrhnuLJm3cBJ30pVuSvyoBuQdgcuAv46IrYMUmZBvk+lb8OGDaMI06y5nPStVJW6Ig9G0g7AD4BPRsSNg5Vzd2TrNk76VqphuyIPJpf/HnBhRHy7hTGaNZ2TvhWpSldkSYdI6geOAc6RtDKv/g5gFjA/d1G+VZIHFbRxwb13rFgVuiIvp8HggRFxMXBxywM0awG39M3MCuKW/gRQO3Db2kVv7mAkZtbt3NI3MyuIk76ZWUGc9M3MCuKkb2ZWECd9M7OCOOmbmRXEXTbNbMzqn/fsrsPdy0l/gnGffTMbipO+mY0LbtA0h8/pm5kVxC19M2up+vP91llu6ZuZFcRJ38ysIE76ZmYFcdI3MyuIk76ZWUHce8fMms49drqXW/pmZgVx0reiSZot6S5JqyWd3GD5LEk3S9osaW7dshMk3Z3/Tmhf1GajVynpVzgwPijpTkkrJF0taZ+aZT4wrCtJmgScDRwFzACOkzSjrtg6YD7wjbp1dwE+DbwSmAl8WtLOrY7ZbKyGTfoVD4xbgJ6IOAC4DDg9r+sDw7rZTGB1RKyJiMeBS4A5tQUiYm1ErAC21K3758BVEbExIh4CrgJmtyNos7Go0tKvcmBcGxGP5Jc3AnvlaR8Y1s32BO6ted2f57V6XbOOqZL0R1q5TwSuGMm6khZI6pPUt2HDhgohmTWFGsyLZq7rum3dpkqXzcoHhqTjgR7g9SNZNyIWA4sBenp6qh5045K7snWVfmDvmtd7AetHsO5hdeteV1+opLpt40OVln6lA0PSG4FPAEdHxGMjWdesQ5YD+0maLmlrYB7QW3HdZcCRknbO16mOzPPMulqVpD/sgSHpFcA5pIT/QM0iHxjWtSJiM3ASqU6uAi6NiJWSFko6GkDSIZL6gWOAcyStzOtuBD5LOj6WAwvzPLOuNuzpnYjYLGngwJgELBk4MIC+iOgFzgC2B74tCWBdRBwdERslDRwYUOCB4dM53S0ilgJL6+adUjO9nKc7JtSvuwRY0tIAzZqs0jAMFQ6MNw6xrg8MM7Mu4TtyzcwK4qRvZlYQJ30zs4J4aGUzG3fqO0isXfTmDkUy/rilb2ZWECd9M7OC+PSOmY2K70EZn9zSNzMriJO+mVlBnPTNzAric/pmVonP4U8MbumbmRXESd/MrCBO+mZmBXHSNzMriC/kTmAen8TM6rmlb2ZWECd9M7OCOOmbmRXESd/MrCBO+lY0SbMl3SVptaSTGyzfRtK38vKfSZqW5z9H0gWSbpe0StLH2x272Wg46VuxJE0CzgaOAmYAx0maUVfsROChiNgXOBM4Lc8/BtgmIl4GHAz87cAXglk3c9K3ks0EVkfEmoh4HLgEmFNXZg5wQZ6+DDhckoAAnidpMvBc4HHg9+0J22z03E+/BTww1bixJ3Bvzet+4JWDlYmIzZI2AbuSvgDmAPcD2wEfiIiNLY/YbIzc0reSqcG8qFhmJvAksAcwHfiQpBc+6w2kBZL6JPVt2LBhrPGajZmTvpWsH9i75vVewPrByuRTOTsCG4F3AldGxBMR8QDwE6Cn/g0iYnFE9EREz5QpU1qwC2Yj46RvJVsO7CdpuqStgXlAb12ZXuCEPD0XuCYiAlgHvEHJ84BXAb9oU9xmo+akb8WKiM3AScAyYBVwaUSslLRQ0tG52HnArpJWAx8EBrp1ng1sD9xB+vL4t4hY0dYdMBuFShdyJc0GzgImAedGxKK65bOAfwYOAOZFxGU1y54Ebs8v10XE0Zh1iYhYCiytm3dKzfSjpO6Z9es93Gi+dV5tRwoPMvhswyb9mr7MR5DOby6X1BsRd9YUWwfMBz7cYBN/jIgDmxCrmVlD7jFXXZWW/lN9mQEkDfRlfirpR8TavGxLC2I0M7MmqXJOv1Ff5j1H8B7b5i5rN0p6S6MC7tZmZtYeVZJ+lb7MQ5kaET2kLm7/LOlFz9qYu7WZmbVFldM7VfoyDyoi1ud/10i6DngFcM8IYrQm8QUuM6vS0q/Sl7khSTtL2iZP7wa8hpprAWZm1l7DJv0qfZklHSKpn9SF7RxJK/Pq+wN9km4DrgUW1fX6MTOzNqrUT79CX+blpNM+9ev9FHjZGGM0M7Mm8R25ZmYFcdI3MyuIk76ZWUGc9M3MCuKkb2ZWECd9M7OCOOmbmRXESd/MrCBO+mZmBal0R66ZlckPJ5l43NI3MyuIk76ZWUF8eqdQ9T/bSx1fX9Js4CxgEnBuRCyqW74NcCFwMPAgcGzN40EPAM4BdgC2AIfkB6mbdS239K1YkiYBZwNHATOA4yTNqCt2IvBQROwLnAmcltedDFwMvCciXgIcBjzRptDNRs1J30o2E1gdEWsi4nHgEmBOXZk5wAV5+jLgcEkCjgRWRMRtABHxYEQ82aa4zUbNSd9Ktidwb83r/jyvYZn8QKFNwK7Ai4GQtEzSzZI+2oZ4zcbM5/StZGowLyqWmQy8FjgEeAS4WtJNEXH1M1aWFgALAKZOnTrmgM3Gyi19K1k/sHfN672A9YOVyefxdwQ25vnXR8RvI+IR0pPlDqp/g4hYHBE9EdEzZcqUFuyC2ci4pW8lWw7sJ2k6cB8wD3hnXZle4ATgBmAucE1EhKRlwEclbQc8DryedKF33PMNWRObk74VKyI2SzoJWEbqsrkkIlZKWgj0RUQvcB5wkaTVpBb+vLzuQ5K+RPriCGBpRDhbdhl3TX42J30rWkQsJZ2aqZ13Ss30o8Axg6x7Manbptm44aTfBP45bGbjhS/kmpkVxEnfzKwgPr1jZsWoPRVb6kVdt/TNzAripG9mVpBKSV/SbEl3SVot6eQGy2fl8Uc2S5pbt+wESXfnvxOaFbiZmY3csEm/4vCz64D5wDfq1t0F+DTwStKIhp+WtPPYwzYzs9Go0tIfdvjZiFgbEStID5Ko9efAVRGxMSIeAq4CZjchbjMzG4UqSb/K8LOtWNfMzJqsSpfNKsPPjmldDz9rZu1W6rg8VVr6VYafHdO6Hn7WzKw9qiT9p4aflbQ1aZTB3orbXwYcKWnnfAH3yDzPzMw6YNiknx8RNzD87Crg0oHhZyUdDSDpEEn9pNEIz5G0Mq+7Efgs6YtjObAwzzMzsw6oNAxDheFnl5NO3TRadwmwZAwxmplZk3jsHbPCeWjwsngYBjOzgrilb4BHHzQrhVv6ZmYFcdI3MyuIk76ZWUF8Tn+U3ONhYpA0GzgLmAScGxGL6pZvA1wIHAw8CBwbEWtrlk8F7gROjYgvtCtus9FyS9+KVXHY8BOBhyJiX+BM4LS65WcCV7Q6VrNmcdK3kg07bHh+fUGevgw4XJIAJL0FWAOsbFO8ZmPmpG8lqzL091Nl8pAkm4BdJT0P+BjwmTbEadY0TvpWsipDfw9W5jPAmRHx8JBvIC2Q1Cepb8OGDaMM06x5fCHXSlZl6O+BMv2SJgM7AhtJjwCdK+l0YCdgi6RHI+IrtStHxGJgMUBPT0/V51CYtYyTvpXsqWHDgftIw4a/s65ML3ACcAMwF7gmIgJ43UABSacCD9cnfLNu5KRvxYqIzZIGhg2fBCwZGDYc6IuIXuA84CJJq0kt/Hmdi9hs7Jz0rWgVhg1/lPSciKG2cWpLgjNrAV/INTMriJO+mVlBnPTNzAripG9mVhAnfTOzgjjpm5kVxF02KyppKOX6ffXjE60EpTwy1EnfrEAlNWLsmXx6x8ysIE76ZmYFcdI3MyuIk76ZWUGc9M3MClIp6UuaLekuSaslndxg+TaSvpWX/0zStDx/mqQ/Sro1/321ueGbmdlIDNtlU9Ik4GzgCNJThJZL6o2IO2uKnQg8FBH7SpoHnAYcm5fdExEHNjluM7OWmcj3qlRp6c8EVkfEmoh4HLgEmFNXZg5wQZ6+DDhcUqNni5qZWQdVSfp7AvfWvO7P8xqWiYjNwCZg17xsuqRbJF0v6XWYmVnHVLkjt1GLvf4Bz4OVuR+YGhEPSjoYuFzSSyLi989YWVoALACYOnVqhZDMbCR8B64NqNLS7wf2rnm9F7B+sDKSJgM7Ahsj4rGIeBAgIm4C7gFeXP8GEbE4InoiomfKlCkj3wszM6ukStJfDuwnabqkrUkPhu6tK9MLnJCn5wLXRERImpIvBCPphcB+wJrmhG5mZiM17OmdiNgs6SRgGTAJWBIRKyUtBPoiohc4D7hI0mpgI+mLAWAWsFDSZuBJ4D0RsbEVO2JmZsOrNMpmRCwFltbNO6Vm+lHgmAbrfQf4zhhjtA6byEPOSpoNnEVq0JwbEYvqlm8DXAgcDDwIHBsRayUdASwCtgYeBz4SEde0NXizUfAduVasmntQjgJmAMdJmlFX7Kl7UIAzSfegAPwW+MuIeBnp1OZF7YnabGyc9K1ko74HJSJuiYiBDg0rgW3zrwKzruakbyUb6z0oA94O3BIRj9W/gaQFkvok9W3YsKFpgZuNlpO+lWws96CkhdJLSKd8/rbRG7g7snUbPy7RSjaSe1D6a+9BAZC0F/A94F0RcU/rw7VOmUidGZz0bUQm2EBUT92DAtxH6mr8zroyA/eg3MAz70HZCfgB8PGI+EkbY7YOG+/HgJN+Hd+uXo4x3oNyErAv8ClJn8rzjoyIB9q7F2Yj46RvRRvDPSifAz7X8gDNmswXcs3MCuKWvtkE5VOV1ohb+mZmBXHSNzMriJO+mVlBnPTNzAripG9mVhD33rExmUi3p4937q1jVbilb2ZWECd9M7OCOOmbmRXESd/MrCC+kGtmNgbjrTODW/pmZgUpvqXvbm5m1izj4QErxSd9a57xUOEnGjdabKR8esfMrCBu6ZuNI27Zjy/deJG3yKTvA6c9urHCj0eurxNPJ0+FFpH0fdDYeOL6aq1UKelLmg2cBUwCzo2IRXXLtwEuBA4GHgSOjYi1ednHgROBJ4F/iIhlTYvexo1uvcjrum3doJ2/iodN+pImAWcDRwD9wHJJvRFxZ02xE4GHImJfSfOA04BjJc0A5gEvAfYAfiTpxRHxZLN3xGykXLetnbrlF1yVlv5MYHVErAGQdAkwB6g9MOYAp+bpy4CvSFKef0lEPAb8StLqvL0bmhO+jVddcr6/a+p2tyQE67yh6kIzjpUqSX9P4N6a1/3AKwcrExGbJW0Cds3zb6xbd8/RBtvqD8M6oxUJr2J96Jq6bdYuVZK+GsyLimWqrIukBcCC/PJhSXdViOuZ2zhtpGt0xG7AbzsdRJt0dF+HqQ/7DBRrsKyddbtb6oPj6K4YYJA4KtbrIVVJ+v3A3jWv9wLWD1KmX9JkYEdgY8V1iYjFwOIqAY9nkvoioqfTcbTDONnXjtbtbvmMHEd3xdDqOKrckbsc2E/SdElbky5e9daV6QVOyNNzgWsiIvL8eZK2kTQd2A/4eXNCNxsz120rzrAt/Xwe8yRgGalb25KIWClpIdAXEb3AecBF+WLWRtLBQy53KenC2Gbg7927wbqF67aVSKnRYu0gaUH+uT/hlbSvo9Utn5Hj6K4YWh2Hk76ZWUE8yqaZWUGc9FtA0t6SrpW0StJKSe/L83eRdJWku/O/O3c61mYYYn9PlXSfpFvz35s6HWunSPpA/mzukPRNSdtKOl/Sr2o+nwPbEMf7cgwrJb0/z2t7vRwkjpbXF0lLJD0g6Y6aeQ33X8mXJa2WtELSQR2K4zBJm2o+l1PG9OYR4b8m/wG7Awfl6ecDvwRmAKcDJ+f5JwOndTrWFu/vqcCHOx1fp/9IN239Cnhufn0pMB84H5jbxjheCtwBbEfqxPEjUq+jttbLIeJoeX0BZgEHAXfUzGu4/8CbgCtI92S8CvhZh+I4DPh+s97bLf0WiIj7I+LmPP0HYBXpwJ8DXJCLXQC8pTMRNtcQ+2tPmww8N/f1344GffrbYH/gxoh4JCI2A9cDb6X99XKwOFouIn5M6oVVa7D9nwNcGMmNwE6Sdu9AHE3lpN9ikqYBrwB+BrwgIu6HlCiBP+lcZK1Rt78AJ+WfxksmyumskYqI+4AvAOuA+4FNEfHDvPjz+fM5U2lEz1a6A5glaVdJ25FasnvT/no5WBzQmfoy2P43GqajlY2Zof4fDpV0m6QrJL1kLG/ipN9CkrYHvgO8PyJ+3+l4Wq3B/v4r8CLgQFKy+2IHw+uYnLzmANNJI3I+T9LxwMeBPwMOAXYBPtbKOCJiFWmU0KuAK4HbSPcYtNUQcXRbfak01EYb3AzsExEvB/4fcPlYNuak3yKSnkNKgF+PiO/m2b8Z+HmY/32gU/E1W6P9jYjfRMSTEbEF+BppFMoSvRH4VURsiIgngO8Cr86nxSLSSJ3/Rhs+n4g4LyIOiohZpNMLd9OBetkojg7Wl8H2v9JQG62OIyJ+HxEP5+mlwHMk7TbaN3HSbwFJIt3JuSoivlSzqPaW/hOAf293bK0w2P7Wnf98K+lnfYnWAa+StF3+rA4HVtUc4CKdv2355yPpT/K/U4G3Ad+kA/WyURwdrC+D7X8v8K7ci+dVpNNy97c7Dkn/I9cRJM0k5e0HR/0urbpKXvIf8FrSz8AVwK35702kIXmvJrWurgZ26XSsLd7fi4Db8/xeYPdOx9rBz+gzwC9IiewiYBvgmvz53AFcDGzfhjj+kzR0xG3A4Xle2+vlIHG0vL6QvuTuB54gteRPHGz/Sad3zgbuyXH1dCiOk4CV+bO6kfQrcdTv7TtyzcwK4tM7ZmYFcdI3MyuIk76ZWUGc9M3MCuKkb2ZWECd9M7OCTLikL2knSe+teX2YpO+PcBvzJe0xyve/TtKQDzTO2//KaLbfjTTI0Mp52TF53pbBPpdh1n+5pBsk3S7pPyTt0I596kadrtsVtz9s/R9PJO0s6Xt5PKCfS3pp3fJJkm4Z7P9B0j6Srs7rXydpr5plUyX9MNf7O/O4VS034ZI+sBPw3mFLDW0+aYwUq2Yz8KGI2J80BO3fS5qRl91BuuPyx6Nc/1zScLMvA74HfKQVOzBOuG633/8Fbo2IA4B3AWfVLX8faVTZwXyBNFLnAcBC4J9qll0InJHr/UzaNCzLREz6i4AX5YcNnJHnbS/pMkm/kPT1mluaD5Z0vaSbJC2TtLukuUAP8PW8jedKOkXScqWHPiweWH8Ix0v6aS4/2Pghe0i6UumBCacPzJT0r5L6cov3MzXzF+XWwApJX6jfmNIDKC7ILYe1kt4m6fTcQr5SaWychvuc5//vvI+3SfqO0uiHKD3o48t5f9bkz+cZYoihlSNiVUTcNdSHNdT6wJ/y9BfGVcDbh9rWBNexui1pf0k/r3k9TdKKQeI8JreKfynpdTXl/1PSzfnv1Xn+7pJ+nOO5Y6B83XuvlfSPSr/4+iQdlPfpHknvqSn3kbwvK+qOncvz57BS0oKa+Q9L+nyu8zdKekGDfZlBujuWiPgFMG2gnFKr/c2khslgnlofuJY08B65UTM5Iq7K2344Ih4ZYjvN0+rbrdv9B0zjmQ8mOAzYRBosaSvgBtKwAc8BfgpMyeWOBZbk6euoueWamtvSSbeK/+UQ738d8LU8Pas2lpoy84E1wI7AtsCvgb1r3wuYlLd1AGkExrt4+pnGOzXY5qnA/8/79XLgEeCovOx7pLFdhtrnXWu29Tng/+Tp84Fv589uBrC6wue/Dtihwecy7G3s9evneOfk6Q8Cf+h0HSu4bt8KvDBPfwz45CD1/4t5+k3Aj/L0dsC2eXo/oC9Pfwj4RE2df36Dba4F/i5Pn0kapuH5wBTggTz/SGAxaeiErYDvA7Pqjqnnkn557ppfx8D+kh5g0mh//hH4Up6eSfpVenB+fRlwMEM85AT4BvC+PP22/J67ko7H75MG37sFOAOY1I56NJky/Dwi+gEk3Uo6eH5HeoLPVblxM4k0FkYj/1PSR0kVdxfSOBj/McT7fRPSgxIk7SBpp4j4XV2ZqyNiU47pTmAf0tjd78itkcmkJ1LNII1R8ihwrqQfkCpLI1dExBOSbs/7c2Wef3ve5z8dYp9fKulzpFMI2wPLarZ7eaSRD+8cpDVE3o8xDSU9yPrvBr6s9Ii4XuDxkW53gmtn3b4UeAfpF8ex+a+RgVFlb8rxQPoi+orSIyGfBF6c5y8HluRfopdHxK2DbLM3/3s7aYyiPwB/kPSopJ1ISf9IUgKFVIf3I/1K/AdJAw9p2TvPf5BUlwaOpZuAIxq87yLgrPzZ3p63v1nSX5C+cG6SdNggMQN8OO/3/BzLfaQvjsnA60jPnlgHfIvUGDxviG01RSlJ/7Ga6SdJ+y1gZUQcOtSKkrYF/oXUOrpX0qmk1vlQ6gc0ajTA0bNikjSdVEkOiYiHJJ1Pah1tVjpNdDgwjzQA0xsG22ZEbJH0ROTmBbCF4ff5fOAtEXFbrqCHDRLrYD//Gw0lXdlg60f6SX1kLvNi0s9pe1o76/a3gG9L+i4QEXH3MDENxAPwAeA3pF+hW5EaMQMNo1mk/9eLJJ0RERcOsc0tPHOfa+v2P0XEOXX7eBhpaOtDI+IRSdfV7GPtMVIb61Ny4+Ov87ZEeuzlr0jH4dFKz/HdFthB0sURcXzd+utJLfyBRs3bI2KTpH7glohYk5ddTrqe1fKkPxHP6f+B9NNvOHcBUyQdCinp6Okn0tRuY6CC/Db/pz3rnHYDx+ZtvpY0HOumirHvAPwXsCm3qI/K29ke2DHSWNrvJz1kYjSG2ufnA/fn5PtXI9loPhgaDSU95vX19BC8WwGfBL460u1PIB2t2xFxDyk5for0BTASOwL351+M/4v06wNJ+5BazF8j1YHRPnx8GfDuvB9I2jPXnR2Bh3LC/zNSYq1MqcfU1vnl3wA/jjS+/ccjYq+ImEb6ArimPuHn9XfLdRfSQ3OW5OnlwM6SpuTXbyD9om+5CZf0I+JB4Cf5otAZQ5R7nFTJT5N0G+l85avz4vOBr+afdI+RHuhwO+mJNcsrhPGQpJ+SEtSJI4j9NtLPx5WkyvGTvOj5wPfzhbPrSa2mERtmnz9FesThVaQhgEfiNaQD+Q35gtytuQWEpLfmVs2hwA8kLcvz95C0dLj1geMk/TLHtJ70sJEidUnd/hZwPOlUz0j8C3CCpBtJp3b+K88/DLhV0i2ki/T1vWMqifT4yW8AN+TTm5eRjpsrSb+iVwCfJQ1NPBL7Aysl/YLUCHvfMOWRtFDS0fnlYcBduQ6/APh8jvdJ0q/6q3O8Iv1ftJyHVjYzK8iEa+mbmdngSrmQ23SSziadlqh1VkQUe/rBJgbX7YnNp3fMzAri0ztmZgVx0jczK4iTvplZQZz0zcwK4qRvZlaQ/wbAdTzA3ygaDAAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 2 Axes>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcJHV9//HXe3aXvXdnj1lg2WVHDpHDJOrGM1EiGAkxmEdCIvxE0RBJ4i9RY4zB45cQo7+Y3y8xmqiJRAlgAI2YCBEwEOOCJBwuiNwgLAt7z+zsNcfuzO7MJ39U9dIMc/R0Vx9V834+HvOY6q7qqk9197zn29+u+pYiAjMzy7+2ZhdgZmbZcKCbmRWEA93MrCAc6GZmBeFANzMrCAe6mVlBONALRNLDkk5v8DYl6R8l7ZZ0zxjz3yXpjgkef7OkC+tU2+mSNtdj3RVs++2SbplgftNqmypJ6yT95jjzjpXUJ2nGZMtOcZuXSvqnWtcz3TjQJyBpo6T96Rt2RxpcC5pdF4CkKyR9svy+iDg1ItY1uJSfAd4ErIqIV071wRHxCxFxZfZlNVdEXB0RP1+6LSkkndCIbWcVqpWIiGcjYkFEDDdiezYxB/rkfikiFgAvB34a+PjoBdJWasOey1JrqEWsATZGRH+zCzGb7hzoFYqILcDNwGlwuBX0KUn/BQwAx0laKekGSbskPSnpPaXHpx8hr5P0dUm9ku6T9JNl809O17kn7To5p2zeFZL+TtJNkvqBi4C3Ax9OPz38W7rcRklnptOzJX1W0tb057OSZqfzTpe0WdIfSOqStE3Su8fb9/H2S9JFwJeB16R1/On4q9DfStor6TFJZ5TNeF5rUtJ7JD2aPkePSHq5pD+U9M1RK/xbSZ9Np5emn562pl0/35pgP74pqVvS05LeVzbvlZLWS9qXfhr7zDjruE3Sr6bTP5O2vM9Ob58p6f50+nBXk6Tb04f/KH2e3la2vjFfA0mLJV2V1vqMpI+XGg2juyMkdaZ1zJT0KeBngc+n2/r8OPvxDUnb09fkdkmnjrVcmeMl3ZMuf72kpaO3PcY2jpb0gKQPle3TV9J93SLpk6qwcTJRvenfxxck3Zi+b+6WdHzZ/JdIujV9/z4u6dcr2WYuRYR/xvkBNgJnptOrgYeBP0tvrwOeBU4FZgKzgNuALwJzgJ8CuoEz0uUvBQ4C56bLfgh4Op2eBTwJfBQ4Angj0AuclD72CmAv8DqSf8Jz0vs+OUG9nwDuAlYAHcB/l9V+OnAoXWYWcDbJP6Ul4zwPE+3Xu4A7JngO35Vu6/fTbb0t3ZelZc/jb6bTvwZsIfkkJOAEkk8ARwP9QHu63EygC3hFevtG4OvAknQbbyjbz83pdBtwL/DH6XN8HLABeHM6/07gHen0AuDV4+zPJ4C/Tac/CjwF/EXZvM+N9bwAAZxQdnvC1wC4CrgeWAh0Ak8AF5W9l/6pbF2d6fpnjn5OJ3hdfiNd92zgs8D9Eyy7Ln1dTgPmA98sbX+8bZfVfHHZer4FfCldxwrgHuC3xtnm6H0ct16Sv4VdwCvT98bVwNfSefOBTcC703kvB3YCpzY7X+qSWc0uoJV/SAKyD9gDPEMSanPTeeuAT5QtuxoYBhaW3ffnwBXp9KXAXWXz2oBtJK2pnwW2A21l868FLk2nrwCuGlXbFUwc6E8BZ5fNezNJ1wgkYbK/9EeY3tfFGCFWwX69i8kDfSugsvvu4bnwXMdzgf7vwPvHWc/NwHvS6bcAj6TTRwMjjPHPiOcH+quAZ0fN/wjwj+n07cCfAssneU+cATyQTn+HJLzuSm/fBvzKWM8LYwf6mK8BMAMYBE4pm/dbwLqy91JNgT5qn9rTxy8eZ/464NNlt08BhtI6x9r2Z0jei+eXPebIdJ/mlt13PvC9cbb5vH2cqF6Sv4Uvl80/G3gsnX4b8P1Rj/8S8CeVPj95+nGXy+R+OSLaI2JNRLw3IvaXzdtUNr0S2BURvWX3PQMcM9byETECbE4ftxLYlN436WMrtDJdR/n6Vpbd7omIQ2W3B0hapmOtZ7L9msyWSP+SxqmlZDXJP6KxXAlckE5fAHy17DG7ImL3JDWsAVYq6dLaI2kPSQv7yHT+RcCLgcck/UDSW8ZZz53AiyUdSfJp5SpgtaTlJC3E28d53FjGew2Wk3yKGP36TeU5H5ekGZI+LekpSftIwpd0u+Mpf/89Q/KpYrzl307Sor+u7L416WO2lT3/XyJpqWdR7/ay6fL38hrgVaNe97cDR0223TxyoNemPKS2AkslLSy771iSN3bJ6tJE2h+6Kn3cVpJQaJvgsaOHxZxsmMytJG/m8vVtneQx461nsv2azDGSVEEtm4Djx7gfko/rPyHpNJIW+tVlj1kqqX2SGjYBT6f/nEs/CyPibICI+HFEnE8SMH8BXCdp/uiVRMQASdfN+4GHImKIpDvrg8BTEbFzkjoqsZOke27061d6zvuBeWXzRofTZO+N/wW8FTgTWEzSyoakm2s8q8umj03rG29fL03nXVPWR76JpIW+vOz5XxQRk/XdV1tvySbgtlGv+4KI+J0KHps7DvSMRMQmkj/sP5c0R9JPkLT6ri5b7BWSfiX9AukDJG/wu4C7Sf5IPyxplpJjyX8J+NoEm9xB0g88nmuBj0vqSFuPfwxM+bjeCvdrMiuA96X79mvAycBNYyz3ZeBDkl6hxAmS1qR1HCBp8V0D3BMRz6b3byPpjvmipCXpNl4/xrrvAfZJ+iNJc9NW32mSfhpA0gWSOtJPSXvSx4x3KN5twO+mvyHpZii/PZbJXq/DIjkE8J+BT0lamD4HH+S51+9+4PVKjgFfTNJ1NJVtLSR57/WQ/GP4vxWUdYGkUyTNI+n3vy7GP1TxIMn3IfOBr0pqS1+nW4C/krRIUpuk4yW9oYJtV1NvybdJPlG9I31vzJL005JOnsI6csOBnq3zSVoPW4F/Jemnu7Vs/vUkfXq7gXeQ9LceTFt55wC/QNKy+SLwzoh4bIJtfQU4Jf0YOdZRHZ8E1gMPAA8C96X31WO/JnM3cCLJvn0KODciekYvFBHfSOdfQ/Kl8LeApWWLXAm8lOe6W0reQRIij5H0Q39gjHUPk/yT/CmSL6N3kvwDWZwuchbwsKQ+4HPAeek/kbHcRhIyt49zeyyXAlemr1clR1n8Hsk/+Q3AHSTPyeXpvtxK8iXwAySfFr496rGfA85VcsTP34yx7qtIuk22AI+QNCom81WSvurtJF+Ov2+ihdP39K+Q/DO/PP30+U6SrqRHSP4GriP5DmQy1dRbqqMX+HngPJL373aST2CzK11Hnuj5XZtWL5IuJflS7ILJlrWxSTqWJLSPioh9za7HrNW4hW65kLbwPkhyOJrD3GwMLzgZwKzVpF9O7iD52H1Wk8sxa1nucjEzKwh3uZiZFURDu1yWL18enZ2djdykmVnu3XvvvTsjomOy5Roa6J2dnaxfv76RmzQzyz1Jz0y+lLtczMwKw4FuZlYQDnQzs4JwoJuZFYQD3cysIBzoZmYF4UA3MyuI6RnoETAyMvlyZmY5Mv0CfWgA/uHn4LMvhd0VHatvZpYL0y/Qf3QNbP0h7NsMd3ym2dWYmWVm+gX6YzfB0uPhtF+Fx26EkfGuomVmli/TK9APDcLG78OL3wwnnQ393bD1/mZXZWaWiekV6N2Pw/AQrFoLa16b3Lf5B82tycwsI9Mr0Hc8lPw+8qWwaCUsOCrpTzczK4DpFejbH4KZc2HZ8cntlS9zoJtZYUyvQO95EpadAG0zktsrXgK7NsDwoebWZWaWgekV6Ls3wpI1z91edgKMHIQ9Ph7dzPJv+gT6yEgS3Es6n7tv2YnJ756nmlKSmVmWJg10SZdL6pL00BjzPiQpJC2vT3kZ6tsBhw6MCvQTkt89P25KSWZmWaqkhX4FcNboOyWtBt4EPJtxTfWxJy2zvazLZd5SmNOe9K2bmeXcpIEeEbcDu8aY9dfAh4HIuqi66N2W/F509HP3SbD8RNjpFrqZ5V9VfeiSzgG2RMSPKlj2YknrJa3v7u6uZnPZ6NuR/F5w1PPvb18Dezc1vh4zs4xNOdAlzQM+BvxxJctHxGURsTYi1nZ0dEx1c9np3Q5tM2Hesuffv3gV7N3i4XTNLPeqaaEfD7wI+JGkjcAq4D5JR034qGbr2wELjoS2Ubu8eFVy6GJ/V3PqMjPLyMypPiAiHgRWlG6nob42InZmWFf2ercngT7a4tXJ7z2bYGFr/08yM5tIJYctXgvcCZwkabOki+pfVh30bh87sNvTQHc/upnlXCVHuZwfEUdHxKyIWBURXxk1v7PlW+cAfeO10FcB8Klrb21wQWZm2ZoeZ4oeGoKBHlh49AvnzVnMvpjLSvU0vi4zswxNj0DvTw+XXDhGCx3YGss5Rq3/IcPMbCLTI9AH0rCeN/YIBV3RzgrtbmBBZmbZmyaBnnanjD4GPdXFElZoTwMLMjPL3jQJ9GTkgjP+/sExZ3dFOx3s9clFZpZr0yrQd8fCMWfviCXM0jDsH2vIGjOzfJgegZ4G9V7mjzm7K9qTidIAXmZmOTQ9An2gB+YsZpgZY85+LtB3NLAoM7NsTZNA3zXmF6Kdl9wIQBdpoPdtb2RVZmaZmiaB3gNzl447uyuWJBPucjGzHJsegb5/7BZ6ySBHsDfmucvFzHJtegT6wK7kcnM8180yWlcscZeLmeXaNAr08VvokH4x6ha6meVY8QP94AE42A9zl0y4WBftyRC7ZmY5VfxAP7A3+T23fcLFuqM9GcQr8nHNazOz0aZPoM+ZONB7YhEc2g9D/Q0oyswse1O+BF3uDO5Lfs9eBBwExv5itIdFyUR/N8xe0KDizMyyMw1a6OkoinMWT7jYzigFusdFN7N8mgaBXupyGTvQS631nkjnly6GYWaWM5VcJPpySV2SHiq77/9LekzSA5L+VdLEHdTNdDjQF024WE+UdbmYmeVQJS30K4CzRt13K3BaRPwE8ATwkYzrys6BtA99ki6X5/Whm5nl0KSBHhG3A7tG3XdLRBxKb94FrKpDbdk4sBfaZsKseRMuNsgRcMRC96GbWW5l0Yf+G8DN482UdLGk9ZLWd3c3ofV7YG9yhIs0+bLzl7uFbma5VVOgS/oYcAi4erxlIuKyiFgbEWs7Ojpq2Vx1BvdN2t1y2PwOB7qZ5VbVx6FLuhB4C3BGRAufXnlg79QCfffGupZjZlYvVbXQJZ0F/BFwTkQMZFtSxg7snfQIl8Pc5WJmOVbJYYvXAncCJ0naLOki4PPAQuBWSfdL+vs611m9A1PschnYCSMj9a3JzKwOJu1yiYjzx7j7K3WopT6m2uUSI7B/N8yfeLhdM7NWMz3OFJ1daaAvT36728XMcqjYgT58KBkLfSotdHCgm1kuFTvQBys7S/QwB7qZ5VixA/3wSIuVHuVSCnSfLWpm+VPwQJ9iC33eUkBuoZtZLhU80NORFmdX2EJvm5FcTNqBbmY5VOwrFg31AfCLl/2Ih2NfZY/x6f9mllPFbqEPJoE+wJzKHzN/ufvQzSyXih3oQ70A9EVlgd55yY1JoA840M0sfwoe6P0A9DO38sfMcwvdzPKp2IGedrns54jKHzN/eXK44/DBOhVlZlYfxQ70oT44YgExld0snf4/0FOfmszM6mRaBPqUzCuN5+JuFzPLl2IH+mAfHDF/ao/xAF1mllPFDvShPpg9xRZ66fR/d7mYWc4UPND74YiFU3tM2uVy6dduq0NBZmb1U+wzRQd7+e6WGVN7zNwlDIdYpgrPLDUzaxEFb6H30T+Vs0SBzo/ezC4WsgwHupnlSyXXFL1cUpekh8ruWyrpVkk/Tn8vqW+ZVRrqp7/Cs0TL7YpFLFVvHQoyM6ufSlroVwBnjbrvEuC7EXEi8N30dusZnHoLHZJAX6a9yVAAZmY5MWmgR8TtwK5Rd78VuDKdvhL45Yzrqt3ICBzsn9pp/6keFrEUt9DNLF+q7UM/MiK2AaS/V2RXUkYOpuO4xOwpP7QnFvpLUTPLnbp/KSrpYknrJa3v7m7gyTrpOC7VtNB3xSLa1c9MDmVdlZlZ3VQb6DskHQ2Q/u4ab8GIuCwi1kbE2o6Ojio3V4X04haVDp1brofkCkdL3O1iZjlSbaDfAFyYTl8IXJ9NORkaquLiFqmeSAJ9mY90MbMcqeSwxWuBO4GTJG2WdBHwaeBNkn4MvCm93VrSLpe+KrtcAJZpb6YlmZnV06RnikbE+ePMOiPjWrJVurhFFV0uO9Mul2XucjGzHCnumaJDpS9FqzkOPRn/ZamPdDGzHCluoA8mretqWuh7WODxXMwsd4ob6NVcTzQVtHk8FzPLnQIHeukol6mfWASl0/8d6GaWH8UN9MFemDmXYaY4fG4qGaDLgW5m+VHcQB/qn/rVisp4PBczy5sCB3oV1xMt0xMLWe7j0M0sR4ob6IN9U7/8XJldsYjFGoDhgxkWZWZWP8UN9GouEF2mNJ6LLxZtZnlR7ECvqcslDfT+Bo4QaWZWg+IG+mAfHFF9C33X4UDfmVFBZmb1VdxAr/Eol53ucjGznClwoNfaQk+/UHWXi5nlRDEDPaLmQC+N5+IuFzPLi2IG+sH9ECM1dbmUxnNhwIFuZvlQzEBPx3GppYUO6RejbqGbWU4UM9DToXNrDfQeB7qZ5UgxAz0dOreWLheAXSxyl4uZ5UZBAz3pcrngqw/XtJqeWOijXMwsN2oKdEm/L+lhSQ9JulbS1C8PVA/pBaKruVpRuZ5YDAf2ejwXM8uFqgNd0jHA+4C1EXEaMAM4L6vCapK20PuquFpRuV2kx6L75CIzy4Fau1xmAnMlzQTmAVtrLykDpasVRXVXKyrxeC5mlidVB3pEbAH+EngW2AbsjYhbRi8n6WJJ6yWt7+5uUDAOZtNC7/F4LmaWI7V0uSwB3gq8CFgJzJd0wejlIuKyiFgbEWs7Ojqqr3Qq0qNcBqixD93juZhZjtTS5XIm8HREdEfEQeBfgNdmU1aNhnoZjJkcZGZNq/F4LmaWJ7UE+rPAqyXNkyTgDODRbMqq0WAf/TW2zsHjuZhZvtTSh343cB1wH/Bguq7LMqqrNkP99Edt/efg8VzMLF9q6pOIiD8B/iSjWrIzlE0LHZLxXDrcQjezHCjmmaKDvZkFusdzMbO8KGagD/XXfJZoyS4W+UtRM8uFggZ6dl0uPeE+dDPLh2IG+mAf/TWeVFRyeDyXQ0OZrM/MrF6KGehDffRl1uXi8VzMLB+KF+jp9URrPUu05PDp/+52MbMWV7xAPzQII4cya6F7PBczy4viBXo6jktmfeg40M0sHwoY6Mn1RAeobejcksPjubjLxcxaXPECvTR0bgan/kMyngtqcwvdzFpe8QL9cJdLNn3oQRvdIwu45nv3ZrI+M7N6KWCgJ10uWZ0pCsl4LsvUm9n6zMzqoXiBXrpAdEZfikJypMtS7ctsfWZm9VC8QE+7XPoy6nKBZDyXZTjQzay1FTDQSxeIzi7Qd8YilmtvZuszM6uH4gX6YNqHnmELvSvaWaT9MDSQ2TrNzLJWvEAf6udgzGCQWZmtspv2ZKJve2brNDPLWgEDvS89qUiZrXJHLEkmendktk4zs6wVL9AH++jL8AgXgK5SoLuFbmYtrKZAl9Qu6TpJj0l6VNJrsiqsakO9mR6DDkkfOsCl1/xnpus1M8tSTReJBj4HfCcizpV0BDAvg5pqcvtDG1mobFvou1nAwZjBCu3JdL1mZlmqOtAlLQJeD7wLICKGgKZf1meB9mc2dG5J0EY3ix3oZtbSaulyOQ7oBv5R0g8lfVnS/NELSbpY0npJ67u763+x5XkMZnqWaElXtLOC3Zmv18wsK7UE+kzg5cDfRcTLgH7gktELRcRlEbE2ItZ2dHTUsLnKLNB++jMaOrdcdyyhwy10M2thtQT6ZmBzRNyd3r6OJOCbaj4H6M9o6NxyO6LdXS5m1tKqDvSI2A5sknRSetcZwCOZVFWD+ezP9CzRkq5Ykoy4eKjpXxOYmY2p1qNcfg+4Oj3CZQPw7tpLqsGhQY7QcGYXtyjXVTpbtL8LFq/KfP1mZrWqKdAj4n5gbUa11O7w0Ln1aKGngd67w4FuZi2pWGeKDmU/MFfJ4UD32aJm1qKKFegZX0+03OHT/3sd6GbWmooV6EP163LpYRHDIT53/ffpvOTGzNdvZlarYgV6HVvow8ygh8WswIcumllrKlag17EPHdKzRX0supm1qGIFeh0uEF1ueyzhKO2qy7rNzGpVrEAfKnW51KeFvi2WcbR66rJuM7NaFSvQ69xC3xbLWKo+5jBYl/WbmdWiWIE+1MtgzORgzSfAjm1rLAPgaHe7mFkLKlagD/bV7QtRgO0sBXC3i5m1pGIF+lA/A3UM9FILfaUD3cxaUMECva8ux6CXbI+0hY4D3cxaTyEC/fCZm4O9de1yGWIW3bGIo9Xjs0XNrOUUItAPG+qjv06HLJZsi2Ws9JeiZtaCihXog3301emQxZJtscwnF5lZSypWoDekhb7UX4qaWUsqVqA3qIW+SAPMZ39dt2NmNlXFCfQIGKrvl6KQBDr4WHQzaz3FCfSD+yFG6K/jYYsAW9JAP8aBbmYtpuZAlzRD0g8lfTuLgqpWGpirzi30TbECgNXqqut2zMymKosW+vuBRzNYT20G07HQ6/ylaDeLGYxZDnQzazk1BbqkVcAvAl/OppwaDO4DoJd5dd1M0Mam6GC1uuu6HTOzqaq1hf5Z4MPAyHgLSLpY0npJ67u76xiCBxoT6ACbooNj1eWzRc2spVQd6JLeAnRFxL0TLRcRl0XE2ohY29HRUe3mJldqoUcjAn2Fu1zMrOXU0kJ/HXCOpI3A14A3SvqnTKqqxoG9AOyr83HoAM/GChZrgEX01X1bZmaVqjrQI+IjEbEqIjqB84D/jIgLMqtsqtIul30xv+6b2hTJJ41j3Uo3sxZSnOPQ0y6Xep8pCuWHLvqLUTNrHZlcqy0i1gHrslhX1Q7sg1nzGT4wo+6bKgW6W+hm1koK1ELfC3MWNWRTvcxjdyzwF6Nm1lKKE+gH9sLsxgQ6PHfooplZqyhQoO+DOYsbtrln4kg6tb1h2zMzm0xxAn1wX8O6XACeipWs0s5kUDAzsxaQyZeirWDD5m08HLMbtr2nRlbSNjNg1wY48tSGbdfMbDyFaaEv1EBDzhIt2RArk4mdTzRsm2ZmEylMoC9igH0NGMelZEMclUzs/HHDtmlmNpFCBPpshpitQ+xrYAt9P3PYHMsd6GbWMgoR6AvT63s2YqTFchtGjnaXi5m1jGIEugaAxoy0WO6pWEnf1sfovKS5F2syM4OCBPoi+gHobcA4LuWeipUs0AGOZHdDt2tmNpZCBPpCpV0uTWihA5zQtqWh2zUzG0sxAp2ky2Uf9R86t9yPR1YB8BJtauh2zczGUohAX6yky6WRR7kA7GQxO6KdU9qeaeh2zczGUohAX5JeOWg3Cxq+7UdG1nCKHOhm1nyFCPTF6mMwZrGfxp36X/JorOEEbYFDgw3ftplZuUIE+hL60ta5Gr7tR0bWMEvD0P1Yw7dtZlauGIGuXnZH47tbAB6JNcnE9gebsn0zs5JCBHq7+tgTC5uy7Y1xFAMx24FuZk1XdaBLWi3pe5IelfSwpPdnWdhEOi+58XnTz3W5NN4IbTwWq7n7znVN2b6ZWUktLfRDwB9ExMnAq4H/LemUbMqamqSF3thj0Ms9OPIiTtPTMHyoaTWYmVUd6BGxLSLuS6d7gUeBY7IqbAqV0E4fe2hOlwvAfSMnMl+D0PVI02owM8ukD11SJ/Ay4O4x5l0sab2k9d3d3Vls7nkWsJ9ZGm7al6IA98aLk4lNL9h9M7OGqTnQJS0Avgl8ICL2jZ4fEZdFxNqIWNvR0VHr5l6gXclJRXua1IcOsDk66Ip22PyDptVgZlZToEuaRRLmV0fEv2RT0tS0l84SbdJRLglx38iJsOmeJtZgZtNdLUe5CPgK8GhEfCa7kqZmqXoBmtrlAnDvyImw+2nYt62pdZjZ9FVLC/11wDuAN0q6P/05O6O6KtbBXgC6aW/0pp/nv0dOSyaevq2pdZjZ9DWz2gdGxB0041z7UTq0B4DuWNzUOh6JY2HuUtiwDn7yvKbWYmbTU+7PFF2hPfTGXPYzp6l1BG1w3BuSQI9oai1mNj3lPtA7tKfprfPDjjsderdB9+PNrsTMpqHcB/oK7Wl6/3nJa76RPp1P3NzcQsxsWsp9oC9nL93RGoG+jWWw8mXw6L81uxQzm4ZyH+gt1eUCcPIvwZZ7Ye/m5w0iZmZWb/kO9KEBFmk/XbGk2ZUc9sabFiUTbqWbWYPlO9D7dgDQTeu00DfESh4a6YT7r2l2KWY2zeQ70PdtAWBHC7XQAb4+fDpsf4BTtbHZpZjZNJLvQN/zLACbIvtBv2px/fBrYcZs3jbje80uxcymkdwH+kiIrbG82ZU8zz4WwGm/yrkzbudll1zrL0fNrCHyHei7n2EHSxhiVrMreaHXvZ95GuTCmbc0uxIzmybyHeh7nmVzi7XOD1vxEm4ZfgXvnvEd2ul1K93M6i7fgb7zCZ4eObrZVYyp85Ib+ctDv84C9vOBmd9sdjlmNg3kN9AHdkF/F0/EqmZXMq4nYjVXD5/JBTP+g5dqQ7PLMbOCy2+gdz8GwJPRhOtST8FfHvo1umjnc7M+D4N9zS7HzAosN4HeecmNz++H3vEwAI+PrG5SRZXZxwI+ePC9dGoH/Mt7YGS42SWZWUHlJtBf4Nm7YOHRbGNpsyuZ1F0jp3DpoXfC4zfBt94LwwebXZKZFVBOAz2SQD/21bTARZMqctXwm+HnPg4PfA2uPhd6tx+e5yNgzCwLuQz0k/Us7NvMx+5vrVP+J/WGP4RzPs/+p/4bvvgauOcf4NBQs6sys4KoKdAlnSXpcUlPSrokq6Imc+6M2zkYM7hp+JWN2mQmOi+5kc5/Xspbhj7F3f0r4KYP0fVnJ/KRmVfD5vXuijGzmlR9kWhJM4AvAG8CNgM/kHRDRDySVXGjzWAYNtzGBTNu5YaR17KbRfXaVF09FcfwtqH/w8+0PcSFM27hN2Z8B758I/0xm4fiRbxq7atg2QmwYAVPBHtNAAAEj0lEQVTMWw7zlsLshTBzNsyc89zvtlkgJT9mNu1VHejAK4EnI2IDgKSvAW8Fsg/073yUx2dfxmwdhKtgaxzJpw+en/lmGkvcMfJS7hh5KUvYx2vaHuFVbY9yatszyVjq+3dVsco2oBTwSm4fnk5/m9WDGxWTe9tX4fg31nUTiiqvUC/pXOCsiPjN9PY7gFdFxO+OWu5i4OL05klAtVdQXg7srPKxeeV9nh68z9NDLfu8JmLyYWVraaGP9S/5Bf8dIuIy4LIatpNsTFofEWtrXU+eeJ+nB+/z9NCIfa7lS9HNQPlZPauArbWVY2Zm1aol0H8AnCjpRZKOAM4DbsimLDMzm6qqu1wi4pCk3wX+HZgBXB4RD2dW2QvV3G2TQ97n6cH7PD3UfZ+r/lLUzMxaSy7PFDUzsxdyoJuZFUTLBfpkwwlImi3p6+n8uyV1Nr7KbFWwzx+U9IikByR9V9KaZtSZpUqHjZB0rqSQlOtD3CrZX0m/nr7OD0u6ptE1Zq2C9/Wxkr4n6Yfpe/vsZtSZJUmXS+qS9NA48yXpb9Ln5AFJL8+0gIhomR+SL1efAo4DjgB+BJwyapn3An+fTp8HfL3ZdTdgn38OmJdO/8502Od0uYXA7cBdwNpm113n1/hE4IfAkvT2imbX3YB9vgz4nXT6FGBjs+vOYL9fD7wceGic+WcDN5Ocx/Nq4O4st99qLfTDwwlExBBQGk6g3FuBK9Pp64AzpFyfdzzpPkfE9yJiIL15F8kx/3lWyesM8GfA/wMONLK4Oqhkf98DfCEidgNERFeDa8xaJfsccHhApsUU4DyWiLgdmGjcjrcCV0XiLqBdUmYXRm61QD8G2FR2e3N635jLRMQhYC+wrCHV1Ucl+1zuIpL/8Hk26T5LehmwOiK+3cjC6qSS1/jFwIsl/ZekuySd1bDq6qOSfb4UuEDSZuAm4PcaU1pTTfXvfUpqOfW/HioZTqCiIQdypOL9kXQBsBZ4Q10rqr8J91lSG/DXwLsaVVCdVfIazyTpdjmd5BPY9yWdFhF76lxbvVSyz+cDV0TEX0l6DfDVdJ9H6l9e09Q1v1qthV7JcAKHl5E0k+SjWhVDE7aMioZQkHQm8DHgnIgYbFBt9TLZPi8ETgPWSdpI0td4Q46/GK30fX19RByMiKdJBrE7sUH11UMl+3wR8M8AEXEnMIdkAKsiq+uQKa0W6JUMJ3ADcGE6fS7wn5F+25BTk+5z2v3wJZIwz3vfKkyyzxGxNyKWR0RnRHSSfG9wTkSsb065Navkff0tki+/kbScpAtmQ0OrzFYl+/wscAaApJNJAr27oVU23g3AO9OjXV4N7I2IbZmtvdnfCo/zLfATJN+Qfyy97xMkf9CQvOjfAJ4E7gGOa3bNDdjn/wB2APenPzc0u+Z67/OoZdeR46NcKnyNBXyG5HoCDwLnNbvmBuzzKcB/kRwBcz/w882uOYN9vhbYBhwkaY1fBPw28Ntlr/MX0ufkwazf1z7138ysIFqty8XMzKrkQDczKwgHuplZQTjQzcwKwoFuZlYQDnQzs4JwoJuZFcT/AFLrCVJurufkAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "fig,ax=plt.subplots(1,2)\n", | |
| "#ax[0].set_xlim(15,30)\n", | |
| "#ax[1].set_xlim(80,115)\n", | |
| "j=ax[0].hist(theta_b,density=True,bins=40)\n", | |
| "j=ax[1].hist(theta_v,density=True,bins=40)\n", | |
| "ax[1].set_xlabel('theta_v has mean '+str(np.round(np.mean(theta_v),2)))\n", | |
| "ax[0].set_xlabel('theta_b has mean '+str(np.round(np.mean(theta_b),2)))\n", | |
| "plt.show()\n", | |
| "fig,ax=plt.subplots(1)\n", | |
| "s=ax.hist(1-p_nolane,bins=50,density=True)\n", | |
| "ax.plot(np.linspace(0,1,1000),beta.pdf(np.linspace(0,1,1000),b=81,a=6.6))\n", | |
| "ax.set_title('Proportion of bicycles without a bike lane')\n", | |
| "plt.show()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 92, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEWCAYAAABG030jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHjRJREFUeJzt3XmYHVWZx/HvzyRA2Je0KGujomyOARpEeUREVBBBZxQFYdwn6qDouDuOI6Kj6DDj4KBiRHFlkUVkwEFwiYBCoMMeFtnCFiDNJgSQ9Z0/zulJ5ebevtVLdffp/D7P00/fW/dU1Vvn1n3vqVNV5yoiMDOzcjxrogMwM7PhceI2MyuME7eZWWGcuM3MCuPEbWZWGCduM7PCTNnELekYSZ+vPP+gpHskLZW0gaRdJd2Qn79pImMdqdZtnMokvULS9RMdRzuSFkrafYjX50l6X0PrXiRpzyaWPVlMZP1OVtMnOoCRkLQI2BB4CngauAb4MTA3Ip4BiIgPVMrPAP4T2CUirsjTDgeOjoijxjf6sVPdxqlGUgBbRsSNABFxPvCiiY2qvYjYdvCxpMOAF0TEwRMX0dQyUfWbvyx+GhGbNL2u4Sq5xb1vRKwFbA4cAXwa+H6HshsCqwELK9M2b3lem6Qiv/BGajy3d2WrW7MRiYji/oBFwJ4t03YGngG2y89/CHwZeCHwCBDAUuB3wE257GN52qrAOqTEfxdwZ553Wl7Wu4A/At8A7ge+nKe/B7gWeAD4NbB5JZ4APgDckF//FqDK6/+Q532YdMSwQ56+EXAqMADcAhw6RD38sBLL7sAdwMeBJXk73j3EvPOArwIXA38Bfgmsn1/rzfG/F7gNOC9P34/0Zfdgnn/rlvfks3lbHgCOA1Zr2d4bc/2dAWzUUleH5Lq6BTgvT3skvz9vG9y+yjxb5xgezDHt11Iv3wLOyvU7H3h+fk35fVySt/tK8j7TUj+vAq6qPP8NcHHl+QXAm6r7I7AX8ATwZI77ikpdf4m0Dz0MnAPM6vC+rAecmd//B/LjTep8FkifgQtzndwFHA2sMox9suP+3LLOs4EPtUy7Avi7yVi/DL3fBqkFv9xnCliDlB+eyetaSmWfnei/CQ9gREG3Sdx5+m3AB6tvQH7cm9+g6Z2WAZwOfDe/Yc8mJbT359feReqW+TCpe2km8CZSIto6T/sX4E8tO8SZwLrAZqQP4l75tf1JXw475R39BaQjgGcBC4B/BVYBngfcDLyuQz1Ut3H3HOPhwAzg9cCjwHod5p2XY9gub/OppMPCan39OL82k2VfgK/Jy/9U3v5VKvV5NbApsH7+EA3GtgdwL7AD6Uvyv8lfBpW6OjfPN7PDB2p3cuLO678R+OdcT3uQPrAvqtTL/aRENh34GXBifu11uY7XzXW/NfDcNvWzGumDOysv425gMbBWro/HgA1a9yXgsMF6bKnrm3IdzszPj+jwvmwAvBlYPa/rZOD0Op8FYEdglxxvLykJf7TmPjnk/tyyzncAf6w834aUFFedbPVL9/22beJu3ecm21/JXSXtLCZ9+IdF0obA3qSd/JGIWEJqNRxQXXZE/HdEPBURjwHvB74aEddGxFPAV4DZkjavzHNERDwYEbcBvwdm5+nvA74eEZdEcmNE3EpK5D0RcXhEPBERNwPfa4ljKE8Ch0fEkxHxK1IrYah+4Z9ExNUR8QjweeCtkqZVXj8s18djpFbvWRFxbkQ8CRxJ+pC8vFL+6Ii4PSLuB/4NODBPPwj4QURcGhGPk1rmL5PUW5n3qxFxf15XN7sAa5Lq94mI+B0pIR1YKXNaRFyc35ufsazunyQlh61Irc1rI+Ku1hVExF+BfmA3oI/UcrwA2DWv/4aIuK9GrIOOi4g/5+37eSWe1vXeFxGnRsSjEfEwqR5fWWcFEbEgIi7K++giUkOkdd5O+2Sd/XnQL1peO4hU348z+eq3zn5bnKmWuDcmtbSGa3PSt/Fdkh6U9CBpp392pcztbeY5qlL+flILY+NKmbsrjx8lJRtIrdKbOsSx0eAy83L/mdRHX8d9+UPXbp3tVLfpVlIdzOrw+ka5DACRTgLfzvLb27q8jTrMuxS4b4h5u9kIuD3HUF1f17rPSf5oUjfBPZLmSlq7w3r+QGp17ZYfzyMlwlfm58PRaV9YjqTVJX1X0q2SHiJ1G63b8oXalqQXSjpT0t153q+w/Ps5VBx19mcA8hfKWSxrUBxA+nKcjPVbZ78tzpRJ3JJ2Ir0ZF4xg9tuBx0n9Yuvmv7WjcjabdEjVOs/7K+XXjYiZEfGnmut7fofpt7Qsc62IeP0ItqmOTSuPNyO1lu6tTKtu82LShxsAScrz3znE8hZ3mHcNUpdAdd7hDFO5GNhUUnX/3axleR1FxDcjYkdgW9Kh9Cc7FG1NLH+ge2IZ7XCbHycdJb00ItbO64aURLv5DnAd6WqctUlf+nXmg+HvzycAB0p6GakF+/vBFyZZ/Xbbbx8ldUsNes4o1jVuik/cktaW9AbgRFLf11XDXUY+lDsH+I+8vGdJer6koQ5RjwE+K2nbHMc6kvavucpjgU9I2lHJC/Jh58XAQ5I+LWmmpGmStstfSk04WNI2klYn9Y2fEhFPdyj7c2AfSa/Ol1d+nPRlV/1gHyJpE0nrk5LGSXn68cC7Jc2WtCqpJTg/H853cg+pj7+d+aR+y09JmpEv29qXtA8MSdJOkl6at+ER4K+kS0rb+RMpie5MOnG2kJQEXkpqCXeKu7flS2U41iL17z6Y6/ELw5z3IWCppK2ADw5j3uHuz78i1cXhwEmDRz+TsH677beXA2/Pn7W9WL5r6R5gA0nr1FzXuCk5cf+PpIdJLYXPka7TfvcolvcO0omuwasiTgGe26lwRPwC+BpwYj4svZrUT95VRJxM6rs8nnRS7XTSFR1PkxLQbNLVFfeSknxTO85PSCdj7iadLDp0iJivBw4mnVi8N8e5b0Q8USl2POkL8Ob89+U8729Jfeinkq52eD7d++0PA36UD93f2hLLE6QrBfbOsXwbeEdEXNdtg4G1SecNHiAdQt9H6vdst82PAJcCCyvbeSFwaz4P0s7J+f99ki6tEU+r/yK1YO8FLiJdwVHXJ4C3k/ap77Hsi7Or4e7PuT/7NNLVHsdXXppU9Vtjv/1InvYgqa/+9Mq815GOLG7O++FGTBKKmLRHA9YgSfNIRyjHjtHyFgHvi4jfjMXyzKyzklvcZmYrJSduM7PCuKvEzKwwbnGbmRWmkQF9Zs2aFb29vU0s2sxsSlqwYMG9EdFTp2wjibu3t5f+/v4mFm1mNiVJurV7qcRdJWZmhXHiNjMrjBO3mVlhnLjNzArjxG1mVhgnbjOzwjhxm5kVxonbzKwwTtxmZoVp5M5Js5L0fuasWuUWHbFPw5GY1eMWt5lZYdziNqupXcvcrXCbCG5xm5kVxonbzKwwtRK3pH+StFDS1ZJOkLRa04GZmVl7XRO3pI2BQ4G+iNgOmAYc0HRgZmbWXt2ukunATEnTgdWBxc2FZGZmQ+mauCPiTuBI4DbgLuAvEXFOazlJcyT1S+ofGBgY+0jNzAyo11WyHvBGYAtgI2ANSQe3louIuRHRFxF9PT21fjbNzMxGoE5XyZ7ALRExEBFPAqcBL282LDMz66RO4r4N2EXS6pIEvBq4ttmwzMysk653TkbEfEmnAJcCTwGXAXObDsysVL7D0ppW65b3iPgC8IWGYzEzsxp856SZWWE8yJTZKNQdEtZsLDlxm40D93vbWHJXiZlZYZy4zcwK464Sm7LcPWFTlVvcZmaFceI2MyuME7eZWWGcuM3MCuOTk2YTpNPNOz6Bat04cVtxfLWIrezcVWJmVhgnbjOzwjhxm5kVxn3ctlLxaH42FdT5seAXSbq88veQpI+OR3BmZraiOj9ddj0wG0DSNOBO4BcNx2UGuIVs1s5w+7hfDdwUEbc2EYyZmXU33MR9AHBCuxckzZHUL6l/YGBg9JGZmVlbtRO3pFWA/YCT270eEXMjoi8i+np6esYqPjMzazGcFvfewKURcU9TwZiZWXfDSdwH0qGbxMzMxk+txC1pdeA1wGnNhmNmZt3UugEnIh4FNmg4FjPDg2hZd75z0qYEX+9tKxOPVWJmVhgnbjOzwjhxm5kVxonbzKwwTtxmZoVx4jYzK4wTt5lZYZy4zcwK48RtZlYYJ24zs8I4cZuZFcaJ28ysME7cZmaF8eiAZgXwUK9WVfeHFNaVdIqk6yRdK+llTQdmZmbt1W1xHwWcHRFvyT8avHqDMZmZ2RC6Jm5JawO7Ae8CiIgngCeaDcvMzDqp01XyPGAAOE7SZZKOlbRGw3GZmVkHdbpKpgM7AB+OiPmSjgI+A3y+WkjSHGAOwGabbTbWcdoU45NtZiNXp8V9B3BHRMzPz08hJfLlRMTciOiLiL6enp6xjNHMzCq6Ju6IuBu4XdKL8qRXA9c0GpWZmXVU96qSDwM/y1eU3Ay8u7mQzMxsKLUSd0RcDvQ1HIuZmdXgW97NzArjW97NCuUrc1ZebnGbmRXGidvMrDBO3GZmhXHiNjMrjBO3mVlhfFWJTRrtrpIwsxW5xW1mVhgnbjOzwjhxm5kVxonbzKwwTtxmZoVx4jYzK4wTt5lZYZy4zcwKU+sGHEmLgIeBp4GnIsI/qmBmNkGGc+fkqyLi3sYiMTOzWtxVYmZWmLqJO4BzJC2QNKddAUlzJPVL6h8YGBi7CM3MbDl1E/euEbEDsDdwiKTdWgtExNyI6IuIvp6enjEN0szMlqmVuCNicf6/BPgFsHOTQZmZWWddE7ekNSStNfgYeC1wddOBmZlZe3WuKtkQ+IWkwfLHR8TZjUZlZmYddU3cEXEz8JJxiMXMzGrwL+CYTSHtfkVo0RH7TEAk1iRfx21mVhi3uM1WQm6Zl80tbjOzwjhxm5kVxonbzKwwTtxmZoVx4jYzK4wTt5lZYXw5oNkU1+7SPyubW9xmZoVx4jYzK4wTt5lZYZy4zcwK48RtZlYYJ24zs8LUTtySpkm6TNKZTQZkZmZDG06L+yPAtU0FYmZm9dRK3JI2AfYBjm02HDMz66Zui/u/gE8Bz3QqIGmOpH5J/QMDA2MSnJmZrahr4pb0BmBJRCwYqlxEzI2Ivojo6+npGbMAzcxseXVa3LsC+0laBJwI7CHpp41GZWZmHXVN3BHx2YjYJCJ6gQOA30XEwY1HZmZmbfk6bjOzwgxrWNeImAfMayQSMzOrxeNxW6M8FrTZ2HNXiZlZYZy4zcwK48RtZlYYJ24zs8I4cZuZFcaJ28ysME7cZmaFceI2MyuMb8CxEfGNNVNPu/d00RH7TEAk1o1b3GZmhXHiNjMrjBO3mVlhnLjNzArjxG1mVhgnbjOzwtT5seDVJF0s6QpJCyV9cTwCMzOz9upcx/04sEdELJU0A7hA0v9GxEUNx2ZmZm10TdwREcDS/HRG/osmgzIzs85q9XFLmibpcmAJcG5EzG9TZo6kfkn9AwMDYx2nmZlltRJ3RDwdEbOBTYCdJW3XpszciOiLiL6enp6xjtPMzLJhXVUSEQ+SfuV9r0aiMTOzrupcVdIjad38eCawJ3Bd04GZmVl7da4qeS7wI0nTSIn+5xFxZrNhmZlZJ3WuKrkS2H4cYrFJykO4mk0uvnPSzKwwTtxmZoXxL+CY2bD4l3ImnlvcZmaFceI2MyuME7eZWWGcuM3MCuPEbWZWGCduM7PCOHGbmRXGidvMrDBO3GZmhfGdk2bWkQcYm5zc4jYzK4wTt5lZYZy4zcwK07WPW9KmwI+B5wDPAHMj4qimA7OJ4T5Ns8mvzsnJp4CPR8SlktYCFkg6NyKuaTg2MzNro2tXSUTcFRGX5scPA9cCGzcdmJmZtTesPm5JvaTfn5zf5rU5kvol9Q8MDIxNdGZmtoLaiVvSmsCpwEcj4qHW1yNibkT0RURfT0/PWMZoZmYVtRK3pBmkpP2ziDit2ZDMzGwoXRO3JAHfB66NiP9sPiQzMxtKnRb3rsDfA3tIujz/vb7huMzMrIOulwNGxAWAxiEWMzOrwXdOmpkVxqMDrsR8l6RZmdziNjMrjFvcKwm3rs2mDiduM2tEu8bCoiP2mYBIph4nbjMbNR/RjS/3cZuZFcYtbjMbN+4+GRtucZuZFcYt7inI/Y1mU5tb3GZmhXHiNjMrjBO3mVlhnLjNzArjxG1mVhgnbjOzwtT56bIfSFoi6erxCMjMzIZWp8X9Q2CvhuMwM7Oa6vx02XmSepsPxcxWRp1uGPOt8J2NWR+3pDmS+iX1DwwMjNVizcysxZgl7oiYGxF9EdHX09MzVos1M7MWvqrEzKwwTtxmZoXpenJS0gnA7sAsSXcAX4iI7zcdmNXjkQDNVj51rio5cDwCMTOzetxVYmZWGP+QwiTlLhBb2flnzjpzi9vMrDBO3GZmhXHiNjMrjBO3mVlhfHLSzIrhE5aJW9xmZoVxi3ucucVgNrZWxs+UE7eZTTlTPZk7cU8CvtnGzIbDidvMVlqltsx9ctLMrDBucTfIXSBmk8dU+jw6cQ9TqYdWZjZ1uKvEzKwwtRK3pL0kXS/pRkmfaTooMzPrrM5Pl00DvgW8BrgDuETSGRFxTRMBjUdXRN11TKU+MTOrp4Tu0Dp93DsDN0bEzQCSTgTeCDSSuNupm0AnqnKd4M2mtsmWg+ok7o2B2yvP7wBe2lpI0hxgTn66VNL1ow9vePS1YRWfBdw7wnknynIxF6C0eMExj5fSYq4V7yjzyOZ1C9ZJ3GozLVaYEDEXmFt3xRNNUn9E9E10HMNRWsylxQuOebyUFvNki7fOyck7gE0rzzcBFjcTjpmZdVMncV8CbClpC0mrAAcAZzQblpmZddK1qyQinpL0IeDXwDTgBxGxsPHImldMt05FaTGXFi845vFSWsyTKl5FrNBdbWZmk5jvnDQzK4wTt5lZYaZ04pa0vqRzJd2Q/6/Xodw7c5kbJL2zMn1evtX/8vz37IbiHHJIAUmrSjopvz5fUm/ltc/m6ddLel0T8Y1lzJJ6JT1WqdNjJlHMu0m6VNJTkt7S8lrbfWQSx/t0pY7H7WKCGjF/TNI1kq6U9FtJm1deG/c6HoOYJ6SeiYgp+wd8HfhMfvwZ4GttyqwP3Jz/r5cfr5dfmwf0NRzjNOAm4HnAKsAVwDYtZf4ROCY/PgA4KT/eJpdfFdgiL2faONTraGLuBa6egH2hTsy9wN8APwbeUmcfmYzx5teWTtI6fhWwen78wcp+Me51PNqYJ6qeI2Jqt7hJt+b/KD/+EfCmNmVeB5wbEfdHxAPAucBe4xQfVIYUiIgngMEhBaqq23EK8GpJytNPjIjHI+IW4Ma8vMkc80TpGnNELIqIK4FnWuadiH1kNPFOlDox/z4iHs1PLyLdFwIT9zkcTcwTZqon7g0j4i6A/L9dV0e7W/o3rjw/Lh8Gfb6hxNNt/cuViYingL8AG9SctwmjiRlgC0mXSfqDpFc0HWxrPNlw6moi6nm061xNUr+kiyS1a7A0Ybgxvxf43xHOO1ZGEzNMTD2X/0MKkn4DPKfNS5+ru4g20wavkTwoIu6UtBZwKvD3pMPSsVRnSIFOZWoNR9CA0cR8F7BZRNwnaUfgdEnbRsRDYx1kzXiannekRrvOzSJisaTnAb+TdFVE3DRGsXVSO2ZJBwN9wCuHO+8YG03MMDH1XH6LOyL2jIjt2vz9ErhH0nMB8v8lbRbR8Zb+iLgz/38YOJ5muiHqDCnw/2UkTQfWAe6vOW8TRhxz7ta5DyAiFpD6F1/YeMSjq6uJqOdRrTMiBvfhm0nnarYfy+A6qBWzpD1JDav9IuLx4czbgNHEPFH1POVPTv47y5+c/HqbMusDt5BOiKyXH69POhqZlcvMIPXTfqCBGKeTTsRswbKTI9u2lDmE5U/0/Tw/3pblT07ezPicnBxNzD2DMZJOCN0JrD8ZYq6U/SErnpxcYR+ZxPGuB6yaH88CbqDlhNsE7hfbk76st2yZPu51PAYxT0g9R8SUT9wbAL/NFfrbwR2BdLhzbKXce0gn9m4E3p2nrQEsAK4EFgJHNZUUgdcDf847x+fytMNJ3+4AqwEn5/guBp5Xmfdzeb7rgb3HsW5HFDPw5lyfVwCXAvtOoph3IrXAHgHuAxYOtY9M1niBlwNX5Tq+CnjvJKrj3wD3AJfnvzMmso5HE/NE1rNveTczK0zxfdxmZisbJ24zs8I4cZuZFcaJ28ysME7cZmaFceIeY5XRwhZKuiKPLPas/FqfpG/mx6tK+k0u+zZJr8jzXC5p5sRuRXtKoyVOmh9MHS0l38yjwl0paYcO5f5N0u2SlrZM7zg63xjHubR7qbbz7S7pzA6vdR1VUunnCufn0fpOUvrpwurrb5EUg/uEpBmSfiTpKknXSvpspWzbEfgknV8ZXW+xpNNHsq0rGyfusfdYRMyOiG2B15CuEf0CQET0R8Shudz2wIxc9iTgIODI/PyxbivJScfv3+jsDWyZ/+YA3+lQ7n9of9fsbcC7SHfVFkPSNqSborYlDeT0bUnT2hT9GvCNiNgSeIA0TsfgMtYCDgXmV8rvT7oh5cXAjsD7lYbxnQZ8i1Tf2wAH5hiIiFfkfX42cCFw2thu7dTkD36DImIJKSF8KCfa3SWdqTSu90+B2bml8X7grcC/SvoZgKRPSroktwS/mKf15pbMt0k3r2wq6bWSLswtv5MlrZnLLpL0xTz9Kklb5elrSjouT7tS0pvz9LbLaWN/SRdL+rPyAFE5rvPzvJdKenme/lxJ5+VtvFptBpTKcX4lr7tf0g6Sfi3pJkkfqJRboT7y9NMlLchHK3Mq05fmlvIVSgMAbdhmW94I/DiSi4B1lYdIaHkfL4o8WFnL9EVRY3S+4caYW7oX5u39UodlriHprDzv1ZLelqfvJek6SRcAf9chpK6jSkoSsAfpjmFYcXTNL5GGTf5rtUqANZSGOJgJPAE8RI0R+PIXwR6AW9w1OHE3LNIYBs+iMjJhTujvA87PrY3vAmcAn4yIgyS9ltQK3BmYDewoabc8+4tIyWZ70h1z/wLsGRE7AP3AxyqrvzdP/w7wiTzt88BfIuLFEfE3pIFxZnVZTtX0iNgZ+Cj5SII0Bsxr8rxvA76Zp78d+HVuTb2EdNdZO7dHxMuA88m3bwO7kO5eo0t9vCcidiTdDXuopMERCNcALoqIlwDnAf/QZr3jNSLdcGM8CvhOROwE3N1hmXsBiyPiJRGxHXC2pNWA7wH7Aq+g/eBrUG+7NwAejDSy43JlJG0PbBoRrd0wp5D2ybtIRyNHRsT9Ndf3t8Bvo/nBxqaE4kcHLMRwh4N9bf67LD9fk5S4bgNuza1DSMltG+CPqYHEKqTDzUGDh50LWNb62pN0mAxARDwg6Q1dllNVXWZvfjwDOFrSbOBplg0adQnwA0kzgNMjolPiHvzlkKuANSMN6vWwpL9KWneI+jiPlAj/Nk/fNE+/j9TaG0wsC0jdVq3Ga0S64ca4K2loAICfkLosWl0FHCnpa8CZEXF+rv9bIuIGAEk/JR3xtRrx6I65e+4bpC6iVjuT3v+NSON4nK80emed9R0IHNumnLXhxN0wpeEenya1SreuOxvw1dwSry6rl9SiqZY7NyIO7LCcwVHMnmbZey3af0iHWk63Zf4TaSyHl5COLv4KEBHn5ZbxPsBPJP17RLQbFndwmc9UHg8+n07n+tid9EX0soh4VNI80hgpAE/GsvEcqrFWNT4i3ShiHPILJCL+rDQs7uuBr0o6h/QFWOeLp85230vqOpqeW92DZdYCtgPm5S/55wBnSNqPdIR1dkQ8CSyR9EfSUcbtQ60vH4HsTGp1Ww3uKmmQpB7gGODoyge0jl8D76n0V2+s9r93eRGwq6QX5HKrS+o2ROo5wIcqMa43wuVUrQPcFRHPkMYsn5aXszmwJCK+B3wfaHvVRg2d6mMd4IGcELciHYEMxxnAO5TsQupCWqEve5RGEuMfWXZUdFC7ApI2Ah6NiJ8CR5Lq9jrSj1Q8Pxfr9EV8BnCA0pVNW5COAC6uFsj76+9J3VYA7wR+GRF/iYhZEdEbEb2kfWe/iOgnHRHuketzjbyt15GOvLbMffer5G2r/j7j/qSjhmp/uQ3BiXvszcwn4xaSRhU7B/hil3mWExHnkK5UuFDSVaS+w7XalBsgHbKeIOlK0odoqy6L/zKwXj6hdQXwqhEup+rbwDslXUTqJhk8KtgduFzSZaRD/6OGscz/N0R9nA1MzzF/Kcc9HL8iDel5I6lv+B8HX5B0eeXx1yXdAawu6Q5Jh+XpO+Xp+wPfze95q5HE+BHgEEmXkBJ/Oy8GLs5xfg74ck58c4Cz8snJW9vNGBELgZ8D1+T4DomIp/M2/Sp/KQB8GviYpBtJfd7f7xL3t0jdWFeTkvVxEXFlbrF/iPQFfC1piN9qXR0AnNBl2Vbh0QHNzArjFreZWWGcuM3MCuPEbWZWGCduM7PCOHGbmRXGidvMrDBO3GZmhfk/rJxbWrw6fU4AAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "fig,ax=plt.subplots(1)\n", | |
| "j=ax.hist(p_nolane-p_lane,bins=60,density=True)\n", | |
| "ax.set_title('Difference in proportions with a lane vs without')\n", | |
| "x=np.linspace(0,1,1000)\n", | |
| "ax.set_xlabel('Difference has mean '+str(np.round(np.mean(p_nolane-p_lane),2))+' and sd '+str(np.round(np.sqrt(np.var(p_nolane-p_lane)),5)))\n", | |
| "\n", | |
| "plt.show()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 3", | |
| "language": "python", | |
| "name": "python3" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython3", | |
| "version": "3.6.5" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 2 | |
| } |