diff --git a/lectures/02-prob-review.pdf b/lectures/02-prob-review.pdf
new file mode 100644
index 0000000..6f85a47
Binary files /dev/null and b/lectures/02-prob-review.pdf differ
diff --git a/lectures/02_binomial_polls.ipynb b/lectures/02_binomial_polls.ipynb
new file mode 100644
index 0000000..4a8d072
--- /dev/null
+++ b/lectures/02_binomial_polls.ipynb
@@ -0,0 +1 @@
+{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"02_binomial_polls.ipynb","provenance":[]},"kernelspec":{"name":"python3","display_name":"Python 3"}},"cells":[{"cell_type":"code","metadata":{"id":"AFSOME3-nOx9","colab_type":"code","outputId":"3201fda1-3490-4150-e84a-78b06534aa45","executionInfo":{"status":"ok","timestamp":1570651515404,"user_tz":-120,"elapsed":1305,"user":{"displayName":"Jan Chorowski","photoUrl":"","userId":"18100051538790809279"}},"colab":{"base_uri":"https://localhost:8080/","height":34}},"source":["%pylab inline"],"execution_count":1,"outputs":[{"output_type":"stream","text":["Populating the interactive namespace from numpy and matplotlib\n"],"name":"stdout"}]},{"cell_type":"code","metadata":{"id":"kdkCEcShnOyM","colab_type":"code","colab":{}},"source":["import scipy.stats"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"sfPy5oy0nOyT","colab_type":"text"},"source":["# Introduction\n","\n","During the first lecture we have seen that the goal of machine learning is to train (learn/fit) a \n","**model** on a **dataset** such that we will be able to answer several questions about the data using\n","the model. Some useful questions are:\n","1. Predict a target $y$ for a new input $x$: predict what is in an image, recognize some audio sample, tell if the stock price will go up or down...\n","2. Generate new sample $x$ similar to those form the training dataset. Alternatively, given part of a generate the other part (e.g. given half of an image generate the other half).\n","\n","Historically, similar questions were considered by statisticians. In fact, machine learning is very similar to statistics. Some people claim that there is very little difference between the two, and a tongue-in-cheek definition of machine learning is \"statistics without checking for assumptions\", to which ML practitioners reply that they are at least able to solve problems that are too complex for a through and formal statistical analysis. \n","\n","For a more in-depth discussion I highly recommend the [\"Two cultures\"](https://projecteuclid.org/euclid.ss/1009213726) essay by Leo Breiman.\n","\n","Due to the similarity of the two fields we will today explore a few examples of statistical inference. Some of the resulting concepts (maximum likelihood, interpreting the outputs of a model as probabilities) will be used through the semester."]},{"cell_type":"markdown","metadata":{"id":"rimfq0JznOyU","colab_type":"text"},"source":["# Statistical Inference\n","\n","Consider the polling problem:\n","1. There exists **a population** of individuals (e.g. voters).\n","2. The individuals have a voting preference (party A or B).\n","3. We want the fraction $\\phi$ of voters that prefer A.\n","4. But we don't want to ask everyone (which means holding an election)!\n","\n","Instead we want to conduct a poll (choose a **sample** of people \n","and get their mean preference $\\bar\\phi$).\n","\n","Questions:\n","1. How are $\\phi$ and $\\bar\\phi$ related?\n","2. What is our error?\n","3. How many persons do we need to ask to achieve a desired error?"]},{"cell_type":"markdown","metadata":{"id":"9HATFBl5nOyV","colab_type":"text"},"source":["# Polling\n","\n","Suppose there is a large population of individuals, that support either candidate A or candidate B. We want to establish the fraction of supporters of A in the population $\\phi$.\n","\n","We will conduct an opinion poll asking about the support for each party. We will choose randomly a certain number of people, $n$, and ask them about their candidates.\n","\n","We want to use the results of the poll to establish:\n","1.  an estimate of the true population parameter $\\phi$\n","2.  our confidence about the interval"]},{"cell_type":"markdown","metadata":{"id":"TOAlBdaxTy1H","colab_type":"text"},"source":["## Sampling Model\n","First, we define a formal model of sampling. We will assume that the population is much bigger than the small sample. Thus we will assume a *sampling with replacement* model: each person is selected independently at random from the full population. We call such a sample IID (Independent Identically Distributed).\n","\n","Having the sampling model we establish that the number of supporters of A in the sample follows a *binomial distribution* with:\n","* poll size == $n$ == number of samples,\n","* fraction of A's supporters == $\\phi$ == success rate.\n","\n","For the binomial distribution with $n$ trials and probability of success $\\phi$ the expected number of successes is $n\\phi$ and the variance is $n\\phi(1-\\phi)$. \n","\n","Alternatively, the *fraction* of successes in the sample has the expected value $\\phi$ and variance $\\frac{\\phi(1-\\phi)}{n}$. \n","\n","Lets plot the PMF (Probability Mass Function) of the number of successes."]},{"cell_type":"code","metadata":{"id":"xSzyj7sEnOyW","colab_type":"code","outputId":"c202963f-b26f-4369-e3ed-035e0c67b898","executionInfo":{"status":"ok","timestamp":1570654638108,"user_tz":-120,"elapsed":633,"user":{"displayName":"Jan Chorowski","photoUrl":"","userId":"18100051538790809279"}},"colab":{"base_uri":"https://localhost:8080/","height":187}},"source":["# Poll variability check: draw samples form the binomial distribution\n","\n","n = 50\n","phi = 0.55\n","\n","# Simulate a few polls\n","for _ in range(10):\n","    sample = random.rand(n)<phi\n","    print (\"Drawn %d samples. Observed success rate: %.2f (true rate: %.2f)\" % \n","        (n, 1.0*sample.sum()/n, phi))"],"execution_count":77,"outputs":[{"output_type":"stream","text":["Drawn 50 samples. Observed success rate: 0.68 (true rate: 0.55)\n","Drawn 50 samples. Observed success rate: 0.62 (true rate: 0.55)\n","Drawn 50 samples. Observed success rate: 0.48 (true rate: 0.55)\n","Drawn 50 samples. Observed success rate: 0.56 (true rate: 0.55)\n","Drawn 50 samples. Observed success rate: 0.52 (true rate: 0.55)\n","Drawn 50 samples. Observed success rate: 0.50 (true rate: 0.55)\n","Drawn 50 samples. Observed success rate: 0.56 (true rate: 0.55)\n","Drawn 50 samples. Observed success rate: 0.46 (true rate: 0.55)\n","Drawn 50 samples. Observed success rate: 0.46 (true rate: 0.55)\n","Drawn 50 samples. Observed success rate: 0.60 (true rate: 0.55)\n"],"name":"stdout"}]},{"cell_type":"code","metadata":{"id":"06jlcRMYnOyd","colab_type":"code","outputId":"5bd53e24-8c5c-4794-bd79-f99ba28728fa","executionInfo":{"status":"ok","timestamp":1570654664011,"user_tz":-120,"elapsed":2315,"user":{"displayName":"Jan Chorowski","photoUrl":"","userId":"18100051538790809279"}},"colab":{"base_uri":"https://localhost:8080/","height":286}},"source":["# model parameters\n","n = 10\n","phi = 0.55\n","\n","# the binomial distribution\n","model = scipy.stats.binom(n=n, p=phi)\n","x = arange(n+1)\n","\n","# plot the PMF - probability mass function\n","stem(x, model.pmf(x), 'b', label='Binomial PMF')\n","\n","# plot the normal approximation\n","mu = phi * n\n","stdev = sqrt(phi*(1-phi) * n)\n","model_norm = scipy.stats.norm(mu, stdev)\n","x_cont = linspace(x[0], x[-1], 1000)\n","plot(x_cont, model_norm.pdf(x_cont), 'r', label='Norm approx.')\n","\n","axvline(mu, *xlim(), color='g', label='Mean')\n","\n","legend(loc='upper left')"],"execution_count":78,"outputs":[{"output_type":"execute_result","data":{"text/plain":["<matplotlib.legend.Legend at 0x7fec65ede048>"]},"metadata":{"tags":[]},"execution_count":78},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcjeX/x/HXZ8ZsyL4ba4mGsU6y\nJKKyhkqylD2UaBHx/VakviXaKUuy/CKSNKlESiT0DY2vXZYsMyq7imHGzPX74zqmMUZzmDNzn+Xz\nfDzOY865z32f8z6Wz1znuq/7usQYg1JKqcAQ5HQApZRSuUeLvlJKBRAt+kopFUC06CulVADRoq+U\nUgFEi75SSgUQLfpKKRVAtOgrpVQA0aKvlFIBJI/TATIqVqyYqVixotMxlPKoncd2AlC1aFWHkyh/\ntWHDhqPGmOJZ7ed1Rb9ixYqsX7/e6RhKeVSzmc0AWNFrhaM5lP8Skf3u7KfdO0opFUC06CulVADR\noq+UUgHE6/r0M5OcnEx8fDxnz551OorKRHh4OJGRkYSEhDgdRSmVBZ8o+vHx8VxzzTVUrFgREXE6\njkrHGMOxY8eIj4+nUqVKTsdRSmXBJ7p3zp49S9GiRbXgeyERoWjRovotTCkf4RNFH9CC78X070Yp\n3+FW0ReRViKyU0R2i8iITJ5/QkS2icgmEflGRCqkey5FRDa6bos8GV4pr5aUBGvXwrRpsH+/vU2a\nBN9+C6dPO51OBagsi76IBANvA62BKKCriERl2C0OiDHG1AQWAOPSPZdojKnturX3UO5cJyIMHTo0\n7fErr7zC6NGjnQukvNf330PPnlCkCDRqBA8+CPv22dvDD0Pz5lC0KLRvD0uXQmqq04lVAHGnpV8f\n2G2M2WuMSQLmAR3S72CM+dYYc8b18Acg0rMxnRcWFsbChQs5evToVR1//vx5Dye6esYYUrXQeF5c\nHNxxBzRpArGx0K0bLFwIv/xitzVpAgcPwuLF8NBDsG4dtGoF9erBd985nV4FCHeKflngYLrH8a5t\nl9MX+DLd43ARWS8iP4hIx6vI6BXy5MlD//79ef311y95bt++fTRv3pyaNWvSokULDhw4AECvXr0Y\nOHAgN910E8OHD2f06NH07NmTJk2aUKFCBRYuXMjw4cOJjo6mVatWJCcnX/La7777LjfeeCO1atXi\nnnvu4cyZMxe9dkxMDNdffz2ff/45ADNnzqRDhw40a9aMKlWq8Nxzz6VlrFq1Kj169KBGjRocPHiQ\nuXPnEh0dTY0aNXjqqacA2L9/P1WqVOHo0aOkpqbSpEkTvvrqqxz5M/UbSUnw7LNw443w00/w+uvw\n668wdSrcdRdUrAhBQfYWGQmtW9t99u+HWbPg+HFo2hT69tVuH5XjPDpkU0TuB2KApuk2VzDGJIhI\nZWC5iGw2xuzJcFx/oD9A+fLl//lNHnsMNm70ZGyoXRveeCPL3QYNGkTNmjUZPnz4RdsHDx5Mz549\n6dmzJ9OnT2fIkCHExsYCdrjpmjVrCA4OZvTo0ezZs4dvv/2Wbdu20bBhQz7++GPGjRvHXXfdxRdf\nfEHHjhf/Xrz77rt58MEHAXj66ad57733GDx4MGAL+Y8//siePXu49dZb2b17NwA//vgjW7ZsIW/e\nvNx44420bduWYsWKsWvXLmbNmkWDBg04dOgQTz31FBs2bKBw4cLccccdxMbG0rFjR5566ikeeugh\n6tevT1RUFHfccUe2/4j91u+/28K+di306AFvvgmFCrl3bGioPaZTJ3j+eXj5ZVizxn47uOGGnM2t\nApY7Lf0EoFy6x5GubRcRkduAfwPtjTHnLmw3xiS4fu4FVgB1Mh5rjJlqjIkxxsQUL57lJHGOKVCg\nAD169OCtt966aPvatWvp1q0bAA888ADff/992nP33nsvwcHBaY9bt25NSEgI0dHRpKSk0KpVKwCi\no6PZt2/fJe+5ZcsWmjRpQnR0NHPmzGHr1q1pz3Xu3JmgoCCqVKlC5cqV2bFjBwC33347RYsWJSIi\ngrvvvjstT4UKFWjQoAEA69ato1mzZhQvXpw8efLQvXt3vnN1MfTr148//viDyZMn88orr2T3j81/\nbd0K9evbRsi8ebbV7m7BTy9vXnjpJVi2zLb6GzWCVas8n1cp3GvprwOqiEglbLHvAnRLv4OI1AGm\nAK2MMYfTbS8MnDHGnBORYkBjLj7Je+XcaJHnpMcee4y6devSu3dvt/bPly/fRY/DwsIACAoKIiQk\nJG24Y1BQUKb9/r169SI2NpZatWoxc+ZMVqxYkfZcxqGSFx5fbnvGLJdz5swZ4uPjAfjrr7+45ppr\n3DouoGzaBC1aQEiILdD16mX/NVu0gB9+sN0/t99uzwu4GgVKeUqWLX1jzHngEWApsB2Yb4zZKiJj\nROTCaJzxQH7gowxDM28A1ovI/4BvgbHGmG0e/xS5qEiRInTu3Jn33nsvbVujRo2YN28eAHPmzKFJ\nkyYee78///yT0qVLk5yczJw5cy567qOPPiI1NZU9e/awd+9eqla1c7UvW7aM48ePk5iYSGxsLI0b\nN77kdevXr8/KlSs5evQoKSkpzJ07l6ZNba/cU089Rffu3RkzZkxa15JKZ9MmOwInPNyegPVEwb+g\nUiVYvRqiomy3Ubpf8kp5glt9+saYxcDiDNueTXf/tssctwaIzk5AbzR06FAmTpyY9njChAn07t2b\n8ePHU7x4cWbMmOGx93r++ee56aabKF68ODfddBN//vln2nPly5enfv36aV0x4eHhgC3o99xzD/Hx\n8dx///3ExMRc0nVUunRpxo4dy6233ooxhrZt29KhQwdWrlzJunXrWL16NcHBwXz88cfMmDGD3r17\n06ZNG6ZNm0aZMmU89vl8zoEDtvUdEWEL8rXXev49ihaFr76yJ3fbtYPly203klIeIMYYpzNcJCYm\nxmRcRGX79u3coCe2LtKrVy/atWtHp06dLto+c+ZM1q9ff9EvpdwQEH9Hp07BzTfbwr9mDVSv7vah\nV7WIyq+/QuPGkJhoh3dG+t1IaOVBIrLBGBOT1X4+Mw2DUo5KTYWuXWHHDvj44ysq+FetdGm+GTeN\n0yf+YHNMM5qPWUxs3CVjKJS6Ij4xy6a61MyZMzPd3qtXL3r16pWrWQLCSy/Bl1/CO+/AbZn2Znpc\nbFwCI/+XRIM7h/HegjEMmv8KI5PtcOGOdf7pUhmlLk+LvlJZWbHCXnzVtSsMHJhrbzt+6U4Sk1P4\n+JeeVCl9kn9tfZPVFWszPl+oFn111bR7R6l/cuQIdO3Kn+UqcnvVblQauZjGY5fnSjfLoZOJACQd\nLsALwf/ih3I1eP6rSYTv2ZXj7638lxZ9pf7JoEGkHj3G/Xc8ya5EwQAJJxMZuXBzjhf+MoUi0u6n\nSjCP3vkk5/KE8vaXr4EXzeWkfIsWfaUu56OP4KOPeLf5A/yvSHmOfx3F8a/tBLOJySmMX7ozR99+\nWMuqRIT8fTX379cUY0zrR6iW8DPoldLqKmnRd5OIcP/996c9Pn/+PMWLF6ddu3YOplI55vBhOw1y\nTAzja9lJZZMOFyDpcIG0XS50v+SUjnXK8tLd0YTlsf9NyxaKoNm/H4K774bRo2H79hx9f+WftOi7\nKV++fGzZsoXERPsffdmyZZQtqyfT/Najj8Iff8CsWZQskj/TXdJ3v+SUjnXKUqd8YRpULsrqEc3p\nWDfSjiDKlw/69IGUlBzPoPyLFv0r0KZNG7744gsA5s6dS9euXdOeO336NH369KF+/frUqVOHTz/9\nFLAzYTZp0oS6detSt25d1qxZA8CKFSto1qwZnTp1olq1anTv3h1vu1AuYC1fbidQGzkSoqIu6WYB\niAgJZljLqs7kK1nSzub5ww8wfbozGZTP8rkhm48teYyNv3l2auXapWrzRqusJ3Lr0qULY8aMoV27\ndmzatIk+ffqwyjUb4n/+8x+aN2/O9OnTOXnyJPXr1+e2226jRIkSLFu2jPDwcHbt2kXXrl25cMVx\nXFwcW7dupUyZMjRu3JjVq1dz8803e/SzqSuUlASPPAKVK4NrjYELwyN7zgni3PlUyhaKYFjLqs4O\nm+zeHd591/5iuuceu0qXUm7wuaLvpJo1a7Jv3z7mzp1LmzZtLnruq6++YtGiRWlTEZ89e5YDBw5Q\npkwZHnnkETZu3EhwcDA///xz2jH169cn0nVpfe3atdm3b58Wfae9+abtK//8czu/jovtZrH3V4xo\n7lC4dERgwgSoWxeeeQbeftvpRMpH+FzRd6dFnpPat2/Pk08+yYoVKzh27FjadmMMH3/8cdpMlxeM\nHj2akiVL8r///Y/U1NS0SdHg72mWAYKDg71qScWAlJAAzz1n165t29bpNFmrWRMGDYKJE6FfP6hz\nyVIVSl1C+/SvUJ8+fRg1ahTR0RdPHtqyZUsmTJiQ1i8fFxcHwKlTpyhdujRBQUG8//77pOiJN+/1\n9NN2/LvDazZckeees7NyDh4Mek5IuUGL/hWKjIxkyJAhl2x/5plnSE5OpmbNmlSvXp1nnnkGgIcf\nfphZs2ZRq1YtduzY4fZCJiqXbdpkV74aPNjOae8rChWCF16wc/C7Bg8o9U90amXlET7/d9SmjR0N\ns2cPFC6c6S7NmtmfV7OuyVVNrezu+54/DzVq2H7+zZshj8/12ioP0KmVlXLX8uV2Bs1//euyBd+r\n5cljZwHdsQM8uICP8k9a9FVgS02F4cOhfHk7VNNXdexoF1QfNQrOnHE6jfJiWvRVYFuwADZssP3i\n6UZW+RwRGDfOrrb15ptOp1FeTIu+ClwpKXb0S1SUvdjJ1zVubNfUfeUVO4WEUpnQoq8C14IFsG2b\n7RIJ8pP/CqNGwfHjduy+Upnwk3/pSl2h1FQYM8a28jMsLu/TYmK0ta/+kV8W/di4BBqPXU6lEV94\nbJWj4OBgateuTa1atS6aOO3QoUN0yuGisX79+kyvDUhvxYoVmU7zvGLFCgoWLEjt2rW54YYbeO65\n59K2iwjTpk1L23fjxo2ISNpUEr169aJSpUrUrl2b2rVr89Zbb3nwUznsQiv/2Wf9p5V/wahRcOKE\nnaZBqQz8bkBvbFwCIxduJjHZXvl6YZUjyN5i0hEREWzcaCd6W7p0KSNHjmTlypWUKVOGBQsWZD/4\nP4iJiSEmJsvht5fVpEkTPv/8c06fPk3t2rW58847AahRowbz58+nX79+gJ05tFatWhcdO378+Bz/\npZbrUlP/7sv3t88Gf7f2X33VXmxWoEDWx6iA4WdNnL8Xk87JVY7++OMPCrvGc+/bt48aNWoAMHPm\nTO6++25atWpFlSpVGD58eNoxc+fOJTo6mho1avCUa/ZGgPz58zNs2DCqV6/Obbfdxo8//kizZs2o\nXLkyixYtAi5uxf/44480bNiQOnXq0KhRI3budP9z5cuXj3r16rF7924AKlSowNmzZ/n9998xxrBk\nyRJat26dvT8cX5C+lR8cnPX+vkhb++oy/K7op19M2pOrHCUmJlK7dm2qVatGv3790qZZyGjjxo18\n+OGHbN68mQ8//JCDBw9y6NAhnnrqKZYvX87GjRtZt24dsbGxgJ2Hv3nz5mzdupVrrrmGp59+mmXL\nlvHJJ5/w7LPPXvL61apVY9WqVcTFxTFmzBj+9a9/uf0Zjh07xg8//ED16tXTtnXq1ImPPvqINWvW\nULdu3YsmgQMYNmxYWvfO5s2b3X4vr2WMHZ55ww3+2cq/ICbGXmX8xhs6bl9dxO+6d8oUiiAhkwKf\n3VWO0nfvrF27lh49erBly5ZL9mvRogUFCxYEICoqiv3793Ps2DGaNWtG8eLFAejevTvfffcdHTt2\nJDQ0lFatWgEQHR1NWFgYISEhREdHs2/fvkte/9SpU/Ts2ZNdu3YhIiQnJ2eZfdWqVdSpU4egoCBG\njBhB9erVWeG6pr9z587cd9997Nixg65du6adq7jA77p3liyxUxXMmuW/rfwLRoyAW26xV+kOGuR0\nGuUl/K6lnxurHDVs2JCjR49y5MiRS5670umSQ0JCEBEAgoKC0o4PCgrK9NhnnnmGW2+9lS1btvDZ\nZ59x9uzZLPM2adKEuLg4NmzYwMCBAy96rlSpUoSEhLBs2TJatGiR5Wv5vJdfhnLlIN2qZ37r5puh\nYUM7kken7VYuPtfS/+3FFzm3fcdln68DzP7rHDvC/49UY7hm7WnKFYmg2L4w9l/mmLAbqlHqCrpJ\nduzYQUpKCkWLFuWMG1+d69evz5AhQzh69CiFCxdm7ty5DB482O33S+/UqVNpa/POnDnzql4jozFj\nxnD48GGC/b3l+9//wsqV8NprEBLidJqcJ2JX/+rYEebPh27dnE6kvIDPFX13FMsfRqmi9v515Qt5\n5DUv9OmDXTBl1qxZbhfJ0qVLM3bsWG699VaMMbRt25YOHTpcVY7hw4fTs2dPXnjhBdp6aKGPRo0a\neeR1vN64cXYqYtdopYBw5532/MW4cfbbjetbpQpgxpgsb0ArYCewGxiRyfNPANuATcA3QIV0z/UE\ndrluPbN6r3r16pmMtm3bdsk25V28/u9o505jRIz597+v+iWaNrW3qzp2RlPTdMbVHZyd9zXGGDNj\nhjFgzJdfZuNFlLcD1hs36nmWffoiEgy8DbQGooCuIhKVYbc4IMYYUxNYAIxzHVsEGAXcBNQHRomI\nD85dq3zeq69CaKgdtx5ounWDyEgYO9bpJMoLuHMitz6w2xiz1xiTBMwDLuqbMMZ8a4y50Ln9AxDp\nut8SWGaMOW6MOQEsw35rUCr3/PabHa3TuzeULOl0mtwXGgqPP27PZ2zY4HQa5TB3in5Z4GC6x/Gu\nbZfTF/jyKo+9LONlK3ypv3n9383EiZCUBEOHOp3EOX37Qv78Ou2y8uyQTRG5H4gBxl/hcf1FZL2I\nrM9sGGR4eDjHjh3z/uISgIwxHDt2jHBvnYv+7FmYMgXat4frrnM6jXMKFoQ+fWDePDvnvgpY7oze\nSQDKpXsc6dp2ERG5Dfg30NQYcy7dsc0yHLsi47HGmKnAVLBr5GZ8PjIykvj4+EzHxSvnhYeHExkZ\nmfWOTvjgAzh6FB591OkkzhsyxE7L8M478PzzTqdRDnGn6K8DqohIJWwR7wJcNOBXROoAU4BWxpjD\n6Z5aCryY7uTtHcDIKw0ZEhJCpUqVrvQwFeiMsd0ZNWv+vbp4ILv2WvuNZ/Jkux5wRPauUle+Kcvu\nHWPMeeARbAHfDsw3xmwVkTEi0t6123ggP/CRiGwUkUWuY48Dz2N/cawDxri2KZXzVq6ETZtsC1fH\np1uPPWa/+XzwgdNJlEPcujjLGLMYWJxh27Pp7t/2D8dOB6ZfbUClrtqbb0LRonolanpNm0KtWnYi\ntj599JdhAPK7uXeUAmDvXvj0UxgwQLsx0hOxrf0tW+Cbb5xOoxygRV/5p4kT7SyaDz/sdBLv06UL\nlCgBr7/udBLlAC36yv/8+Se8956dL7/s1a+W5rfCw+Ghh+DLL2HPHqfTqFymRV/5n1mz7KLgOkzz\n8vr3t9+EJk92OonKZVr0lX8xBt5+G268ERo0cDqN9ypTBu66C6ZPh8TsrSqnfIsWfeVfvvsOduzQ\nvnx3PPwwHD8OH37odBKVi7ToK/8yZYqdM79zZ6eTeL+mTSEqyn4zUgFDi77yH4cPw4IF0KMH5M3r\ndBrvJ2Jb++vXw48/Op1G5RIt+sp/zJwJycl2bL5yzwMP2Nk333nH6SQql2jRV/4hNRWmToVbbrFd\nFso9BQrYb0bz5tnpGZTf06Kv/MM339gx59rKv3IPPwznzsGMGU4nUblAi77yD5MnQ7FicM89Tifx\nPdWr25O6kyZBSorTaVQO06KvfN+hQ3aend69ISzM6TS+adAg+OUXWLLE6SQqh2nRV75v+nTbQu3f\n3+kkvqtjRyhdWodvBgAt+sq3paTYE7i33RbYyyFmV0gI9OtnW/r79jmdRuUgLfrKty1ZAgcPwsCB\nTifxff362bH706Y5nUTlIC36yrdNngylStllAFX2lC8PrVvb7rLkZKfTqByiRV/5rgMHYPFi6NvX\ndk+o7BswAH79FT7/3OkkKoe4tVyiUl5p2jQ7q+aDDzqdxH+0bg2RkTBlCrEV6zN+6U4OnUykTKEI\nhrWsSsc6uj6Br9OWvvJNycm26LduDRUqOJ3Gf+TJA/36Yb76ionvLSPhZCIGSDiZyMiFm4mNS3A6\nocomLfrKN332me2G0BO4nte3L6kId21YzPGvozj+tZ3WIjE5hfFLdzocTmWXFn3lm6ZMgXLloE0b\np5P4n8hIll97I503LyP19wiSDhdIe+rQSV1wxddp0Ve+Z88e+OorO8QwONjpNH5paeMOFD99kjsT\nv7xoe5lCEQ4lUp6iRV/5nqlTbbHv29fpJH7r5ke6c6hACfr9NSttW0RIMMNaVnUwlfIELfrKZ8TG\nJdDshSUcnTiFldUaEnvY6UT+q2NMeU5170GLcyupnLyXsoUieOnuaB294we06CufEBuXwMiFm4n+\ncTnFzpzivajbdTRJDrvh6cdJIZinQxawekRzLfh+Qou+8gnjl+4kMTmFzt9+xy+h5VhVqY6OJslp\nZcqwpuidtP5tBiQlOZ1GeYgWfeUTDp1M5NpjB2ny13+ZFtEbI0Fp21XO+azMAAonH4FPPnE6ivIQ\nLfrKJ5QpFEG3jUtIIoT/y9f1ou0q56wrfAe/hle0J8+VX9Cir3zCU03L02nLN8RGtOVIcHFAR5Pk\nBiNBfF7qQVi+HHbtcjqO8gC3ir6ItBKRnSKyW0RGZPL8LSLyk4icF5FOGZ5LEZGNrtsiTwVXgaX9\nrjUUPPsXswr2AdDRJLnoy9J97PQM2tr3C1lOuCYiwcDbwO1APLBORBYZY7al2+0A0At4MpOXSDTG\n1PZAVhXIJk+GqlU5XbIdDURYMaK504kCxvHQUtChA8ycCS+8oEtS+jh3Wvr1gd3GmL3GmCRgHtAh\n/Q7GmH3GmE1Aag5kVIFu0yZYu9ZO+yvidJrANGAAHD0KCxc6nURlkztFvyxwMN3jeNc2d4WLyHoR\n+UFEOl5ROqXAzrMTFgY9ezqdJHC1aAGVK9tvXMqn5caJ3ArGmBigG/CGiFybcQcR6e/6xbD+yJEj\nuRBJ+Yy//oL334f77oMiRZxOE7iCguzC8999B9u3O51GZYM7RT8BKJfucaRrm1uMMQmun3uBFUCd\nTPaZaoyJMcbEFC9e3N2XVoFg7lz480/bvaCc1bu3XaFMT+j6NHeK/jqgiohUEpFQoAvg1igcESks\nImGu+8WAxsC2fz5KqXQmT4boaGjY0OkkqkQJuOcee0I3US+K81VZFn1jzHngEWApsB2Yb4zZKiJj\nRKQ9gIjcKCLxwL3AFBHZ6jr8BmC9iPwP+BYYm2HUj1KXt349/PSTXShFT+B6hwED4ORJmD/f6STq\nKrm1Rq4xZjGwOMO2Z9PdX4ft9sl43BogOpsZVaCaPBny5oXu3Z1Ooi5o2hSqVrUn1/XEuk/SK3KV\ndzp1yvbnd+sGBQs6nUZdIGK/ea1dC//7n9Np1FXQoq+80+zZcOaMroHrjXr0sENop0xxOom6Clr0\nlfcxBiZNgnr17E15lyJF7BDa2bPtkFrlU7ToK++zejVs3QoPPeR0EnU5AwfaobRz5zqdRF0hLfrK\n+0yaZPvxu3RxOom6nAYN7FDaSZPsNzPlM7ToK+9y5AgsWGD7jfPlczqNupwLJ3Tj4uzQWuUztOgr\n7zLDtTSfXoHr/bp3t0Nq9YSuT9Gir7xHaqotILfcAtWrO51GZaVgQTukdu5ce8GW8gla9JX3WLYM\n9u7VE7i+ZOBAO7R29mynkyg3adFX3mPyZCheHO66y+kkyl0XhtVOmaIndH2EFn3lHeLjYdEi6NtX\nV2byNQMHwpYtsGaN00mUG7ToK+8wbZptKfbv73QSdaW6dIECBXSBFR+hRV85LzkZ3n0XWrWCSpWc\nTqOuVP78cP/98NFHcOyY02lUFrToK+d9/jkcOqTz7PiyAQPg3DmYNcvpJCoLWvSV8yZNgnLloG1b\np5Ooq1WzJjRqZLt49ISuV9Oir5y1e7cdqvnggxAc7HQalR0DBsCuXbB8udNJ1D/Qoq+cNWUK5MkD\n/fo5nURlV+fOULQoTJzodBL1D7ToK+ecOQPvvQcdO0Lp0k6nUdkVHm6/sS1aBPv3O51GXYYWfeWc\nDz6AEydg8GCnkyhPuXA19aRJzuZQl6VFXznDGJgwwZ4AbNLE6TTKU8qXhw4d7BDcxESn06hMaNFX\nzli1CjZtsq18EafTKE8aPBiOH4d585xOojKhRV85Y8IEKFzYztKo/EuzZnaW1AkTdPimF9Kir3Lf\nwYPwySd2xE7evE6nUZ4mAo88YhdYWbvW6TQqAy36KvdduIDn4YedTqJyyv332/n2J0xwOonKQIu+\nyl1nz8LUqXDnnVCxotNpVE7Jnx9697ZLX/76q9NpVDpa9FXu+vBDOHpUh2kGgkGDICVFl1P0Mlr0\nVe4xBt56C6KioHlzp9OonHbdddC6tS36SUlOp1EuWvRV7lm7Fn76yZ7k02GageGRR+C332w3j/IK\nWvRV7nnrLXty74EHnE6ickvLllC1Krz2mg7f9BJa9FXu2L/ftvYefNCe5FOBISgIHn8cNmywF+Qp\nx7lV9EWklYjsFJHdIjIik+dvEZGfROS8iHTK8FxPEdnluvX0VHDlY95803bpDBnidBKV2x54wM6+\n+eqrTidRuFH0RSQYeBtoDUQBXUUkKsNuB4BewAcZji0CjAJuAuoDo0SkcPZjK59y6pRdA7dzZ7tY\nigosefPaazI++wx+/tnpNAHPnZZ+fWC3MWavMSYJmAd0SL+DMWafMWYTkJrh2JbAMmPMcWPMCWAZ\n0MoDuZUvefdd+PNPGDrU6STKKYMGQUgIvP6600kCnjtFvyxwMN3jeNc2d2TnWOUPkpNt106zZlC3\nrtNplFNKlrRX6c6aZa/TUI7xihO5ItJfRNaLyPojR444HUd50oIFEB8PTzzhdBLltCeesNMtT57s\ndJKA5k7RTwDSd8RGura5w61jjTFTjTExxpiY4sWLu/nSyusZY0/eVa2qi54rO/Nmq1Z2OcWzZ51O\nE7DcKfrrgCoiUklEQoEuwCIA80BYAAAUsElEQVQ3X38pcIeIFHadwL3DtU0Fgu++s0P1Hn/cDt1T\nauhQ+P13u2qackSW/xONMeeBR7DFejsw3xizVUTGiEh7ABG5UUTigXuBKSKy1XXsceB57C+OdcAY\n1zYVCF59FYoVgx49nE6ivEWLFna1NL1YyzF53NnJGLMYWJxh27Pp7q/Ddt1kdux0YHo2MipftGWL\nHaI3ahRERDidRnkLEdva79kTvvgC2rVzOlHA0e/cKmeMHQv58ulsmupSXbtChQrwn/9oa98BWvSV\n5+3dC3PnwsCB9kpMpdILCYHhw+GHH2DlSqfTBBwt+srzxo2DPHl0mKa6vN697dj9F190OknA0aKv\nPOvQIZgxw/6nLlPG6TTKW0VE2EbBsmWwbp3TaQKKFn3lWa+/DufP26/vSv2TgQOhUCF46SWnkwQU\nLfrKc44fh0mT7Im6ypWdTqO8XYEC9kT/J5/Atm1OpwkYWvSV50yYAKdPw4hLZt9WKnNDhthZOF9+\n2ekkAUOLvvKMU6fgjTegfXuoUcPpNMpXFCsGAwbAnDmwe7fTaQKCFn3lGW+8ASdPwujRTidRvmb4\ncAgNheefdzpJQNCir7LvxAl7Wf1dd0GdOk6nUb6mVCm7yMrs2bBzp9Np/J4WfZV9r70Gf/yhrXx1\n9YYPh/BwGDOG2LgEGo9dTqURX9B47HJi49yd1Fe5w625d5S6rGPHbNfOvffaibSUuholSsDgwZhx\n45h2zS0kFLJTeSWcTGTkws0AdKyj6y95grb01RVL3xL7v06DMadP24nVlMqOJ58kMSScgSvncPzr\nKI5/bZfiTkxOYfxS7fbxFC366orExiUwcuFmEk4mUvjMKe5Z/QmLb7iF2KRCTkdTvq5YMabXa0+7\nHauoEh9P0uECaU8dOpnoYDD/okVfXZHxS3eSmJzC8a+j6PPhN4SfT+K1hl20JaY84vPbuvBHaF5G\n/XHxVbplCun03J6iRV9dkQstrtKHTvLg4dl8XKM5e4qV05aY8oiBd93I9EadaJ/4JQ3P/ReAiJBg\nhrWs6nAy/6FFX12RCy2u0adeJFWCeO3m+y/arlR2dKxTlmv/8zS/BpfixZOjKVswnJfujtaTuB6k\nRV9dkWEtq1L32C90PbOAifkH8FuBYtoSUx51Z6MqzL72eRomrWN11VNa8D1Mi766Ih3rlGXylo84\nHlSYVwoMoWyhCG2JKY9bUqoX+/LeACNHQnKy03H8ihZ9dWWWLaPED9/xQaVnuOG6yqwe0VwLvvK4\nFMnDlEovw88/w3vvOR3Hr2jRV+5LSbFXTlasyKdlHnY6jfJza4u2gyZN7JXef/3ldBy/oUVfuW/a\nNNi4EcaOJTkozOk0yt+JwPjx8PvvuqyiB2nRV+45cQL+/W+45Rbo3NnpNCpQ3HQTPPAAvPqqTr3s\nIVr0lXtGj7aF/803bQtMqdzy8st26uUnnnA6iV/Qoq+ytmULvP029O8PtWs7nUYFmtKl4dln4bPP\n4MsvnU7j87Toq39mDDz2GFxzjS5yoZzz6KNw/fX2Z1KS02l8mhZ99c8++gi++QbGjLFL2ynlhNBQ\n27W4axe8/rrTaXyaFn11eSdO2IWr69WDhx5yOo0KdK1aQYcO8NxzsHev02l8lhZ9dXkjRsCRIzB1\nKuTR9XaUF5g40f5bHDjQdj2qK6ZFX2Xu++9tsX/sMahb1+k0SlmRkfDSS7BsGcyZ43Qan+RW0ReR\nViKyU0R2i8iITJ4PE5EPXc//V0QqurZXFJFEEdnouk32bHyVI86dsyN1KlSwX6WV8iYPPQQNG8Lj\nj8PRo06n8TlZFn0RCQbeBloDUUBXEYnKsFtf4IQx5jrgdeDldM/tMcbUdt0Geii3ykljx8L27fDO\nO5A/v9NplLpYUJD9FnrqlI7dvwrutPTrA7uNMXuNMUnAPKBDhn06ALNc9xcALUT0Ch6ftGEDvPAC\ndOsGbdo4nUapzNWoYc85vf8+fPqp02l8ijtFvyxwMN3jeNe2TPcxxpwHTgFFXc9VEpE4EVkpIk2y\nmVflpLNnoUcPKFHCnjBTyps9/bS9WPDBB+HwYafT+IycPpH7K1DeGFMHeAL4QEQKZNxJRPqLyHoR\nWX/kyJEcjqQu65lnYNs2O5Vt4cJOp1Hqn4WG2pb+qVP2HJSO5nGLO0U/ASiX7nGka1um+4hIHqAg\ncMwYc84YcwzAGLMB2ANcn/ENjDFTjTExxpiY4sWLX/mnUNm3apWd1GrAADseWilfUKOGnYHz009h\n5kyn0/gEd4r+OqCKiFQSkVCgC7Aowz6LgJ6u+52A5cYYIyLFXSeCEZHKQBVAr6rwNsePQ/fuUKkS\nvPKK02mUujKPPw5Nm9opGvSirSxlWfRdffSPAEuB7cB8Y8xWERkjIu1du70HFBWR3dhunAvDOm8B\nNonIRuwJ3oHGmOOe/hAqG4yB3r3ht99g3jwdraN8T1AQzJplf953nx1yrC7LrcssjTGLgcUZtj2b\n7v5Z4N5MjvsY+DibGVVOmjABFi2y85nceKPTaZS6OhUqwIwZcPfdMGwYvPWW04m8ll6RG8g2bIAn\nn4T27e1XY6V82V132SvIJ0yABQucTuO1tOgHqiNHbKuoVCmYPl0XRlH+4eWXoX596NvXzsipLqFF\nPxAlJ8O999qxzQsXQtGiWR+jlC8IDYX58+2kbB062OGc6iJa9APRE0/AypXw7rsQE+N0GqU8q0IF\n272za5e9sjwlxelEXkWLfqCZNs1ebTt0KNx/v9NplMoZt95qT+YuXgwjRzqdxqvoJOmB5Msv7Tzk\nLVvaSdWU8mcPPWTXdx4/HqpVgz59nE7kFbToB4r1620/fs2adglEXRRFBYI33oDdu+00DSVKQLt2\nTidynHbvBIK9e6FtW7vG7Rdf2EXOlQoEISG2f79OHejcGdascTqR47To+7v4eLj9djh/HpYsgdKl\nnU6kVO665hrb2ImMtC39rVudTuQoLfr+7NdfoXlzOyZ/8WLbr6lUICpRApYuhfBwaNHCziYboLTo\n+6vff7cF/9Ah28K/6SanEynlrEqVYPlyeyFi8+YBW/i16Puj+Hg7ZO3AAdvCb9TI6URKeYdq1eDb\nb/8u/AHY1aNF39/8/DM0bmwL/+LFcMstTidSyrukL/xNmgTcyV0dt+fDYuMSGL90J4dOJlKmUAQv\nlE/i1sd62CdXrrQjFpRSl6pWzRb7O+6A226DDz+EO+90OlWu0Ja+j4qNS2Dkws0knEzEANV/XE79\nXndxJk8YfP+9FnylslKpEqxeDdWr2xk6p0xxOlGu0Ja+jxq/dCeJySkcX3YDQ397h1GHXmVj6et5\nttcLLLr+khUplVKZKVHCdvV07gwDB/LL16vpVas7B/46T5lCEQxrWZWOdco6ndKjtKXvow6dTCTf\nuTNM2z6UUYde5ZOoZtzX9SU2p+Z1OppSviV/fvjsM37u+RCVFrzPuClDKXr6BAknExm5cDOxcRmX\nBPdtWvR9VLPTB/l81qPck7iIpws+w+PthnIuJIwyhSKcjqaU7wkOpne1Tgy5cxjRCXv4fPIT3PxL\nHInJKYxfutPpdB6lRd/XpKbCW28xbcqjRJxPomXxWF4p8CiIEBESzLCWVZ1OqJRPOnQykUVRTbml\nxFKOU4TZ859hxLfTOXLsD6ejeZQWfV+yaxc0awaPPkrwHbfzU+y3rMt/MwBlC0Xw0t3Rftf/qFRu\nufAteUtodRqXXMbs2q0Z+ONCvpjzJKxb53A6z9Gi7wvOn4dXXrEzZG7ebJc3/Owz2jSPpk75wjSo\nXJTVI5prwVcqG4a1rEpESDAAiUF5ebrlIB7uPIrI839BgwZ28aHTpx1OmX1a9L3d0qVQqxYMG2bH\nFG/dCr1765q2SnlYxzpleenuaMLy2LJYtlAEd4x4kIhdO+3UzK+/DlFRMHcuGONw2qunRd9bbdtm\nZwRs1QrOnYNPPoHYWChTxulkSvmtjnXKXvrtuWBBmDQJvvsOihSxSzA2aGDH+PsgLfreZutW6NIF\natSw/8jGjbPbOnbU1r1STmrSxC5GNGOGnebk5puhTRufm8ZBi743MMZeRXvvvRAdbef+HjHCLn4y\nbBiEhTmdUCkFEBwMvXrZOa5efNGe4G3c2E5wuGSJHV3n5bToOykx0bYa6tWzrYivv7aLOO/bZ/9B\nFSvmdEKlVGby5fv7/+prr8HOndC6NVStah+fOOF0wsvSop/bUlPtZd99+kDJkvZncrKd9yM+Hv7z\nHyha1OmUSil35MsHjz8Ov/wCH3wApUrB0KH23FvXrvDZZ5CU5HTKi+jcO7nh3Dk76+Vnn9mTsfHx\ndgm3Tp2gZ087/bH21yvlu8LCbJHv2hU2bbKNuPnzYd48KFwY7rnHDsxo0cJO++AgLfrZlHF642Et\nq9Kxdhk7+mblSrtSz9Kl8NdfEBFh16sdPx7at4e8Ok+OUn6nZk14+2144w3bZfvBB3bq5mnTIDTU\nXmDZurVt7NWqZc8T5CIt+tlwYXrj0D9O0vD3PdT4fQ95Z/3Mud+2E3byuN2pXDno3t3O1d28uS38\nSin/FxJii3vr1raL5/vv7SCNL76wXUJgv/E3asS2a2sx5VwxVuQvT/5SxXN0dk+3ir6ItALeBIKB\nacaYsRmeDwP+D6gHHAPuM8bscz03EugLpABDjDFLPZY+tx0/bs/au26FY1fyVfwuyp36PW2XAwVL\n8s21MbR5pIv9jV6xonbdKBXoQkNto695c3j1VTh4EFatglWr+GPZt0QtXcqbwPbiFWndZyIjF24G\nyJHCn2XRF5Fg4G3gdiAeWCcii4wx6VcV7gucMMZcJyJdgJeB+0QkCugCVAfKAF+LyPXGmBRPf5BM\nu1nc/QNLSrIF/dgxOHzY9rmnvyUk2LP0x479fUyePJQrUJKNpa9nEgPYGFqLQ13PcTKiAAL80rut\npz+iUspflCtnL/Lq1o3WY5fz529HqLg4iTBjT/pemN3TkaIP1Ad2G2P2AojIPKADkL7odwBGu+4v\nACaKiLi2zzPGnAN+EZHdrtdb65n4VmxcAs/O30DNvZuonnyWiOSz/LQqlmurFSa6UB47X8bp07Zf\n/fRpO5zqQpE/dsxuz0zhwhAZaW8xMXD99X/fKlXigVdXkXAykd8+aABAqYgfAHR6Y6WU2w6dTMSE\n5+ers646wg9p23OCO0W/LHAw3eN44KbL7WOMOS8ip4Ciru0/ZDjW47+6xi/dSZ7TfzF7/jMXP7HY\n/jBBQaTmyYNx3VJDQ0kJCyM1LIzUyEhSw8Ls49AwUsLDScmbl5S8eTF50v3xnD4DcRvtzeWtv86x\n9+hpzobNBiB01R8EBQmVi+Vj/wMzPf0xM/W468vH/gdy5e0C9n2z+95ny+5wHdsjV983u/TvOee9\nfuAk586nkBQ2m52pVZjtKq851Xj0ihO5ItIf6A9Qvnz5Kz7+0MlEgsLz8/5N7SiReILzwcGkBOXh\nfHAw9a4rAUE5czlCsfz2StmfT6SSkmoIyxNMuSIRadtzg1OjvwLtfZ18b/3M/v3e5YpEsPfoaSQk\nheBg272Tk2tjuFP0E4By6R5HurZltk+8iOQBCmJP6LpzLMaYqcBUgJiYmCuevq5MoQgSTibyTLOB\nF20vWyiC1SOaX+nLXZEK2LPXTqmg7+sT7x0+s5l9jbH/l6vvm13695w773cwLoE3XOcky+bw2rzu\nFP11QBURqYQt2F2Abhn2WQT0xPbVdwKWG2OMiCwCPhCR17AncqsAP3oq/AXDWlZl5MLNJCb/fX5Y\nV5FSSvmKjnXK5tp6GFkWfVcf/SPAUuyQzenGmK0iMgZYb4xZBLwHvO86UXsc+4sB137zsSd9zwOD\ncmLkzoU/rKsevaOUUgFCjJctBhATE2PWr1/vdAylPKqZq3tnRa8VjuZQ/ktENhhjYrLaTydcU0qp\nAKJFXymlAogWfaWUCiBa9JVSKoBo0VdKqQDidaN3ROQIsD8bL1EMOOqhOL4i0D5zoH1e0M8cKLLz\nmSsYY4pntZPXFf3sEpH17gxb8ieB9pkD7fOCfuZAkRufWbt3lFIqgGjRV0qpAOKPRX+q0wEcEGif\nOdA+L+hnDhQ5/pn9rk9fKaXU5fljS18ppdRl+E3RF5FWIrJTRHaLyAin8+Q0ESknIt+KyDYR2Soi\njzqdKbeISLCIxInI505nyQ0iUkhEFojIDhHZLiINnc6U00Tkcde/6y0iMldEwp3O5GkiMl1EDovI\nlnTbiojIMhHZ5fpZ2NPv6xdFP93i7a2BKKCra1F2f3YeGGqMiQIaAIMC4DNf8Ciw3ekQuehNYIkx\nphpQCz//7CJSFhgCxBhjamCndO/ibKocMRNolWHbCOAbY0wV4BvXY4/yi6JPusXbjTFJwIXF2/2W\nMeZXY8xPrvt/YguB3y8gICKRQFtgmtNZcoOIFARuwa5ZgTEmyRhz0tlUuSIPEOFaiS8vcMjhPB5n\njPkOu/5Ieh2AWa77s4COnn5ffyn6mS3e7vcF8AIRqQjUAf7rbJJc8QYwHEh1OkguqQQcAWa4urSm\niUg+p0PlJGNMAvAKcAD4FThljPnK2VS5pqQx5lfX/d+Akp5+A38p+gFLRPIDHwOPGWP+cDpPThKR\ndsBhY8wGp7PkojxAXWCSMaYOcJoc+MrvTVz92B2wv/DKAPlE5H5nU+U+Y4dWenx4pb8UfbcWYPc3\nIhKCLfhzjDELnc6TCxoD7UVkH7YLr7mIzHY2Uo6LB+KNMRe+xS3A/hLwZ7cBvxhjjhhjkoGFQCOH\nM+WW30WkNIDr52FPv4G/FP20xdtFJBR70meRw5lylIgItp93uzHmNafz5AZjzEhjTKQxpiL273i5\nMcavW4DGmN+AgyJS1bWpBXbNaX92AGggInld/85b4Ocnr9NZBPR03e8JfOrpN8hyYXRfcLnF2x2O\nldMaAw8Am0Vko2vbv4wxix3MpHLGYGCOq0GzF+jtcJ4cZYz5r4gsAH7CjlKLww+vzhWRuUAzoJiI\nxAOjgLHAfBHpi51tuLPH31evyFVKqcDhL907Siml3KBFXymlAogWfaWUCiBa9JVSKoBo0VdKqQCi\nRV8ppQKIFn2llAogWvSVUiqA/D+LPErx1DwAGQAAAABJRU5ErkJggg==\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[]}}]},{"cell_type":"markdown","metadata":{"id":"p2xTO4CxnOyj","colab_type":"text"},"source":["## Parameter Estimation\n","\n","In Statistics and Machine Learning we only have access to the sample. The goal is to learn something useful about the unknown population. Here we are interested in the true heads probability $\\phi$.\n","\n","The MLE (Maximum Likelihood Estimator) for $\\phi$ is just the sample mean $\\bar\\phi$. However, how precise it is? We want the (sample dependent) confidence interval around the sample mean, such that in 95% of experiments (samples taken), the true unknown population parameter $\\phi$ is in the confidence interval.\n","\n","Formally we want to find $\\bar\\phi$ and $\\epsilon$ such that $P(\\bar\\phi-\\epsilon \\leq \\phi \\leq \\bar\\phi + \\epsilon) > 0.95$ or, equivalently, such that $P(|\\phi-\\bar\\phi| \\leq \\epsilon) > 0.95$.\n","\n","Note: from the sampling model we know that for a large enough sample (>15 persons) the random variable denoting the sample mean $\\bar\\phi$ is approximately normally distributed with mean $\\phi$ and standard deviation $\\sigma = \\sqrt{(\\phi(1-\\phi)/n)}$. However we do not know $\\phi$. When designing the experiment, we can take the worse value, which is 0.5. Alternatively, we can plug for $\\phi$ the estimated sample mean $\\bar\\phi$. Note: we are being too optimistic here, but the error will be small.\n","\n","For a standard normal random variable (mean 0 and standard deviation 1) 96% of samples fall within the range $\\pm 1.96$. \n","\n","Therefore the confidence interval is approximately $\\bar\\phi \\pm 1.96\\sqrt{\\frac{\\bar\\phi(1-\\bar\\phi)}{n}}$.\n"]},{"cell_type":"code","metadata":{"id":"j0sweWStnOyl","colab_type":"code","outputId":"f0f5999e-4f7e-4447-a64e-94b24255780b","executionInfo":{"status":"ok","timestamp":1570654946326,"user_tz":-120,"elapsed":3665,"user":{"displayName":"Jan Chorowski","photoUrl":"","userId":"18100051538790809279"}},"colab":{"base_uri":"https://localhost:8080/","height":336}},"source":["phi=0.55\n","n=100\n","n_experiments=1000\n","\n","samples = rand(n_experiments, n)<phi\n","phi_bar = samples.mean(1)\n","\n","hist(phi_bar, bins=20, label='observed $\\\\bar\\\\phi$')\n","axvline([phi], color='r', label='true $\\\\phi$')\n","title('Histgram of sample means $\\\\bar\\\\phi$ form %d experiments.\\n'\n","      'Model: %d trials, %.2f prob of success'%(n_experiments,n,phi))\n","legend()\n","xlim(phi-0.15, phi+0.15)\n","\n","confidence_intervals = zeros((n_experiments, 2))\n","confidence_intervals[:,0] = phi_bar - 1.96*np.sqrt(phi_bar*(1-phi_bar)/n)\n","confidence_intervals[:,1] = phi_bar + 1.96*np.sqrt(phi_bar*(1-phi_bar)/n)\n","\n","#note: this also works, can you exmplain how the formula works in numpy?\n","confidence_intervals2 = phi_bar[:,None] + [-1.96, 1.96] * np.sqrt(phi_bar*(1-phi_bar)/n).reshape(-1,1)\n","assert np.abs(confidence_intervals-confidence_intervals2).max()==0\n","\n","good_experiments = (confidence_intervals[:,0]<phi) & (confidence_intervals[:,1]>phi)\n","\n","print (\"Average confidence interval is phi_bar +-%.3f\" \n","       % ((confidence_intervals[:,1]-confidence_intervals[:,0]).mean()/2.0,))\n","\n","print (\"Out of %d experiments, the true phi fell into the confidence interval %d times.\"\n","       % (n_experiments, good_experiments.sum()))"],"execution_count":96,"outputs":[{"output_type":"stream","text":["Average confidence interval is phi_bar +-0.097\n","Out of 1000 experiments, the true phi fell into the confidence interval 961 times.\n"],"name":"stdout"},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAYIAAAEdCAYAAAABymAfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu8VXWd//HXW0BJBUXAG5BgiYp4\nI0Tn52UsMkkd0VKhyRGLxkrLcZxUujzSmvxF8zMzJ3WG8q6pDKlpWg1p3iosQFAUL4jHOIiCeAm8\n4O3z+2N9Dy4O+5yzz9l7n33OWe/n47EfZ+/vd10+37XWWZ+9Lnt9FRGYmVlxbVLvAMzMrL6cCMzM\nCs6JwMys4JwIzMwKzonAzKzgnAjMzArOicDMrOCcCMzMCs6JwMys4JwIzCog6XZJa9Pr7nrHY9YR\n8iMmzMyKrcccEUh6VNKh9Y6jliTtKmmBpDWSTq93PKVIapD08XrH0dV0h3VXNEXYZ5SrWySCUjsX\nSSdLeqDpc0TsERH3tHc63czZwO8jol9EXFzvYCwj6VJJ325jsLquO0lfkTRX0jpJV5Wo30bSLZJe\nk/SspH8sp66c+q6qnH1GLXTF/VC3SARdgaTe9Y4B2Al4tN5B2EYOAh5oY5gOrbsqbnfPAd8Drmih\n/hLgLWA74LPAZZL2KKOunPoupYv8L3ctEdHlX0AD8PFmZScDD5QaBjgHWA6sAZ4AxgPXAu8BbwBr\ngbPTsGOAh9Kw/wPcBHwvN81zgIeBdUBvYBrwdBr+MeDYErGelcZ5Dbic7B/k12mc3wEDWmjn7sA9\nwCtkO42jc3V3A+8Cb6b4R5YYf6N2p/IWY25vvGn4r6fpvAxcCfRtYT3sCPwCWAU8A5zexjpuTxwt\nTruM9n4tzefVtL77trUMS8S7Cdm3/AYggJXAmS0Mu9G6a2NdN7Dxdteu5dPKcv4ecFWzsi3IduQj\nc2XXAtNbq2tr3Bbm39p6+xDwEjAmN+wq4NC2tr22trVWlunHc/XV2v4aaGEbo8R+iDK3uZruYzt7\nhh0Ksh2JANgVWAbsmMqHAx8qNR1gU+BZ4F+APsCn0kadTwQLgGHAB1LZ8Wkj2ASYlDaaHZrFMSdt\nREPIdhDzgX2BvmQ7hXNLtLEPsAT4RorrY2nD2DU3zD3AF1pYRq21u8WY2xtvGn5RWibbAH9oWl7N\n1sMmwDzg26k9OwNLgcNbWcdlxdHWtMto759T/TbAYuBLbS3DEvF+E7gXOJrsn3cc2c5+WAvDr193\nba1rSm937VpPrfwvlUoE+wKvNyv7GnB7a3VtjVti3m1uE8A/k+3oNwd+C1zQ1rZX5nRbWqb5RFCt\n7a+BFraxEvMte5ur5as7nRq6VdIrTS/g0haGexfYDBglqU9ENETE0y0MewDZN4OLI+LtiLiZbAXm\nXRwRyyLiDYCI+J+IeC4i3ouIm4CnyHYCef8ZES9ExHLgfuDBiHgoIt4EbiHbuErFsiXZN6m3IuJu\n4FfAZ1pZJmW1u4yY2xvvT9IyeQk4v4UY9wMGR8R3U3uWAj8FJrfShnLjaHXaZbT34lT/EtnObp+2\nlmGepH5kiWAK2Y7loYj4M9k/9G6ttK9JOet6g+2uncunvbYE/tas7FWgXxt1bY3bXJvbRET8lCxJ\nPgjsQLac80pte+Vua6WWaV5Vtr/cvEptY821Z39VM90pERwTEVs3vYBTSw0UEUuAM4DzgJWSbpS0\nYwvT3BFYHikVJ8uaDbPBZ0knpbs/mhLSaGBQs3FeyL1/o8TnLVuIZVlEvJcre5bs20mbWmt3GTG3\nN978Mnk2xd7cTsCOzZL3N8i+cbWk3DhanXYZ7X0+9/71pum2Y9v5GPBERDQAewMPSdoEGED2TbIt\n5azr5tshdGy7KsdaoH+zsv5kRymt1bU1bnPlbhM/JVtn/xkR65rVldr2yp1uqWWaV5XtLym5jTXX\nzv1VzXSnRFC2iPh5RBxEtsIC+EFTVbNBVwBDJClXNqz55JreSNqJbCP9CjAwJaRFgKjcc8CwtENp\n8kGyc4dlKdXuGsWcX0YfTLE3twx4Jp+8I7tj5ogK5tvmtCttbyvbTt72ZOeyIfum9xBwCNk34YfL\nmE0567ozf+DzJNBb0i65sr3Jrl20VtfWuM21uU1I2hK4iOwc/XmStmk2jVLbXrnbWrWWaaXb9gZx\nlLnN1VSPSwTpfu2PSdqM7OLcG2QXZyDL8DvnBv8T2aHZVyT1ljSRjU/z5G1BtqJWpXl9juybSzU8\nSPbN4WxJfdL9zf8A3FjOyK20uxYxnyZpaPon/SbZxbDm/gyskXSOpA9I6iVptKT9Kpx3W9PucHvb\n2HbyHgc+ImnnNO2XgJ+Q3YBQzs6monXdEWn77gv0AnpJ6tt090xEvAbcDHxX0haSDgQmAte2VtfW\nuCXCKGeb+DEwNyK+ANwB/FezaZTa9mq5rZVS6fzW74fasc3VVI9LBGTn26YDL5Idnm1LdqcBwPeB\nb6XDua9FxFtkF4inkt29cSLZudrmh6MARMRjwA/JEsgLwJ5kF6wqlmL5B+CTKfZLgZMi4vEyJ1Gy\n3TWK+efA/5JdIHua7ILdBiLiXeAosm/Mz6S4fgZsVeG8W512he1tbdvJu4/sG+s8sguJ1wM/iIiy\nduRVWNcd8S2yncw0su38jVTW5FTgA2Sntm4AvhwRj5ZRV0490PY2kb6ITQC+nEY5Exgj6bO5yWy0\n7dVyWyulCvNbvx8iu5mhxW1O0q8lfaN60ZfmR0w0I+lB4L8i4sp6x9IVSWogu/vld/WOpd4kTQI+\nFxET6h1LEXjbq52eeETQLpL+XtL26dB5CrAX8Jt6x2Xdwq5kt46adWv+hV32zzyT7NzyUuC4iFhR\n35Csm9gV+GO9gzCrlE8NmZkVXOFPDZmZFZ0TgZlZwTkRmJkVnBOBmVnBORF0E5KGSwqV8Sx1Neu0\np7tRO3qOSsvkwzUOqVuQdJWkjX7c18FpuUe1AnEiqAFlPRC9JWlQs/KH0o5reH0iK03SDElPSHpP\n0skl6v9V0vOS/ibpivRz+Ka64ZJ+L+l1SY+rlZ6Xyt1RRSf3HKXMDyStTq8fSCr5bCJJh6bltDb3\nmpKrv0fSm7m67vo7A/eGVyBOBLXzDLnHCkvak+wZ613RQrLHBMxvXiHpcLLHEowneyjWzsB3coPc\nQPbQtYFkz36ZJWlwR4Io52inRk4BjiF7WNpeZI9/+GIrwz8XEVvmXlc3q/9Krm7XagRYh2WzE+4N\nrzCcCGrnWuCk3OcpwDX5ASRtJekaSauU9fX6LaUnUqYHWV0g6UVJS4EjS4x7uaQVkpZL+p6kXh0J\nNCIuiYi7yB561dwU4PKIeDQiXgb+naxTICSNJOvh7dyIeCMifgE8Any6+UQknULWjeHZ6Zvy7am8\nIT2862HgtfQL7/V9ukoaJ+lP6flQKyT9RNKmpdoh6QhJj6XTGcslfa3MRTAF+GFENKZn0f+wqY21\nJOk8SbMk3ZRini9p71x9qWWzezrqeCWdQju62WQHSZqdpnevsqextjT/o9M0XknT3D2V3w18FPhJ\nWlcjS4x7sqSlaT7PKD0PKLXputxwG5zSVNa/8ZWSnpP0sqRbc8NOVHY66m+SnpY0IZW3uK1L+nBq\n56vpf+WmVC5JP5K0Mk3vEUnVekBkzxOd3BNOEV6830vXE2RdEvYCGnn/MbPD03DXAL8k68RjONkj\nfaemui+RPeWyqTem36dxe6f6W4D/JvtF9LZkT0T8Yqo7mQ17b/sVMK2MuB8ATm5WthCYlPs8KMUx\nEDgWWNxs+J+QPUe+1PSvItebWW5ZtdZz1Ed4vwOh4WS9PZ2RGz+AD6f3K4CD0/sBpC4Py2j3q8D+\nuc9jgTUtDHsoWS92L5Ad9f0I2CJXfw/Zk09fJHvY3aGtzPc84G3gOLJey76Wptmn1LKh7Z7Nrkqf\nDyF7gN6P89tBs3mPJOu57bA03bPTtDfNtaOl3vC2IOuMpmm+OwB75Np0XW7Y4Wy43d5B9sTQAWm+\nf5/Kx6X1cBjZF9QhwG5lbOs3kB2JbkL2AMCDUvnhZA8F3JrsEeS7k+tJ0K8NXz4iqK2mo4LDyHZg\n6583n77RTCZ7QuiayDo5+SHwT2mQE4CL4v3emL6fG3c74AiyHeJrEbGSbIdUsveviDgqIqZ3sA1b\nkv2DNml6369EXVN9qd6pWtNiz1ERMS8i5kTEO2kZ/Tfw9y1M522ynp76R8TLEbHRqa4WlGrjli1c\nJ3ic7KmTO5DtiD8CXJirP4fs9NkQYAZwu6QPtTLveRExKyLeTtPpS5b4muSXTTk9m90REfdF1qHL\nN4G/k9S8jw3Innp5R0TMTvO+gCzZ/J9WYs17Dxgt6QMRsSJKPG20OUk7kD1x9Utp/bwdEfem6qnA\nFSme9yJieUQ8Xsa2/japo5iIeDMiHsiV9yPrMU4RsTj86JgWORHU1rXAP5J9Q7+mWd0gsm9Ez+bK\n8r1U7cjGvTE12SmNu0Lv95D032TflqqteQ9UTe/L6b2qXC32HCVppKRfKV2sBv4vG/cI1+TTZDuN\nZ9Ppgr8rc/6l2rg20lfLvIh4PiIeSzurZ8i+SX86V/9gSuzrIrt28IcUU0vWtz2yHssa2bDHt/yy\naVfPZhGxlqyvhFI9Xu1IbptK01xGGT3iRdYHwSSyo9YVku6QVE4XncOAlyI7xViqrlQXjW1t62eT\nfeP/czrN9fkU491kR6eXkPX8NUNS823VEieCGoqIZ8kO9Y8g67wj70Xe/zbTJN9L1Qo27o2pyTKy\nPhMGxfs9JPWPiD2qGX/yKNlF1CZ7Ay9ExOpUt7OyPnzz9S19O2zpwVatPfDqMrJv4btERH+y0yIl\n7+iJiL9ExESyncStZA8TLEepNpZ7oTRo/f8oaL13tPXrWNn1oaFs2ONbftmU07NZfnpbkp1WLNWD\n3HPktr109DOMMnvEi4jfRsRhZEdGj5P1CgfZ6ab8TRHb594vA7aRtHWJSS4DSh05tbqtp8T8zxGx\nI9kF/kuVbieOiIsj4iPAKLJTYWeV07YiciKovanAx9K3qPUi69xiJnC+pH7pot6ZQNOFtpnA6cp6\nYxpAdudO07gryDrn+KGk/pI2kfQhSS2dMmmVpE2V9V4loI+y3quato1rgKmSRqV/4G+RnYsmIp4k\nO4d9bhrnWLK7bn7Rwqya9xBXjn5k56PXpm+dXy41UGrDZyVtlU51/I1cT0/pguWhLczjGuBMSUOU\n9Rf7b01tLDGfj0raKV2MHEbWqcgvU93Wkg5Py6J3uoB6CK0/1vwjkj6VLqaeQbbTm9PCsOX0bHaE\npIOUXVD/d2BORJQ64poJHClpvKQ+qc3rKONpqpK2Sxd2t0jjrOX9Zb0AOETSByVtRa6TlbTd/pps\nZz0gteGQVH058LkUzyZpXezW1rYu6XhJQ9M0XiZLnO9J2k/S/qltr5HdCNHpPX91G/W+SNETX+Qu\ndjYr782GF4sHkO34V5F98/k2sElu2B8Bq8mOKk5jw4tuW5F9W24kO6f9EDA51Z3MhheLfw18o5V4\n70nTzr8OzdWfSbYT/xtwJbBZrm54Gv8NsovjG7U7N+wuZDuKV4BbW1pWbHix+BCyb5xrgfuB7zZr\nWwAfJrt4+huyncHfgL/w/oXDYalsYAtxCfgPstMoL6X3ytWv5f2L0GeSfWt+Pa2zi4F+qW5wmu+a\n1MY5wGGtLI/zgFlkF0/XpHU4ptRyyJXtAdyb1vljwLG5uqvIunacnWK+DxjRyvyPTdN4NU1zj2bb\nREsXi3fIxfBKGnZUrv6SVL4E+Gc23G63Aa5O29PLwM3N4nk4LYslwOFlbOv/kdbHWrJTS6ek8vFp\nWmvJjr6vB7as976hq778GGrr8SSdSLaTK9XtZN1IOo/sjqcT6x2LFZs7prEeLyKua3sos+LyNQIz\ns4LrEqeGBg0aFMOHD693GGYbeyI9KmjXqjwpwqyq5s2b92JEdOiRLnld4tTQ8OHDmTt3br3DMNvY\noYdmf++5p55RmJUk6dm2h2qbTw2ZmRWcE4GZWcE5EZiZFVyXuEZgZt3H22+/TWNjI2++Weqp5VYL\nffv2ZejQofTp06cm03ciMLN2aWxspF+/fgwfPpzSD2i1aooIVq9eTWNjIyNGjKjJPHxqyMza5c03\n32TgwIFOAp1EEgMHDqzpEZgTgZm1m5NA56r18nYiMDMrOCcCM7OC88Vis5zh0+7Y4PONS1cDMLlZ\neVsaph9ZtZjMas1HBGbW7bzyyitceumlNZ/PLbfcwle/+tWaz6fenAjMrNtpKRFEBO+9V72OyObP\nn8+YMWOqNr2uyqeGzKzbmTZtGk8//TT77LMPu+yyCw8//DD7778/8+bN48477yQiOOqoo1i0aBEA\nF1xwAWvXruW8887juuuu4+KLL+att95i//3359JLL6VXr14bTP/JJ5/ktNNOY86cOQwcOJBXX32V\nM844ox5N7RROBGbWcWecAQsWVHea++wDF13U6iDTp09n0aJFLFiwgIaGBnbeeWeuvvpqDjjgAAAa\nGhpKjrd48WJuuukm/vCHP9CnTx9OPfVUrr/+ek466aT1w6xbt44TTjiBa6+9lokTJ/LHP/6RUaNG\n8aUvfYm+fftWrZldiROBmXV7O+200/ok0Jq77rqLefPmsd9++wHwxhtvsO22224wzOzZs9l7773Z\ncccd6d+/P9tvvz19+/bl3XffrUnsXYETgZl1XBvf3DvLFltsscHn3r17b3CtoOlXuRHBlClT+P73\nv9/itBYuXMiee+7Jww8/zF577cXKlSvp16/fRvPoSXyx2My6nX79+rFmzZoW67fbbjtWrlzJ6tWr\nWbduHb/61a8AGD9+PLNmzWLlypUAvPTSSzz77IZ9u/Tr14/HH3+chQsXstdee3Huuedy2mmn1a4x\nXUCbRwSSrgCOAlZGxOhc+VeB04B3gTsi4uxU/nVgaio/PSJ+W4vAzay4Bg4cyIEHHsjo0aPZfffd\nN6rv06cP3/72txk3bhxDhgxht912A2DUqFF873vf4xOf+ATvvfceffr04ZJLLmGnnXZaP+6JJ57I\nsccey80338yAAQOYPHlyj7+FtJxTQ1cBPwGuaSqQ9FFgIrB3RKyTtG0qHwVMBvYAdgR+J2lkRPTc\nk2tmVhc///nPW60//fTTOf300zcqnzRpEpMmTWpxvG222YZ7772XffbZh9/97ncMGjSo4li7ujZP\nDUXEfcBLzYq/DEyPiHVpmJWpfCJwY0Ssi4hngCXAuCrGa2ZWc+vWrePVV18tRBKAjl8jGAkcLOlB\nSfdK2i+VDwGW5YZrTGUbkXSKpLmS5q5ataqDYZiZVd9mm23GM888U+8wOk1H7xrqDWwDHADsB8yU\ntHN7JhARM4AZAGPHjo0OxmEGbPyMIDMrX0ePCBqBmyPzZ+A9YBCwHBiWG25oKjMzsy6qo4ngVuCj\nAJJGApsCLwK3AZMlbSZpBLAL8OdqBGpmZrVRzu2jNwCHAoMkNQLnAlcAV0haBLwFTImIAB6VNBN4\nDHgHOM13DJmZdW1tJoKI+EwLVSe2MPz5wPmVBGVmZp3Hvyw2Mys4P2vIzKzODj74YNasWUPv3r2Z\nO3dup8/ficDMrM7uv//+us7fp4bMzArOicDMeoyGhgZGjx7d9oCd7LzzzuOCCy5osb7efSP71JCZ\nVaTav+pumH5kVadXiYggIthkk9p+Z65338g+IjCzbunCCy9k9OjRjB49motyHeS88847fPazn2X3\n3XfnuOOO4/XXX+e1117jyCOPZO+992b06NHcdNNNAFx33XWMGzeOffbZhy9+8Yu8++67NDQ0sOuu\nu3LSSScxevRopk6dyiWXXLJ++vlv96XGb3L++eczcuRIDjroIJ544omSbXjyySc57LDDuOiii/jO\nd76zQTs6kxOBmXU78+bN48orr+TBBx9kzpw5/PSnP+Whhx4C4IknnuDUU09l8eLF9O/fn0svvZTf\n/OY37LjjjixcuJBFixYxYcKEDfovXrBgAb169eL6668H4KmnnuLUU0/l0Ucf5fTTT2fmzJnr5z1z\n5kwmTZrU6vjz5s3jxhtvZMGCBdx555385S9/2agNTX0jX3jhhQwePJg5c+bw3e9+d31vap3JicDM\nup0HHniAY489li222IItt9yST33qU+vvvBk2bBgHHnggkHUy88ADD7Dnnnsye/ZszjnnHO6//362\n2mqrDfov3meffbjrrrtYunQpsGEfyPvuuy8rV67kueeeY+HChQwYMIBhw4a1Ov7999/Psccey+ab\nb07//v05+uijN2pDV+ob2dcIzKxHkbTR55EjRzJ//nzuvPNOvvWtbzF+/HgGDBhQsv/ihoaGjfon\nPv7445k1axbPP//8+k5tyun/uDVdqW9kHxGYWbdz8MEHc+utt64//3/LLbdw8MEHA/DXv/6VP/3p\nT0DWi9lBBx3Ec889x+abb86JJ57IWWedxfz588vqv7jJpEmTuPHGG5k1axbHH3880Hr/x4cccgi3\n3norb7zxBmvWrOH222/faJpdqW9kHxGYWbczZswYTj75ZMaNyzpA/MIXvsC+++67/kLvJZdcwuc/\n/3lGjRrFl7/8Ze6//37OOussNtlkE/r06cNll13WYv/F22+//Ubz22OPPVizZg1Dhgxhhx12AFrv\n/3jMmDFMmjSJvffem2233Zb99ttvo2l2pb6RlT00tL7Gjh0b9fhZtfUcteqY5safTwNg8j9Ob9d4\nXekWyGpbvHhxyQ7jrWPK7Ru51HKXNC8ixlYag08NmZnVSVfpG9mJwMysTrpK38htJgJJV0hamTqh\naV73b5JC0qD0WZIulrRE0sOS6vdTOTMzK0s5RwRXAROaF0oaBnwC+Guu+JNk3VPuApwCXFZ5iGZm\nVkttJoKIuA94qUTVj4CzgfzV5onANalT+znA1pJ2qEqkZmZWEx26RiBpIrA8IhY2qxoCLMt9bkxl\npaZxiqS5kuauWrWqI2GYWZ10hbsNi6TWy7vdiUDS5sA3gG9XMuOImBERYyNi7ODBgyuZlJl1or59\n+7J69Wong04SEaxevZq+ffvWbB4d+UHZh4ARwML0U+6hwHxJ44DlwLDcsENTmZn1EEOHDqWxsREf\nyXeevn37MnTo0JpNv92JICIeAbZt+iypARgbES9Kug34iqQbgf2BVyNiRbWCNbP669OnDyNGjKh3\nGFZF5dw+egPwJ2BXSY2SprYy+J3AUmAJ8FPg1KpEaWZmNdPmEUFEfKaN+uG59wHU56lJZmbWIf5l\nsZlZwTkRmJkVnBOBmVnBORGYmRWcO6Yx68Kq1c9CT+4fwSrnIwIzs4JzIjAzKzgnAjOzgnMiMDMr\nOCcCM7OCcyIwMys4JwIzs4JzIjAzKzgnAjOzgnMiMDMrOCcCM7OCK6eHsiskrZS0KFf2/yQ9Lulh\nSbdI2jpX93VJSyQ9IenwWgVuZmbVUc4RwVXAhGZls4HREbEX8CTwdQBJo4DJwB5pnEsl9apatGZm\nVnVtJoKIuA94qVnZ/0bEO+njHGBoej8RuDEi1kXEM2R9F4+rYrxmZlZl1bhG8Hng1+n9EGBZrq4x\nlW1E0imS5kqau2rVqiqEYWZmHVFRIpD0TeAd4Pr2jhsRMyJibESMHTx4cCVhmJlZBTrcMY2kk4Gj\ngPEREal4OTAsN9jQVGZmZl1Uh44IJE0AzgaOjojXc1W3AZMlbSZpBLAL8OfKwzQzs1pp84hA0g3A\nocAgSY3AuWR3CW0GzJYEMCcivhQRj0qaCTxGdsrotIh4t1bBm5lZ5dpMBBHxmRLFl7cy/PnA+ZUE\nZWZmnce/LDYzKzgnAjOzgnMiMDMrOCcCM7OCcyIwMys4JwIzs4JzIjAzKzgnAjOzgnMiMDMrOCcC\nM7OCcyIwMys4JwIzs4JzIjAzKzgnAjOzgnMiMDMruDYTgaQrJK2UtChXto2k2ZKeSn8HpHJJuljS\nEkkPSxpTy+DNzKxy5RwRXAVMaFY2DbgrInYB7kqfAT5J1j3lLsApwGXVCdPMzGqlzUQQEfcBLzUr\nnghcnd5fDRyTK78mMnOArSXtUK1gzcys+trsqrIF20XEivT+eWC79H4IsCw3XGMqW0Ezkk4hO2rg\ngx/8YAfDsHoaPu2OiqfRMP3IKkRiZpWo+GJxRAQQHRhvRkSMjYixgwcPrjQMMzProI4mgheaTvmk\nvytT+XJgWG64oanMzMy6qI4mgtuAKen9FOCXufKT0t1DBwCv5k4hmZlZF9TmNQJJNwCHAoMkNQLn\nAtOBmZKmAs8CJ6TB7wSOAJYArwOfq0HMZmZWRW0mgoj4TAtV40sMG8BplQZlZmadx78sNjMrOCcC\nM7OCcyIwMys4JwIzs4JzIjAzKzgnAjOzgnMiMDMrOCcCM7OCcyIwMys4JwIzs4JzIjAzKzgnAjOz\ngnMiMDMrOCcCM7OCcyIwMys4JwIzs4KrKBFI+ldJj0paJOkGSX0ljZD0oKQlkm6StGm1gjUzs+rr\ncCKQNAQ4HRgbEaOBXsBk4AfAjyLiw8DLwNRqBGpmZrVR6amh3sAHJPUGNgdWAB8DZqX6q4FjKpyH\nmZnVUIcTQUQsBy4A/kqWAF4F5gGvRMQ7abBGYEip8SWdImmupLmrVq3qaBhmZlahNjuvb4mkAcBE\nYATwCvA/wIRyx4+IGcAMgLFjx0ZH47Dubfi0O+odglnhVXJq6OPAMxGxKiLeBm4GDgS2TqeKAIYC\nyyuM0czMaqiSRPBX4ABJm0sSMB54DPg9cFwaZgrwy8pCNDOzWqrkGsGDZBeF5wOPpGnNAM4BzpS0\nBBgIXF6FOM3MrEY6fI0AICLOBc5tVrwUGFfJdM3MrPP4l8VmZgXnRGBmVnAVnRoys9J8W6x1Jz4i\nMDMrOCcCM7OCcyIwMys4JwIzs4JzIjAzKzgnAjOzgnMiMDMrOCcCM7OCcyIwMys4JwIzs4JzIjAz\nKzgnAjOzgqsoEUjaWtIsSY9LWizp7yRtI2m2pKfS3wHVCtbMzKqv0iOCHwO/iYjdgL2BxcA04K6I\n2AW4K302M7MuqsOJQNJWwCGkrigj4q2IeAWYCFydBrsaOKbSIM3MrHYqOSIYAawCrpT0kKSfSdoC\n2C4iVqRhnge2KzWypFMkzZWxJOYwAAAIoklEQVQ0d9WqVRWEYWZmlagkEfQGxgCXRcS+wGs0Ow0U\nEQFEqZEjYkZEjI2IsYMHD64gDDMzq0QliaARaIyIB9PnWWSJ4QVJOwCkvysrC9HMzGqpw4kgIp4H\nlknaNRWNBx4DbgOmpLIpwC8ritDMzGqq0j6LvwpcL2lTYCnwObLkMlPSVOBZ4IQK52FmZjVUUSKI\niAXA2BJV4yuZrpmZdR7/stjMrOCcCMzMCq7SawRmViDDp91R8TQaph9ZhUismnxEYGZWcE4EZmYF\n51ND3Ug1DsvBh+ZmtiEfEZiZFZwTgZlZwfnUkFkBVOu0ovVMPiIwMys4JwIzs4JzIjAzKzgnAjOz\ngnMiMDMrON81ZGadyj+M7Hp8RGBmVnAVJwJJvSQ9JOlX6fMISQ9KWiLpptR7mZmZdVHVOCL4F2Bx\n7vMPgB9FxIeBl4GpVZiHmZnVSEWJQNJQ4EjgZ+mzgI8Bs9IgVwPHVDIPMzOrrUqPCC4CzgbeS58H\nAq9ExDvpcyMwpNSIkk6RNFfS3FWrVlUYhpmZdVSHE4Gko4CVETGvI+NHxIyIGBsRYwcPHtzRMMzM\nrEKV3D56IHC0pCOAvkB/4MfA1pJ6p6OCocDyysM0M7Na6fARQUR8PSKGRsRwYDJwd0R8Fvg9cFwa\nbArwy4qjNDOzmqnF7wjOAc6UtITsmsHlNZiHmZlVSVV+WRwR9wD3pPdLgXHVmK6ZmdWef1lsZlZw\nftZQAbm3KjPL8xGBmVnBORGYmRWcE4GZWcE5EZiZFZwTgZlZwTkRmJkVnBOBmVnBORGYmRWcf1DW\nCfwDLjPrynxEYGZWcE4EZmYF50RgZlZwTgRmZgVXSZ/FwyT9XtJjkh6V9C+pfBtJsyU9lf4OqF64\nZmZWbZUcEbwD/FtEjAIOAE6TNAqYBtwVEbsAd6XPZmbWRVXSZ/GKiJif3q8BFgNDgInA1Wmwq4Fj\nKg3SzMxqpyrXCCQNB/YFHgS2i4gVqep5YLsWxjlF0lxJc1etWlWNMMzMrAMqTgSStgR+AZwREX/L\n10VEAFFqvIiYERFjI2Ls4MGDKw3DzMw6qKJEIKkPWRK4PiJuTsUvSNoh1e8ArKwsRDMzq6VK7hoS\ncDmwOCIuzFXdBkxJ76cAv+x4eGZmVmuVPGvoQOCfgEckLUhl3wCmAzMlTQWeBU6oLEQzM6ulDieC\niHgAUAvV4zs6XTMz61z+ZbGZWcE5EZiZFZwTgZlZwTkRmJkVnBOBmVnBORGYmRWcE4GZWcE5EZiZ\nFZwTgZlZwVXyiAkzs25v+LQ7qjKdhulHVmU69eAjAjOzguuRRwTO8GZm5fMRgZlZwfXIIwIzs87W\nnc9EOBG0olor1sysK/OpITOzgqvZEYGkCcCPgV7AzyJielvj+Bu4mZXL+4vqqckRgaRewCXAJ4FR\nwGckjarFvMzMrDK1OjU0DlgSEUsj4i3gRmBiSwM/svxVZ3czszpRRFR/otJxwISI+EL6/E/A/hHx\nldwwpwCnpI+jgUVVD6TrGAS8WO8gasjt6756ctug57dv14joV+lE6nbXUETMAGYASJobEWPrFUut\nuX3dW09uX09uGxSjfdWYTq1ODS0HhuU+D01lZmbWxdQqEfwF2EXSCEmbApOB22o0LzMzq0BNTg1F\nxDuSvgL8luz20Ssi4tFWRplRizi6ELeve+vJ7evJbQO3ryw1uVhsZmbdh39ZbGZWcE4EZmYFV/NE\nIGmCpCckLZE0rZXhPi0pJI3NlX09jfeEpMNrHWtHdLR9koZLekPSgvT6r86LujxttU3SyZJW5drw\nhVzdFElPpdeUzo28PBW2791ceZe8EaKcbVPSCZIek/SopJ/nyrv9+kvDtNS+br/+JP0o14YnJb2S\nq2vf+ouImr3ILhQ/DewMbAosBEaVGK4fcB8wBxibykal4TcDRqTp9KplvJ3cvuHAonq3oZK2AScD\nPykx7jbA0vR3QHo/oN5tqlb7Ut3aerehCu3bBXioad0A2/aw9VeyfT1l/TUb/qtkN+V0aP3V+oig\n3EdN/DvwA+DNXNlE4MaIWBcRzwBL0vS6kkra19W16zEhzRwOzI6IlyLiZWA2MKFGcXZUJe3rDspp\n3z8Dl6R1RESsTOU9Zf211L7uoL3b52eAG9L7dq+/WieCIcCy3OfGVLaepDHAsIho/rChNsftAipp\nH8AISQ9JulfSwTWMsyPKXf6flvSwpFmSmn5E2CPWXVKqfQB9Jc2VNEfSMTWNtGPKad9IYKSkP6R2\nTGjHuPVWSfugZ6w/ACTtRHbW5O72jtukrh3TSNoEuJDsELzHaaN9K4APRsRqSR8BbpW0R0T8rTNj\nrNDtwA0RsU7SF4GrgY/VOaZqaq19O0XEckk7A3dLeiQinq5bpB3Tm+z0yaFkv/6/T9KedY2oukq2\nLyJeoWesvyaTgVkR8W5HJ1DrI4K2HjXRj+yBc/dIagAOAG5LF1S7w2MqOty+dMprNUBEzCM7Hziy\nU6IuT5vLPyJWR8S69PFnwEfKHbcLqKR9RMTy9HcpcA+wby2D7YBy1kEjcFtEvJ1Ovz5JtuPsEeuP\nltvXU9Zfk8m8f1qoveNmanzBozfZhYoRvH/BY49Whr+H9y+m7sGGF4uX0vUuFlfSvsFN7SG7ILQc\n2KbebWpP24Adcu+PBebE+xerniG7UDUgve8ybatC+wYAm6X3g4CnaOVCXhdu3wTg6lw7lgEDe9D6\na6l9PWL9peF2AxpIPw5OZe1ef53RoCPIMvHTwDdT2XeBo0sMu35HmT5/M433BPDJeq+carYP+DTw\nKLAAmA/8Q73b0t62Ad9PbVgI/B7YLTfu58ku8C8BPlfvtlSzfcD/AR5J5Y8AU+vdlg62T2SnLh9L\n7Zjcw9Zfyfb1lPWXPp8HTC8xbrvWnx8xYWZWcP5lsZlZwTkRmJkVnBOBmVnBORGYmRWcE4GZWcE5\nEZiZFZwTgZlZwf1/YWV+ov/9piMAAAAASUVORK5CYII=\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[]}}]},{"cell_type":"markdown","metadata":{"id":"7D8Pdy6uVJYe","colab_type":"text"},"source":["## Bootstrap estimation of confidence interval"]},{"cell_type":"code","metadata":{"id":"BGIfn9Uf4e6y","colab_type":"code","colab":{"base_uri":"https://localhost:8080/","height":298},"outputId":"8438dc44-99c6-4819-8729-0d9dbb3ab9ff","executionInfo":{"status":"ok","timestamp":1570655010420,"user_tz":-120,"elapsed":1021,"user":{"displayName":"Jan Chorowski","photoUrl":"","userId":"18100051538790809279"}}},"source":["# Here we make a bootstrap analysis of one experiment\n","n_bootstraps = 200\n","exp_id = 1\n","exp0 = samples[exp_id]\n","\n","# sample answers with replacement\n","bootstrap_idx = np.random.randint(low=0, high=n, size=(n_bootstraps, n))\n","exp0_bootstraps = exp0[bootstrap_idx]\n","\n","# compute the mean in each bootstrap sample\n","exp0_bootstrap_means = exp0_bootstraps.mean(1)\n","\n","# Estimate the confidence interval by taking the 2.5 and 97.5 percentile\n","sorted_bootstrap_means = np.sort(exp0_bootstrap_means)\n","bootstrap_conf_low, bootstrap_conf_high = sorted_bootstrap_means[\n","    [int(0.025 * n_bootstraps), int(0.975 * n_bootstraps)]]\n","\n","hist(exp0_bootstrap_means, bins=20, label='bootstrap estims of $\\phi$')\n","axvline(phi, 0, 1, label='$\\\\phi$', color='red')\n","axvline(phi_bar[exp_id], 0, 1, label='$\\\\bar{\\\\phi}$', color='green')\n","axvspan(confidence_intervals[exp_id, 0], confidence_intervals[exp_id, 1], # ymin=0.5, ymax=1.0, \n","        alpha=0.2, label='theoretical 95% conf int', color='green')\n","axvspan(bootstrap_conf_low, bootstrap_conf_high, # ymin=0.0, ymax=0.5, \n","        alpha=0.2, label='bootsrap 95% conf int', color='blue')\n","legend()\n","_ = xlim(phi-0.15, phi+0.15)\n","title('Theoretical and bootstrap confidence intervals')"],"execution_count":100,"outputs":[{"output_type":"execute_result","data":{"text/plain":["Text(0.5, 1.0, 'Theoretical and bootstrap confidence intervals')"]},"metadata":{"tags":[]},"execution_count":100},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8FVX6x/HPExIICEpVREFs9N4R\nEHBRUPSHLCgonUXEFVxcy7JW7F2w0AVEdmFxLaCw66poFBRwkQUEFUGJFFF6qEGB8/tjJuEm3NR7\nQ8p8369XXrn3THvOzOTJ3DNzzzHnHCIiUvTF5HcAIiJyaijhi4gEhBK+iEhAKOGLiASEEr6ISEAo\n4YuIBEQgE76ZjTazv+V3HOGYWTUzO2BmxSJcT4KZDYlWXJlsp7qZOTOLzWB6opl1yus4gsA8081s\nj5l9YWbtzGxdJvO/amaPnsoYM2Nm95jZK/kdR3acqr+fUy3sH2lhZ2YHQt6WAo4Ax/z3N5/6iDJm\nZonAEOfchwDOuU1A6XwNqhAwswTgb865LBOImb0KbHHO3ZfXceWxtsDlwLnOuYN+Wc18jCdHnHOP\nZ3deMxsNXOSc65t3EQVPkbzCd86VTvkBNgHXhJT9/VTFkdFVrxQsheg4nQckhiR7yUAhOqanVJFM\n+NlU3MxeM7P9ZrbWzJqlTDCzKmb2ppntMLONZnZbyLQSZjbWzH7yf8aaWQl/Wgcz22JmfzGzn4Hp\nfvnVZrbSzPaa2edm1sAvnwlUA971m3HuTt9EYmbl/Y/xP/kf5ef65eXMbL4f4x7/9bnZqbiZtTCz\nJX4828zsZTMrHjLdmdkwM1vvzzPOzMyfVszMnjWznWb2A9A1G5tsbmZf+3FON7P4kG3dZGYbzGy3\nmb1jZlVCpl1iZv81syT/9yV++WNAO+Blf7+97Dd3jDGz7Wa2z8y+MrN6ZjYU6APc7c/7rr+ORP84\nrQYOmlmsmY0ys+/9c+JrM+seEstAM/vM31aSmX1rZr/LZB9XNbO3/OOzy8xe9stjzOw+M/vRj/U1\nMzvDn5Zy7AeY2SZ/H9/rT/sD8ArQ2q/HQynnW8g2G5vZCj/+OUB8upjCnoch++NOM1vt129OuuPU\nzV92n7+PuvjlZ5jZVP882mpmj1oGzZEW0pSaRV27APcAvfy6rspqWyHHZ4yZ7QIe8etZL2T7lczs\nsJmdaTn4+zGzi8zsE3+/7PT3beHknCvSP0Ai0Cld2WggGbgKKAY8ASz1p8UAXwIPAMWBC4AfgM7+\n9IeBpcCZQCXgc+ARf1oH4CjwFFACKAk0BrYDLf1tDfBjKhEuPqA64IBY//0CYA5QDogD2vvlFYAe\neE1WZYB/AnND1pOA11QUbp80BVrhNelVB74BRoZMd8B8oCzeP6QdQBd/2jDgW6AqUB74ODTeDPb/\nmpD5PwMe9addBuwEmvj76yXgU39aeWAP0M+P8wb/fYVw9QM6+8etLGBAbeBsf9qrKdtMF9dKP66S\nftl1QBX/HOgFHAxZx0D/2N7uH4deQBJQPkydiwGrgDHAaXiJt60/bTCwAe+8Kg28BcxMd+yn4J07\nDfGaI2uHxLA4ZDsd8JqqwDtXfwyJryfwW8i+zs55+IVf//J458Qwf1oLv66X+/vmHKCWP+1tYJJf\nzzP9ddycwbkwGq8ZLjt1TZ03ZPkMtxVyfEbgnS8lgWnAYyHL3wq8l9O/H2A2cK9f99RjWRh/8j2A\nPK9gxgn/w5D3dYDD/uuWwKZ08/8VmO6//h64KmRaZ7yP2Sl/gL8C8SHTJ+D/QwgpW8eJxJ0mvpA/\nhFjgbOA4UC4b9WwE7Al5n3rCZmPZkcDbIe9d6EkNvA6M8l9/hJ8I/PdXkHXCD53/KuB7//VU4OmQ\naaXxklR1vET/Rbp1LQEGhqsf3j+P7/D+kcWkW+5Vwif8wVnsl5VAN//1QOAnwEKmfwH0C7Nca7x/\nkiftE2Ah8MeQ9zX9Oqf883V4bfSh2+gdEkNGCf/SMPF9zomEn53zsG/ItKeBif7rScCYMHU5Cy9J\nlwwpuwH4OIP9OZqTE35GdU2dNzvb8vdN+r/bTinnmv/+M6B/Tv9+gNeAyaGxFtafILdz/Rzy+hAQ\nb14zynlAFTPbGzK9GLDIf10F70oqxY9+WYodzrnkkPfnAQPMbERIWfF0y2SkKrDbObcn/QQzK4V3\nBdkF7+ofoIyZFXPOHUs/f7plawDPA83wrnBi8a6OQ6XfPyk3kqsAm0Omhe6LjKSfP6XuVYAVKROc\ncwf8j+PncPJ+Tln2nHAbcM595DebjAPOM7O3gDudc/uyGRdm1h/4M14yAq/OFUNm2er8DBCmLqGq\nAj86546GmRbu/InFS2gpMtr3mamSQXwpsnMept9uyrSqwL/CbPM8vE8T28xr8QPvKnhzmHkzkt26\nZmdb6bf7MVDKzFoCv+Al9bchx38/dwOPAF+Y2R7gOefctOxVr2AJcht+RjYDG51zZUN+yjjnrvKn\n/4R38qWo5pelSN/96Ga8j5Wh6yvlnJudwfzply1vZmXDTLsD7+qwpXPudLwrPPCaM7IyAa9Z5mJ/\n2XuyuRzANrwEkKJaNpZJP3/K/kqzL83sNLyP2lvTTwtZdqv/+qT95px70TnXFO8TWw3grozmTV9u\nZufhNS8Mx2s2KovXFBW6X86xkGzDycc+xWagmoW/cRju/DmKl5AisS2D+EJjyuw8zMxm4MIMyo8A\nFUPWebpzrm6ua3FCuL+jrLaVZhk/cb+O90ngBmC+c26/Pznbfz/OuZ+dczc556rgPeU33swuirB+\n+UIJ/2RfAPv9G3olzbtJWc/MmvvTZwP3+TeAKuK19Wf2TP8UYJiZtTTPaWbW1czK+NN/wWvPPYlz\nbhvwb7wTrJyZxZlZyolZBjgM7DWz8sCDOahjGWAfcMDMagG35GDZ14HbzOxcMysHjMrGMrf685fH\nawtNuek1GxhkZo3Mu/H9OLDMOZeId0VZw8xuNO+Gai+8RD7fXzbNfjOz5v4+jsNre0/Gaw47ad4M\nnIaXMHb46xsE1Es3z5l+3ePM7Dq8+wThrny/wEvAT/rHO97M2oTU+XYzO9/MSvt1npPBp4GcWIL3\njyMlvt/jtb2nyOo8zMxUvOP0O/NuOp9jZrX88/N94DkzO92fdqGZtY+wLuAds+pmFgOpfwu52dYs\nvPstffzXKbL992Nm14Xc0N2Dd54cz2j+gkwJPx3/quBqvI9/G/FuKr4CnOHP8iiwHFgNfIXXJJHh\nl1ucc8uBm4CX8U6WDXjtjSmewPsHstfM7gyzin54bbzf4t10G+mXj8W7MbUT7ybyezmo5p3AjcB+\nvESQk6cOpgD/wbspuQLvpmNWZuH9sf6Adw/kUQDnfffgfuBNvAR5IdDbn7YL7zjcAezC+1h9tXNu\np7/OF4Ce/hMWLwKn+7HtwWvK2AU84887Fajj7+O54QJ0zn0NPIeXOH8B6uO1+YZaBlyMt88fA3r6\ncaZf1zHgGuAivMeCt+AlHfBuJM4EPsU7v5LxbjRGxDn3K/B7vHNrt7+9t0KmZ3UeZrbuL4BBeE0g\nScAnnPiU0h+vaehrf71v4N17itQ//d+7zCyl2S/H23LOLcO7AKiCd/GUIid/P82BZeZ9v+cd4E/O\nuR9yVJsCwtI2+YlIOGY2EO8mXtv8jkUkt3SFLyISEEr4IiIBoSYdEZGA0BW+iEhAnNIvXlWsWNFV\nr179VG6yyDn460GOZfK9qsOHinE8069dSXoVft4EwK7K2flKQc7sPOqtu2Js9NcteSPGilEyzutG\nqFgxOO20fA4I+PLLL3c65ypFup5TmvCrV6/O8uXLT+Umi5z3NrxHpVIZH/fPPz6dchUifaQ7WC6/\newAAHzw9I+rrfuxrb9331on+uiVv7EnewyVVLwFgxw7o0iWfAwLMLDvfaM+SmnRERAJCCV9EJCCU\n8EVEAiLfe8v87bff2LJlC8nJyVnPLJz525kUS854uNsWNQ8To3/jYSX/ZmzZW4xjLrv9xIkULfme\n8Lds2UKZMmWoXr06aTv6k3CSkpOIjcn4sB3YX4zYWH23Ij3nHPuT9gK7+HFPvp/2Ivki368Fk5OT\nqVChgpK95Ckzo8wZZYmP0z9DCa58T/iAkr2cEjrPJOgKRMIXEZG8p4QvIhIQunslImmMnLMyx8uM\n7dUoDyKRaNMVfjpvv/02I0ZEPACRiEiBo4SfzooVK2jSpEl+hyEiEnVq0vF999133HrrrSxdupQK\nFSqQlJTEyJEjs15QRKSQKFgJf+RIWJnz9sNMNWoEY8dmOsuRI0e4/vrrmTlzJt26dePzzz+nTp06\nDBs2jPj4+OjGIyKST7Js0jGzqmb2sZl9bWZrzexPfvloM9tqZiv9n6vyPty88cEHH9CwYUOqVKnC\n6aefTuXKlYmPj+fYMXUsLyJFR3au8I8CdzjnVphZGeBLM/vAnzbGOfds1KLJ4ko8r6xatYr69euz\nevVqGjRowPbt2ylTpgynFYSRD0REoiTLhO+c2wZs81/vN7NvgHPyOrBTqUyZMqxevZrY2FgaNGjA\ngw8+yK233prfYYmIRFWOntIxs+pAY2CZXzTczFab2TQzK5fBMkPNbLmZLd+xY0dEweaVvn37sn79\neh5++GEmTJhA+fLl9WimiBQ52U74ZlYaeBMY6ZzbB0wALgQa4X0CeC7ccs65yc65Zs65ZpUqRTwk\nY54oX748n3zyCdWqVWPZsmU89thj6ndFRIqcbCV8M4vDS/Z/d869BeCc+8U5d8w5dxyYArTIuzDz\n3pEjR0hKSqJixYr5HYqISJ7IzlM6BkwFvnHOPR9SfnbIbN2BNdEP79QpUaIEGzduzO8wRETyTHae\n0mkD9AO+MrOUh+TvAW4ws0aAAxKBm/MkQhERiYrsPKWzGAjXoP2v6IcjIiJ5RX3piIgEhBK+iEhA\nKOGLiASEEr6ISEAUrN4y81m7du3Yv38/sbGxLF++PL/DERGJKiX8EIsWLcrvEERE8oyadEREAkIJ\nPx2NaSsiRZUSfjoa01ZEiiq14fs0pq2IFHVK+GhMWym6Rs6J8hjRUqgVqIQ/8r2RrPw5uidoo8qN\nGNsl86ETNaatiASB2vDRmLYiEgwF6go/qyvxvKIxbUUkCApUws8vffv2pXv37rz11luUK1eO3r17\n69FMESlylPA5MaZto0aN+PDDDzXMoYgUSWrD92lMWxEp6pTwfRrTVkSKOiV8EZGAUMIXEQkIJXwR\nkYBQwhcRCQglfBGRgFDCFxEJCCV8EZGAUMIXEQkIJXwRkYAocH3pVB+1IKrrS3yya9bzJCZy9dVX\ns2bNmoi2tXfvXmbNmsUf//jHiOY5VZL27uXNN+YweMjNqWVXXtGRf7//8SnZ/uSJ45k+bQoNGjZi\n0pTpYedZ8O48Fn36CU8+8/wpiUmkKNMVfhTt3buX8ePH53oe5xzHjx/Pi9DCSkpKYtrUKWnKTlWy\nB5g2dTJvvv1uhskeYPWqlTRo2OiUxSRSlCnh+44ePUqfPn2oXbs2PXv25NChQwA8//zz1KtXj3r1\n6jF27In++sOVjxo1iu+//55GjRpx1113cfDgQbp27UrDhg2pV68ec+bMOWmexMREatasSf/+/alX\nrx6bN2/m2muvpWnTptStW5fJkycD3qeQWrVqcdPAm2jSoAl9evdJjTG91+fM5vLL2tGhbUv+PHI4\nx44d4+DBg/S+vjvt27SkbetmvP3WGzz80P0kbvyBDm1b8uD99wBw3jmVANj044+0at6I4bcMpUXT\nBtx80yA+SfiIqzpfRvMm9Vnx5X8Bwq43vfEvv0jb1s1o27oZE8e/DMAdt4/gx8SN9LquOxPGvXTS\nMhs2rKfHtVczceI4nnnq8dTlRCT3smzSMbOqwGvAWYADJjvnXjCz8sAcoDqQCFzvnNuTd6HmrXXr\n1jF16lTatGnD4MGDGT9+PB07dmT69OksW7YM5xwtW7akffv2HD9+PGz5k08+yZo1a1i50hum8c03\n36RKlSosWOA1UyUlJdGyZcs08yQmJrJ+/XpmzJhBq1atAJg2bRrly5fn8OHDNG/enB49eqTG+ML4\nF2jbti23DL2FKROn8Kc//ylNPb777lvmvvUG//rPR8TFxXHXHX/ijdf/QalSp1G58tn84/W3AdiX\nlETTps359puvSVi8LOw+2fjD90x79W+8WHsinTq25c1/zmHBewv597/mM+a5Z5g563U++vCDk9Yb\nauXKFcyeNZP/fPgJzjk6d2rPJW3a8tyYl/joww+Y++6/qVAhbQ+lR44cYcjAvoyfNJV+N17Pv97/\niDYtmzJw8BCNMSwSgexc4R8F7nDO1QFaAbeaWR1gFLDQOXcxsNB/X2hVrVqVNm3aAN6AKIsXL2bx\n4sV0796d0047jdKlS/P73/+eRYsWZVieXv369fnggw/4y1/+wqJFizjjjDPCbvu8885LTfYAL774\nIg0bNqRVq1Zs3ryZ9evXp8bY6hJvvl439GLJ50tOWteiTz9m1ar/cXnHtnRo25JPP0kgMXEjtevW\n5ZOPP+KhB+9jyeefcXoGsYSqdl516tStR0xMDLVq1aFd+46YGXXq1mPzph8BslzvsiVLuKrrNan7\n6uqr/4+lSz7PdLsJHy+kbr0GVD77bMqUKcNZZ1WmRHwJjTEsEqEsE75zbptzboX/ej/wDXAO0A2Y\n4c82A7g2r4I8Fcws0/e5UaNGDVasWEH9+vW57777ePjhh8POFzp2bkJCAh9++CFLlixh1apVNG7c\nmOTk5PAxhQnROUfvG/qSsHgZCYuXsWz5Kv7y1/u46KKL+ejTz6ldpy6PP/oQzzz1eJbxlyhRIvV1\nTExM6vsYi+Gon3xzs96srF3zFbXr1mXtmjXUqVufHTu2U7q0xhgWiVSO2vDNrDrQGFgGnOWc2+ZP\n+hmvyafQ2rRpE0uWeFfMs2bNom3btrRr1465c+dy6NAhDh48yNtvv027du0yLC9Tpgz79+9PXedP\nP/1EqVKl6Nu3L3fddRcrVqw4aZ70kpKSKFeuHKVKleLbb79l6dKlaWL8YukXAPxzzj+55JJLTlq+\n3aUdeWfe2+zYsR2APXt2s3nTJrZt+4mSJUtxfa8bGH7bSFavWknpMqU5kEks2RFuvaFatb6Efy+Y\nn7qvFix4l1atT447VOnSZVj/3TrWrvmKOvXq8dQTj/KHkCeJRCR3sv1YppmVBt4ERjrn9oVebTrn\nnJm5DJYbCgwFqFatWpbbyc5jlHmhZs2ajBs3jsGDB1OnTh1uueUWSpUqxcCBA2nRogUAQ4YMoXHj\nxgAZlrdp04Z69epx5ZVX0qlTJ+666y5iYmKIi4tjwoQJVKhQIc086QdL79KlCxMnTqR27drUrFkz\nTVNPzZo1eWXSKwwfNpxatWsx5OYhYepRm3vue5Drul/D8eOO2LhYnn52LPv2JTH6/nuJiTFi4+J4\n5vkXKF++Ai1ataZt62b8rtMVPPRIzq/Ov/l67UnrDdWwUWN639iXK353KQB9+w3M8qmb63r15t0+\nc1nw7juULVuW7j16ctPNt+Q4NhFJy5wLm6fTzmQWB8wH/uOce94vWwd0cM5tM7OzgQTnXM3M1tOs\nWTO3fPnyNGXffPMNtWvXzm38gZHyXYHPln9GbEzG/6cP7C9GbGzWx7Qw6NC2JW/OW3DSTd1I/PD9\netZtj0tTdvndAwD44OkZ4RaJyGNfe+u+t070150dI+eszHqmKBjbq+g8OrsneQ+XVPU+he7YAV26\n5HNAgJl96ZxrFul6smzSMe9SfirwTUqy970DDPBfDwDmRRqMSIojR46wb9++qCZ7kaDLTht+G6Af\ncJmZrfR/rgKeBC43s/VAJ/+95JHq1atH/E3gwqREiRKsWP1NfochUqRk2YbvnFtM2OdBAPhddMMR\nEZG8om/aiogEhBK+iEhAKOGLiASEEr6ISEAo4YuIBESBGwBl6eal7D2yN2rrK1uiLK2qtspwevoB\nSRISEnj22WeZP39+1GLISkJCAsWLF0/tKmHixImUKlWK/v3753hdF5xXkU1bd5xUPmnCOGa+Nh3n\nHP36D2LYH4cD8NQTjzLztelU9J93v/eBh7j8ii4sW7qEu/78J+KKxzF56gwuvPAikvbu5Q+D+vH6\nm/OIicm7a4V5c9/iyccf4cwzz2Le/PdSy7dt+4l7/nIn01+blenyY557mtvvuDvP4hMprApcwt97\nZC+VSlWK2vp2HDo5+aXZnj8gSV6PQHX06FFiY8Pv7oSEBEqXLp2a8IcNGxbVbX/z9Vpmvjad9xd+\nSvHixbm+Rzeu6HIlF1xwobe9P45g+IiRaZYZ//IL/OOfb7Fp0yZenfYKjzz2JM89+xQj/3xXniZ7\ngL/PfJUxL4w7qc+ds8+ukmWyBxj7/DNK+CJhBL5JJ/2AJAAHDhygZ8+e1KpViz59+pDS/cSXX35J\n+/btadq0KZ07d2bbNq/vuJUrV9KqVSsaNGhA9+7d2bPHGxagQ4cOjBw5kmbNmvHCCy+wY8cOevTo\nQfPmzWnevDmfffYZiYmJTJw4kTFjxtCoUSMWLVrE6NGjefbZZwHYsGEDnTp1omHDhjRp0oSNP2zk\nwIEDdO3clTYt29CiSQvmv5P5p5HvvltH06bNKFWqFLGxsVzSpi3z3838i9GxcXEcOnyYw4cPERcX\nx8aNP7B16xbatrs0w2VWrFjOlVd0pH2bllx+WTv2799PcnIyI/44lHaXNKdju1Ys+vQTAGb/fSYD\n+vbm+h7/R/Mm9Rn9wL0APPPU4yxbuoQ/jbgldVCWFJt+/JG2rZtluvzDo+/n8OHDdGjbkptvGpRp\nHUWCpsBd4Z9q6QctSUhI4H//+x9r166lSpUqtGnThs8++4yWLVsyYsQI5s2bR6VKlZgzZw733nsv\n06ZNo3///rz00ku0b9+eBx54gIceeih1FKxff/2VlP6DbrzxRm6//Xbatm3Lpk2b6Ny5M9988w3D\nhg2jdOnS3HnnnQAsXLgwNb4+ffowatQounfvTnJyMnsO7SE+Pp7Z/5zN6aefzs6dO7ms3WV0vaZr\nhl06165dh8ceGc3u3buIjy/Jhx/8h0aNmqROnzp5Iq/PnkWjxo15+LEnKVu2HCNvv5Nbhw0hPr4k\nEya9wgP338M99z2Y4X789ddfuWlQf6ZMf40mTZqxf98+SpYsyaQJ4zAzFn3+X9Z/t46ev7+GZctX\nA7Dmq9V8/OkSipcoQatmDblp6C3c9Zd7WPTpJzz06OM0btw002MXbvkHRj/C1CkTMxzURSTIAp/w\nw2nRogXnnnsuAI0aNSIxMZGyZcuyZs0aLr/8cgCOHTvG2WefTVJSEnv37qV9+/YADBgwgOuuuy51\nXb169Up9/eGHH/L111+nvt+3bx8HDhzIMI79+/ezdetWunfvDkB8fDylKIU75hh9/2g+W/wZMTEx\n/PTTT2z/ZTtnVQ7fQ3WNmrW47U9/pmf3ayhV6jTq1W9AsWLFABj0h5u48+6/YmY88dhDPHDvKF4c\nN4n6DRrynw+9q/HPP1vMWWdVxjnHHwb1Iy4ujocffYIzzzyxvQ3rv+OsypVp0sS7Ai9z+ukALFv6\nOUOGej1dXlyjJlWrVuP7Dd6ALu3ad0gdMKVmrVps3ryJc/z9nh2RLi8SNEr4YYQO/FGsWDGOHj2K\nc466deum9pmfIindkH7phQ7acfz4cZYuXRrxMH1zZs9h586dLF66mLi4OOrUqJM6SEpG+vYfSN/+\nAwF49OEHqFLlHIA0Sbtf/8Hc2LtHmuWcczz/7FNMmTaDUXffweiHHmPTph+ZMmk8997/UET1SDPA\nir+fM7Mt6TC/HTvO5t2H2H3wV35zxdi82xvX99djsG3PATbvPsRxR2p51fKlIopRpCgJfBt+VgOS\npKhZsyY7duxITfi//fYba9eu5YwzzqBcuXKpQxzOnDkz9Wo/vSuuuIKXXjoxYHdKM1JGMZQpU4Zz\nzz2XuXPnAl4PkocOHSIpKYlKlSoRFxfHJwmfsOnHTVnGnzIgypbNm5n/7jv06Ol98vj5522p8yyY\n/w61atdJs9yc2X+n0+WdKVeuPIcPHSImJoaYmBgOHTqcZr6LLq7BLz//zIoVXvPV/v37OXr0KK1a\nt+HNf/4D8AYm37J5MxddXCPLeCMRGxvLb7/9lqfbECmMCtwVftkSZbN8sian68tM+gFJunYNPwBL\n8eLFeeONN7jttttISkri6NGjjBw5krp16zJjxgyGDRvGoUOHuOCCC5g+fXrYdbz44ovceuutNGjQ\ngKNHj3LppZcyceJErrnmGnr27Mm8efPS/EMA7x/IzTffzAMPPEBcXBzT/jaNXjf04vrfX0+LJi1o\n0rQJNWpmnUAH9b+R3bt3Excbx9PPjuGMst5+eeiB+1izZjWGUbVaNZ4be2L7hw4dYvasv/HG2+8C\ncMutt9H7+u7ExcUx6ZVXT9o/U6a/xl/vvoPkw8nEl4znzbkLGDxkKHf9+TbaXdKc2GLFeGn85DRX\n9nmhR+/+9LrqUmrVbcDM12bm6bZECpNsDYASLRoAJXJJyUmBGQAlvZRmmpxI36SjAVDyhgZAyVun\nbAAUEREpGpTwRUQCQglfRCQglPBFRAKiwD2lIxJN6W/07jn4KyPnrE1TNm+79+W3lBucRekGpEgo\nXeGLiAREgbvCX7oU9kavd2TKloVWGfeOTGJiIldffTVr1qyJeFtjx45l6NChlCp16r7dOe6lcbw6\n7VWccwwaPIgBg24Dima3x1PHjeEPt96eZ/GJFHUFLuHv3QuVotc7Mjui9x2uLI0dO5a+fftmO+Fn\n1mVydqxdu5ZXp73KJ599QvHixbn26mu5tH1XLq5xAVD0uj2eOmGsEr5IBNSkg5d4+/TpQ+3atenZ\nsyeHDnntvgsXLqRx48bUr1+fwYMHc+TIkQzLX3zxRX766Sc6duxIx44dOXbsGAMHDqRevXrUr1+f\nMWPGACd3mfzuu+/SsmVLGjduTKdOnfjll18AGD16NP369aN169ZcfPHFTJky5aS41327juYtmqd2\ne9z20rYsmD8307oW5m6PjyQfplfXDtwz8ubMD6iIhFXgrvDzw7p165g6dSpt2rRh8ODBjB8/nuHD\nhzNw4EAWLlxIjRo16N+/PxN/2fBxAAANX0lEQVQmTGDYsGFhy0eOHMnzzz/Pxx9/TMWKFfnyyy/Z\nunVralPR3pB2qtAuk/fs2cPSpUsxM1555RWefvppnnvuOQBWr17N0qVLOXjwII0bN6Zr166cVv5E\nZ2x16tTh4QceZteuXZQsWZL333ufevWLbrfHUyZPZM6ChJwdXBFJpSt8oGrVqrRp0waAvn37snjx\nYtatW8f5559PjRpePzUDBgzg008/zbA8vQsuuIAffviBESNG8N5773G6310wpO0yecuWLXTu3Jn6\n9evzzDPPsHbtiSdIunXrRsmSJalYsSIdO3bkiy++SLONWrVrcfudt9OtazeuveZa6jeon6bb4+Ur\n15KweClnVa7MA/eOAkjt9nje/PdITExM0+3xsKGD2b79lzTbCNftcWxsLMuWfk7P628AMu72OD4+\nPrXb4pyIdHkRCU8JH04aOCSjgURyoly5cqxatYoOHTowceJEhgwZkjottMvkESNGMHz4cL766ism\nTZqUppvj7MQ1YNAAFi9dzPsL36dcuXJccOHFgNftcbFixYiJiaFf/8GsWPFlmuVSuj2+8+5RPPPU\n44x+6DH69R/ElEnjI657Trs9jvbyIhKemnSATZs2sWTJElq3bs2sWbNo27YtNWvWJDExkQ0bNnDR\nRReldnucUTmc6Oa4YsWK7Ny5k+LFi9OjRw9q1qxJ3759w247KSmJc87x+qafMSNtB1vz5s3jr3/9\nKwcPHiQhIYEnn3zypOW3b9/OmWeeyeZNm5k3dx7z/+192vj5521Urnw2EL1uj5s0acb+/fspWbJk\narfHl7bvkKbb49Wr8q6zrpRuj+Pi4rKeOQI57XBMz+1LYVHgEn7ZstF9sqZs5r0jA15f9+PGjWPw\n4MHUqVOHW265hfj4eKZPn851113H0aNHad68OcOGDaNEiRJhywGGDh1Kly5dqFKlCmPHjmXQoEEc\nP34cgCeeeCLstkePHs11111HuXLluOyyy9i4cWPqtAYNGtCxY0d27tzJ/fffT5UqVUhKTjvgSp/e\nfdi9azdxcXE8/8LznHFGWcAV+W6PHx87KU+3JVIUqXvkAmr06NFpxrlNoe6RI7N9y0ZunLY+Tdm8\nfz8CQLcr78/VOjO7wlf3yIWPukcWEZFCr8A16Yhn9OjR+R2CiBQxWV7hm9k0M9tuZmtCykab2VYz\nW+n/XBVJEKeyWUmCyzmHzjQJsuw06bwKhGvFGuOca+T//Cu3AcTHx7Nr1y4lfclTzjmOJu9n487k\nrGcWKaKybNJxzn1qZtXzKoBzzz2XLVu2sONUdnpTiB3+7TDFYoplOD052cjjLnHyzZ6Dv+Z6WQds\n3JnMMx/8FL2ARAqZSNrwh5tZf2A5cIdzbk+4mcxsKDAUoFq1aidNj4uL4/zzz48gjGB5b8N7VCqV\nce9yX3x8OuUqFM0vKqXvx15Ecia314ITgAuBRsA24LmMZnTOTXbONXPONasUzW4wRUQkR3KV8J1z\nvzjnjjnnjgNTgBbRDUtERKItVwnfzM4OedsdiHz0EBERyVNZtuGb2WygA1DRzLYADwIdzKwR3r2w\nREAdlIuIFHDZeUrnhjDFU/MgFhERyUNF9AE+ERFJTwlfRCQglPBFRAJCCV9EJCCU8EVEAkIJX0Qk\nIJTwRUQCQglfRCQglPBFRAJCCV9EJCCU8EVEAkIJX0QkIJTwRUQCQglfRCQglPBFRAJCCV9EJCCU\n8EVEAkIJX0QkIJTwRUQCQglfRCQglPBFRAJCCV9EJCCU8EVEAkIJX0QkIJTwRUQCQglfRCQgYvM7\nABHJnpFzVuZ3CFLI6QpfRCQglPBFRAJCCV9EJCCU8EVEAiLLhG9m08xsu5mtCSkrb2YfmNl6/3e5\nvA1TREQilZ0r/FeBLunKRgELnXMXAwv99yIiUoBlmfCdc58Cu9MVdwNm+K9nANdGOS4REYmy3D6H\nf5Zzbpv/+mfgrIxmNLOhwFCAatWq5XJzUpDp+XCRwiHim7bOOQe4TKZPds41c841q1SpUqSbExGR\nXMptwv/FzM4G8H9vj15IIiKSF3Kb8N8BBvivBwDzohOOiIjklew8ljkbWALUNLMtZvYH4EngcjNb\nD3Ty34uISAGW5U1b59wNGUz6XZRjERGRPKRv2oqIBIQSvohIQKg//CIup8/Ij+3VKI8ikVBB/+5C\nbuqvczNyusIXEQkIJXwRkYBQwhcRCQglfBGRgFDCFxEJCCV8EZGAUMIXEQkIPYcvIhEL+vcKCgtd\n4YuIBIQSvohIQCjhi4gEhBK+iEhAKOGLiASEEr6ISEAo4YuIBISew5c09Dy1SNGlK3wRkYBQwhcR\nCQglfBGRgFDCFxEJCCV8EZGAUMIXEQkIJXwRkYDQc/giEcrsuwtJlQ9kOY/IqaIrfBGRgFDCFxEJ\nCCV8EZGAUMIXEQmIiG7amlkisB84Bhx1zjWLRlAiIhJ90XhKp6NzbmcU1iMiInlITToiIgER6RW+\nA943MwdMcs5NTj+DmQ0FhgJUq1Ytws0VHdVHLcjxMolPds2DSEQkKCK9wm/rnGsCXAncamaXpp/B\nOTfZOdfMOdesUqVKEW5ORERyK6KE75zb6v/eDrwNtIhGUCIiEn25TvhmdpqZlUl5DVwBrIlWYCIi\nEl2RtOGfBbxtZinrmeWcey8qUYmISNTlOuE7534AGkYxFhERyUN6LFNEJCCU8EVEAkL94YtIoXCq\nxhR4sNt5OV6msHyvRlf4IiIBoYQvIhIQSvgiIgGhhC8iEhBK+CIiAaGELyISEEr4IiIBoYQvIhIQ\nSvgiIgGhhC8iEhBK+CIiAaGELyISEEr4IiIBoYQvIhIQSvgiIgGh/vALkeqjFnDIllOMshnOk/x9\nNYqdlnwKoxIpWh6Yt4Z4dwyAYweLUyphZz5HFD26whcRCQglfBGRgFDCFxEJCCV8EZGAUMIXEQkI\nJXwRkYBQwhcRCQglfBGRgFDCFxEJCCV8EZGAUMIXEQkIJXwRkYCIKOGbWRczW2dmG8xsVLSCEhGR\n6Mt1wjezYsA44EqgDnCDmdWJVmAiIhJdkVzhtwA2OOd+cM79CvwD6BadsEREJNrMOZe7Bc16Al2c\nc0P89/2Als654enmGwoM9d/WA9bkPtwCryJQdDrPPllRrl9RrhuofoVdTedcmUhXkucDoDjnJgOT\nAcxsuXOuWV5vM7+ofoVXUa4bqH6FnZktj8Z6ImnS2QpUDXl/rl8mIiIFUCQJ/7/AxWZ2vpkVB3oD\n70QnLBERibZcN+k4546a2XDgP0AxYJpzbm0Wi03O7fYKCdWv8CrKdQPVr7CLSv1yfdNWREQKF33T\nVkQkIJTwRUQCImoJP7vdLJhZDzNzZtYspOyv/nLrzKxztGKKltzWzcyqm9lhM1vp/0w8dVFnX1b1\nM7OBZrYjpB5DQqYNMLP1/s+AUxt59kRYv2Mh5QXyoYTsnJ9mdr2ZfW1ma81sVkh5gT5+Edat0B87\nMxsTUofvzGxvyLScHzvnXMQ/eDdtvwcuAIoDq4A6YeYrA3wKLAWa+WV1/PlLAOf76ykWjbgKQN2q\nA2vyuw6R1g8YCLwcZtnywA/+73L+63L5Xado1c+fdiC/6xCF+l0M/C/l2ABnFobjF0ndisqxSzf/\nCLyHY3J97KJ1hZ/dbhYeAZ4CkkPKugH/cM4dcc5tBDb46ysoIqlbYRBJFxmdgQ+cc7udc3uAD4Au\neRRnbhX1LkCyU7+bgHH+McI5t90vL+jHL5K6FQY5PTdvAGb7r3N17KKV8M8BNoe83+KXpTKzJkBV\n59yCnC6bzyKpG8D5ZvY/M/vEzNrlYZy5ld3938PMVpvZG2aW8oW7gn7sILL6AcSb2XIzW2pm1+Zp\npLmTnfrVAGqY2Wd+PbrkYNn8FEndoGgcOwDM7Dy8FpCPcrpsqDzvWgHAzGKA5/E+OhcpWdRtG1DN\nObfLzJoCc82srnNu36mMMQreBWY7546Y2c3ADOCyfI4pmjKr33nOua1mdgHwkZl95Zz7Pt8izZ1Y\nvKaPDnjfiP/UzOrna0TRE7Zuzrm9FI1jl6I38IZz7lgkK4nWFX5W3SyUwes4LcHMEoFWwDv+zc2C\n3kVDruvmN1PtAnDOfYnXXlfjlESdfVnuf+fcLufcEf/tK0DT7C5bAERSP5xzW/3fPwAJQOO8DDYX\nsnMMtgDvOOd+85tNv8NLkgX9+EVSt6Jy7FL05kRzTk6XPSFKNx9i8W4anM+Jmw91M5k/gRM3NuuS\n9qbtDxSsm7aR1K1SSl3wbsxsBcrnd51yWj/g7JDX3YGl7sSNo414N43K+a+LUv3KASX81xWB9WRy\nU60A168LMCOkHpuBCgX9+EVYtyJx7Pz5agGJ+F+U9ctydeyiGfxVeP9dvwfu9cseBv4vzLypSdF/\nf6+/3Drgyvw+ENGqG9ADWAusBFYA1+R3XXJTP+AJvx6rgI+BWiHLDsa70b4BGJTfdYlm/YBLgK/8\n8q+AP+R3XXJZP8Nrdvzar0fvwnL8clu3onLs/PejgSfDLJvjY6euFUREAkLftBURCQglfBGRgFDC\nFxEJCCV8EZGAUMIXEQkIJXwRkYBQwhcRCYj/B4f0gf1S3kquAAAAAElFTkSuQmCC\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[]}}]},{"cell_type":"markdown","metadata":{"id":"V8jrPKUmnOyr","colab_type":"text"},"source":["## Practical conclusions about polls\n","Practical outcome: in the worst case ($\\phi=0.5$) the 95% confidence interval is $\\pm 1.96\\sqrt{\\frac{0.5(1-0.5)}{n}} = \\pm \\frac{0.975}{\\sqrt{n}}$. To get the usually acceptable polling error of 3 percentage points, one needs to sample 1056 persons. Polling companies typically ask between 1000-3000 persons."]},{"cell_type":"markdown","metadata":{"id":"vvv3EZBfnOys","colab_type":"text"},"source":["Questions:\n","1. How critical is the IID sampling assumption?\n","2. What do you think is a larger problem: approximating the PDF with a Gaussian distribution, or people lying in the questionnaire?"]},{"cell_type":"markdown","metadata":{"collapsed":true,"id":"0ISDMIYznOyu","colab_type":"text"},"source":["# Bayesian reasoning\n","\n","We will treat $\\phi$ - the unknown fraction of A supporters in the population as a random variable. Its probability distribution will express *our subjective* uncertainty about its value.\n","\n","We will need to start with a *prior* assumption about our belief of $\\phi$. For convenience we will choose a *conjugate prior*, the Beta distribution, because the formula for its PDF is similar to the formula for the likelihood."]},{"cell_type":"code","metadata":{"id":"kYPNSAgGnOyw","colab_type":"code","outputId":"3ba215b5-74bb-48d9-9b57-fefc132be3fd","executionInfo":{"status":"ok","timestamp":1570655099685,"user_tz":-120,"elapsed":634,"user":{"displayName":"Jan Chorowski","photoUrl":"","userId":"18100051538790809279"}},"colab":{"base_uri":"https://localhost:8080/","height":298}},"source":["support = linspace(0,1,512)\n","\n","A=1\n","B=1\n","\n","plot(support, scipy.stats.beta.pdf(support, A,B))\n","title(\"Prior: Beta(%.1f, %.1f) distribution\" %(A,B))"],"execution_count":101,"outputs":[{"output_type":"execute_result","data":{"text/plain":["Text(0.5, 1.0, 'Prior: Beta(1.0, 1.0) distribution')"]},"metadata":{"tags":[]},"execution_count":101},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFgFJREFUeJzt3X+QZWV95/H3R4fBRMARZnT5JSO7\nYBgNiWyLqEGIplwgQSIxCuIPKHepFc3WZuNmMcZCcd2sMVoWpdFghUzwB6LGUGPERVclRCMubRDC\nYMBBUWaGZUbJoCNmEfzuH+c0Xtruvre773TbPO9XVRf3nvOc83yfe29/zjnPuT2kqpAkteERy12A\nJGnpGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9BuQ5FNJXrHcdYxLknVJ/inJzy13LQuV5NQkl89z\nm41J/nv/+Pgkt4yxngc/I0nOTvKFMe77rCSfHtf+tDiG/gqU5PYkP0yyO8ldfRjsM1v7qjq5qv5y\nD9VyYpIf97XsTrItyZvmsf2DQTYP5wMbq+qH/T5elOTvk9yb5OoR+nxJkm8l+UGSK5LsP2Ktq5N8\nrH/9K8mJQ9rvn+Sv+36+leQlU+uq6hPAk5McPUrf01XV31XVk0ao+Y1JPjDC/sbyGUmyvn9tVg3s\n+4NV9bzF7lvjYeivXKdW1T7AMcAE8IfTG6Sz4Pd48Bd3iO1VtU9fz68Ar0zymwvtd0hNewOvAAaD\n7G7gncD/HGH7JwN/BrwMeDxwL/Cn8yjhC8BLgf87Qtt3A/f1/ZwFvKfvf8plwLnz6HvsFvsZ0crj\nm73CVdU24FPAUwCSXJ3kLUm+SBdoh/fL/n2//hFJ/rA/89yR5NIkj+nXTZ2lvTLJt4HPLaCebwJ/\nD2yYWpbkF5J8JsndSW5J8qJ++bl0Yfj7/VXCJ/rl5ye5Lcn3k9yc5AUDXTwd2FVVWwf6/N9V9RFg\n+wglngV8oqquqardwBuA05PsO8LY7quqd1bVF4AH5mqb5NHAbwFvqKrd/Tab6A42U64Gfn2OfTw1\nyT/0r8PlwKMG1p2YZOvA8//WX2V9v3+Nn5vkJOAPgBf3r+8Nfds5PyM/2WXeleSefirtuQMrbk/y\nawPPB68mrun/u6vv8xnTp4uSPDPJdf2+r0vyzIF1Vyd5c5Iv9mP5dJK1c73Wmh9Df4VLcihwCnD9\nwOKX0Z1B7gt8a9omZ/c/vwocDuwDvGtamxOAo4B/1/dx4+DUxJB6jgCeBVzbP3808BngQ8DjgDOA\nP02yoaouBj4I/HF/pXBqv5vbgOOBxwBvAj6Q5MB+3S8Ci5nLfjJww9STqrqN7mz8yEXscyZHAvdX\n1a0Dy27o+5/yNWB9kv2mb5xkNXAF8H5gf+CjdAeRn5LkScBrgKdV1b5079vtVfW/gP8BXN6/vr80\nsNlcnxHoDq63AWuBC4CPjzgN9uz+v2v6Pr80rdb9gU8CFwEHAO8APpnkgIFmLwHOofu8rAZeO0K/\nGpGhv3JdkWQX3XTD39L9ck/ZWFWbq+r+qvrRtO3OAt5RVd/oz3RfB5wxbSrnjVX1g6k586o6uqo+\nNEctByXZleR7wK3Al/u6AH6DLoD+oq/neuCvgN+ebWdV9dGq2l5VP66qy4GvA8f2q9cA35+jlmH2\nAe6ZtuweuvAbp32A7w3pZ2oca2bY/jhgL+CdVfWjqvoYcN0sfT0A7A1sSLJXVd3eH8zmMtdnBGDH\nQN+X0x1oZ70qmYdfB75eVe/v+74M+Cfg1IE2f1FVt/afv48AvzyGftUz9Feu36yqNVV1WFWdNxXQ\nvTvm2O4gHnpm9y1gFd288yjbz2R7X8t+dAH2Q2DqpuBhwNP7g8Ku/kB1FvCvZttZkpcn+epA+6fQ\nnXEC/DOLC+jdwPQz6/1Y3IFkof1MjWPXDNsfBGyrh/6LiDOdkVNVW4D/DLwR2JHkw0kOGlLfsPd4\npr6H7XMU0z9/U/s+eOD54P2Se+kOoBoTQ//haa5/OnU7XRBPeQJwP3DXiNvP3XHVPXRTOVNnbncA\nf9sfFKZ+9qmqV83UV5LDgPfRTVccUFVrgJuA9E1uZHFTMZuBB6c5khxOd5Z866xbLMytwKp+umvK\nL/X9TzmK7ipo+hUBwJ3AwUkysOwJs3VWVR+qql+he28LeOvUqtk2GVL/TH1P3TP5AfDzA+sGD+DD\n9jv98ze1721DttOYGPrtuQz43SRPTPc1z6k53/vHsfN+n2fwk3D7G+DIJC9Lslf/87QkR/Xr76K7\ntzDl0XTBsbPf3zn0N6l7/wdYk+TBM8Mkj0zyKLorlkckeVSSvWYp8YPAqem+5/5o4ELg41X1/X5f\nG5NsnGN8e/d9Aazu+8r0dlX1A+DjwIVJHp3kWcBpdHP0U06guwk/ky/RHYz/U/+anc5Pprim1/Sk\nJM9J982mf6G70vpxv/ouuvsG8/1df9xA379Nd4C6sl/3Vbopwb2STAAvHNhuZ9/34czsSrrPw0uS\nrEryYrqb/n8zz/q0QIZ+ey6hC55rgG/ShcTvzLVBks1JzpqjyUH9NzV2012q7083hUMfps+jOxBs\np7t0fyvd2TXAn9PNRe9KckVV3Qy8nS707qK7cfvFqY6q6j5gI93XJqe8jC7o3kN3A/iHdFcLU/Xv\nTnJ8v/1m4D/Shf8OuimW8wb2dehgfzO4pd//wcBV/ePD+n7+IMlgiJ8H/Fzfz2XAq/r+p5xJ9/XR\nn9KP83S6m+53Ay+mO4jMZG+6r6t+h+71fRzdvRrobgADfDfJP8wxrum+DBzR7/MtwAur6rv9ujcA\n/5puqu1NdFd2U3Xf27f/Yv+eHjdtXN+lu8/ze8B3gd8HfqOqvjOP2rQI8X+iopUmyTrg74CnTruX\nsdj9rqb7hs3Rs9zcHJskpwIvq6oX7cl+pOkMfUlqiNM7ktQQQ1+SGmLoS1JDRv0HtZbM2rVra/36\n9ctdhiStKF/5yle+U1XrhrX7mQv99evXMzk5udxlSNKKkmTGv9iezukdSWqIoS9JDTH0Jakhhr4k\nNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JD\nDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWrI0NBPckmSHUlummV9klyU\nZEuSG5McM239fkm2JnnXuIqWJC3MKGf6G4GT5lh/MnBE/3Mu8J5p698MXLOQ4iRJ4zU09KvqGuDu\nOZqcBlxanWuBNUkOBEjyb4HHA58eR7GSpMUZx5z+wcAdA8+3AgcneQTwduC1w3aQ5Nwkk0kmd+7c\nOYaSJEkz2ZM3cs8DrqyqrcMaVtXFVTVRVRPr1q3bgyVJUttWjWEf24BDB54f0i97BnB8kvOAfYDV\nSXZX1flj6FOStADjCP1NwGuSfBh4OnBPVd0JnDXVIMnZwISBL0nLa2joJ7kMOBFYm2QrcAGwF0BV\nvRe4EjgF2ALcC5yzp4qVJC3O0NCvqjOHrC/g1UPabKT76qckaRn5F7mS1BBDX5IaYuhLUkMMfUlq\niKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY\n+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEv\nSQ0x9CWpIUNDP8klSXYkuWmW9UlyUZItSW5Mcky//JeTfCnJ5n75i8ddvCRpfkY5098InDTH+pOB\nI/qfc4H39MvvBV5eVU/ut39nkjULL1WStFirhjWoqmuSrJ+jyWnApVVVwLVJ1iQ5sKpuHdjH9iQ7\ngHXArkXWLElaoHHM6R8M3DHwfGu/7EFJjgVWA7eNoT9J0gLt8Ru5SQ4E3g+cU1U/nqXNuUkmk0zu\n3LlzT5ckSc0aR+hvAw4deH5Iv4wk+wGfBF5fVdfOtoOquriqJqpqYt26dWMoSZI0k3GE/ibg5f23\neI4D7qmqO5OsBv6abr7/Y2PoR5K0SENv5Ca5DDgRWJtkK3ABsBdAVb0XuBI4BdhC942dc/pNXwQ8\nGzggydn9srOr6qtjrF+SNA+jfHvnzCHrC3j1DMs/AHxg4aVJksbNv8iVpIYY+pLUEENfkhpi6EtS\nQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE\n0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9\nSWqIoS9JDTH0JakhQ0M/ySVJdiS5aZb1SXJRki1JbkxyzMC6VyT5ev/zinEWLkmav1HO9DcCJ82x\n/mTgiP7nXOA9AEn2By4Ang4cC1yQ5LGLKVaStDirhjWoqmuSrJ+jyWnApVVVwLVJ1iQ5EDgR+ExV\n3Q2Q5DN0B4/LFlv0bN70ic3cvP17e2r3krRHbThoPy449cl7tI9xzOkfDNwx8Hxrv2y25T8lyblJ\nJpNM7ty5cwwlSZJmMvRMfylU1cXAxQATExO10P3s6SOkJK104zjT3wYcOvD8kH7ZbMslSctkHKG/\nCXh5/y2e44B7qupO4CrgeUke29/AfV6/TJK0TIZO7yS5jO6m7NokW+m+kbMXQFW9F7gSOAXYAtwL\nnNOvuzvJm4Hr+l1dOHVTV5K0PEb59s6ZQ9YX8OpZ1l0CXLKw0iRJ4+Zf5EpSQwx9SWqIoS9JDTH0\nJakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+S\nGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakh\nhr4kNcTQl6SGjBT6SU5KckuSLUnOn2H9YUk+m+TGJFcnOWRg3R8n2Zzka0kuSpJxDkCSNLqhoZ/k\nkcC7gZOBDcCZSTZMa/YnwKVVdTRwIfBH/bbPBJ4FHA08BXgacMLYqpckzcsoZ/rHAluq6htVdR/w\nYeC0aW02AJ/rH39+YH0BjwJWA3sDewF3LbZoSdLCjBL6BwN3DDzf2i8bdANwev/4BcC+SQ6oqi/R\nHQTu7H+uqqqvLa5kSdJCjetG7muBE5JcTzd9sw14IMm/AY4CDqE7UDwnyfHTN05ybpLJJJM7d+4c\nU0mSpOlGCf1twKEDzw/plz2oqrZX1elV9VTg9f2yXXRn/ddW1e6q2g18CnjG9A6q6uKqmqiqiXXr\n1i1wKJKkYUYJ/euAI5I8Mclq4Axg02CDJGuTTO3rdcAl/eNv010BrEqyF91VgNM7krRMhoZ+Vd0P\nvAa4ii6wP1JVm5NcmOT5fbMTgVuS3Ao8HnhLv/xjwG3AP9LN+99QVZ8Y7xAkSaNKVS13DQ8xMTFR\nk5OTy12GJK0oSb5SVRPD2vkXuZLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS\n1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kN\nMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNWSk0E9yUpJbkmxJ\ncv4M6w9L8tkkNya5OskhA+uekOTTSb6W5OYk68dXviRpPoaGfpJHAu8GTgY2AGcm2TCt2Z8Al1bV\n0cCFwB8NrLsUeFtVHQUcC+wYR+GSpPkb5Uz/WGBLVX2jqu4DPgycNq3NBuBz/ePPT63vDw6rquoz\nAFW1u6ruHUvlkqR5GyX0DwbuGHi+tV826Abg9P7xC4B9kxwAHAnsSvLxJNcneVt/5fAQSc5NMplk\ncufOnfMfhSRpJOO6kfta4IQk1wMnANuAB4BVwPH9+qcBhwNnT9+4qi6uqomqmli3bt2YSpIkTTdK\n6G8DDh14fki/7EFVtb2qTq+qpwKv75ftorsq+Go/NXQ/cAVwzFgqlyTN2yihfx1wRJInJlkNnAFs\nGmyQZG2SqX29DrhkYNs1SaZO358D3Lz4siVJCzE09Psz9NcAVwFfAz5SVZuTXJjk+X2zE4FbktwK\nPB54S7/tA3RTO59N8o9AgPeNfRSSpJGkqpa7hoeYmJioycnJ5S5DklaUJF+pqolh7fyLXElqiKEv\nSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLU\nEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkNSVctdw0Mk2Ql8\naxG7WAt8Z0zlrBStjbm18YJjbsVixnxYVa0b1uhnLvQXK8lkVU0sdx1LqbUxtzZecMytWIoxO70j\nSQ0x9CWpIQ/H0L94uQtYBq2NubXxgmNuxR4f88NuTl+SNLuH45m+JGkWhr4kNWRFhn6Sk5LckmRL\nkvNnWL93ksv79V9Osn7pqxyvEcb8X5LcnOTGJJ9Ncthy1DlOw8Y80O63klSSFf/1vlHGnORF/Xu9\nOcmHlrrGcRvhs/2EJJ9Pcn3/+T5lOeoclySXJNmR5KZZ1ifJRf3rcWOSY8ZaQFWtqB/gkcBtwOHA\nauAGYMO0NucB7+0fnwFcvtx1L8GYfxX4+f7xq1oYc99uX+Aa4FpgYrnrXoL3+QjgeuCx/fPHLXfd\nSzDmi4FX9Y83ALcvd92LHPOzgWOAm2ZZfwrwKSDAccCXx9n/SjzTPxbYUlXfqKr7gA8Dp01rcxrw\nl/3jjwHPTZIlrHHcho65qj5fVff2T68FDlniGsdtlPcZ4M3AW4F/Wcri9pBRxvwfgHdX1T8DVNWO\nJa5x3EYZcwH79Y8fA2xfwvrGrqquAe6eo8lpwKXVuRZYk+TAcfW/EkP/YOCOgedb+2Uztqmq+4F7\ngAOWpLo9Y5QxD3ol3ZnCSjZ0zP1l76FV9cmlLGwPGuV9PhI4MskXk1yb5KQlq27PGGXMbwRemmQr\ncCXwO0tT2rKZ7+/7vKwa1470syHJS4EJ4ITlrmVPSvII4B3A2ctcylJbRTfFcyLd1dw1SX6xqnYt\na1V71pnAxqp6e5JnAO9P8pSq+vFyF7YSrcQz/W3AoQPPD+mXzdgmySq6S8LvLkl1e8YoYybJrwGv\nB55fVf9viWrbU4aNeV/gKcDVSW6nm/vctMJv5o7yPm8FNlXVj6rqm8CtdAeBlWqUMb8S+AhAVX0J\neBTdP0z2cDXS7/tCrcTQvw44IskTk6ymu1G7aVqbTcAr+scvBD5X/R2SFWromJM8FfgzusBf6fO8\nMGTMVXVPVa2tqvVVtZ7uPsbzq2pyecodi1E+21fQneWTZC3ddM83lrLIMRtlzN8GnguQ5Ci60N+5\npFUurU3Ay/tv8RwH3FNVd45r5ytueqeq7k/yGuAqujv/l1TV5iQXApNVtQn4c7pLwC10N0zOWL6K\nF2/EMb8N2Af4aH/P+ttV9fxlK3qRRhzzw8qIY74KeF6Sm4EHgP9aVSv2KnbEMf8e8L4kv0t3U/fs\nlXwSl+QyugP32v4+xQXAXgBV9V66+xanAFuAe4Fzxtr/Cn7tJEnztBKndyRJC2ToS1JDDH1Jaoih\nL0kNMfQlqSGGviQ1xNCXpIb8f9MZWi5yEP9LAAAAAElFTkSuQmCC\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[]}}]},{"cell_type":"markdown","metadata":{"id":"_U-CfKE7nOy3","colab_type":"text"},"source":["Then we will collect samples, and after each sample update our belief about $p$."]},{"cell_type":"code","metadata":{"id":"z2MbXyxPnOy4","colab_type":"code","colab":{}},"source":["n_successes = 0\n","n_failures = 0\n","phi = 0.6"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"LBaHphhMnOy8","colab_type":"code","outputId":"6330ae4f-441d-4b96-a9fd-d0073b085cd6","executionInfo":{"status":"ok","timestamp":1570655389663,"user_tz":-120,"elapsed":881,"user":{"displayName":"Jan Chorowski","photoUrl":"","userId":"18100051538790809279"}},"colab":{"base_uri":"https://localhost:8080/","height":343}},"source":["for _ in range(10):\n","    if rand() < phi:\n","        n_successes += 1\n","    else:\n","        n_failures +=1\n","\n","plot(support, scipy.stats.beta.pdf(support, A+n_successes, B+n_failures), label='posterior')\n","axvline(phi, color='r', label='True $\\\\phi$')\n","conf_int_low, conf_int_high = scipy.stats.beta.ppf((0.025,0.975), A+n_successes, B+n_failures)\n","axvspan(conf_int_low, conf_int_high, alpha=0.2, label='95% conf int')\n","title(\"Posterior after seeing %d successes and %d failures\\n\"\n","      \"Prior pseudo-counts: A=%.1f, B=%.1f\\n\"\n","      \"MAP estimate: %f, MLE estimate: %f\\n\"\n","      \"conf_int: (%f, %f)\"% (n_successes, n_failures, A, B, \n","                             1.0*(A+n_successes-1)/(A+n_successes+B+n_failures-2),\n","                             1.0*n_successes/(n_successes+n_failures),\n","                             conf_int_low, conf_int_high))\n","legend()"],"execution_count":140,"outputs":[{"output_type":"execute_result","data":{"text/plain":["<matplotlib.legend.Legend at 0x7fec6332a978>"]},"metadata":{"tags":[]},"execution_count":140},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAXcAAAE1CAYAAAAcUKCZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXl8VcX5/99P9gCBsIc97PsiIqKI\noigibnW3LmBrtVZ/2tZa6/LVoqV7rUuxblUpLnXfxRVBWUVQQPZN9kBCgJCQhGzz+2PmwsnlJrmB\nm5yb5Hm/XveVmzNz5jznnLnPmfPMzGfEGIOiKIpSv4jx2wBFURQl8qhzVxRFqYeoc1cURamHqHNX\nFEWph6hzVxRFqYeoc1cURamHqHOvJUQkT0S61dKxRorIOnfMH9XGMSOBiDwpIvf5bYcSWURkqohM\njlBZvUVkiYjkishtVeTt7H4Dse7/WSLys0jYUReod85dRDaJSIG7qbtcxWpyDOWli4gRkbhjscsY\n08QYs/FYyqgGDwJT3DHfqSuV2hhzkzHmDzV5DBGZ4O5n1F+PhoKItBaRl0UkR0T2ishLlWS/E5hp\njEkxxjxWWbnGmC3uN1AaWYvrBvXOuTvON8Y0AYYCw4D/88uQY30oHOX+XYAVx3LcIBtiI1WWn4hI\nc+AeInhtlIjwFrAT6Ay0Af5RSd6I1u1wEEvd85XGmHr1ATYBZ3r+/zvwgfveHngP2AOsB27w5BsO\nLAL2A7uAf7rtWwAD5LnPSW77T4FVwF7gE6CLpywD3AKsA37wbOvhvjcDpgFZwGbswyfGpV0HzAUe\nBrKBySHOcTgwH9gHZABTgASXtgEoAwqcvX8GSoFC9/8Ul68P8Jm7FmuAyz3lTwWeAKYDB7zX05Pn\nOmAjkAv8AFztSavs2lR13Mnu+2hgG/AbINOd5088eVsC77v79Q0wGZhTRd14ErgZmAX8rJJ8Ic8N\nmAS86MmX7u5rnPu/BfA8sMOd+zuevBcCS5y9G4BxnrrwrDu/7e48Yl1aD+BLIAfYDbzqtourH5mu\nvO+BAS4tEesct2Dr8ZNAsktrBXyArTd7gNm4ehfiGjwKbHXlLwZGedImAa9h63Au1tkO86QfB3zr\n0l4FXiFEPXZ5x2J/s7Fh/La/oHxd7gWcC3zn7NwKTKrk/hy672Hcy1nAH7G/xQJ3L6p9r3z1hX4b\nEPET8jh3oJOreH9w/38F/BtIAoZgnesZLm0+cK373gQYEeqme36o64G+QBzWOc/zpBusA2vh+WF5\nnfs04F0gxZW/FrjepV0HlAC3urKTQ5zj8cAIl56OdaS/CnUNgiu1+7+x+yH8xJVxnKuQ/Vz6VFdJ\nR2Lf7pKCjt8Y+2Pq7f5vB/Sv6tqEeVyvcy/BhpjigfFAPtDcpb/iPo2Afq7cCp07hx/eMcHXoxrn\nNonKHcKHWGfW3Nl8mufYOcBZ7vgdgD4u7W3gKXfcNsBC4Ocu7X/AvYF7AJzitp+NdbipWEffF2jn\n0h7GNmBaYOvX+8CfXdqfsc4+3n1GAVLBdbgG+wCNwz5gdwbqgbsOhe6exLpyF7i0BGyD5dfuGJcC\nxVTs3O/HNgBexDZmvglctwryl7t3rp4MdNdoEPaB9qMK7s+hfcO4l7OwD8j+7hrEH8298tUX+m1A\nxE/IOrY8bOtkM9aZJ2MdfSmQ4sn7Z2Cq+/4V8ADQKqi8cjfdbfsI54zd/zFYx9PF/W9wDw1PHoN9\nuscCRTiH5tJ+Dsxy368DtlTznH8FvB10DSpz7lcAs4PKeAr4vfs+FZhWyfEau+t7CUEPn8quTZjH\n9Tr3gqDrnol9qMViHUZvT1qFLXeXfxGHH9jlrkc1zm0SFTgE7EOgDPfwCXGOD4fY3hY46D0O8GNs\nTBlsI+BpoGPQfmdgGwQj8LS8sY7+ANDds+0kDr89PohtVPQ4it/VXmCw5zp87knrBxS476di31zE\nkz6Pip370+4aXo91oFe669+qgvwV3juX/kjgWnPszv3BY71Xfn7qXhwpPH5kjEk1xnQxxtxsjCnA\nhmT2GGNyPfk2Y1tRYCtXL2C1iHwjIudVUn4X4FER2ScigVdc8ZQFtiUZilbYSry5Ajsq2xcAEekl\nIh+IyE4R2Q/8yZUbLl2AEwP2u3O4GkgLxwZjzAGso74JyBCRD0Wkj6fsiq5NOMf1km2MKfH8n499\nq2qNdaheGyu7ZjcDy4wxCyrJE865VUYnbP3aW0HahhDbu2DrQobnejyFbRWC7TwUYKGIrBCRnzob\nv8CG4h4HMkXkaRFpir0ujYDFnvI+dtvBhijXA5+KyEYRuauikxGRO0Rklevk3IcNSXjr2E7P93wg\nyfUPtQe2G+fxHN66HkwBsMkY86wxptgY8wr2Xo6sZB+vnSeKyEwRyRKRHOx9q85voTK8deqo7pWf\n1FfnHoodQAsRSfFs64yNnWGMWWeM+TH2Zv0VeENEGmOf5sFsxb6OpXo+ycaYeZ48ofYDG4YoxlaW\nI+yoYt8ATwCrgZ7GmKbYTkKpJH9weVuBL4Psb2KM+UW4NhhjPjHGnIVtsa4GnvGUXdG1Cee44ZCF\nDdl09GzrVEn+McBF7mG4EzgZeEhEplTz3A5gnWeA4IdhCxFJDVHkVqB7BdsPYlupgevR1BjT39mx\n0xhzgzGmPfbt7t8i0sOlPWaMOR7bau4F/BZbtwqwYaRAec2MHVyAMSbXGPMbY0w34ALgdhEZE2yU\niIzCOqvLsW8iqdiwUmV1LEAG0EFEvHk7V5J/GUfWtarqv5eXsWGoTsaYZtiwUzh2VnYvQ9lx1PfK\nLxqMczfGbMW+Hv5ZRJJEZBC2tf4igIhcIyKtjTFl2NdCsK/ZWe6vd4z6k8DdItLf7dtMRC4L045S\nbGfUH0UkRUS6ALcH7AiTFGxcOM+1KqtyjruC7P8A6CUi14pIvPucICJ9wzm4iLQVkQvdw+8gNgxW\n5pIruzbHdNwA7hq+BUwSkUbuGkyoZJfrsHHpIe6zCBuCu7ea57YEONWNn24G3O2xKQMbkvq3iDR3\n53aqS34W+ImIjBGRGBHpICJ93D6fYh80TV1adxE5zdlymYgEHmB7sc6mzF2zE0UkHuukCoEyV3ef\nAR4WkTaujA4icrb7fp6I9HCONwcbpgycm5cU7MMzC4gTkfuBppVcXy/z3b63uWtwMbbPoSLeBpqL\nyEQRiRWRS7EP7blhHi8F+8ZUKCLDgavC3K/CexmKo71XYdpSIzQY5+74MTa2tgNbqX5vjPncpY0D\nVohIHnakwJXGmAJjTD6u19y9jo0wxryNbd2/4sIiy4FzqmHHrdgf5UZgDrb18Vw19r8DW4lzsT/m\nV6vI/yhwqdgxxI+50NRYbHxzB/YV+6/YkRbhEIN9IO3Ahl1Owz1gKrs2ETiul/+HDRXsBF7Admgd\nDJXRGLPPtax2GmN2Yvs89htjcqp5bp9hr/UybIfmB0H7Xot9K1uN7R/4ldtvIbYT+WGsU/2Sw29u\nE7CdkCuxTuEN7BsDwAnA165Ovgf80ti5Ek2x930vNuSRjQ25APwOG3pZ4K7/50Bvl9bT/Z+HdcL/\nNsbMDHENPsGGc9a68gupIlQYwBhTBFyMfaDuwYa43qok/x7sW8Qd2GtzF3ChMWZ3OMfDhtweFJFc\nbOfsa2HaWdW9DMXR3CvfkPKhMUWpm4jIX4E0Y8xEv21RlGigobXclXqCiPQRkUFiGY4Nsb3tt12K\nEi0c0+xJRfGRFGwopj22T+Eh7DA/RVHQsIyiKEq9RMMyiqIo9RB17nUEEflIROpVZ6FYdUZfxwIr\nSn1FnbtPSDWliY0x5xhj/lubNirlifTDyHUGbxSRlcdQxuUiMk9E8kVkVhj5rxKRzSJyQETeEZEW\nYR5ntIiUufqaJyLbReSBo7T5DyLyvYiUiMikKvKKiPxVRLLd569unL5SBerc/aVKaWJXuY/6Pskx\nSg4rNcqp2BnR3UTkhKMsYw9WT+UvVWV0E8uewo7Hb4uVDfh3NY61w80obgKcAlwvR7cYzHrsDNgP\nw8h7I/AjYDBWGOx87AxQpQrUuUcBxpjt2NmNA+DQijF/FJG52B9gN/EsuOFmx/2fa4Flisg0N8vO\nu7jI9SKyBSuTWg7XCtsmIveIyG73FnG1J328iKwUu9rNdhG5w5N2ntiVcPa5FuMgT1q5lq0ErcAj\nIr8VkQwR2SFB2htiZ7JOE6sRstmdX4X1U0T6i8hnIrLHvfnc47Ynisgj7hg73PdEl3adiMwJKueQ\nzc7ex8XqyeSKyNci0t2lfeV2WeparleISCuxGj/7nB2zq/kgnogd4TPdfa82xpjPjTGvYSddVcXV\nwPvGmK+MMXnAfcDFUl6SI9zj/oCd8d3vKPb9rzHmI+wkvKqYCDxkjNnmficPYSdIKVWgzj0KEJFO\nWPnU7zybr8W2WlI4UnjpOvc5HSsr0AQrJOXlNOyU+7MrOGwaVmCpA/YH9LSIBGYyPovVh0nBPnC+\ncHYeh51J+3OsHOxTwHsB51nFOY7DzkI8CztT8sygLP/Czjjt5myfgJ3VGaqsFOxMy4+xQyF7ADNc\n8r1YtcQh2NbecKq3WMuVWGmC5tgW5h8BjDEBKYHBrvX6KlYKdxtWmKstVuPHOBv/LSIVtopFpBFW\nDvcl97lSRBI86f8Wj8Ba0GdZNc7HS39gaeAfY8wG7GzdXtUtSER6YsW9Fni2LavE5uq8IVRos/ve\n/yjLaliEkorUT81/qECa2KXNwiM36tkWkCudAdzsSeuNnfYe0Hc3QLdKjj0aq//R2LPtNeA+930L\n1oE3DdrvCZw2vmfbGg7rlhs8crKUl/B9DviLJ61XID9VyCCHsP/HwHcVpG0Axnv+PxurOgj2gTgn\nKP8hm529//GkjQdWh8rr/j8WCd1rcNotWP3vHOCiY6hPP6voennyzABuCtq2HRgdRvmjsVop+7C6\nRgYrK5BwDDa/iGdxjQrylOK0793/Pd2xQ+rQ6+fwR1vu/hJKmjhAZVoe7TlSMjgO23oMZ3+AvcbK\n23rLaO++X4J1bJtF5EsROclt7wL8RspL9nby7FcZ7YNs8tpfqQyyWAnVQEfeKCqW0A0cJ7iccOwL\nECxlW9n6u2FL6IZgIvCaMabEGFMIvMlRhmaqQR5HCoA1JbzwCNiYe6qxSqSpWAXKmu7kD7a5KZBn\nnKdXKkade/RSWeXdwZGSwSXYmZrh7A9Wia9xUBk7AIwx3xhjLsR29r3DYTGmrcAfTXnJ3kbGmP+5\n9HwqllHNoLwsr1cGtlIZZGNMf+M68owxs50dXpVLL6GuTSAeXU7mVUQq0pEPCxOmhG4wYtUDzwCu\nkcMyxJcC40WklcvzpOeBFvw52jVEV2BDVQE7umFF29ZWtyBjRddexnZwBspbUYnNT0bCZvdd18AN\nA3XudZP/Ab8Wka5ih0/+CbtmY0kV+wXzgIgkuNbwecDr7v+rRaSZMaYY+woekC59BrhJrNysiEhj\nETnX0yG3BLhKrHTrOGzsPMBrwHUi0s/Fm38fSDDVl0H+AGgnIr9yHagpInKi59r8n4i0do7yfk85\nS4H+IjJERJKwq/FUh3LSyRK+hG4w12Idam8OyxD3wsbvfwxgjLnJ80AL/hyKObtrnYR9c4sRK2cd\nX8FxXwLOF5FR7sH+IPCWcQvYuA7lqeFcCFfvrsTjaIMewsGfmzz7xjubY7CSwklS8SLs07APzQ4i\n0h7bzxGWjQ0ev+NCDfVD0FJ4QWmzCFpKjPIx9xis09qKjdu+yOG1RdMJWhYwRPmjsY7kXmyreQuH\n149NwHZU7uXw4tOnePYd57YFFud+Hbd0IXY45wrsa35AhneyZ9+7sGGPHdhFtL3x7ubuPLLced1P\nBYs3u/wDsDHkva7Mu9z2JOAxZ1uG+57k2S9wzluxce/gmPvk4Ovk+f8mV+Y+7EIWv3b38YC7nvd5\n8j4JPFmB7auBW0NsvxNYVM16dJ07B+9nqic9j/KLW1/l7vcBbH9BC0/aDDyLxoeoM2UcXig+GzuU\n8Wj6G6aGsPk6lzYKG3YJ5BXgb9ghn3vcd423h/FRbZkGiIiMxq4f2bGqvErDwI3UWQoMMvaNTanj\n6AQXRVEwdpGNaq2IpUQ3GnNXFEWph2hYRlEUpR6iLXdFUZR6iDr3KMCNDx7ttx1K/cWNm7/PbzuU\n2qNeOHexwldFgQkgnu3fiRWGSg/aPsltPzFo+3UiUuomXewXK5B1XoRtLSemBYfGB8+K5HHcsUaL\nyLYIlzlGRFaLlZid6cakV7XPae56Tw7a/ms3iWe/iDwnHo0asQJoM91xVovImbWxbxXncZ07j4eD\ntl/otk/1HN9ICEVOV/eKgyb47Avn+OEiIQTSjB03/4dIHsdzvEhLIVd6/yrYp4VY0blgYbgK66ub\nI/Gcqwc7ReT22ti3tqgXzt3xA24CCICIDKT8bMnAdsGKUu1xf4OZb6ykaSpWQOs1EWleIxbXMdzD\n8y2smmALYBHwahX7xAOPAl8HbT8bO+59DHZGaTesYFeA/2GF1Fpix6a/ISKta2HfqtgAXB7kuCdS\nvVmer5ryE3xSq7FvQ6DC+1cJfwVWeTeEUV8nYbVqumBF+O4UO/mupvetHfweaB+JD3Yiyf8B33i2\n/QNbMQyQ7tl+KlYT42rsRIwET9p1eISlgMZu/2EVHPen2Aq1F/gE6OK2C/AwkImdCPQ9dtLNjdhp\n9kXYiSDve+w/032fhJ0Y9CJ2MtD32NmLd7vytgJjPTb8xNmQC2zEqjkGbC+g/MST9tgH+l1YJ5WN\nnRnaIszrfCMwL+j6FOARdgqxz13YiSdTKT9B6GXgT57/xwA73fdewEHc5Ci3bTZO9Kqm9g3j/K8D\n5mAneZ3rtrXATqL6O27yEJVMJHP398Vq1O0RWGndfdhx6KOD7Nno7v0P2DrdFyjEzpbNA/a5vIeu\nP4cnsd3p6lQGVjN9PPYhtQe4x3Oc4cB8Dk9cm4L73QBfuXM94I53hdt+HnbG8j5n/6Awz7fS+1fB\nPic7+35C+d9vpfUVO5nO+1v6A/BKTe9bW5/61HJfADQVkb5ipzJfSejp6xOB9zmsl3J+iDyBRS5+\nhq2w60KkX4iVeL0YK/k6G9viABiLfYj0wsrYXg5kG2Oexk4B/5uxLbaQx3Y2vYCdtfkd9sERgxXS\nehArtRsgE/tDaoqt3A+LyFBjRcHOwbPAgjFmB3Ar9od8GtbZ7wUe95zXMhG5qgK7giVjD2AfEiEl\nWN2r6E+dzZWW5b63FZGWLm2jcdPiPen9a3jfcJnG4be+K7EzPQ9WY/+wEJEO2Fmgk7EPkTuAN8VK\nKzTGzr49x1hp5pOBJcaYVdiZtPNN5W8FadjZvB2ws4Gfwc7YPR47S/Q+Eenq8pZiZ+O2Ak7CPhBv\nhtBSyFKFNLRULodc1f0Lvkax2IfN/+NIPaUK66t7G29HxXLCNbJvBedcI9Qn5w7WIU7AaoavwglP\nBRCraXIZ8LKxs/De4MjQzAgXA92JDfNcZKxIUjA3AX82xqwyVtPlT8AQ59CKsTrsfbDDTVcZYzKq\ncR6zjTGfuHJfxz48/uJsfgVIF5FUAGPMh8aYDcbyJfAp9sdZETcB9xq7+MFBbEvy0kCYwRgzyBjz\ncgX7NsFqqHjJcecaisewU/Lzwigr8D0ljOPU1L7h8jYwWuwCKROwzr46XC7ltc5nVpDvGmC6MWa6\nMabMGPMZ9hV/vEsvAwaISLIxJsMYUx1BrWKsCFygTrUCHjVWDG0FsBIn2GWMWWyMWWCsguUmrLM+\nrYJywbZcnzLGfG2MKTV2eciD2LcQjFVAvbmCfatbx24DvjbGLK5mWU08/4c6Tk3tW2vUR+d+FfZ1\nNdQP7iKseuJ09/9LwDlB8bwFxqodtjLGjDDGfF7BsboAj8ph6ds92HBMB2PMF9jWxONApog8LSLB\nUquV4VV3LAB2GyuuFfgfXAUTkXNEZIHYlYD2YX/45TqWQ9j9tsfuVdiWWdtK9gkQtmSsiJyPfbWu\nKNYYSsoVV1ZVx6mpfcPCWGnmD7GhwJbGmLnh7ut4zZRX1jy9gnxdgMukvMTyKUA71xq8Aqd3I3b1\nqD7VsCE7RJ0KrneBOtZL7IpTO0VkP7YhU1UdO1pp6OrUsfZY537vUZSV5/k/1HFqat9ao145d2PM\nZmzscTy2QyOYidgKu0WszOrrWB3xisIQlbEVG9/2/kiTjTHznC2PGWOOxy5D1gv4bcDMozhWSNxr\n7pvY/oW27hV8OvYhU9GxtmJf5b12Jxm7hFlVBEvGNga6E1qCdQwwTA5L2l4B/EpE3g1Vlvu+yxiT\n7dK6Sfnl37xSrzW1b3WYhlUorEi5MhJsBV4IuleNjTF/AXBvd2dhQwSrsaEViGAdczzhyu9prJb7\nPRyuYxXZXZk0dGVUdf+8DMee+0pXxx4Fhrs6F0sl9dUYsxfbf1CRnHCN7BvG+UeMeuXcHdcDZ5jy\nC1EE4pdjsPHpgMzqYGwve6hRM1XxJHC32EWHA2uAXua+nyBWFjce29FUyGEp2HKyscdIAlaPOwso\nEZFzsPH+ALuAli584LX7j4GhWS5+e2GYx3sbGwa4RKxk6/3AMmPM6hB578M+1ALX+j2s8wksnTcN\nu8ByPxdi+j+clKsxZi22M+73YuVgL8IujvxmTe7rrscsEZkUxrX4Ehv++1cleRKdDYFPdX9vL2Il\nes8WJ+0rdnhrRxFpK3YIZmNsyCOP8nWso3iW7TtGUrADA/Lc28EvgtKD63RV0tAVEsb98/IRtvM6\nUMfux/ZRDXFvJVXV12lYeejm7rxu4HBdqMl9a4fq9sBG44cK5HOxwmgGWwHuAhaHyNMeG38cQIhl\n2Ko47rXY0Sz7sa2V59z2McAy7A9uNzb808Sl9eTwKIJ3gu0naDQFdq3RTSHOqaP7/xbsj2sfNiz1\nCuVHpTyHHRWzj8OjZW7HLo+Xi+3o8Y4eWQFcXck5n4ltxRVgZYjTPWmVydxO9drltt3ubN8PPA8k\netLSXfkFztYza2nfDcBZFZxDhfUD2+k51XP8YElb467dJFff8oI+bSoo90Tsg2QP9iH+IXYBknZu\ne467t7NwyxRiH/ofun12B19/jpQyPvQ78WybA1zjvp/q7nkeduDAg5QflVJOCtltq0wausJ6UtX9\nw44IWhHu/aHy+pqI/X3sd/Xh9trYt7Y+qi2jKA6xKyS9Zow52W9bFOVYUeeuKIpSD6mPMXdFUZQG\njzp3RVGUeog6d0VRlHqIOndFUZR6iDr3Bowbg/y8iOwVkYWV5BslImtq2JazReSdmjyGEl2IyK0i\n8le/7aivqHNv2JyCnYjT0RgzvKJMxpjZxpje4RQoR68h/0fgL55ywtb0FpHLRWSeyzsrRHqsiEwW\nkR0ikitW5z/VpT0p5bXVD4pIqKnuPUWkUEReDNp+q4j8IFbXe5GInOJJS3Tl7xIrD/G+m0wXFlJN\nTXAR+aWz5YCIrBKRXm67iMi9IrLF2fmKeOQwROQfIrLOXZvVIjIhqNwhIrLY2bFYRIZ40lJF5L8i\nkuk+k4L23SQiBZ7r+6kn+RngahFpE+41UcJHnXvDpgt2gtSBKnPWICJyAtDMGLPAs7k6mt57gEfw\nPByCeACrmngSVuPjWuysYYxdxOKQtro77ushyngcOynHa/eJ7piXYtU/n8Xq9sS6LL90xxzEYQXO\nyma0esuulia4iPwMOzv7XKzExnnYCXRgZ2BfC4x0diQH2XEAq0TaDCvR8aiInOzKTcCqXr6IVSn9\nL/CuZ/brw9h1E9KxcgDXishPKM/5nmt8aAa1MaYQO8v0aGaIK1VR27Om9HN0H6zw0lvYWYrZwBS3\nPQY7hX4zVv53GtZRwuGZkhOBLdgf+70u7XrK634/UMmxR1N+RuMmrPzsMuwMyVex8rEhNeTDOLf7\ngf94/q+2prfL8zNgVtC25s6O7mHY0Rg7a/e0oO1XYiWiJ1F+9vAVwMKg/Q1W2AusJsvfPOnnAmvC\nvN9ha4K7OrAVGFNBWW8Av/X8f7K7940qyP8e8Bv3fSxWXVU86VuAce77buAET9o9WFVTb105Yva4\nJ/1qYKbfv6/6+NGWex3AtQQ/wDrwdKwG9ysu+Tr3OR2r79EEq0jp5RSgN1YW4X4R6WuMeZbyut+/\nr6ZZl2OnmHfFtkyvMxVoyIvIKVL5UnIDsdPMA1RL07sKBmKVQC8VKyi1VkRuqSDvJdiH51eBDS58\n8SBWsiCYj4BYsRoqsVjt+iVYuWiwLfmRItJerNz01W6fcKiOJnhH9xkgIltdaOYBKa9lI0HfE7FS\nGOUQkWTgBA6LXPXH6qJ4ZzsuC7IjuOwBQcW+JHYJvE9FZHBQ2irKC3ApEeKINR6VqGQ49nX6t8Zq\nvIPV/gDrMP5pjNkIICJ3A8uDXo0fMFamdqmILMX+mMotSXYUPGbs4h+IyPtY4aaQGGPmYJctrIhU\nysuhVqSHHXa82kNHbLihF/ZB1BOYISJrjdVH9zIRmBbkyP4APGuM2SZyhBBiLlbQag7Wqe3DKm4G\n9l+HbVFvx74hfY9dVCIcmmAfNF4q0gTv6P6OxT7MUrG6/tuwce2PscvAvYYNDf3O5T9iGUqs7stS\n7AIxATsq0yb/GLhLRCZiZaN/GlTu1cC32OvzS+ATEeljjAk87HOx90eJMNpyrxt0AjZ7HLuX9tgW\nfYDN2Ie2V599p+d7PocXGzgWIlnmXso7rUjqYQe0yh80xhQYY5Zh33rGezOJSGds+GmaZ9sQrABU\nuQWxPVyPVbnsjxXrugb4QKzOONg4fSK236AxNqwWbsu9OtcgcI5/M8bsM4cX1Aic43PYvoRZ2BZ5\nYHGQch3fIvJ3bKv7cs8Dqio7bnPHX4eNzf/PW64xZq677vnGmD9jH4DexWRSOPLhoUQAde51g61A\nZym/KHOAHdiO0QCdsWGIXSHy1gZHI1a0DNuyDlAdTe9wyg62K5SN1wJzA29AjtHYMFhA//8O4BIR\n+dalDwE+MMasNXalpI+xCogne9KnGmP2GLvq1b+weuOVLXQRoDqa4Guw6/KGPEdn2++NMenGmI6u\njO14VioTkQewIbWxxpj9QXYMkvKvLYMCdrhzu9oYk2aM6Y/1KRUOq3V2ecvqS/nl6pQIoc69brAQ\n6zT+IlYbO0lERrq0/wG/FpHv8pEvAAAgAElEQVSuItIEu0rOqxW08muDUBryVTEdz7Jtpnqa3oGh\njknYN5YYt0+8K2sDtjP2Xjc0sS+2g/SDoGIm4NF1dzyNdagBvfAnsVK6Z7v0b4BzRaSbG254FvYh\ntdyTPkGs1n88dt3RHcaY3c7uqSISfMwAYWuCG2PysZ3ad4pIilh1yxsD5ygiLUSku7OxH/BP7JtM\nmUu/G7tgzZnmyEVLZmFDSre56xcIK33h9u0uIi3dPTjHHXeyS+ssIiNFJMHdk99iV3Dyrlx1GuG/\nzSjVwe8eXf2E98G2yN/BjpTZjY15g31A349t3Wfhhqy5tHRsSynOU84s4Gfu+3WEoV9P6NEyXo3t\nSZQfRRKsIT8KyKviGN8AJ3r+TydMTW93HsHa6VM96R2wseE8YCN2BS3vsU/CDgdMqcLG4PMUbGfr\nFmyYYhVwrSe9JVbLP9NdiznAcE/6DOCGSo4XtnY+NlTyirNjq6sTAdXXXu4a5mPDdsHa44bDC34E\nPvd40o8DFjs7vgWO86Rdjn17zMc+kM/2pPXHvjkdcPVhBjDMk56EDeG09fv3VR8/KvmrRAUiMha4\n2RjzI79tqQ3cOPGlwCBjF6lucIjIrUAnY8ydfttSH1HnriiKUg/RmLsCgIjcI+Wn4Qc+Gg9VlDqI\nttwVRVHqIb5NYmrVqpVJT0/36/BRQ0FRqd8mKFFC4sZ1ABzsdsTE0RojOSG26kxKVLF48eLdxpiK\ndJYO4ZtzT09PZ9GiRX4dPmr4fpvO31AsXS87F4AfXv+w1o45sKNODq1riMjmqnNpzF1RFKVeos5d\nURSlHqLOXVEUpR6iqpCKopSjuLiYbdu2UVhY6LcpDZqkpCQ6duxIfHz8Ue2vzl1RlHJs27aNlJQU\n0tPTCSFzrNQCxhiys7PZtm0bXbt2PaoyNCyjKEo5CgsLadmypTp2HxERWrZseUxvT+rcFUU5AnXs\n/nOs90DDMkrUUFpm2LInnw2ZeeQUFlNUUkZqo3jSmibRt11TkuJ1wo2ihIs6d8V3tu3NZ/r3GXy1\nbjc5BaEFEuNihEEdUzl3YBrHd2lBbIy2LJXKeeedd+jVqxf9+vWr1n7vvfceK1eu5K677qohy2oH\nde6Kb+zLL2Lags3MWLWLGBFO7NaSE7u2oFebFFo2SSAuRthXUMyWPfks2bqPr9Zm8YcPV9GlRSN+\nNqobQzpVtiyr0tB55513OO+886rl3EtKSrjgggu44IILqrVPXFz0udLos0hpEMzfmM3jM9dz4GAJ\n5w1qz2XHdyS1UcIR+Vo1SaRVk0SGdm7OhBFdmLchm2kLNnHfu8s5o3cbfn5aNxolaDWub2zatIlx\n48Zx/PHH8+2339K/f3+mTZvG/PnzueOOOygpKeGEE07giSeeIDExkbvuuov33nuPuLg4xo4dy8UX\nX8x7773Hl19+yeTJk3nzTbuI1y233EJWVhaNGjXimWeeoU+fPlx33XUkJSXx3XffMXLkSAYNGsSi\nRYuYMmUKmzZt4qc//Sm7d++mdevWPP/883Tu3PmIff75z3/6fMWORH8VSq1SZgwvLtjM64u30b11\nY/500UA6t2gU1r5xsTGc2qs1I7q15LVFW3l98VZW7dzPveP70qVl4xq2vGHywPsrWLljf9UZq0G/\n9k35/fn9q8y3Zs0ann32WUaOHMlPf/pT/vnPf/LUU08xY8YMevXqxYQJE3jiiSe49tprefvtt1m9\nejUiwr59+0hNTeWCCy7gvPPO49JLLwVgzJgxPPnkk/Ts2ZOvv/6am2++mS+++AKwwz/nzZtHbGws\nU6dOPWTDrbfeysSJE5k4cSLPPfcct912G++8884R+0QjOlpGqTVKSsv4+ydreH3xNs7un8bfLx0c\ntmP3khAXwzUjuvCniwZSWFzK795cxtKt+2rAYsVPOnXqxMiRdqnga665hhkzZtC1a1d69bJrqU+c\nOJGvvvqKZs2akZSUxPXXX89bb71Fo0ZH1qm8vDzmzZvHZZddxpAhQ/j5z39ORkbGofTLLrsspJOe\nP38+V111FQDXXnstc+bMqXKfaEFb7kqtUFJaxt8+WcP8jdn85OR0LjquwzEP9erfvhn/uHQwD3yw\nkknvr+Ce8X05Ib1FhCxWgLBa2DVFcP1ITU0lOzt4/W6Ii4tj4cKFzJgxgzfeeIMpU6YcapEHKCsr\nIzU1lSVLloQ8VuPG1X/zO5p9ahNtuSs1TpkxPPz5WuZvzOaGUd24eGjHiI2jbtM0ib9eMoj0Vo35\n0/RVLNq0JyLlKv6zZcsW5s+fD8DLL7/MsGHD2LRpE+vXrwfghRde4LTTTiMvL4+cnBzGjx/Pww8/\nzNKlSwFISUkhNzcXgKZNm9K1a1def/11wM4ADeSrjJNPPplXXnkFgJdeeolRo0ZF/DxrCnXuSo3z\n4oLNfLVuNxNPSueCwe0jXn6TxDj+cMEA0ls25s8fr2btrtyIH0OpfXr37s3jjz9O37592bt3L7/+\n9a95/vnnueyyyxg4cCAxMTHcdNNN5Obmct555zFo0CBOOeWUQ52bV155JX//+9857rjj2LBhAy+9\n9BLPPvssgwcPpn///rz77rtV2vCvf/2L559/nkGDBvHCCy/w6KOP1vRpRwzfltkbNmyY0cU66v9i\nHV+s3sXDn6/j7P5p3DK6e43OfNyXX8RvXl9KUWkZ/7h0MG2bJtXYsWqCaFmsY9WqVfTt27fWbAjF\npk2bOO+881i+fLmvdvhNqHshIouNMcOq2ldb7kqNsTn7AI/P2sCgDs246dRuNT6lPbVRApPO709x\naRkPfLCS/KKSGj2eokQz6tyVGqGwuJS/frKGRvGx/GZsb+Jia6eqdWrRiLvP6cv2vflMmbkeXQC+\nbpKent7gW+3Hijp3pUZ4ZvZGtu3J5/azetGi8ZGTk2qSwR1TuebELsxet5uPV+ys1WMrSrSgzl2J\nOIs37+XTlbu4ZGhHjuvc3BcbLjm+I0M7p/LM7I1szMrzxQZF8RN17kpEyS8q4fFZ6+nYPJmrTuzs\nmx0xItx+Vm9SkuJ56LO1FJeW+WaLoviBOnclorwwfzO7cw9y2xk9ia+lOHtFNEuO59bTe7BlTz7/\nW7jFV1sUpbZR565EjNU79/Ph9xmcN6gdfds19dscAIalt2BMnza8+e021un4d6UBUaVzF5FOIjJT\nRFaKyAoR+WWIPKNFJEdElrjP/TVjrhKtlBnD019tpHnjBK4dke63OeX42ahuNG+UwCMz1lGi4Rml\ngRBOy70E+I0xph8wArhFREIJJM82xgxxnwcjaqUS9cxcncm6zDwmnpROckJ0iSk1SYzj5tHd2bIn\nn/eW7vDbHKUKsrOzGTJkCEOGDCEtLY0OHToc+r+oqKhGjvn2229z66231kjZflGlcJgxJgPIcN9z\nRWQV0AFYWcO2KXWEgqJSps3fTK+2TRjdu7Xf5oRkeFe7EMjLC7cwqmdrWqck+m2SUgEtW7Y8JPA1\nadIkmjRpwh133HFEPmMMxhhiYo49uvztt98ydOjQYy4nmqjWVRGRdOA44OsQySeJyFIR+UhEQkrJ\niciNIrJIRBZlZWVV21glOnnj223syS/ihlHdiInihZVvGNUNA/xnzka/TVGOkk2bNtG7d28mTJjA\ngAEDmD17NgMGDDiU/o9//INJkyYB8OKLLzJ8+PBDEr+lpaVHlLd27VrOOussHnnkER544AEeeeSR\n2jqVGidsyV8RaQK8CfzKGBOs3v8t0MUYkyci44F3gJ7BZRhjngaeBqstc9RWK1HD3gNFvLtkO6f2\nbEWftOjoRK2Itk2TuGJYJ15YsJlvt+xlqE9j8OsUv/oVVCCTe9QMGQLH4ETXrVvHf//7X0aMGMGm\nTZtC5lm1ahWvvvoqc+fOJT4+nptvvpmXXnqJCRMmHMpz8OBBLr/8cl544QUuvPBC5s2bR79+/bjp\npptISqpbukShCMu5i0g81rG/ZIx5Kzjd6+yNMdNF5N8i0soYsztypirRyGuLt1JcWsbVJ3bx25Sw\nuOi4Dny2chfPz/2BwR1TdaHtOkiXLl0YMWJEpXlmzJjB4sWLOeGEEwAoKCigTZs25fJ89tlnDB48\nmPbt29O0aVPS0tJISkoK2cKvi1Tp3MWqPT0LrDLGhFwoUETSgF3GGCMiw7HhniNV9ZV6RWZuIR8v\n38mZfdvSPjXZb3PCIj42hoknp/PXj1fzxepdnNUvzW+TopsoDFN4F8mIi4ujrOzwCKjCwkLAxuMn\nTpzIn//85wrLWbp0KQMHDmTZsmUMGjSIzMxMUlJSon4RjnAJJ+Y+ErgWOMMz1HG8iNwkIje5PJcC\ny0VkKfAYcKVRxaZ6zyvfbAXgyhP8m4l6NIzs3pLebVN4ccEWCovrRyutodK2bVsyMzPJzs7m4MGD\nfPDBB4BdL/WNN94gMzMTgD179rB58+Zy+6akpLB69WqWLl3KoEGD+P3vf88tt9xS6+dQU4QzWmYO\nUOm7qzFmCjAlUkYp0c/OnEJmrNrFeYPa17mRJyLC9ad05c43l/H2d9v58fC69XBSDhMfH8/999/P\n8OHD6dChA3369AGgX79+TJ48mbFjx1JWVkZ8fDyPP/44XbocDh9ec801XHTRRbz11ls0b96cK6+8\nsl4Nh9TFOnymri7W8e9Z6/ls5S7+M2EYLZvULece4C8frWLxlr08dc2wWleuDIUu1uEPQ4YM4fPP\nP6dVq1Z+m3IEuliHUqvsPVDE56t2MaZv2zrr2AEmnJROUUkZry/e6rcpik8cPHiQnJycqHTsx4o6\nd6XavLt0O6VlhkuGdvDblGOifWoyZ/Zty8fLd5KVe9BvcxQfSExM5IcffvDbjBpBnbtSLfIOljD9\n+52c0qMV7ZrVjREylXHFsE4A2npX6h3q3JVq8eH3GRQUl3Lp8R39NiUitGmaxFn92vLZyl3s2l/o\ntzmKEjHUuSthc7CklPeX7mBYl+Z0bdXEb3MixhXDOiECry7S1rtSf1DnroTNl2uzyCko5uLj6nas\nPZiWTRI5Z0A7ZqzaRUZOgd/mKEpEUOeuhIUxhveX7iC9ZSMGdDhy+Fxd59KhHYmNEd5cvM1vUxQl\nIqhzV8Ji+fYcNmXnc/7g9kgUKz8eLc0bJ3Bm37bMWJ1Jdp6OnFHqPurclbB4b9kOUpLiOK1XdOq1\nR4KLh3akzBjeWaILenj5fltORD/h8OijjzJgwAD69+9fToZ30qRJ5RbvmD59OgBz585l0KBBDBs2\njHXr1gGwb9++QzNUa5LXX3+dvn37cvrpp5fbvmPHDi699NIq9//Tn/5UI3apc1eqZOf+Qr7euIdx\n/dNIjIuuVZYiSVrTJE7t2ZqPV2SQW1jstzkNluXLl/PMM8+wcOFCli5dygcffMD69esPpf/6179m\nyZIlLFmyhPHjxwPw0EMPMX36dB555BGefPJJACZPnsw999wTkcU8KuPZZ5/lmWeeYebMmeW2t2/f\nnjfeeKPK/dW5K77x4bIMRGD8wHZ+m1LjXHp8RwqLy/hgWYbfpjRYVq1axYknnkijRo2Ii4vjtNNO\n4623jlAaL0d8fDz5+fnk5+cTHx/Phg0b2Lp1K6NHj65wn2+++YaTTz6ZwYMHM3z4cHJzcyksLOQn\nP/kJAwcO5LjjjjvksKdOncrFF1/MuHHj6NmzJ3feeScADz74IHPmzOH666/nt7/9bbnyN23adGgh\nkYr2v+uuuygoKGDIkCFcffXVR3vJQhL2Yh1Kw6SwuJTPVu5kZI9WtKrDUgPh0qVlY07s2oL3l+7g\nR0M6RN16sA2BAQMGcO+995KdnU1ycjLTp09n2LDDUipTpkxh2rRpDBs2jIceeojmzZtz9913M2HC\nBJKTk3nhhRe44447mDx5coXHKCoq4oorruDVV1/lhBNOYP/+/SQnJ/Poo48iInz//fesXr2asWPH\nsnbtWgCWLFnCd999R2JiIr179+bWW2/l/vvv54svvuAf//hHORtDEWr/v/zlL0yZMuXQsoKRRFvu\nSqXMWb+bA0WljB9Q/1vtAS49viO5B0v4ZOVOv01pkPTt25ff/e53jB07lnHjxjFkyBBiY+1D9he/\n+AUbNmxgyZIltGvXjt/85jeAFf9asGABM2fOZOPGjbRr1w5jDFdccQXXXHMNu3btKneMNWvW0K5d\nu0OLeTRt2pS4uDjmzJnDNddcA0CfPn3o0qXLIec+ZswYmjVrRlJSEv369TtCQrgqjnX/6qLOXamU\nT1bspENqMv3bR/cSepGkT1pTBnZoxjvfbaektGY745TQXH/99SxevJivvvqK5s2b06tXL8Dqt8fG\nxhITE8MNN9zAwoULy+1njGHy5Mncd999PPDAA/ztb3/jhhtu4LHHHjtmmxITD7+5xsbGUlJSUqv7\nVxd17kqFbNp9gNU7cxnXP61eDn+sjIuO60D2gSLmrNeVIv0gsMjGli1beOutt7jqqqsAyMg43Bfy\n9ttvl1scG2DatGmMHz+eFi1akJ+fT0xMDDExMeTn55fL17t3bzIyMvjmm28AyM3NpaSkhFGjRvHS\nSy8BdvHsLVu20Lt37xo7T7D9BcXFke/A15i7UiGfrNxJXIxwep82VWeuZxzfpTkdUpN5d+kOTuvV\nusE93LyE0nyvaS655BKys7MPLbKRmpoKwJ133smSJUsQEdLT03nqqacO7ZOfn8/UqVP59NNPAbj9\n9tsZP348CQkJvPzyy+XKT0hI4NVXX+XWW2+loKCA5ORkPv/8c26++WZ+8YtfMHDgQOLi4pg6dWq5\nFndNcOONNzJo0CCGDh166MESCXSxDp+J1sU6CotLuW7qQo7v3ILfnl2zLZdo5aPlGfx71gb+cvFA\n+reveQeni3UowehiHUrEmbdhNwcOljKuf1u/TfGN03u3ISUxjnd1UpNSB1HnroTk4xW7aN8sqV7q\nyIRLUnws4waksWBjNjtzVA5YqVuoc1eOYHP2AVZl7OfsBtiRGsy5A9sRGyO8v6xhtd79CtcqhznW\ne6DOXTmCT1fuIi5GGNO34YZkArRsksionq34bOUuDhys2aFr0UJSUhLZ2dnq4H3EGEN2djZJSUlH\nXYaOllHKUVJaxpdrsxjetQXNkuP9NicquHBIB2auyeLTlTu56Lj6sQJVZXTs2JFt27aRlZXltykN\nmqSkJDp2PPr6ps5dKcfiLXvJKShmTAMc/lgR3Vs3YWCHZry/LIMLBncgNqZ+h6ri4+Pp2rWr32Yo\nx4iGZZRyzFiVSbPkeIZ2bu63KVHF+YPbk5V7kIU/ZPttiqKEhTp35RD7C4r5ZtMeRvdqTVysVg0v\nw9Nb0KpJIh9+r2qRSt2gyl+wiHQSkZkislJEVojIL0PkERF5TETWi8gyERlaM+YqNcnsdVmUlBnG\n9NWQTDCxMcI5A9JYui2HrXvyq95BUXwmnOZZCfAbY0w/YARwi4j0C8pzDtDTfW4EnoiolUqtMGN1\nJl1bNaZrqyZ+mxKVjO3XlrgYYbq23pU6QJXO3RiTYYz51n3PBVYBHYKyXQhMM5YFQKqINByN2HrA\n1j35rMvM4wztSK2Q1EYJnNKzFTNWZ5Jf1DCGRSp1l2oFVkUkHTgO+DooqQOw1fP/No58AChRzIzV\nmcQI9XqN1Ehw7sB2FBSXMmuNDhNUopuwnbuINAHeBH5ljNl/NAcTkRtFZJGILNIxtNFDaZlh5ppM\nju/SnOaNEvw2J6rp3TaF7q0b8+H3GTrJR4lqwnLuIhKPdewvGWNCLWa4Hejk+b+j21YOY8zTxphh\nxphhrVtrCzFa+H57DnsOFHFGH52RWhUiwrkD27FlTz7Lt0enoqeiQHijZQR4FlhljPlnBdneAya4\nUTMjgBxjjPY61RG+WptFcnwsJ6Tr2PZwOLVXa1IS43RYpBLVhDNDdSRwLfC9iARWcb0H6AxgjHkS\nmA6MB9YD+cBPIm+qUhMUl5Yxb8NuTureksQ4XQw6HBLjYjmzX1veXbKd7LyDtGwAC4crdY8qnbsx\nZg5Q6XxrY4OPt0TKKKX2WLR5LweKSjmtp4bJqsP4Ae1457vtfLRiJ9ec2MVvcxTlCHQaYgPny7VZ\nNEuOZ3CnVL9NqVOkNUvi+C7N+XTFTl1EW4lK1Lk3YPKLSvjmhz2c0qNVvRfDqgnGDUhjb34xCzft\n8dsURTkCde4NmAUb91BUWsapOrb9qBjWpQUtGyfw8fKdfpuiKEegzr0B89W6LNqkJNInLcVvU+ok\nsTHC2H5t+W7rPl2GT4k61Lk3UHIKivluy15G9WxNTANfSu9YGNs/jRiBT1Zo612JLtS5N1Dmrt9N\nmVG5gWOlVZNETkhvweerdlGsHatKFKHOvYHy5dosOrVoRHrLRn6bUucZ1z+NfQXFfP2Ddqwq0YM6\n9wZIZm4hKzP2c1rPVoiGZI6Z4zo3p3VKIh8v1xmrSvSgzr0BMnvdbgAdJRMhYmOEs/u1Zem2HHbs\nK/DbHEUB1Lk3SL5am0Wvtk1o1yzZb1PqDWf2basdq0pUoc69gbF9bwEbdx9glMoNRJSWTRI5sWtL\n7VhVogZ17g2MORtsSGZk91Y+W1L/GNc/jf2FJczfkO23KYqizr2hMXf9bvqkpdA6RZUMI82Qzqm0\nbZqooRklKlDn3oDYsa+AH3YfYGQPbbXXBDEinN0vjWXbc9i+VztWFX9R596AmLNeQzI1zZl92xIb\nI3ysrXfFZ9S5NyDmrt9N77YakqlJmjdOYETXFsxYvYuiEu1YVfxDnXsDYcc+O0rmFA3J1DjjBrQj\nt7CEea7zWlH8QJ17A2GuC8mc3KOlz5bUfwZ1bEa7ZkkamlF8RZ17A2HOBhuSaZOS5Lcp9Z4YEcb2\nS2PFjv1s3ZvvtzlKA0WdewNgx74CNmZpSKY2ObNvG+JihE90IQ/FJ9S5NwA0JFP7pDZKYES3lnyx\nOlM7VhVfUOfeANCQjD+M659G7sES5mrHquID6tzrORk5GpLxi4GuY1VnrCp+oM69njNHQzK+ESPC\nuP62Y3XLHu1YVWoXde71nMDEJQ3J+MMZfVzHqrbelVpGnXs9JiOngA1ZBxiprXbfSG2UwEndbcfq\nwZJSv81RGhDq3OsxqiUTHYzrn0bewRLmrlcpYKX2qNK5i8hzIpIpIssrSB8tIjkissR97o+8mcrR\nMHf9bnq1bUKbphqS8ZOBHZrRXmesKrVMOC33qcC4KvLMNsYMcZ8Hj90s5VjZmVPIBh0lExWICGf3\nT2NVxn42Zx/w2xylgVClczfGfAXsqQVblAiiIZnoYkzfttqxqtQqkYq5nyQiS0XkIxHpX1EmEblR\nRBaJyKKsrKwIHVoJhYZkootmyfGc3L0VX6zRjlWldoiEc/8W6GKMGQz8C3inoozGmKeNMcOMMcNa\nt9YFmmuKnfsLWZ+Vp632KGNc/7YcOFjKnHU6Y1WpeY7ZuRtj9htj8tz36UC8iKhX8ZF5gZCMxtuj\nigEdmtEhNVlDM0qtcMzOXUTSRETc9+GuTB3z5SNz1u+mR5smtNWQTFQhbsbqqp252rGq1DjhDIX8\nHzAf6C0i20TkehG5SURuclkuBZaLyFLgMeBKY4ypOZOVysjcX8i6zDwdJROlnO5mrH6sUsBKDRNX\nVQZjzI+rSJ8CTImYRcoxEVAg1Hh7dNIsOZ6RPVoxc00mE09OJyk+1m+TlHqKzlCtZ8xdn0331o1J\na6YhmWhlXP80DhSVHhquqig1gTr3ekRW7kHW7MrVVnuU0799Uzo2T9bQjFKjqHOvR8zboKNk6gKB\nGatrduXyw27tWFVqBnXu9Yi5G7Lp2qox7VOT/TZFqYIzerchPlZnrCo1hzr3ekJ23kFWZexnZHeV\n960LNE2OZ2R327FaWKwzVpXIo869njBvg51aoCGZusO4AWnkF+mMVaVmUOdeT5i7YTddWjSiY/NG\nfpuihEm/dk3p1DxZpYCVGkGdez1gz4EiVu7Yr632Okb5jtU8v81R6hnq3OsB8zfsxqAhmbrIGX1s\nx+pHOixSiTDq3OsBczdk06l5Mp1baEimrpGSFM8pPVoxa00WZSraoUQQde51nL35RazYkaOt9jrM\n2f3TKCguJaegyG9TlHqEOvc6zoKN2ZQZ1ZKpy/Rr15ROLRqx54A6dyVyqHOv48xdv5sOqcl0aakh\nmbqKiDB+QBoFRaUUFOmYdyUyqHOvw+QUFPP9dhuScZL6Sh3ljD5tiBEhW1vvSoRQ516HCYRkTumh\ns1LrOo0S4mjeKIGcgmL25auDV44dde51mDnrd9OuWRLpLRv7bYoSAVo0ScAYw6crd/ltilIPUOde\nR9lfUMyybfsY2V1DMvWFxLgYmiTG8dHyDEp1XKRyjKhzr6N8/YMbJaNDIOsVLZoksjuviK9/0GWI\nlWNDnXsdZc76bNo2TaR7aw3J1CeaJsXRJiWRD5dl+G2KUsdR514HySkoZum2fZyio2TqJeMHtmPZ\n9hw2Z+tCHsrRo869DjJ/QzalZYZRPVv7bYpSA5zVty0JsTF8+L223pWjR517HWT2+iw6pCbTrZWG\nZOojTZPjObVXK75YnUnewRK/zVHqKOrc6xh7DhTx/bYcRvXUkEx95tyB7TlYUsYXq3VYpHJ0qHOv\nY8xdb+V9NSRTv+nRpgl90lL4cFkGZUaHRSrVR517HWP2+t2kt2yk8r4NgHMHtmNHTiHfbt7rtylK\nHUSdex0iM7eQVRn7tdXeQBjZoxUtGyfw7tIdfpui1EHUudchAgspj+qpE5caAvGxMZw3qD1Ltu7j\nh906LFKpHlU6dxF5TkQyRWR5BekiIo+JyHoRWSYiQyNvpgI2JNOjTRPaNUv22xSllhjXP43EuBje\nXbLdb1OUOkY4LfepwLhK0s8BerrPjcATx26WEsyOfQWsz8zjVG21NyiaJMVxZt+2fLk2SxfzUKpF\nlc7dGPMVsKeSLBcC04xlAZAqIu0iZaBimb3ehmRO6aHx9obGBYPbU1pmmK6TmpRqEImYewdgq+f/\nbW7bEYjIjSKySEQWZT6GIZkAABOASURBVGVlReDQDYc567Lo264prVMS/TZFqWXapyYzvGsLpi/P\noLBYV2pSwqNWO1SNMU8bY4YZY4a1bq0t0HDZnH2ATdn5GpJpwPxoSAdyC0uYuSbTb1OUOkIknPt2\noJPn/45umxIhZq/fTYzoItgNmf7tm9KjdRPeXbJDJzUpYREJ5/4eMMGNmhkB5BhjNDgYIYwxfLkm\ni0EdU2neOMFvcxSfEBEuHNKe7fsKWLSpsi4wRbGEMxTyf8B8oLeIbBOR60XkJhG5yWWZDmwE1gPP\nADfXmLUNkFU7c9m5v5DTe2sYq6FzSo9WtE5J5I1v9cVYqZq4qjIYY35cRboBbomYRUo5Zq3JJCEu\nhhHddBHshk5cbAwXDenA07M3smJHDv3bN/PbJCWK0RmqUUxxaRmz1+3mpG4taZRQ5XNYaQCc1a8t\nzZLjeX3xNr9NUaIcde5RzKLNe8k7WMLpvdv4bYoSJSTFx3L+4PYs3ryXH3bn+W2OEsWoc49iZq7O\nJLVRPEM6pfptihJFnDugHcnxsbyhrXelEtS5Ryl5hSV8s2kPp/ZsTWyMLsqhHKZJUhznDEhjzvrd\n7NhX4Lc5SpSizj1KmbN+NyVlRkMySkguHNKB2Bjhre905IwSGnXuUcrMNZl0ap5M99a6TqpyJC0a\nJzCmT1tmrNqlgmJKSNS5RyE79xeyMmM/p/duo+ukKhVyydCOlBnD29p6V0Kgzj0K+dLph5ymE5eU\nSkhrlsRpvVozfXkGe7X1rgShzj3KMMbwxepMBrRvSpuUJL/NUaKcK0/oTElpGW99pyNnlPKoc48y\nVuzYz46cQs7q19ZvU5Q6QPvUZEb3bsP05Tu19a6UQ517lPHZyl00SojlZFWAVMLkyhM6UVJaxhvf\nautdOYw69yjiwMES5mzYzak9W5MUH+u3OUodoV2zZM7o04aPl+/UkTPKIdS5RxFfrcuiqKRMQzJK\ntbliWGdKjeGNxVurzqw0CNS5RxGfrdxFestG9GzTxG9TlDpGWrMk23pfsZOs3IN+m6NEAerco4RN\nuw+wLjOPs/q11bHtylFx5bBOGAP/W7jFb1OUKECde5Tw2apdxMUIo3up3IBydLRpmsR5g9oxY/Uu\nNmcf8NscxWfUuUcBxaVlzFydyYhuLWmaHO+3OUod5rLjO5EUH8u0+Zv9NkXxGXXuUcCCjdnkHizR\njlTlmGmaHM+lQzuycNMeVuzI8dscxUfUuUcBH36fQVrTJNVtVyLC+YPb06JxAlPnbcKugqk0RNS5\n+8ym3QdYsWM/5wxII0Y7UpUIkBQfy1XDO7N6Zy4LNmb7bY7iE+rcfWb68gziY4Uz+2pIRokcZ/Zt\nS6fmyTw/bxPFpWV+m6P4gDp3H8ktLGbWmixG9WytHalKRImNEa4/pRsZOYW8t3SH3+YoPqDO3Ufe\n/m47BcWlnDuwnd+mKPWQ47s0Z3h6C179ZqvKEjRA1Ln7hDGGafM306NNE3q1TfHbHKWecv0pXSku\nLeO/8zb5bYpSy6hz94n5G7JZn5nHuQO01a7UHO1Tk/nRkA58sSaT1Tv3+22OUouoc/eJ/8z5gRaN\nExjVS6V9lZrlsmEdadEogae/2kiZDo1sMKhz94H1mbl8sTqTa0d0ITFOpX2VmqVRQhzXjUxnXWYe\nn6zY6bc5Si0RlnMXkXEiskZE1ovIXSHSrxORLBFZ4j4/i7yp9Ydn5/xAQlwM157UxW9TlAbC6F6t\nGdihGf+dt0k7VxsIVTp3EYkFHgfOAfoBPxaRfiGyvmqMGeI+/4mwnfWG3XkHefPb7VwytAOtmiT6\nbY7SQBARbhndg6LSMp6evdFvc5RaIJyW+3BgvTFmozGmCHgFuLBmzaq/vDB/M0UlZVx/Sje/TVEa\nGB2aJ3PFsE7MXb+bhT/s8dscpYYJx7l3ALzLu2xz24K5RESWicgbItIpVEEicqOILBKRRVlZWUdh\nbt2msLiUFxZs5ow+beihC3IoPnDx0I50atGIJ77cQH5Rid/mKDVIpDpU3wfSjTGDgM+A/4bKZIx5\n2hgzzBgzrHXr1hE6dN3h9UV2MskNo7TVrvhDfGwMt57eg+y8gzw3d5Pf5ig1SDjOfTvgbYl3dNsO\nYYzJNsYE1vb6D3B8ZMyrPxwsKeXfszYwrEtzRnRr4bc5SgOmb7umXHRcBz5ZsZOZazL9NkepIcJx\n7t8APUWkq4gkAFcC73kziIh3Js4FwKrImVg/eHPxdjJyCrltTE9dRk/xnatP7ELnFo343RvL2Jev\no2fqI1U6d2NMCfD/gE+wTvs1Y8wKEXlQRC5w2W4TkRUishS4DbiupgyuixSVlPH4zPUM6ZTKqJ46\naUnxn4S4GG4/qxd7DhRx37sr/DZHqQHiwslkjJkOTA/adr/n+93A3ZE1rf7w9nfb2L6vgMk/GqCt\ndiVq6N66Cb8c05OHPlvLGX1ac9FxHf02SYkgOkO1hikuLWPKzPUM6tiM0b0bXieyEt38YnR3hqe3\n4N63l7M+M89vc5QIos69hnlz8Ta27ingtjM01q5EH3GxMTz24+NIio/l/738LYXFpX6bpEQIde41\nSEFRKQ9/vpahnVMZ07eN3+YoSkjSmiXxz8sHs3pnLg+8v9Jvc5QIoc69Bnlu7g/s2n+Qu8f31Va7\nEtWM7t2Gm07rzv8WbuH1RVur3kGJetS51xB7DhTx5KwNnNm3LSek67h2Jfq5Y2wvRvZoyb1vL2fx\n5r1+m6McI+rca4jHZqzjQFEJvxvX229TFCUs4mJjmPLjobRLTeKmFxeTkVPgt0nKMaDOvQZYvXM/\nLyzYzJXDO9NTl9BT6hDNGyfwzIRh5B8s4cZpiyko0g7Wuoo69whjjOH+d1aQkhTHb8dqq12pe/Rq\nm8KjVx7Hih053PLytxSXlvltknIUqHOPMO8u2cHCTXu48+w+NG+c4Lc5inJUnNmvLQ9eOIAvVmdy\n91vfY3R5vjpHWDNUlfDIyS/mj9NX/f/27jw4yvoM4Pj32c0mIRByGDBcSUCIHAGFBKQ6IhapgU6h\n9QIUxSlF8arHdFR0pvUYp6W2aqWMgBQrZVpBvNKK2kEFlAoV5AiHXJEjnAFJCEdCyD79410wxYQs\nYbNvdvN8Znb23bxvkufJu3ned3/v7/396NMxiVH9ax312JiIMXZgJgePVvLSwi2ktYrj8WHd3Q7J\nnAcr7iH07Psb+PbYSWaN64/XY10fTeR7cEg3SsormbZ4G/E+Dw9dl+12SCZIVtxD5NOvDzB/ZTH3\nXXsJvTsmuR2OMSEhIjwzMoeKKj8vLdyCX+Hh6+xu60hgxT0Eyk5UMentQrIvbsUvh3RzOxxjQsrr\nEX5/Ux884nTxVVUeGZptBb6Js+J+gVSVJ94ppORoJTPuyCUuxut2SMaEnNcjTL6xDx4RpnyylcPH\nT/L0iBxrfmzCrLhfoDnLd/L+2r08lt+dPh2T3Q7HmEbj8Qi/vaE3yS19TF9cxL6ySqaM6UuLWDuh\naYqsK+QFWLe7jGf/uYHBl7bh7kE2L6qJfh6PMGlYD54e0YuPv97PmFeXsf9IhdthmVpYcW+gkvJK\nJs5ZSWrLWF645XI89vHUNCPjrsxi2thcNu8v58cvf87yokNuh2TOYsW9ASqqqpkwewUHj1Yy/fZc\nUu1mJdMMXd8rnXfvu4rW8THcOnM5ry4pwu+3m52aCivu58nvVx6Zt5o1xaW8NKovl3WydnbTfGVf\nnMh791/F0B4X89yCjdw+azm7S23AsabAivt58PuVx95ay4LCfTwxrAf5Oeluh2SM6xLjfbwyth+/\nvaE3q3eWkv/iEuat2GVDFrjMinuQThf2N1cW8+CQbkywC6jGnCEijBmQwYcPDaJH+9Y8On8tt0z/\ngvV7ytwOrdmy4h6EiqpqHnhj1ZnC/vBQuwXbmNp0Sk3gjQkDmXxjb7aVHOMnUz7nyXcKrUeNC6yf\nez0OlFcwYfZK1haXMmlYd+6+5hK3QzKmSfN4hFH9M8jv1Y4XF25mzrIdzF9ZzG1XZDJxcBfaJsa7\nHWKzYMX9HD7fcpBfvbmGshNVvHJbrrWxG3MekhJ8PDWiFz+/qjNTPtnC619sZ87yHfz08vaMuzKL\nXu1tDKbGZMW9FsdPnuL5jzbx2tLtdGnTkpnj8sjpYG9EYxoi46IEnr/5Mu69tiszPyvi7a92M29F\nMXmZKdzQryPDe6eTnGDdiUPNinsN1X7lrZXF/OHfmzhQXsmdV2bxWH53u73amBDonNaS537Wm0ev\n787cFTuZ++UunninkN8UrOOa7Lbk56QzKDvNmm1CxIo7zpn6W1/t5rWl31BUcoy+Gcm8MrYfuZmp\nbodmTNRJSvBx16BLmHB1F9bvOcJ7q3dTsGYPCzfuB6Bnu9YMym5DXmYKl2ckk9YqzuWII1NQxV1E\n8oE/AV5gpqr+7qz1ccBsIBc4BIxS1e2hDTW0yiuq+M+2Q3xQuJePNx6gvPIUfTomMfXWfgzvnW7D\nmRrTyESEnA5J5HRIYtKwHmzcd4TFm0tYtKmEmZ8VMW2x008+IzWBnA6t6do2ka5tW9GtbSs6p7Uk\n3mefqM+l3uIuIl5gKjAUKAa+FJECVd1QY7PxwGFV7Soio4HJwKjGCPh8VJ6qpvR4FYePn2T34RN8\nc/AY20qOsWZXKV/vO4JfITnBR35OOqMHdKJfRooVdWNc4PEIvdon0at9EvcO7sqJk9Ws21PGqp2H\nWbWzlI17y/lw3T5qjm6QkuCjXVIL2iXFk54Uz0Wt4khu4SPp9CPBR2J8DPExXuJ9XuJiPGeem8NY\nUMGcuQ8AtqpqEYCIvAGMBGoW95HAU4Hl+cCfRUS0EW5RW7TpAM/+awN+ddrIq/2KX///udqvVFUr\nJ6qqv/f9KQk+erRrzf0/7MYVnVMZ0DkVn9e6+xvTlLSI9dI/K5X+Wd81jVZUVbP90DG27D/KjkPH\n2Hekgr2lFewtq+CrnYcpPVFFsBUn1ushNsaDR5wDi1cEEcEjztj1HhGk5jJAHceDcx0m6jpZHN2/\nE7+4unFvhAymuHcAdtV4XQxcUdc2qnpKRMqAi4CDNTcSkbuAuwAyMjIaFHBivI9L0xPxiOAN7BRP\nzWcPeEXweT0kJ/hIToglJSGW9KQ4uqS1IsUG+TImIsX7vHRPb0339Na1rq/2K+UVVZSd+O5RXnGK\niqpqKk/5v/dcWeXHr1rj4dyJXutyHUeNcx5LzrEyHNcRwnpBVVVnADMA8vLyGnRWn5uZQm5mbkjj\nMsZEPq9HSE6ItW6VAcG0R+wGOtV43THwtVq3EZEYIAnnwqoxxhgXBFPcvwS6iUhnEYkFRgMFZ21T\nAIwLLN8EfNIY7e3GGGOCU2+zTKAN/X7gI5yukLNUdb2IPAOsUNUC4C/A30RkK/AtzgHAGGOMS4Jq\nc1fVBcCCs7726xrLFcDNoQ3NGGNMQ1kfQGOMiUJW3I0xJgpZcTfGmChkxd0YY6KQuNVjUURKgB0N\n/PY0zrr7tRmwnJsHy7l5uJCcM1W1TX0buVbcL4SIrFDVPLfjCCfLuXmwnJuHcORszTLGGBOFrLgb\nY0wUitTiPsPtAFxgOTcPlnPz0Og5R2SbuzHGmHOL1DN3Y4wx52DF3RhjolCTLu4iki8im0Rkq4g8\nXsv6OBGZG1i/XESywh9laAWR8yMiskFE1orIxyKS6UacoVRfzjW2u1FEVEQivttcMDmLyC2Bfb1e\nRP4e7hhDLYj3doaIfCoiqwLv7+FuxBkqIjJLRA6IyLo61ouIvBz4e6wVkX4hDUBVm+QDZ3jhbUAX\nIBZYA/Q8a5t7gWmB5dHAXLfjDkPO1wIJgeV7mkPOge0SgSXAMiDP7bjDsJ+7AauAlMDrtm7HHYac\nZwD3BJZ7AtvdjvsCcx4E9APW1bF+OPABzjSsA4Hlofz9TfnM/czE3Kp6Ejg9MXdNI4HXA8vzgSFS\n14y0kaHenFX1U1U9Hni5DGdmrEgWzH4GeBaYDFSEM7hGEkzOE4CpqnoYQFUPhDnGUAsmZwVOT5Ca\nBOwJY3whp6pLcOa3qMtIYLY6lgHJItIuVL+/KRf32ibm7lDXNqp6Cjg9MXekCibnmsbjHPkjWb05\nBz6udlLV98MZWCMKZj9nA9kislRElolIftiiaxzB5PwUMFZEinHmj3ggPKG55nz/389LWCfINqEj\nImOBPOAat2NpTCLiAV4A7nQ5lHCLwWmaGYzz6WyJiPRW1VJXo2pcY4C/quofReQHOLO75aiq3+3A\nIlFTPnNvjhNzB5MzInId8CQwQlUrwxRbY6kv50QgB1gkIttx2iYLIvyiajD7uRgoUNUqVf0G2IxT\n7CNVMDmPB+YBqOoXQDzOAFvRKqj/94ZqysW9OU7MXW/OItIXmI5T2CO9HRbqyVlVy1Q1TVWzVDUL\n5zrDCFVd4U64IRHMe/tdnLN2RCQNp5mmKJxBhlgwOe8EhgCISA+c4l4S1ijDqwC4I9BrZiBQpqp7\nQ/bT3b6iXM/V5uE4ZyzbgCcDX3sG558bnJ3/JrAV+C/Qxe2Yw5DzQmA/sDrwKHA75sbO+axtFxHh\nvWWC3M+C0xy1ASgERrsdcxhy7gksxelJsxr4kdsxX2C+/wD2AlU4n8TGAxOBiTX28dTA36Mw1O9r\nG37AGGOiUFNuljHGGNNAVtyNMSYKWXE3xpgoZMXdGGOikBV3Y4yJQlbcjTEmCllxN8aYKPQ/dUMC\nSMdhXxwAAAAASUVORK5CYII=\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[]}}]},{"cell_type":"markdown","metadata":{"id":"_OSBnPk_nOzC","colab_type":"text"},"source":["Please note: in the Bayesian framework we treat the quantities we want to estimate as random variables. \n","\n","We need to define our prior beliefs about them. In the example, the prior was a Beta distribution.\n","\n","After seeing the data we update our belief about the world. In the example, this is vary easy - we keep running counts of the number of failures and successes observed. We update them seeing the data. The prior conveniently can be treated as *pseudo-counts*. \n","\n","To summarize the distribution over the parameter, we typically take its mode (the most likely value), calling the approach MAP (Maximum a Posteriori)."]},{"cell_type":"code","metadata":{"id":"xvDuX8GznOzE","colab_type":"code","colab":{}},"source":[""],"execution_count":0,"outputs":[]}]}
\ No newline at end of file
diff --git a/lectures/02_entropy.ipynb b/lectures/02_entropy.ipynb
new file mode 100644
index 0000000..1e5b28f
--- /dev/null
+++ b/lectures/02_entropy.ipynb
@@ -0,0 +1 @@
+{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"02_entropy.ipynb","provenance":[],"collapsed_sections":[]},"kernelspec":{"name":"python3","display_name":"Python 3"}},"cells":[{"cell_type":"markdown","metadata":{"id":"j2l0x3S2XOXs","colab_type":"text"},"source":["# An intuitive introduction to the Entropy\n","\n","Let $X$ be a discrete random variable (RV) taking values from set $\\mathcal{X}$ with probability mass function $P(X)$.\n","\n","*Definition* the entropy $H(X)$ of the discrete random variable $X$ is\n","\\begin{equation}\n","H(X) = \\sum_{x\\in\\mathcal{X}}P(X)\\log \\frac{1}{P(X)} = -\\sum_{x\\in\\mathcal{X}}P(X)\\log P(X).\n","\\end{equation}\n","\n","How to make sense out of this definition? We'll, rather informally, argue below that the entrpy of an RV provides a lower bound on the amount of information provided by the RV, which we'll dfine as the average number of bits required to transmit the value the RV has taken.\n","\n","As a motivating example consider asking your friend for advice. The probabilities of his answers are given in the table below:\n","\n","| $x$      | $P(x)$ |\n","|----------|--------|\n","| OK       | $1/2$  |\n","| Average  | $1/4$  |\n","| Bad      | $1/8$  |\n","| Terrible | $1/8$  |\n","\n","To transmit the answer of your friend you must introduce an *encoding*, e.g.:\n","\n","| $x$      | $P(x)$ | Code 1 |\n","|----------|--------|--------|\n","| OK       | $1/2$  | 00     |\n","| Average  | $1/4$  | 01     |\n","| Bad      | $1/8$  | 10     |\n","| Terrible | $1/8$  | 11     |\n","\n","Under this encoding, we spend 2 bits per answer.\n","\n","However, we could also consider a variable length code, that uses shorter codewords for more frequent symbols:\n","\n","| $x$      | $P(x)$ | Code 2 |\n","|----------|--------|--------|\n","| OK       | $1/2$  | 0      |\n","| Average  | $1/4$  | 10     |\n","| Bad      | $1/8$  | 110    |\n","| Terrible | $1/8$  | 111    |\n","\n","Under this encoding the average number of bits to encode an answer is:\n","\\begin{equation}\n","\\mathbb{E}[L] = \\frac{1}{2} \\cdot 1 + \\frac{1}{4} \\cdot 2 + \\frac{1}{8} \\cdot 3 + \\frac{1}{8} \\cdot 3 = \\frac{7}{8}\n","\\end{equation}\n","\n","Thus, the new code is more efficient. Is it the best we can do?"]},{"cell_type":"markdown","metadata":{"id":"uOpiKOWvn92s","colab_type":"text"},"source":["### The code space\n","\n","We'll now try formalize the coding task, i.e. the assignment of code lengths to possible values of the RV. \n","\n","Let's first observe an important property of our code: in a variable length coding, no codeword can be the prefix of another one. Otherwise, decoding is not deterministic. Therefore, whenever a value is assigned a symbol of length $L$, $1/2^L$ of the code space is reserved and not available to other codes.\n","\n","This can be visualised as a code space. Below, we indicate the codes assigned to symbols in the example and  grey-out codes that are not available because the shorter codes are used:\n","\n","<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAMsAAACjCAYAAAAkXOoEAAAP30lEQVR4Xu2dz2sd1xXHv0q8URGtINmpJMUEVOjCISkmkGBENv0Tuoi6KsiQGFGtWihKI5uuslARNsE2WtStTCNMF9lVi2Ar1sLGXjRkk0WEswgkJTH2ImgRJypHmvtyNe89vZk536s3T/f7wEjWu/fMnXPPZ84599eMQR9pQBqopIExALuVSg6p0Nzc3JCurMsehQZu3LiBhw8ffg7gBwDfO38yZJTbYDLfMV3swbK7215ezp49i5dffvko+q3RNe7fv59F+z799FPcunWro6MXXngBZ86cwYkTJxrpLVQymdevX/8FgKcBPHXIz8O+K9etU7bKdZcEi6ub9yvnAMujR49w8+ZNzMzMYHJyEk+ePMHm5iYmJiZw+vRplxZNf1euXLGHdus/8izOLsoFlk8++QSvvPJKx5OYp3n8+LFgcdoPtbrCMJ86GTCbZxEsyll8lphRGCZYBEtrYEmVQMuzuLu4I0A5i1OXLGNMmUB7RxMVhu0biWBpCSwpwxwGLClh1miY0whD9RwS/NRPbi8s1hcpw0TBIlgqayA1LJUbMqSCgoWkeHkW/6Qfw7OE7rx79+7er97JyCBPk5IkUEyMYGkHLHEY9uKLLwoWoo3TROUCS8oEmulZWDP38iw0RH4UlAMsqRNowcIxTA0dO/XImGdxNuHQ6uz2ybNksEQ/5dAn88nNBkew8DSahWdJvcRcsDQ3SI2GNdddV01GzpJ6HkOwNO9wwdJcd4KFqDsTpTCMp9BswrCUa6943ZGnJM3gk/pdYZhPkSzPknKARLD4+rhTW7D4FMmAJfUAiWDx9TEdllGZISeprSOGBUvKMFawkHqd4VlGaYacpDbBwlakNn/5Ncp4cvtb0V8Co32ph97lWUgWwPIspOZ0iWEYYyw0xRJ47zyQYNnvoSyGjlOBwpzHSLkEXrBwLECwOPXI9iwpFioyYEk5QKIwzGmEoXpuYVgbYUk9QCJYBEsjDaSApVFDjrCSYCEpW57Fp0h2mJhiAEKw+PqYOilJakpPMWxjTOFZvDlLOQzTHvyUFuWQLc/iUJ5WHfuUV6qt0TCnOnPxLEFNKTyfwjCnEeY6GkZSW0fMKMCcHSxra2uYnZ09kGssLy9jfHzc1f+5hWEuZfWofAxhMYNatiPlotu9DGABwA6A6eIdkOcAfBOVeRXAbQBvA7jQRM+UMMxc8/z8PFZWVjA9PY2dnR0sLCxgamoKi4uLTdqVbYLvUtbxhyWAYnca4Cj/7bkesBhA7wN4C8BWUx3TYFldXcXS0lLHk5in2d7eFixNe4ZU75h5FnvyTkWgBC0FYL4AsF6ChQKKXUiwOI2SbYzO5nRVZ7dvyAm+wbINYK2Hnt4AcLIEy7MMjxKuNRKwsA1I8tqlgRoJfh1YLgK4BOCUJ0+JNTUSsDAm1VKZB/vJzW4nq30t2YNfFRbLTwyS1wB8zfIugsVpnSxjdDajb3VG+1q0B79qzmKwvBuFazYSZnUtVItHyGqpnQaLRsNq6f3ICrNgacke/KajYabvfqBV7gsKLHY1zbNU1vmRFjxmsJjums6zPFN4Ghs6Ht48S8re16SkT7vHEBafQhy1aZ7F0YZDq+YCS8oE2jtAoj34+yYqWJyUs57cKbftChZnJxfVBYtTjyxYUibQDFhSwlxjnsXZW77qgsWnP8op9anDHC8spqKUYaJgcRphqJ5DzpIaFlJXJBMjWEiqFSzteLV36E7twc/gnZIkdrvE5JCzlMMw7cFPZU1Oubl4lpQJNCNnCd045FXHTmvyVVeC79MfJcFPnUALFmcna+iYo0BGGMZpSW8p7PbJs2SUs6RIUJlPbjY4goWn0WzCsDafUs/rzm5JgoWn3WxgSZmgyrM0N0iDWfMszfV3oCZ7NCxFzC1Ymne2YGmuu66aowAL8XazFCXPQur2UYAlB8+itWEZLtFXGFb/KdaiPfj1G0+soQTfqUz2aJOzOV3VGe1LvdBTYRip1xWG+RQpWHz6i2tn51l4qtuXxDBGdptieYz2ybPsa1SwOC2VYYxxE9q4wkCwCBYnJlzP0uYVBoJFsLQKltCYNo7WaTRMsAiWGhrQPItylhrm0rsoO2dJ4VncN5lYgIaOSQpmDx2TmtURMwqwMFcYpBiAECwkqxQsPkWyYE45ACFYfH3cqS1YfIpkwZJyAEKw+PpYsJz2HYEUFChYSIaoSUm/ItnG6G/RQQns9qUYgJBnIfV6bmEYSW0jNQAhWEi9Llh8ipRn8ekvrq21YU5dso3R2Zyu6uz2KQzL6Cikthtj29uXPSzsDspJ3tzcXE63S7/XBw8eYGNj43MAPwD4vvgZ/x7+NuhnkzpVr7cEYLf1YdjY2Bh25fkaGykrDEu1Nmxzc9Ne3vs8gKcAPD3gJ6tMuE4VeVb2vHWAYGlshvsVcxiA0Krj/b4WLIJloAa0n0WwDDSSKgVy8Swp33mpeZYqllahjHKWCko6pAgjZ5FnkWfxWWFRm+VZUiXQgoXSzcpZGGpkwJI6gfbuZ5FnkWdhsEIZDUttjAxYUr7GTzkLxRSBHHKWtsNiXZkyTBQsgqWyBlLDUrkhQyooWEiKl2fxbQJjJPhxV2oPfouXkwiWdsCiPfiawXf7vxxGw2IlZb/quM0LFXPwLKkTaO9omGDR0LHbq5gAhmehNKSPEHbOIs+SQc6ytraG2dnZjkmZkS8vL2N8fNxlq4LFpb69V3ZcuXLlJwCWC0kLAHZ8UtPUzmLVsT0N5+fnsbKygunpaezs7GBhYQFTU1NYXFx0aVawuNQXYHkNwO8A/AzAOzat45OapnY2sKyurmJpaanjSczTbG9vC5aadpUoDHsbwDaAk8XPNQDPALhYgueNosyHAG5HTb8MIHgke/p9BeBNAG8B+BrA+wBOFeUNzK3i9+noO2uDfS4AeLWXfMEiz1ILFzYsW1tbuHbt2r8KKJ4tPExs+AaRwWPxsoVqtwD8sQDBjD78/aOinMEyVcBj92ZbglcLb2UQ2PcGnX1M7j8A/LuQ/QWA9QIgA+2AfMFCgKWWtanwAQ18+eWX+OCDD4JnsNwl9iZm3BaeGTzPFUD9DcAfAJwD8E0hzAAIUMW/l7Vt3uqvAP5cfBFfK3gtg8VCwS75goUAC3Nols0SazQs1dqwq1ev4t69e+XbtpGYcij26yIEC0/+EFaFuqGOwWJhWgi17Pvgfc4C+E/kWcpgvh55li75gkWwDOQz1RaCb7/9FpcuXcJnn332yyipj72JjYoFT2H5jEFgOUj5yR/fQwxLyEns+98Wdft5FrtugKWn/Gxg0WjYQCb6Fki10PPOnTv4+OOPzbOYHYZPObE3I7bk+2ERGlleYwl7yCnK5cuwrACYL2C0UMvCunLOYl4sDIuWc5aO/Cxg2cvkNM/SmJYUsHz33XdYX1/H5OSk5SwxLNbOOO8wYzVjtrDKRqrsUx6tCiFYqBuHYQbGP4t6Foa9BMASeZNVdTRsT342sDS2lAEVc5hnSQFLUGsxKVmGJVV3leXGXuiwgYG9eoLF2S1sWFIsgfcOQBxjWGLPEs/V9LQKwdISWFIugRcszk4uqgsWpx7ZnoU96ccYOk41GmaqH3IYVqv3BUstdXUXzgEWu+tU8yyCxWmAcfVc9rOEe07hWYjdkUSU9uCT1CpYfIpkhGFxC1IMQAgWXx93agsWnyJZsKQcgBAsvj4WLHq1N8mCeGKU4Dt1mUuCnzKnkmdxGmGonlsYRlJbRwwrDBMsCWbwL1zYX77j3a4rWDjYCBaOHk0KLQyLFyqeP39esPD6yCVJsLjUd6AyDZYglbW3XZ6F08mChaNHqmcRLLxOYUoSLDxtjoRn4d0uX9Lc3BxfaEYSHzx4gI2Njc8BxO+xb/JO+yZ1vi+uO6iuHXqxOxKw6HjZ5vSwhrZTrQ3b3Ny0jXnPA6jyTnpWGXu3vcmqIs/KnlcY1twGOzVzGNrWquP97pZncQKTCyx6tbdgcaKi1/iddi7L0RJ9twn+KKDtT+62t4+RsxzjbcW1LJUehtW6eoXCbTdGVvvafPpMC2CxwyT2kuzo89/iLLAmh4jHRyFVsMJEOUvlK1csyDLGiperXYzRvraf8t8SWMLxrKGP4jPAwjGuVftPsFTVFLMcC5Y2n/LfgtGwXscU2cks8cmR8Ukt1sX9Tsu341vt3/+K88gqm4PCsMqq6l0wB1jszlPNs1RM8HvBEk6QtJW7dpZxldPy7bC+ANW7gsVp/HWrjwIsde/pqMtX2M9SN2c57LR8u72BB+r10oE8i9MyRgEW77lhsYqGtAe/XxgWn2NszaxyWr6VU87itPtG1XOBZch78HvBEsCwlxGFt3tZHw46LV+wNLJ0QiUWLKNyyn+Ko5oqhmHl0bCQe4RX4VU9LT8cNG6QWQ5T+aMwrLKq0iX4JrnN8yzxnQ8RlvI8izUrPj2/6mn5Nj+zAcB+Chan/deqzvAstS5YszBjBr8FsNS86zTFs/MsOiPAZ0hD8iy+RpNqZwOLzgjgWIxg2d3laDKBFHaYozMCfJ0kWARLYwtiw9y4IX0qKmfhaTSbMCyoLIVn4XUHX1LbzwjIeg8+u7vZT+4UsLT5jACWZ0m1NizrPfiCha0BnzwGLC1YdexTAqm2wjCnItmez9mcruosWLQHP8EefHZns41RYVj9HmrB5q/6jU5QIzvPwtYhG2Z2++RZeBoVLE5dsmFhrzBoOyxbW1u4du2a7Vz8TY+usL/bmq+624aDqFcBvA7ANojZyuUPAWyVrmMLMn8P4C8Adg4zB8HSElhSrTBoOyylnZJmuOU9Ks4e6lQXLCxNNpXD9izsnIoFy82bNzEzM4PJyUk8efIENuQ7MTEB8rlhvWAp762PvU3wCtZ9PwXwXrSxy/6/DsC8S/AsXwF4s+hr2/diJ8PEnuUlALcjW7gMYCF4HHmWppQU9XKAxW411TzLAM/Sa+9JvPc+3tNi4VX5/+UwbKowfoPC5FiI92wRhtly/b8DsP0xJitsLvsoLOUXLBnA4rzF5NWjzV9lz2KwXCxOcQnng5mBnyy8Ra8TXuITX8qwhA1kcb0Ylj8BOBflSAd2aAoWpynk4lmCmhLvwe8Fiz3xy8n/2xEscXJeTtb7Jfj9YDHPcqpkEp0NZoJFsFTSwBHtwa/iWeL2luEYBEtdz3JAN4Klkqn0L5SbZ0m8RL+fZwn75eNDKkKOUsez1MlZukJAwSJYamngiGGxtpVHw0IIFr6rA8uvitNf4nOSDxsNi/f489/PUkvzFQqzn9wVLlmrSNvbxxg6jhWSGJZauj/qwvIsTo0LFp8CKx7f6rsIqbZgcSpSsPgUKFh8+js4AjE2hjZvrhIsvs4WLD79CRai/pSz8JS5F4bxxElS2zSgPfgY9I57+35QGXudxa7Boo80MEwN/LzGO+mrvLe+Shl7t72Vq1p27+jY/wNE+8URRoxnFgAAAABJRU5ErkJggg==\" />\n","\n","We can observe, that the length 1 code for \"OK\" uses $1/2$ of the available codes, the langth 2 code for \"Average\" uses $1/4$ and the two length 3 codes for \"Bad\" and \"Terrible\" each use $1/8$ of the code space.\n","\n","In general, a code of length $L$ uses $1/2^L$ of the code space. Equivalently, assignign a fraction $f$ of the code space to a symbol makes it use a symbol of length $L=\\log_2(1/f)$."]},{"cell_type":"markdown","metadata":{"id":"T_jOnTuVo0jj","colab_type":"text"},"source":["Assuming that we assign a fraction of a bit to a symbol, our optmal coding problem can be formulated as partitioning the code space into four regions (one for each value of the RV) such that the average length of the code is minised. \n","\n","Formally, let $p_1, p_2, p_3, p_4$ be the proibabilities asigned to the 4 symbols and let $f_1, f_2, f_3, f_4$ be the coding space fractoins assigned to them.\n","\n","We want to:\n","\\begin{align}\n","\\text{minimize } &p_1 \\log_2 \\frac{1}{f_1} +  p_2 \\log_2 \\frac{1}{f_2} + p_3 \\log_2 \\frac{1}{f_3} + p_4 \\log_2 \\frac{1}{f_4} \\\\\n","\\text{subject to: } & f_1 + f_2 + f_3 + f_4 = 1\n","\\end{align}\n","\n","For simplicity, we will solve this problem for the case of only two symbols:\n","\\begin{align}\n","\\text{minimize } &p_1 \\log_2 \\frac{1}{f_1} +  p_2 \\log_2 \\frac{1}{f_2} \\\\\n","\\text{subject to: } & f_1 + f_2 = 1\n","\\end{align}\n","\n","Notice first that $p_2 = 1-p_1$ and likewise $f_2 = 1-f_1$. Then our minimization objective becomes\n","\\begin{equation}\n","\\text{minimize } C = p_1 \\log_2 \\frac{1}{f_1} +  (1-p_1) \\log_2 \\frac{1}{1-f_1} \n","\\end{equation}\n","\n","To get the minimum over $f_1$ we compute the derivative of the expression with respect to $f_1$ and set it to zero:\n","\\begin{equation}\n","\\frac{\\partial C}{\\partial f_1} = \\frac{p_1}{\\log 2}\\frac{-1}{f_1} + \\frac{1 - p_1}{\\log 2}\\frac{1}{1 - f_1}\n","\\end{equation}\n","\n","Multiplying both sides by $\\log 2 f_1 (1- f_1)$ we obtain:\n","\\begin{align}\n","p_1(1-f_1) &= (1-p_1)f_1 \\\\\n","p_1 - p_1f_1 & =f_1 - p_1f_1 \\\\\n","f_1 &= p_1\n","\\end{align}\n","\n","Thus the optimal fraction of code space allocated to symbol 1 is $p_1$, the probability assigned to this symbol and the optimal code length is $\\log_2(\\frac{1}{p_1})$!\n","\n","We now see, that the entropy \n","\\begin{equation}\n","H(X) = \\sum_{x\\in\\mathcal{X}}P(X)\\log \\frac{1}{P(X)}\n","\\end{equation}\n","is simply the average code length!\n","\n","### A note about logarithm basis\n","\n","It is custommary to copute the entropy using natural logarithms, which gives its value in \"nats\". IF $\\log_2$ were used, the entrpy has units of bits and corresponds and lowerbounds the average amount of bits needed to transmit a value of the RV."]},{"cell_type":"markdown","metadata":{"id":"KOBFV5X6aeHg","colab_type":"text"},"source":["# Further Reading\n","1. Chris Olah \"Visual Information theory\": https://colah.github.io/posts/2015-09-Visual-Information/\n","2. JA Thomas ad TM Cover, \"Elements of Information Theory\", chapter 2"]}]}
\ No newline at end of file