{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Think Bayes\n", "\n", "This notebook presents example code and exercise solutions for Think Bayes.\n", "\n", "Copyright 2018 Allen B. Downey\n", "\n", "MIT License: https://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Configure Jupyter so figures appear in the notebook\n", "%matplotlib inline\n", "\n", "# Configure Jupyter to display the assigned value after an assignment\n", "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n", "\n", "# import classes from thinkbayes2\n", "from thinkbayes2 import Suite, Joint\n", "import thinkplot\n", "\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Lincoln index problem\n", "\n", "A few years ago my occasional correspondent John D. Cook wrote an excellent\n", "blog post about the Lincoln index, which is a way to estimate the\n", "number of errors in a document (or program) by comparing results from\n", "two independent testers. \n", "\n", "http://www.johndcook.com/blog/2010/07/13/lincoln-index/\n", "\n", "Here's his presentation of the problem:\n", "\n", ">\"Suppose you have a tester who finds 20 bugs in your program. You\n", "want to estimate how many bugs are really in the program. You know\n", "there are at least 20 bugs, and if you have supreme confidence in your\n", "tester, you may suppose there are around 20 bugs. But maybe your\n", "tester isn't very good. Maybe there are hundreds of bugs. How can you\n", "have any idea how many bugs there are? There's no way to know with one\n", "tester. But if you have two testers, you can get a good idea, even if\n", "you don't know how skilled the testers are.\"\n", "\n", "Then he presents the Lincoln index, an estimator \"described by\n", "Frederick Charles Lincoln in 1930,\" where Wikpedia's use of\n", "\"described\" is a hint that the index is another example of Stigler's\n", "law of eponymy.\n", "\n", ">\"Suppose two testers independently search for bugs. Let `k1` be the\n", "number of errors the first tester finds and `k2` the number of errors\n", "the second tester finds. Let `c` be the number of errors both testers\n", "find. The Lincoln Index estimates the total number of errors as `k1 * k2 / c`\"\n", "\n", "I changed his notation to be consistent with mine.\n", "\n", "So if the first tester finds 20 bugs, the second finds 15, and they\n", "find 3 in common, we estimate that there are about 100 bugs.\n", "\n", "Of course, whenever I see something like this, the idea that pops into\n", "my head is that there must be a (better) Bayesian solution! And there\n", "is." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def choose(n, k, d={}):\n", " \"\"\"The binomial coefficient \"n choose k\".\n", "\n", " Args:\n", " n: number of trials\n", " k: number of successes\n", " d: map from (n,k) tuples to cached results\n", "\n", " Returns:\n", " int\n", " \"\"\"\n", " if k == 0:\n", " return 1\n", " if n == 0:\n", " return 0\n", "\n", " try:\n", " return d[n, k]\n", " except KeyError:\n", " res = choose(n-1, k) + choose(n-1, k-1)\n", " d[n, k] = res\n", " return res" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def binom(k, n, p):\n", " \"\"\"Computes the rest of the binomial PMF.\n", "\n", " k: number of hits\n", " n: number of attempts\n", " p: probability of a hit\n", " \"\"\"\n", " return p**k * (1-p)**(n-k)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "class Lincoln(Suite, Joint):\n", " \"\"\"Represents hypotheses about the number of errors.\"\"\"\n", "\n", " def Likelihood(self, data, hypo):\n", " \"\"\"Computes the likelihood of the data under the hypothesis.\n", "\n", " hypo: n, p1, p2\n", " data: k1, k2, c\n", " \"\"\"\n", " n, p1, p2 = hypo\n", " k1, k2, c = data\n", "\n", " part1 = choose(n, k1) * binom(k1, n, p1)\n", " part2 = choose(k1, c) * choose(n-k1, k2-c) * binom(k2, n, p2)\n", " return part1 * part2" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "from itertools import product\n", "\n", "ns = range(32, 350)\n", "ps = np.linspace(0, 1, 31)\n", "hypos = product(ns, ps, ps)\n", "\n", "suite = Lincoln(hypos);" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8.754332814652824e-06" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = 20, 15, 3\n", "suite.Update(data)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "post mean n 105.4574266270967\n", "MAP n 72\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4FNeV8P/vUWtf0cYqQGK1WWwMGLyAd8c4dkySwW9wNmdez+s3GTvrZGbsyTMejzOeX5zJMs4vzuLESZxljAl2EmJjk3jFeAHEYkBmE4uRWAUSkkC7dN4/qtRqyVpRV1e3dD7Po0dVt29Vn6KB0/fWrXtFVTHGGGOiTZzfARhjjDHdsQRljDEmKlmCMsYYE5UsQRljjIlKlqCMMcZEJUtQxhhjopIlKGOMMVHJEpQxxpioZAnKGGNMVIr3O4BIyMvL08LCQr/DMMYYA2zevPmUqub3VW9YJKjCwkKKi4v9DsMYYwwgIu/3p5518RljjIlKlqCMMcZEJUtQxhhjotKwuAdljDFDTXNzM+Xl5TQ0NPgdSo+Sk5MpKCggISHhvI63BGWMMTGovLycjIwMCgsLERG/w/kAVeX06dOUl5dTVFR0XuewLj5jjIlBDQ0N5ObmRmVyAhARcnNzB9XCswQ1BLS2ttHS0up3GMaYCIvW5NRusPFZF18M23voBE89v4mS/Udpa21j3Khsblo0g5uunEkgYN89jDGxzRJUjFqzbge/eOZNNKSs/EQVTzzzJuu37Oef7voQIzJSfYvPGGMGy75mx6DXN+3liS7JKdSeg8d54AerOVNbF9G4jDEmnKwFFWOOn6rhpyvfCO5PnTiSL336OrIzU3nu9R2sfKGYNlWOnDzDt5/4Cw/d+xHi4wM+RmyM8drffPknnp37mUc/3+Nrhw4d4uabb2bRokW89dZbjBs3jj/96U+kpKSE5b2tBRVjfvHMmzQ2NQMwNj+LB75wK2NHjiAlOZHbb5rHV+68gfbbknsOHue3f97gX7DGmCFv37593HPPPZSUlDBixAieeeaZsJ3bElQMKSk9yub3nDkWBfjyZ64nNSWxU50rL5nMp2+7LLj/3GvbeW//sUiGaYwZRoqKipgzZw4A8+bN49ChQ2E7t3XxxZBVa7cEtxfPn8qUiSO7rbf0uovZsfcI23aXocDPfv8G3/2nZcTF2fcRY4ai3rrhvJaUlBTcDgQC1NfXh+3c9j9WjCg7XsX2veWA03r6xM2X9lhXRPjC8qtJTHC+fxw+VskrG/ZEIkxjjAkbS1AxYu36kuD2pbMLGZ2X2Wv9vOx0PnbDnOD+757bSF19k2fxGWNMuFmCigEtLa28sXlfcP/mxbP6ddzS6y4md0QaADVn6/nTK9s8ic8YMzwVFhayc+fO4P7Xv/51HnzwwbCd3xJUDNi2p5yzdY0A5I5IY/a0cf06LikxgU/dujC4//y6nZyrb/QkRmOMCTdLUDFg/ebS4PbieVMHNL/V4nlTGJufBUB9QxMvvFHSxxHGGBMdPE1QIrJERPaISKmI3NfN60ki8rT7+gYRKXTLc0XkVRE5KyI/DKmfKiLPi8huESkRkW95GX80aGlpZXPJ+8H9RXOnDOj4uLg4Pn7j3OD+c69tDz5HZYyJbao9zScTHQYbn2cJSkQCwGPAzcAM4A4RmdGl2l1AlapOAb4PPOKWNwD/Cny9m1N/R1UvAC4BrhSRm72IP1rsPnicugZncEPuiDQKx+UO+ByL500hLzsdgNpzDfz1rV1hjdEYE3nJycmcPn06apNU+3pQycnJ530OL5+DWgCUquoBABFZASwF3gupsxR40N1eBfxQRERVzwHrRaRTc0FV64BX3e0mEdkCFHh4Db7bXHI4uD1/5vktTBYfH+Cj18/h56vWA7Bm3U4+fNUsey7KmBhWUFBAeXk5FRUVfofSo/YVdc+XlwlqHFAWsl8OLOypjqq2iEg1kAuc6uvkIjIC+AjwaA+v3w3cDTBhwoSBxh41tu7qSFDzZp7/dVy3cDpPPb+Jc/WNnDhdw5ZdZcyfOTEcIRpjfJCQkHDeK9XGCi+/Qnf3Vb9rW7Q/dT54YpF44CngB+0ttA+cRPVxVZ2vqvPz8/P7DDYanamto+x4FQCBQByzpo4973MlJSZw/WUXBPdfWLezl9rGGOM/LxNUOTA+ZL8AONpTHTfpZAGV/Tj348A+Vf3vMMQZtXbu7fjjml44iqTEhEGdb8nimcFvBNt2l3Hk5JlBnc8YY7zkZYLaBEwVkSIRSQSWA6u71FkN3OluLwNe0T7u+InIf+Aksq+EOd6os7P0SHB71tT+PfvUm1G5mcwL6dZ78Q1rRRljopdnCUpVW4B7gbXALmClqpaIyEMicptb7QkgV0RKga8BwaHoInII+B7wOREpF5EZIlIAfANnVOAWEdkmIn/n1TX4rWRfRwtqMN17oT589ezg9uub9tHU3BKW8xpjTLh5Opu5qq4B1nQpeyBkuwG4vYdjC3s47cCHscWgmrP1HK2oBpz7T1N7mLl8oC6aNo6RORmcrKzlXH0jG7cfYtG8gT1bZYwxkWDjjKPU3vdPBrcnFeQFZyYfLBHh2oXTg/svv7M7LOc1xphwswQVpfYcOB7cnl44OqznvnbB9GAzdMfeck5W1ob1/MYYEw6WoKLUnkMngtvTJ40K67nzczK4+AJngKUCr9paUcaYKGQJKgqpKvvLOp4OnzYxvAkK4LqQZ6Je3bAnaqdLMcYMX5agotDRimoaGp0JXbMyUoJrOoXTglmFpKc6SzVXVNWyfe+RPo4wxpjIsgQVhfYf7hggMXl8/nnNv9eXhIQAV82fGty3wRLGmGhjCSoK7T/cMRXhpPHeTdMUOvXRxu0HqW+wJeGNMdHDElQUOlDecf9psocJqnBcHhPG5ADQ3NLKhu0HPXsvY4wZKEtQUUZVOXTkdHB/UkGep++3eF5HN9+64n2evpcxxgyEJagoU1l9LrhAYWpyoicDJEItDplFYvuecqpq6jx9P2OM6S9LUFHm8LGq4Pb4MTmeDJAIlZ+TwYzJYwDnmag3t5R6+n7GGNNflqCizOFjHauNTBiTHZH3tG4+Y0w0sgQVZUIT1PjRORF5z8vnTCIQcP4q7C+r4KitE2WMiQKWoKJMWacWVGQSVEZaMnMv7FhO/o3N1s1njPGfJagooqrBJd4hcgkKYHHIQ7tvbN5nUx8ZY3xnCSqKnDhdG1xAMDM9hayMlIi996WzJpKc5Cwpf6yimv2HK/o4whhjvGUJKor4MUCiXWJCPJddPCm4v26zDZYwxvjLElQUOezD/adQoXPzrd9SSltbW8RjMMaYdpagoojfCWr21LGMyEgFoLq2np37jkY8BmOMaWcJKor4MYIvVFxcHFdc0tHNZ6P5jDF+sgQVJVSVYxXVwf2C0ZG9B9Uu9KHdd949QHNzqy9xGGOMJagocarqLM0tTjLISEsmLSXJlzimThzJyJwMAOoamti6u8yXOIwxxhJUlAhtPY3Jz/ItDhFh0dyOCWTX29x8xhifWIKKEtGSoAAWhcxwvmnHoeDy88YYE0meJigRWSIie0SkVETu6+b1JBF52n19g4gUuuW5IvKqiJwVkR92OWaeiOxwj/mBeD3dd4REU4KaMCaHglHOPbCm5haKd77vazzGmOHJswQlIgHgMeBmYAZwh4jM6FLtLqBKVacA3wceccsbgH8Fvt7NqX8M3A1MdX+WhD/6yDt+qia47XeCEpFOrSjr5jPG+MHLFtQCoFRVD6hqE7ACWNqlzlLgSXd7FXC9iIiqnlPV9TiJKkhExgCZqvq2OpPF/Rr4qIfXEDGdWlB5/iYooNN9qC27DnO2rtHHaIwxw5GXCWocEDoErNwt67aOqrYA1UBuH+cs7+OcAIjI3SJSLCLFFRXRPa9cW1sbx05FTxdfewyTx+cD0NraxobtB3yOyBgz3HiZoLq7N9R1iuz+1Dmv+qr6uKrOV9X5+fn5vZzSf6fOnKO11ZlWKDM9hdSURJ8jcoQ+E2UP7RpjIs3LBFUOjA/ZLwC6zp0TrCMi8UAWUEnPyt3z9HbOmBNNAyRCXXHJpOA3gp17j1BVU+drPMaY4cXLBLUJmCoiRSKSCCwHVnepsxq4091eBryivSxEpKrHgFoRucwdvfdZ4E/hDz2yjp2MzgSVOyKdCyePAZxm6ltb9/sbkDFmWPEsQbn3lO4F1gK7gJWqWiIiD4nIbW61J4BcESkFvgYEh6KLyCHge8DnRKQ8ZATgF4CfA6XAfuAFr64hUqK1BQWdu/lsNJ8xJpLivTy5qq4B1nQpeyBkuwG4vYdjC3soLwZmhS9K/0Vzgrp8ziR+tmo9bW1t7D10gpOVtcGpkIwxxks2k0QUOFZxJrg9NsoSVEZaMnMu6Ljtt94GSxhjIsQSlM9UlROVtcH90XmZPkbTvdBnot6wlXaNMRFiCcpnVTV1wSHmGWnJpCRHxxDzUAtmF5IQHwCcRRXLjlf5HJExZjiwBOWzk6c7Wk/Rem8nJTmR+bMKg/vrrRVljIkAS1A+O1nZMQdffpQmKIArL5kc3F6/pZRengYwxpiwsATls5OVZ4Pb0dqCApg3c0Kw+/H4qRr2H47u6aOMMbHPEpTPTp7uaEGNzI3eBJWYEM/Ci4qC+/ZMlDHGa5agfFYR0oKK5i4+6Dya782t+62bzxjjKUtQPgu9BxXNXXwAF00bR2Z6CgCV1ed4b/8xnyMyxgxllqB8pKpUVMXGPSiAQCCOyy+eFNy3Z6KMMV6yBOWjyupznZ6BSk5K8Dmivi0OWWn37W0HaGlp9TEaY8xQZgnKRxUxMoIv1AWTRpM7Ig2As3WNbN97xOeIjDFDlSUoH8XKM1ChRIQrL7Gpj4wx3rME5aNYeQaqq9Buvg3bD9HU3OJjNMaYocoSlI9i5RmorooK8oKzrjc2NVNc8r7PERljhiJLUD6KpWegQokIV4a0ot60JTiMMR6wBOWjWHoGqqvQh3aL3zvMufpGH6MxxgxFlqB8EmvPQHVVMCqbwnF5ALS0tLJx+yF/AzLGDDmWoHxSfbY++AxUanJiTDwD1dWiuR0znL+51br5jDHhZQnKJ5VnzgW3258rijWh3Xzv7i6nurbex2iMMUONJSifnK7uSFA5WbGZoPJzMpheNBqANlXe3nbA54iMMUOJJSifhLagcmK0BQWdn4myJTiMMeFkCconldWhXXzpPkYyOFfMmUycCAC7DhzjVMjAD2OMGQxLUD4J7eLLjdEuPoCsjBRmTxsX3H+9eK+P0RhjhhJPE5SILBGRPSJSKiL3dfN6kog87b6+QUQKQ1673y3fIyI3hZR/VURKRGSniDwlIsleXoNXhkoXH8A1C6YFt1/bsMcWMjTGhIVnCUpEAsBjwM3ADOAOEZnRpdpdQJWqTgG+DzziHjsDWA7MBJYAPxKRgIiMA74EzFfVWUDArRdzKodICwpg4UVFwWHyRyuq2XvohM8RGWOGAi9bUAuAUlU9oKpNwApgaZc6S4En3e1VwPUiIm75ClVtVNWDQKl7PoB4IEVE4oFU4KiH1+CZ02difxRfu6TEBK68pOOZqFc37vExGmPMUOFlghoHlIXsl7tl3dZR1RagGsjt6VhVPQJ8BzgMHAOqVfUv3b25iNwtIsUiUlxRURGGywmfxqZm6hqaAGeV2sz0mOyl7OS6hRcEt9dv2W8znBtjBs3LBCXdlHW9OdFTnW7LRSQbp3VVBIwF0kTk0929uao+rqrzVXV+fn7+AML2XqfWU2YaIt1dbmyZXjSK0XmZANQ3NNnUR8aYQfMyQZUD40P2C/hgd1ywjttllwVU9nLsDcBBVa1Q1WbgWeAKT6L3UOj9p1gfINFORLhmwfTgvnXzGWMGy8sEtQmYKiJFIpKIM5hhdZc6q4E73e1lwCvqDAFbDSx3R/kVAVOBjThde5eJSKp7r+p6YJeH1+CJyiEwi0R3rrl0WrDp++7uMk6fsWeijDHnz7ME5d5TuhdYi5NEVqpqiYg8JCK3udWeAHJFpBT4GnCfe2wJsBJ4D3gRuEdVW1V1A85gii3ADjf+x726Bq+cqho6I/hC5edkMHPqWMDpp319ky0Hb4w5f/FenlxV1wBrupQ9ELLdANzew7EPAw93U/5vwL+FN9LIqqoZel187a5beAE79zk9ua9u2M3HbpgzJO6xGWMiz2aS8EGnmcyHUAsKnGeikhLtmShjzOBZgvJB6DRH2VmpPkYSfslJnZ+Jeunt3T5GY4yJZZagfDBUB0m0u/GKC4Pb67eU2nLwxpjzYgkqwtra2qiqrgvux+pihb2ZOnEkE8bkANDU3MIbxbYMhzFm4CxBRdiZ2nra3MlU01OTSEzwdJyKL0SEGy7vaEX99e1dNoGsMWbALEFFWOUQmoOvN1dfOo2E+AAAh46cYv/h6JpuyhgT/SxBRVhlzdDu3muXnprEFSGDJf76dsw9T22M8ZklqAgbLi0ogA9d0bG6yhubS6l3J8g1xpj+sAQVYUNxHr6eTC8aRcGobMCZwf2NzTZYwhjTf5agImyoLPXeHyLSaci5dfMZYwbCElSEDacuPnAGS8S7gyUOlFXYYAljTL9Zgoqw0Bm+h/IgiXYZaclcMWdScP/F9SU+RmOMiSWWoCLs9BCfRaI7N105M7j9xuZ91Jyt9zEaY0ys6DVBicivQrbv7KWq6Yf6hiYaGpsBiI8PkJEW+0u998f0olFMGu+satzc0mrz8xlj+qWvFtTFIdtf9jKQ4aBT6ykzddgsQyEifHjxrOD+i+t30tra5mNExphY0FeCsvlpwqjTAIlhcP8p1JVzJ5OZngLA6TPn2LjjkL8BGWOiXl8JqkBEfiAi/3/IdvAnEgEOJUN9FvPeJCbEc2PI/HwvvLHTx2iMMbGgr5lK/zFku9jLQIaD0C6+vBHpPkbij5sWzeAPL22lTZWS0qO8f/Q0E8fm+h2WMSZK9ZqgVPXJSAUyHAy3Z6C6yh2RzsKLJ/H2tv0ArFm3ky8sv9rnqIwx0arXBCUiq3t7XVVvC284Q9twmuaoJ7dcNSuYoF7ftJdPf2ThsBnNaIwZmL66+C4HyoCngA3A8Bh25pHTZ4bPNEc9uWDSaArH5XHoyCmaW1r561u7+PiNl/gdljEmCvU1SGI08C/ALOBR4EbglKq+rqqvex3cUGMtKGfI+a1Xzw7ur1m3g+bmVh8jMsZEq14TlKq2quqLqnoncBlQCrwmIl+MSHRDSGtrG2dC1oLKyUz1MRp/LZo7hWz3+qtq6li3ea/PERljolGfUx2JSJKIfBz4LXAP8APgWa8DG2rO1NYFHyrLTE8JTqA6HCUkBLj1mouC+396+V1bEt4Y8wF9TXX0JPAWMBf4d1W9VFW/qapH+nNyEVkiIntEpFRE7uvm9SQRedp9fYOIFIa8dr9bvkdEbgopHyEiq0Rkt4jsEpHL+3mtvjo9zEfwdXXjFReSkpwIwJGTZ9i0832fIzLGRJu+WlCfAabhTHP0tojUuD+1IlLT24EiEgAeA24GZgB3iMiMLtXuAqpUdQrwfeAR99gZwHJgJrAE+JF7PnDuhb2oqhfgTMUUE4sM2QCJztJSkrjpyo6/Dn98eZuP0RhjolFf96DiVDUj5CfT/clQ1cw+zr0AKFXVA6raBKwAlnapsxRof9ZqFXC9OBPULQVWqGqjqh7Eufe1QEQygauAJ9z4mlT1zEAu2C+dB0gM3/tPoW65ejaBgPNXcM/B4+w+cNzniIwx0aSvLr5kEfmKiPxQRO4Wkb6GpYcahzNEvV25W9ZtHVVtAaqB3F6OnQRUAL8Uka0i8nMR6bY54sZbLCLFFRX+L5JXNYynOepJTlYaV8+fFty3VpQxJlRfXXxPAvOBHcCHge8O4NzdPTPV9U54T3V6Ko/HuR/2Y1W9BDgHfODeFoCqPq6q81V1fn5+fv+j9kinpd6H6RDz7iy9vmPC/E07D1F2vMrHaIwx0aSvBDVDVT+tqj8FlgGLB3DucmB8yH4BcLSnOm7rLAuo7OXYcqBcVTe45atwElbUC+3iy860BNWuYFQ2l84qDO7/4aWt/gVjjIkqfSWo5vYNtwtuIDYBU0WkSEQScQY9dJ06aTXQvhDiMuAVdcYbrwaWu6P8ioCpwEZVPQ6Uich095jrgfcGGJcvQufhsxZUZ6EzSbxRvI9jFdU+RmOMiRZ9LlgYOnIPuKi/o/jchHYvsBZnpN1KVS0RkYdEpH0OvyeAXBEpBb6G212nqiXASpzk8yJwj6q2TzfwReB3IrIdmAP850AvOtJUldPVHQ/p5g7Dmcx7M61wFBdNKwCgTZVVf9nic0TGmGjQ12zmg3qaVFXXAGu6lD0Qst0A3N7DsQ8DD3dTvg3nvljMqGtoorHJaYwmxAdIS0n0OaLo87+WzGP73nIA1m3ay7IPzWVMfpbPURlj/NTnTBJm8Co7tZ7Shs1S7wNx4eQxzJ7mDPJsU+WZv1orypjhzhJUBAznlXQH4n8t6WgYv75xL8dP9dqLbIwZ4ixBRUCnhQptgESPZkwew6ypY4H2e1GbfY7IGOMnS1AR0OkZKGtB9Sq0FfXahj32XJQxw5glqAg4feZscNu6+Ho3c8pYLp7ujOhTYMXzG/0NyBjjG0tQEWBdfAPzqVsXBrff2X6Q0vdP+hiNMcYvlqAiwLr4BmbyhHwunzM5uP+756wVZcxwZAkqAmwU38DdcculxLnD8bfvLWf7nnKfIzLGRJolKI+1tLRSU1sPODPgZg/jpd4HYtzIEVx32QXB/V+vfsdW3TVmmLEE5bGqmpCl3jOG91LvA3X7TfNIcP+8Dpaf4rWNe32OyBgTSZagPGbde+cvLzudpdd1LMfxu+c20NDY3MsRxpihxBKUx2yAxOB87IZLgt2iVTV1/MEWNTRm2LAE5bHOy2zYLOYDlZyUwCdvWRDc/9PL2zhVdbaXI4wxQ4UlKI916uKzZ6DOyzULplE4Lg+A5pZWfvPnd3yOyBgTCZagPGZdfIMXFxfH337s8uD++s2llJR2XZzZGDPUWILymM0iER6zpo7r9PDuz37/Bi0trb0cYYyJdZagPGaj+MLncx+9nKTEBADKjlfx/LqdPkdkjPGSJSgPqSqnQ1tQWfaQ7mDkZafziZs7Zjt/+oViGzBhzBBmCcpDZ+saaXa7oZISE0hNtqXeB+uWq2YxfnQ2AI1Nzfzy2Td9jsgY4xVLUB6q7DRAItWWeg+D+PgA/+f2xcH9d7YfZMP2gz5GZIzxiiUoD522ARKemDllLNcunB7cf3zlG5yta/QxImOMFyxBeaiqxgZIeOVzH72CERnOPb0ztXU8+ce3fY7IGBNulqA8FNqCsmegwis9NYm/W7YouP/Kht28a0tyGDOkeJqgRGSJiOwRkVIRua+b15NE5Gn39Q0iUhjy2v1u+R4RuanLcQER2Soiz3kZ/2CF3oPKtgQVdpfPmcRlFxUF93/81OvU1Tf5GJExJpw8S1AiEgAeA24GZgB3iMiMLtXuAqpUdQrwfeAR99gZwHJgJrAE+JF7vnZfBnZ5FXu4VJ6pC25bF583/u72xaSlJAFQUVXLL/5go/qMGSq8bEEtAEpV9YCqNgErgKVd6iwFnnS3VwHXizPUbSmwQlUbVfUgUOqeDxEpAG4Bfu5h7GHRaZojGyThiezMVO4OGdX36oY9vPPuAR8jMsaEi5cJahxQFrJf7pZ1W0dVW4BqILePY/8b+CegLfwhh1elJaiIWDRvClfOnRLc/8nT66iqqevlCGNMLPAyQXX30E/XNbt7qtNtuYjcCpxU1c19vrnI3SJSLCLFFRUVfUcbZs3NrdSc7VjqvX3EmfHG3bcvDn4JqD3XwI+ees2WiDcmxnmZoMqB8SH7BUDXKaiDdUQkHsgCKns59krgNhE5hNNleJ2I/La7N1fVx1V1vqrOz8/PH/zVDFBVbcc3+BGZqQQCNmDSS+mpSdz7yWuD+1veO8yfX9vuY0TGmMHy8n/NTcBUESkSkUScQQ+ru9RZDdzpbi8DXlHna+9qYLk7yq8ImApsVNX7VbVAVQvd872iqp/28BrO2+mQOeJsgERkXDS9gNuu7Vgi/jerN7D30AkfIzLGDIZnCcq9p3QvsBZnxN1KVS0RkYdE5Da32hNAroiUAl8D7nOPLQFWAu8BLwL3qGpMra0Q+gxUfratpBspn7p1AVMmjASgra2N7/3qJWrPNfgclTHmfHja76Sqa1R1mqpOVtWH3bIHVHW1u92gqrer6hRVXaCqB0KOfdg9brqqvtDNuV9T1Vu9jH8wTp3paEHlWoKKmPj4AP/wtzcGJ+atqKrlsf+x+1HGxCK7MeKR06EJaoQlqEgamZPBvZ/quB+1aechnvnrVh8jMsacD0tQHgldp8iGmEfewouKuPXqi4L7K57fSHHJ+z5GZIwZKEtQHgm9B5VnLShffOa2hcycMhZwnl3471+/zJGTZ/wNyhjTb5agPNLpHpS1oHwRHx/gHz53I3nuPcD6hiYe+dmLnKu3pTmMiQWWoDzQ3NxKdW3HQ7o2zNw/WRkp/PNdN5EQ70zleOTkGb79xFpaWmJqUKgxw5IlKA9U1nSexdwe0vXXpPH5nR7i3bnvKI/ZTBPGRD37n9MDNkAi+iyaN4U7blkQ3F9XvI8VLxT7GJExpi+WoDxQGbpQoQ2QiBp/c+Ml3HD5hcH9VWs38/I7Ub9qizHDliUoD1SEtKBsBF/0EBHuvn0xl1zYMc3jT1asY+OOQ/4FZYzpkSUoD3R6SDfbuviiSSAQxz987kYKx+UB0KbKd375F7buKuvjSGNMpFmC8sBp6+KLainJiXzj/97M6LxMAFpb23jk5y+yc98RnyMzxoSyBOWB0GegbKLY6JSTlcaD93yE/OwMAJpbWvnPx19k94HjPkdmjGlnCcoDnVtQ1sUXrfJzMvi3e24lO9NZTLKxqZn/+Oka9hy0JGVMNLAEFWZNzS3BlXTjRIL/+ZnoNCY/iwfv/QiZ6SmAM9vEv//oebbvKfc5MmOMJagwC209ZWelEhdnf8TRrmBUNg/ec2swSbW3pGx0nzH+sv89wyz0Id089/6GiX4Tx+byzS/dFpyWqrW1jf96Yi1vFO/zOTJjhi9LUGF2srIRrXRXAAAXsklEQVQmuJ2fYwMkYknBqGwe/spHg6P72lR59Dcv8+dXt9u0SMb4wBJUmJ08XRvcHpWT6WMk5nyMzMngm19ayvjR2YCzTMev/vgWP/v9elpb2/wNzphhxhJUmJ0ITVB51sUXi3Ky0vjml5YyrXBUsGztmyX8fz97gbr6Jh8jM2Z4sQQVZicrOxLUSGtBxayMtGT+/d6PcOXcKcGyrbvK+Majf+z0GRtjvGMJKsxOnu64BzUy11pQsSwxIZ6vfvZ6ln1obrDs8LFK/vG/VrHlvcM+RmbM8GAJKoyamluoqqkDnGegbKLY2Cci3HHLAu795LXBdb3O1jXynz9dw9MvFNvgCWM8ZAkqjELvP+Vlp9tChUPItQun880vdgxDV2Dli8U8/NM1wQezjTHhZf+DhpF17w1t04tG851/XMasqWODZVt3lfG1R37Ptt02G7ox4eZpghKRJSKyR0RKReS+bl5PEpGn3dc3iEhhyGv3u+V7ROQmt2y8iLwqIrtEpEREvuxl/ANlAySGvqyMFB74wq187Po5wbKqmjq++ePn+cWzb9LU3OJjdMYMLZ4lKBEJAI8BNwMzgDtEZEaXancBVao6Bfg+8Ih77AxgOTATWAL8yD1fC/APqnohcBlwTzfn9E2nZ6DyLEENVYFAHJ++7TLuv/vm4PRIAM+/voN/+s4zHCw/5WN0xgwdXragFgClqnpAVZuAFcDSLnWWAk+626uA60VE3PIVqtqoqgeBUmCBqh5T1S0AqloL7ALGeXgNAxLaxTcqx7r4hrr5Myfy/ftuZ96MicGysuNV/NN3nuG3q9+x1pQxg+RlghoHhHbMl/PBZBKso6otQDWQ259j3e7AS4AN3b25iNwtIsUiUlxRUXHeFzEQx0NaUHYPangYkZHK/Xcv4e7bF5MQHwCcKZL+8PI2vvqtlezYa4sgGnO+vExQ0k1Z1zG5PdXp9VgRSQeeAb6iqjXd1EVVH1fV+ao6Pz8/v58hD07nQRLWxTdciAg3LZrJ9/75dmZMHhMsP36qhgcf+zOP/uZlKqvP9XIGY0x3vExQ5cD4kP0C4GhPdUQkHsgCKns7VkQScJLT71T1WU8iPw9n6xqpa3CmwUmIDzAiI6WPI8xQM3bkCB764m18YfnVpCYnBsvXFe/j3v9YwbN/3Upzc6uPERoTW7xMUJuAqSJSJCKJOIMeVnepsxq4091eBryizpOPq4Hl7ii/ImAqsNG9P/UEsEtVv+dh7APW6f5TbiZOqGa4ERFuuPxCHv2XT3DZxZOC5Y1NzfzuuQ185VtP8/a2A/aArzH9EO/ViVW1RUTuBdYCAeAXqloiIg8Bxaq6GifZ/EZESnFaTsvdY0tEZCXwHs7IvXtUtVVEFgGfAXaIyDb3rf5FVdd4dR39dexU5wRlhrecrDT+8X9/iO17yvnFs29SdrwKcLr9vvPLvzB5fD6f+shCLpo2zr7MGNMDGQ7f5ObPn6/FxcWevsfKF4t5+gXnPT5yzUV87mNXePp+Jna0trbx4voSVqzZFOwGbjdzylg+ecsCLpg02qfojIk8EdmsqvP7qmczSYTJ0ZPVwe2xI0f4GImJNoFAHLdcPZvH/vUObrv2YuLd0X4AJaVH+cajf+Q/fvI8JaVHrevPmBCedfENN0dPnglujxtlCcp8UGZ6Cnd+9HJuuXo2q/6ymZff3k2bm5C27ipj664ypk4cycduuIQFswut688Me9aCCgNV5UhIgrIWlOlNXnY6n//E1fzgG8tZNG9Kp2cq9r1/km8/sZYv/+fTvPzOLhqbmn2L0xi/WYIKg6qaOhoanf9IUpMTbYi56Zcx+Vl89bM38Og3lnP9ZRd0mv3+yMkz/Oip1/k/D/yWX/3hLY5VVPdyJmOGJuviC4MjJzq3nqxrxgzEuJEj+Ps7rmH5hy/l+dd38OL6kuAXnnP1jfz5te38+bXtXHLheG5aNJO5F06wpVzMsGAJKgwOH6sMbheMzvYxEhPLcrLS+Mxtl/HxGy/hpbd3s3Z9CSdCnq9rv0+VlZHC1fOncc2C6Uwcm+NjxMZ4yxJUGIQmKPsPwwxWWkoSS6+7mNuuvYgt7x3mxfUlbH3vcHCur+raela/+i6rX32XSePzuebSaVw+Z1JwMUVjhgpLUGHw/tHTwe2JY3N9jMQMJSLCvJkTmTdzIscqqnnp7V28vmkvVTV1wToHyio4UFbBL599kwsnj+GKSyZz2cWTyM5M9TFyY8LDEtQgqSqHj1UF960FZbwwJj+Lz9x2GZ+8ZQHv7innlQ172LjjIK2tbYAzk/J7+4/x3v5jPLFqfTBZzZ85kXxb+sXEKEtQg3TidG1wKHBmegojMuybq/FOIBDH3BkTmDtjAmfrGnlzSylvbdtPyb6jwS7A0GT181XrmTAmh/luS2xa4Uji4myAhYkNlqAG6dCRjtVTrfVkIik9NYmbFs3kpkUzOVNbxzvbDvLWtv28V3q007o2h49VcvhYJc++tJX01CTmzpjAxdMLmD1tHLkj0n2L35i+WIIapNL3Twa3i8bl+RiJGc5GZKSyZPFMliyeSVVNHRvePUhxySG27z0S7AYEZ1mYdcX7WFe8D4Cx+VnMnuYkq1lTx5KRluzXJRjzAZagBmnf4Y4ENbVwpI+RGOPIzuxIVg2NzWzfe4TinYfYXHKYM7V1neoerajmaEU1a98sQYCJ4/KYMXk0F0waw/TCUeRlWwvL+McS1CC0tbVRerhjOfmpEyxBmeiSnJTAgtmFLJhdiKpyoOwUW3eXsXPfEXYfOE5zS8cCiorTZX3oyCnWrNsJQO6INKYXjeYC92fi2JxOk90a4yVLUINQfuJM8In/ERmp9m3TRDURYfKEfCZPyGfZh+bS1NzCnoMn2LH3CNv3lrP/cEVw8tp2p8+c462t+3lr637AGaRRODaXSePzmFSQz6SCPCaOzSUhwZKWCT9LUIOw+8Dx4PbUiSNtiiMTUxIT4pk9bRyzp43jkyygrr6J3QePs+fgcXYfPM6+9ys+MFlta2sb+8sq2F9WAewCIC4ujgljcphUkEfhuFwmjMlhwpgcsmxOSjNIlqAGYfveI8HtGVPG+BiJMYOXmpIYHMIOTjI6fKySXQeOsfvgCfYePEFFVe0Hjmtrawt2DYbKSEtm/Ohsxo92Etb4MdmMH51NRlqyfZkz/WIJ6jypKjv2lgf3L55e4GM0xoRfIBBHUUEeRQV5fPiq2QDUnmvgQPkp9h+u4ED5KQ6WV3D8VE23x9eeawg+jxUqNTmRMflZjM7PYkx+FmPd32Pys2wUoenEEtR5Olh+irN1jQBkZaQwYYw9A2WGvoy0ZC6eXtDpC9m5+kYOlp9if9kpDh+rpOxYJWXHq2hqbun2HHUNTSHdhJ2lpSQxOi+TkbmZ5Genk5edTn5ORnA7PTXJWl/DiCWo87Rhx6Hg9kXTCuwfjRm20lKSmDV1HLOmjguWqSonK2spO17F4aOVlB2v5PCxKo5VVPe6COO5+sYekxdAUmIC+dnp5Oc4CSsnK40RGalkZ6WS7f4ekZFqy5EMEZagzoOqsn7zvuD+wouKfIzGmOgjIozKzWRUbibzZ04MlqsqVTV1HKuo5vipao6drOaY+yzW8VM1Pba62jU2NVN+ooryE1U91hEgIz2F7MxUcrJSGZHpJK+sjBQy05PJSEshMy2ZjPRkMtOSSUqMty+YUcoS1HnY9/7JYL97SnIi82ZO8DkiY2KDiJCTlUZOVhozp4zt9Fp78jp+qoZTVbWcrDzLqapaTlWdpaLyLBVVZ3ttfQXPA9ScrafmbH2nlQZ6khAf6DZxpaUmkZ6SRGpKAqnJSaSnJpGanEhqSiJpKYmkJifaM2EeswR1Hv708rbg9mUXF5GYYH+MxgxWaPKCD46KVVXO1jU6CavqLBWVtZypqaOyps75XX2OM7VOYhqI5pZWTp85x+kz5wYcc2JCPGkpiaSlJAUTV3JSIslJ8SQnJpCcGE9ycqLzOynBKUt2yxMTSEpKICUpwX0t3hJeF/Y/6wAdLD/FO9sPBvdvcUc3GWO8JSJkpCWTkZZMUUHP8162tLRSfbaequo6qmrrnN81ddSea6DmXAO1Z93f5+qpPttAS8hsGgPV1NxCU3NLpzW6BiNOhISEeBITAu5PPAnxzu/2/cSEQEedeOf1pMR44uPjSEpMICE+jsSEeOIDccQHAgQCcSQkBNx9pyw+EEd8fBwBdzshPkB8fMhrgTgCgTjfuz49TVAisgR4FAgAP1fVb3V5PQn4NTAPOA18QlUPua/dD9wFtAJfUtW1/Tmnl87WNfL9J18K7l86q7DXfyjGmMiLjw+QOyK9XzO1qyqNTS0fSFw1ZxuorWukrr6RuoZmztU1UtfQxLn6puB2XX0j2uc7DEybKo1Nzf3qyoyEQJeEFprgPvWRhZ7ff/csQYlIAHgMuBEoBzaJyGpVfS+k2l1AlapOEZHlwCPAJ0RkBrAcmAmMBV4SkWnuMX2dM2xOVZ3l9JmznKtv4siJM6xZt4OTlc6DiokJ8Xz6toVevK0xJkJExOleS0pg5AAXdlRVGhqbnaRV7ySscw1NNDQ6CaahsYX6xmYaG5tpCN1vaqahsZmGphYaGpqc341OWdeppvzW2tpGa2sbjd28Fokk6mULagFQqqoHAERkBbAUCE0mS4EH3e1VwA/FaVMuBVaoaiNwUERK3fPRj3OGze/Xbualt3d1+9rnP3EVBaOyvXhbY0wMEBFSkhNJSU4kL0z/FbS0tNLc0kpTc6vTfdjSSnNzS8e++7u5uZWmlvb9kLL2Oi0ttLa20dzSRmtrKy2tbbS0tNHSvt3aRktLq/PT2kZrm1O3fb+ltY22trZeYw0EvL9f5mWCGgeUheyXA12bHME6qtoiItVArlv+Tpdj2x+y6OucAIjI3cDdABMmnN8ou9TkxA+UpSQn8vd3XM0Vcyaf1zmNMaYn8fEB4uMDpETBhBqqSktLe/Jq7UhqbgLLyfJ+9XAvE1R3d9e6tl97qtNTeXdP33XbJlbVx4HHAebPn39e7eZRuZlMmTCS1OREsrNSmTZxFIvnTyEtJel8TmeMMTFDREhICJBAgOSkBF9i8DJBlQPjQ/YLgKM91CkXkXggC6js49i+zhk27Yu+GWOMiTwv5wPZBEwVkSIRScQZ9LC6S53VwJ3u9jLgFVVVt3y5iCSJSBEwFdjYz3MaY4wZAjxrQbn3lO4F1uIMCf+FqpaIyENAsaquBp4AfuMOgqjESTi49VbiDH5oAe5R1VaA7s7p1TUYY4zxj2iUDWv0wvz587W4uNjvMIwxxgAisllV5/dVz6b8NcYYE5UsQRljjIlKlqCMMcZEJUtQxhhjopIlKGOMMVFpWIziE5EK4H2/4wDygFN+BxFGQ+l67Fqi01C6Fhha1zOYa5moqvl9VRoWCSpaiEhxf4ZWxoqhdD12LdFpKF0LDK3ricS1WBefMcaYqGQJyhhjTFSyBBVZj/sdQJgNpeuxa4lOQ+laYGhdj+fXYvegjDHGRCVrQRljjIlKlqCMMcZEJUtQHhKRQyKyQ0S2iUixW5YjIn8VkX3u72y/4+yOiPxCRE6KyM6Qsm5jF8cPRKRURLaLyFz/Iu9eD9fzoIgccT+fbSLy4ZDX7nevZ4+I3ORP1B8kIuNF5FUR2SUiJSLyZbc8Jj+bXq4nFj+bZBHZKCLvutfy7255kYhscD+bp9217HDXu3vavZYNIlLoZ/yhermWX4nIwZDPZY5b7s3fM1W1H49+gENAXpeybwP3udv3AY/4HWcPsV8FzAV29hU78GHgBUCAy4ANfsffz+t5EPh6N3VnAO8CSUARsB8I+H0NbmxjgLnudgaw1403Jj+bXq4nFj8bAdLd7QRgg/tnvhJY7pb/BPiCu/33wE/c7eXA035fQz+u5VfAsm7qe/L3zFpQkbcUeNLdfhL4qI+x9EhV1+EsIhmqp9iXAr9WxzvACBEZE5lI+6eH6+nJUmCFqjaq6kGgFFjgWXADoKrHVHWLu10L7ALGEaOfTS/X05No/mxUVc+6uwnujwLXAavc8q6fTftntgq4XkQkQuH2qpdr6Yknf88sQXlLgb+IyGYRudstG6Wqx8D5xwmM9C26gesp9nFAWUi9cnr/Tyaa3Ot2SfwipLs1Jq7H7RK6BOfbbcx/Nl2uB2LwsxGRgIhsA04Cf8Vp4Z1R1Ra3Smi8wWtxX68GciMbcc+6Xouqtn8uD7ufy/dFJMkt8+RzsQTlrStVdS5wM3CPiFzld0Ae6e5bXyw8v/BjYDIwBzgGfNctj/rrEZF04BngK6pa01vVbsqi6lqg2+uJyc9GVVtVdQ5QgNOyu7C7au7vmLoWEZkF3A9cAFwK5AD/7Fb35FosQXlIVY+6v08Cf8D5C3uivenr/j7pX4QD1lPs5cD4kHoFwNEIxzZgqnrC/UfYBvyMjq6iqL4eEUnA+c/8d6r6rFscs59Nd9cTq59NO1U9A7yGcz9mhIjEuy+Fxhu8Fvf1LPrfDR0xIdeyxO2SVVVtBH6Jx5+LJSiPiEiaiGS0bwMfAnYCq4E73Wp3An/yJ8Lz0lPsq4HPuiN5LgOq27ubolmXPvKP4Xw+4FzPcneUVREwFdgY6fi6496jeALYparfC3kpJj+bnq4nRj+bfBEZ4W6nADfg3FN7FVjmVuv62bR/ZsuAV9QdceC3Hq5ld8iXIMG5lxb6uYT/71mkR4cMlx9gEs5oo3eBEuAbbnku8DKwz/2d43esPcT/FE7XSjPOt6O7eoodp3n/GE5/+w5gvt/x9/N6fuPGu939BzYmpP433OvZA9zsd/whcS3C6TrZDmxzfz4cq59NL9cTi5/NRcBWN+adwANu+SScJFoK/B5IcsuT3f1S9/VJfl9DP67lFfdz2Qn8lo6Rfp78PbOpjowxxkQl6+IzxhgTlSxBGWOMiUqWoIwxxkQlS1DGGGOikiUoY4wxUckSlDEhRERF5Lsh+18XkQfDdO5ficiyvmsO+n1ud2cHf7VL+TUi8pzX729MuFiCMqazRuDjIpLndyChRCQwgOp3AX+vqtd6FY8xkWAJypjOWoDHga92faFrC0hEzrq/rxGR10VkpYjsFZFvicin3PV0dojI5JDT3CAib7j1bnWPD4jIf4nIJncSzv8bct5XReR/cB5+7BrPHe75d4rII27ZAzgPv/5ERP6rm+vLFJE/iMh7IvITEYkLvRZ3e5mI/Mrdniwi77ixPRRyzWNEZJ04awLtFJHFA/lDNqY/4vuuYsyw8xiwXUS+PYBjLsaZGLQSOAD8XFUXiLMA3xeBr7j1CoGrcSZCfVVEpgCfxZka5lJ3dug3ReQvbv0FwCx1lpYIEpGxwCPAPKAKZ9b8j6rqQyJyHc5aSsXdxLkAZ02l94EXgY/TsRREdx4FHlXVp0Tk8yHlnwTWqurDbusuta8/IGMGylpQxnShzmzavwa+NIDDNqkzkWYjznQv7QlmB05SardSVdtUdR9OIrsAZ57Gz7pLG2zAmbZoqlt/Y9fk5LoUeE1VK9RZquF3OIsy9mWjqh5Q1Vac6Z8W9VH/cpzpeAD+J6R8E/C37v252eqs5WRMWFmCMqZ7/41zLyctpKwF99+MO1lmYshrjSHbbSH7bXTuqeg6t5jizGP2RVWd4/4UqWp7gjvXQ3znu7Bdd+/ftTy5z5M4C0BeBRwBfiMinz3PeIzpkSUoY7qhqpU4S3XfFVJ8CKdLDZwVRBPO49S3i0ice19qEs6Ep2uBL7jLTiAi09wZ8HuzAbhaRPLcLrY7gNf78f4LRKTIvff0CWC9W35CRC50yz8WUv8d4G/c7eXthSIyETipqj/DmY18bj/e25gBsQRlTM++C4SO5vsZTlLYCCyk59ZNb/bgJJIXgM+ragPwc+A9YIuI7AR+Sh/3h9VZyuB+nKUc3gW2qGp/lm55G/gWzmzUB3HWKQO4D3gOZ7bq0GUSvgJ8zb3mMTirvgJcA2wTka04CezRfry3MQNis5kbY3okIqlAvaqqiCwH7lDVpX7HZYYHG8VnjOnNPOCH7j23M8D/9jkeM4xYC8oYY0xUsntQxhhjopIlKGOMMVHJEpQxxpioZAnKGGNMVLIEZYwxJir9P/vraw+x7g/YAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n_marginal = suite.Marginal(0)\n", "\n", "print('post mean n', n_marginal.Mean())\n", "print('MAP n', n_marginal.MaximumLikelihood())\n", "\n", "thinkplot.Pdf(n_marginal, label='n')\n", "thinkplot.decorate(xlabel='Number of bugs',\n", " ylabel='PMF')" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "post mean p1 0.22994325572783006\n", "MAP p1 0.2\n", "post mean p2 0.17521894732003349\n", "MAP p2 0.13333333333333333\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8XPWZ6P/Poxl1Sy6SXGVbcsUVg2VTDCahBEMITghJ4JIEUi4/yGazm+zNTbZld9mabbnJTXJDSQJpEAIpTigGQi8G29i4IvciV1mSJVl9Zp7fH+fM6GhQGUlzRiPpeb9e8/I5Z075zpE8j77lPF9RVYwxxph0kzHUBTDGGGO6YwHKGGNMWrIAZYwxJi1ZgDLGGJOWLEAZY4xJSxagjDHGpCULUMYYY9KSBShjjDFpyQKUMcaYtBQc6gKkQnFxsZaVlQ11MYwxxgCbN28+o6olfe03KgJUWVkZmzZtGupiGGOMAUTkcCL7WROfMcaYtGQByhhjTFqyAGWMMSYtjYo+KGOMGY46OjqoqqqitbV1qIsyIDk5OZSWlpKZmTmg4y1AGWNMmqqqqqKgoICysjJEZKiL0y+qSk1NDVVVVZSXlw/oHNbEZ4wxaaq1tZWioqJhF5wARISioqJB1f4sQI0wkYjNkGzMSDIcg1PUYMtuTXwjhKpy/wsH2XK4jvcvnMhHV5QSyBi+v9jGGONrDUpE1ohIpYjsE5Gvd/P+ahF5W0RCInKzZ/v7RWSr59UqIh9233tQRA563lvm52cYLvafbuKtA7V0hJVntp/i/z6zl5b28FAXyxgzQn33u99lzpw5iAhnzpzx5Rq+BSgRCQDfA64DFgK3isjCuN2OAHcAv/BuVNUXVHWZqi4DrgSagWc8u3w1+r6qbvXrMwwnmw/WdVnfUdXAv/3+Xc40tg1RiYwxI9mqVat47rnnmDlzpm/X8LOJbyWwT1UPAIjII8BaYFd0B1U95L4X6eU8NwNPqWqzf0Ud3lSVzYfq3rP9WF0L/7JuN1+8Zg6zJo4ZgpIZY5Llo3/2A9/O/fi37+rxvUOHDrFmzRouuugitmzZwrx58/jJT37CBRdc4Ft5ovxs4psGHPWsV7nb+usW4OG4bf8sIttE5Fsikj3QAo4Uh840U3uuHYC8rAC3Xz6ToNv/1NAS4j+eqGTjgdqhLKIxZhirrKzkzjvvZNu2bRQWFvL9738/Jdf1M0B110PfryFmIjIFWAKs92z+S+A8YAUwAfhaD8feKSKbRGRTdXV1fy477LztqT0tmzmOy+eX8OXr5pGfHQCgI6zc+/wB/rDlOKo2ys8Y0z/Tp09n1apVAHzyk5/k1VdfTcl1/WziqwKme9ZLgeP9PMfHgd+oakd0g6qecBfbROTHwP/q7kBVvQ+4D6CiomLEfiurapf+p+Vl4wGYP6WAv167gG+v38upeqcf6rebj3OyvpXbLy8jM2BPGBgznPTWDOe3+OHiqRr67ue31EZgroiUi0gWTlPdun6e41bimvfcWhXi3KEPAzuSUNZhq6q2hdMNTgDKzsxgwbTC2HsTC3P4qxsXcN7Ugti2Dftq+a8n99DQ0vGecxljTHeOHDnCG2+8AcDDDz/MZZddlpLr+hagVDUEfBGneW438Kiq7hSRe0TkRgARWSEiVcDHgHtFZGf0eBEpw6mBvRR36p+LyHZgO1AM/JNfn2E48A6OOH/6OLKCXX+k+dlB/vzauayeXxzbtu/UOf513bscr2tJWTmNMcPXggULeOihh1i6dCm1tbXcfffdfOc736G0tJSqqiqWLl3K5z//+aRfV0ZDn0RFRYWO1AkL//axHZw466QSufuq2SwvH9/tfqrO81GPbawi+iPPywrwF9fPY2ZxfqqKa4zph927d7NgwYIhLcOhQ4e44YYb2LFjYI1V3X0GEdmsqhV9HWsdEcPY8bqWWHDKDAiLSwt73FdEuHbpZP7k6jmxWlZze5jHNx5LSVmNMaa/LEANY97Re0umjyU7M9DnMctmjuOrH5wfW6880UhzW8iX8hljhr+ysrIB154GywLUMOYdvVdRPiHh48pL8ikryQMgHFF2VDUkvWzGGDNYFqCGqdMNrRytdQY5BDOEJdPH9uv482eMiy1vOfzeLBTGGDPULEANU97a06LSQnKz+m7e87pgZmeA2lHVQCjcW7YpY4xJPQtQw5R3ePmFZd2P3OvNtPG5FI3JAqClPcyek+eSVjZjjEkGC1DD0JnGNg5VO7lzAxnCMk9tKFEi0qUWZc18xpj+uO2225g/fz6LFy/ms5/9LB0dyX/43wLUMLTl8NnY8nlTC8jPHljGqvM9AWrr4bOWp88Yk7DbbruNd999l+3bt9PS0sIDDzyQ9GvYjLrDUHe59wZi7qQx5GUFaG4PU9fUwdHaFmYU5SWjiMaYJPv8A/4lG3jg8z0/M9vTdBvXX399bJ+VK1dSVVWV9HJZDWqYOdvUzv7TTn+RCANq3osKBjJYOqNz9N+WbuaUMsaY3qbb6Ojo4Kc//Slr1qxJ+nUtQA0zbx8+G0tVNH9KAYW5mYM6n3e4+TtH6gd1LmPMyNTbdBtf+MIXWL16NZdffnnSr2tNfMPM24McvRdvSelYghlCKKIcqWnmTGMbxQWjfg5IY9JOb81wfutpuo1/+Id/oLq6mnvvvdeX61oNahhpaOmg8kQj4DTvXTiI5r2onKwA8z3Tcbxz5GwvextjRqPuptt44IEHWL9+PQ8//DAZGf6EEgtQw8hWT/Pe7IljGJeflZTzLrNmPmNML7qbbuOuu+7i1KlTXHLJJSxbtox77rkn6de1Jr5hxPtwbk/TagzEspnj+PnrR4DO5LF5Axy6bowZeTIyMvjBD37QZVso5H+SaatBDRNNbSHePd4YW7+wbPDNe1Hj87MseawxJu1YgBomth4+SzjitO+VleRRNCa5Axm8o/m2HrZ+KGOMw6bbMH1K9ui9eN60R9ur6i15rDFpYjhneBls2S1ADQMt7WF2eprdKpLY/xRlyWONST85OTnU1NQMyyClqtTU1JCTkzPgc1hP+DCw/Wg9Ibd5b/qEXCYWDvwH3pNo8tjndp4GnOSxC6f1PIW8McZ/paWlVFVVUV1dPdRFGZCcnBxKS0sHfLwFqGFg08Ha2PKFPtSeos73BKith8/yPy6Z8Z4H9IwxqZOZmUl5eflQF2PI+NrEJyJrRKRSRPaJyNe7eX+1iLwtIiERuTnuvbCIbHVf6zzby0XkTRHZKyK/FJHkPAyUpto6wmw/2vls0mCSw/YlmjwWiCWPNcaYoeJbgBKRAPA94DpgIXCriCyM2+0IcAfwi25O0aKqy9zXjZ7t3wS+papzgTrgc0kvfBrZUdVAR9hp3psyLoep43N9u5YljzXGpBM/a1ArgX2qekBV24FHgLXeHVT1kKpuAxIaMiZOe9OVwGPupoeADyevyOmny9QaPjbvRVnyWGNMuvAzQE0DjnrWq9xticoRkU0iskFEokGoCDirqtFHmHs8p4jc6R6/abh2MKoqO6pS07wXtdhNHgtwpKaZmnNtvl/TGGO642eA6q53vT9jJWeoagXwP4D/IyKz+3NOVb1PVStUtaKkpKQfl00fDS0hmtvDAORkZlA6wb/mvajcuOSx9tCuMWao+BmgqoDpnvVS4HiiB6vqcfffA8CLwAXAGWCciERHH/brnMPNqYbW2PLEwpyUjaiz5LHGmHTgZ4DaCMx1R91lAbcA6/o4BgARGS8i2e5yMbAK2KXO02ovANERf7cDv0t6ydPEqfrOADVpbOrmaPL2Q0WTxxpjTKr5FqDcfqIvAuuB3cCjqrpTRO4RkRsBRGSFiFQBHwPuFZGd7uELgE0i8g5OQPo3Vd3lvvc14Csisg+nT+qHfn2GoXa6obP/Z5IPD+f2ZMKYLGYWW/JYY8zQ8vVBXVV9Engybts3PMsbcZrp4o97HVjSwzkP4IwQHPG8NaiJKaxBgTMFx+EzzYDTD7Vy9oSUXt8YYywXXxobqhoUWPJYY8zQswCVplS1S4BKdQ3KkscaY4aaBag0dba5g/aQU2vJywowJsUz3IoIy2baHFHGmKFjASpNxfc/DUXSVm+A2nK4blim/DfGDF8WoNLUUPY/RVnyWGPMULIAlaa6PqSb2v6nqPjkse8csWY+Y0zqWIBKU6frPTWosUNTgwJYUtoZoCpPNA5ZOYwxo49NWJim0qEGFQqFOXPyNJUHTyIihMNhOsIRMgP2d40xxn8WoNKQqlLdMHQ1qMamVp55fRdPvbyDuoZmWmUCHRJg1/4TbD1whhVzJ6a0PMaY0ckCVBqqbWqPTVKYnx0gP0VDzI+dPssTL27n+TffpSMUjm3PpZ0OcmnrCPMfP3+VH/3vDzImb2hqdcaY0cMCVBpKZf+TqrJj73H+8OI2Nu08/J73xxXksaC8lCd21KDA0bNt/Mt9T/F3X/gg2VmZvpbNGDO6WYBKQ10ySPjU/xQKhXlty37WvbCNQ8fOvOf9smnFfOh9S1h1wRwa28LsrHudg8dqaJVM3j14kv/88bN87XPXEgwGfCmfMcZYgEpD3gESftSgmlra+Lvv/p6DVe8NTMsXzuRD71/K4rlTYw8HT8gMsGBGEaFwhKMn62gjyNu7jvC9h1/kS5+8ckgeIjbGjHwWoNKQ3w/pPvXKzi7BKTMY4MqLzuOD71vCtInjuj1m3uQCTje0EQqFaaluIkdDvLxpLwX5OXzmI5dakDLGJJ0FqDTUJc1Rkpv4QqEwT7+yI7Z+wxVLufnaCynI7z0Qzp9SwKt7zjBt0jhasqH16H4AnnhpO4Vjcrn5AxcmtZzGGGMPtKSZSMTfIeZvbD1AXYMzz9O4gjw+deNFfQYngHlTxrhLQmHRBFYuLY+99/ATb7H+1Z3dH2iMMQNkASrN1Da1E4o4Q8wLcoLkZiVvEIKq8vsXt8XW11y+KOFBDkVjsikucKbf6Agra9dczNJ5nXNN3v+rV3hty/6kldUYYyxApZnTPtaeKg+eYv/RagCCwQAfuHRhv46fP6Ugtnygupn//bkPMHt6CQAKfPunf+SdyqqkldcYM7pZgEozfvY//eGl7bHl1cvnMrYgt1/Hz5vcGaAqTzaSm5PF39x1fWxgRTgc4ZsPrGfv4VPJKbAxZlSzAJVm/KpBVdc28uY7B2LrN7xvSb/P4a1B7Tt1jlA4QuGYXL7xhRsoGpcPQFt7B//0gyepOWsz8BpjBscCVJrxqwb19Ks7ibgTDi6eO5WZU4v6fY7iguzYNPBtHRGO1DiDLYrHj+EbX7ghlv7oXHMbv352S5JKbowZrXwNUCKyRkQqRWSfiHy9m/dXi8jbIhISkZs925eJyBsislNEtonIJzzvPSgiB0Vkq/ta5udnSDU/HtJtbevgmdd2xdZveN/SAZ/LW4t61zP9Rumk8fzZp66KrT+34V1q65sGfB1jjPEtQIlIAPgecB2wELhVROJ75Y8AdwC/iNveDHxaVRcBa4D/IyLeJ0i/qqrL3NdWXz7AEAhHlDON7bH1SUmqQb341h6aW53zTi4upGLRzAGfa54nQO2Jmx/qggXTY4MmQqEw655/Z8DXMcYYP2tQK4F9qnpAVduBR4C13h1U9ZCqbgMicdv3qOped/k4cBoo8bGsaaHmXBthd4j5uLxMsjMHP8RcVXnipc6h5devXjKorA/eGtTeU+di5QUQET62Znlsff1ru6hvtGnijTED42eAmgYc9axXudv6RURWAlmA9yGbf3ab/r4lIt1WM0TkThHZJCKbqqur+3vZIXGqPvlJYrfsPsrx6noAcnOyuPKi+YM6X/GYLCZ4+qEOn+najFexaGasf6u9I8QfPM9dGWNMf/gZoLr7M1272dbzCUSmAD8FPqOq0VrWXwLnASuACcDXujtWVe9T1QpVrSgpGR6Vr9NdZtFNTv/TE56h5VdddB65OVmDOp+IMG/ymNh6/DTwIsLN13amPXrylR2ca27DGGP6y88AVQVM96yXAscTPVhECoEngL9R1Q3R7ap6Qh1twI9xmhJHhC7TbIwdfA3q6Mk6tr7rVGIFuP6KxYM+J3Rt5ttz8r3DyS85fxalk8YDzgANb5A0xphE+RmgNgJzRaRcRLKAW4B1iRzo7v8b4Ceq+qu496a4/wrwYWDHe88wPHmHmCcji7m372nl0nImFRUO+pwQ1w91srFLPxQ4taiPfuACTzm209zSjjHG9IdvAUpVQ8AXgfXAbuBRVd0pIveIyI0AIrJCRKqAjwH3ikg04+jHgdXAHd0MJ/+5iGwHtgPFwD/59RlSretDuoOrQTU2tfLiW3ti6x+8ov8P5vakpCCb8fnObLqtnuehvFZdMIfJxU5AbGpp42lLJmuM6Sdfp9tQ1SeBJ+O2fcOzvBGn6S/+uJ8BP+vhnFcmuZhpIRSOcKaxM0CVDHKQxLOv76YjFAac2XEXzp4yqPN5iQjzpxSwYV8t4Aw3Ly/J77JPIJDBTddcwPcffgmAdS+8w/WrF5OTbdPEG2MSY5kk0sSZxnaiLWXj8zPJHsRU6qFQmKdf9c75NLih5d3xNvPFD5SIuqJiHsXjnQEVjU2tPPv67qSWwRgzslmAShOnkjiCb8O2g9ScdYZ/F47J5bIL5wzqfN2Jfx4qEnnvAM1gMMBHrursi/rd81tp7wglvSzGmJHJAlSa6DrN++Ca97zPHq25bBGZSXjgN15JQTbj8pzmupb2cLf9UABXXjyf8YV5ANQ1NPP8hsqkl8UYMzJZgEoTXZLEDiIH355Dp9h7+DTg9ANde1n/5nxKVLQfKnbdk90382VlBll7ZWe6xN/8cQsht2/MGGN6YwEqTSSrBuWd8+ny5XMZV5A3qHL1psv8UD30QwF8YNUCCsc4c0+dqTvHS5v29LivMcZEWYBKE8moQZ2pO8cbWz1zPiVxaHl3uj4P1X0/FEB2ViYf8mRQ//WzWwiHI93ua4wxURag0kBHOEJtk/Mgq4jTvzMQ61/dSSTifPEvnD2F8tLipJWxO5PGdvZDNbeHqarrOTHsdZcvIj/X+VwnzzTw6tv7fC2bMWb4swCVBqob2tDYEPMssoL9/7F0dIR55vXkzPmUKCcvX2LNfLk5WV1m8X38mbdR7VdqRmPMKGMBKg0ko//p3YMnY0lZS8YXsGLxwOd86o9EnoeKun71ktiDusdOn+V1T3OkMcbEswCVBpLR/7RjX2ce3gsXziAjIzU/2nlTOjOb7znZ2GutaExeNh9c3VmLemz9ZqtFGWN6ZAEqDXSZ5n2ANagde4/FlhfNnTroMiVq8tgcCnOdjFnNbWGqanufoPCDVywhK9PZ/8iJWjbuOOx7GY0xw5MFqDTQtYmv/zWotvaO2LNPAIvnpC5AxT8P1Vcz39iCXNZctii2brUoY0xPLEClga5NfP2vQb178FRs2Pb0yeMZW5CbtLIlostAiR4e2PX60PuXEnRzDe4/Ws2u/Sd8K5sxZviyADXE2kMR6po6gIEPMd+xp7N5b/HcaUkrW6K6ZJQ40Xs/FMCEsfm8f+W82Lp39KExxkRZgBpi1Z4pNorGZBEM9P9Hst3b/5TC5r2oKeNyKMhx+pWa2sIc6+V5qKhrV3U2872x9QD1jX0fY4wZXSxADbHTg5xFt6W1nf1HqgFnWvfFKRwgEdXffiiA8tJi5s6cCEA4HOH5N9/1rXzGmOHJAtQQO+UZIDGQ/qfdB04ScZvUZkwtoiB/8FPFD8S8fgYo6FqLevb13TZYwhjThQWoITbYGpR3ePmSIeh/iuqa2fxcQsFm1YWzycvJAuBUTQPvVFb5Vj5jzPBjAWqIDbYGtWNv5wO6qXz+Kd7UcTmMcfuhzrWGOF7X2scRzlQcV150Xmz9mddssIQxppMFqCF2umHgNaimljYOHO3sf1o0Z0oyi9YvTl6+zqwS755oSOi4a1YtiC1v3H6ImrPnkl42Y8zw5GuAEpE1IlIpIvtE5OvdvL9aRN4WkZCI3Bz33u0istd93e7ZvlxEtrvn/I6IiJ+fwU9toXBsiHmGOKP4+mPX/hNEG9LKp5fEsoUPlQVTC2PLO6sSC1Clk8bHRh5GVPnjBhssYYxx+BagRCQAfA+4DlgI3Coi8dO7HgHuAH4Rd+wE4O+Ai4CVwN+JyHj37f8H3AnMdV9rfPoIvqv2NO8VF2T3e4j5jj2dzXupzB7Rk0WlnQGq8mQjoQTnfPrAqs5fi2df321zRRljAH9rUCuBfap6QFXbgUeAtd4dVPWQqm4D4r+RrgWeVdVaVa0DngXWiMgUoFBV31CnF/4nwId9/Ay+OlXv6X8aQA4+b4LYoRheHq+kIDtWC2zriLD/dFNCx128tDw2425tfRObdx3xrYzGmOHDzwA1DTjqWa9ytw3m2Gnu8kDOmXa6JIntZxbzxqZWDh87A0CGCAtnD13/U5SIdKlF7T6WWDNfMBjg6ou9gyV2Jr1sxpjhx88A1V3fUKIPuvR0bMLnFJE7RWSTiGyqrq5O8LKpdXoQNShv/9PsGSXk5vSv/8ovC6d1BqhdxxMLUABXX7og9sPduvsoJ88kfqwxZmTyM0BVAdM966XA8R72TfTYKne5z3Oq6n2qWqGqFSUlJQkXOpUGU4PyPv+UDv1PUedNKSQ6bOVgdRNNbaGEjptUVMiyBc6PXIHnLD+fMaOenwFqIzBXRMpFJAu4BViX4LHrgQ+IyHh3cMQHgPWqegJoFJGL3dF7nwZ+50fhU2Ew02x4n39aPC99WjnH5ASZWZwHgCq8ezyxrBIAH/Bklvjjm5WEQuGkl88YM3z4FqBUNQR8ESfY7AYeVdWdInKPiNwIICIrRKQK+Bhwr4jsdI+tBf4RJ8htBO5xtwHcDTwA7AP2A0/59Rn81Noepr7ZGWIezBAm9GOIeX1jC0dOOLcjEMjgvPLJvpRxoBZNGxtb7k8z3/KFMygalw9Aw7kWNrxzMOllM8YMH0E/T66qTwJPxm37hmd5I12b7Lz7/Qj4UTfbNwGLk1vS1DvtzWJekEUgI/HHuXbu76w9zZkxkZzszKSWbbAWTC3gia3OHE+JDpQAJ9hefckCfvnUJgDWv7aTy5bP8aWMxpj0Z5kkhsipQeTg8z7/tCQNhpfHmz1pDFlB51frdENblylF+nL1JQvIcDuxdu0/wdGTdb6U0RiT/noNUCLyoGf59l52Nf3Upf+pnwMkdnZ5/il9+p+iMgMZXdIe7UowqwQ4kxmuXFIWW3/WBksYM2r1VYM637P8Z34WZLTpMs17P4aY1zU0U3XKqVUEAhnML5+U9LIlw6LSgfVDAXzgss7BEi+8WUlbe0fSymWMGT76ClA2QY9PBlqD2ukZvTe/bBJZmb52Iw7YQk9evt3HG4hEEv9VWjpvGpOLneObW9t5fcuBpJfPGJP++gpQpW5C1v/rWY69UlHAkWqgNagd+zzTu6dh/1PU1PE5jMtzBm80t4U5XNOc8LEiwjWXdubnW2+ZJYwZlfoKUF8FNgObPMvelxmA5rYQja3OA6zBgPQri7n3+aehnKCwLyLCeVM7JzHc1Y/RfABXXjSfgJs8d+/h0xysOpPU8hlj0l+v7UOq+lCqCjKaeCcpLCnIJtEZQ2rOnuNEdT0AmcEA82amZ/9T1KJpY9mwz3lea9exBj64LPF8gYVjcrl02Wxe2bwXcGpRd33iCl/KaYxJT70GKBHpNfODqt6Y3OKMDtUD7H/y1p7ml08iMzOQ1HIl2wJPDWr/qXO0dYTJ7keZr121MBagXt60j0/feAl5uemRc9AY47++etgvwckq/jDwJt0nazX95M3B16/+p73pPbw83rj8LKaNz+VYXQuhiLLn5DmWTB/b94Gu82ZNZvrk8Rw9WUdbewcvb9rLmssX9X2gMWZE6KsPajLwVziZG74NXAOcUdWXVPUlvws3Ug30IV1vgth07n/y6pLdvJ/9UCLSZTLD9a/txJkGzBgzGvQaoFQ1rKpPq+rtwMU4+e9eFJE/TUnpRijvEPOJYxOrQZ2ubeR0rZN4NSszyJwZ6ZmhPd5gAhTAFSvmxYbSHzlRS+XBU0krmzEmvfWZ6khEskXkJuBnwJ8A3wF+7XfBRrLac+2x5eKCxAKU9/mnBbMmEwymd/9T1LzJYwi6eQaP1bVwtqm9jyO6ys/NZnXF3Nj6U6/uSGr5jDHpq69URw8BrwMXAv+gqitU9R9V9Vhvx5mehcIR6luczAgiMD4vsUSv273zPw2T5j2A7MwAsyd1pj3a3Y/pN6LWeDJLvLH1APWNLUkpmzEmvfVVg/oUMA8nzdEbItLgvhpFxKY8HYC65g6i3ShjczMJBvrO16uqXScoTOMHdLvjbebbeay+38eXlxYzd+ZEAMLhCH/c8G7SymaMSV999UFlqGqB51XovgpUtbC3Y033vM17ic4BdfJMAzVnmwDIyc5k9vTh0f8U5Q1Qu483Dmigw3WXd86w8sxru4hEIkkpmzEmffXVxJcjIn8uIt8VkTtFJD0Tvw0jXQJUfmIBypu9fOHsKbEMC8PFzKI88rKdPrP65g6O1fW/ie6SZbMYk+f011XXNfL27qNJLaMxJv309U33EFABbAeuB/7L9xKNcLVN/a9BDdf+p6iMDGHB1MGN5svKDHLVxefF1te/avn5jBnp+gpQC1X1k6p6L3AzcHkKyjSi1ZzzzKSbQIBS1S4j+BbPGV79T1ELBxmgAK65dGHsSfEtu45w8ox1gxozkvUVoGIT8ahqyOeyjAreJr6iMX0PMT9eXU9dg5MJPC8ni/LSIt/K5qeFpZ0BqvJEIx3h/vchTSkZy7IF0wFnHhibzNCYka3PCQu9I/eApTaKb3Bq+jlIYscez/Qac6aSkTG8+p+iSgqyY2mdOsLKvlPnBnSeaz1Dzp97YzftHfZ3kzEjVV+j+AJxI/eCNopv4FS1ax9UAoMkdngGSCwaps17UQu8o/kG2My3fOEMisc7z1Wda27jja02maExI5Wvf46LyBoRqRSRfSLy9W7ezxaRX7rvvykiZe7220Rkq+cVEZFl7nsvuueMvjfRz8+QTM3tYdo6nKat7MwM8rN7zwahql1G8C2ZN7wDVJd+qH5OAx+VkZHRJT/f0zZYwpgRy7cAJSIB4HvAdcBC4FYRWRi32+cEtIfXAAAgAElEQVSAOlWdA3wL+CaAqv5cVZep6jKch4UPqepWz3G3Rd9X1dN+fYZkix9i3tc8UDVnm2JZE3KyM5k5dXj2P0UtmFpA9CMfPtPMudaBNc9dffGC2FD7PYdO2WSGxoxQftagVgL7VPWAqrYDjwBr4/ZZizOUHeAx4Cp577f2rTjTfQx7/X1I94Dni3dWaXHCExumq7zsIOUl+QCowu4B1qLGFuRy8fmzYutWizJmZPIzQE3DmUsqqsrd1u0+7ijBeiC+mvAJ3hugfuw27/1tNwENAPfB4k0isqm6unqgnyGpajz9T0UJ9D/tP9pZ7uGWPaIn3ma+gfZDAVznGSzx8qa9NLW09bK3MWY48jNAdRc44nPc9LqPiFwENKuqN4X1baq6BOeZrMtxmgDfexLV+1S1QlUrSkrS48u93zUoT4CaNb3YlzKl2oIuefkaBjy/03mzJjNjygQA2jtCvPjWnqSUzxiTPvwMUFXAdM96KXC8p33cNEpjgVrP+7cQV3uKZlJX1UbgFzhNicNC/wOUp4lvhNSgZk/MJzvT+bWrOddOdePAaj4i0iXL+fpXbTJDY0YaPwPURmCuiJSLSBZOsFkXt8864HZ3+WbgeXW/ZUQkA/gYTt8V7ragiBS7y5nADcCwmSCoP0PMa+ubONvoPKCbnZXJtInjfC1bqgQDGcyfXBBb31k18Ga+1RVzyc5ypis5dvosO/bG//1jjBnOfAtQbp/SF4H1wG7gUVXdKSL3iMiN7m4/BIpEZB/wFcA7FH01UKWq3gddsoH1IrIN2AocA+736zMkW9c0R71nkfD2P42EARJei0oHP9wcIDcni/etmBdbt8ESxowsvmYnV9UngSfjtn3Ds9yKU0vq7tgXcaaZ925rApYnvaApEApHONvcOVHhuPzeJyrs2rw3MvqforyJYyuPNxKOKIGMgQXgay9bxPrXnMD01vZD1NY3MWFsflLKaYwZWsMzb84wdDZuosLMPqbMODACR/BFTRmXw3g3QDe3h9l7sv+z7EbNnDqBhbOnABCJRHjujd1JKaMxZuhZgEqR/k6z0aWJb4QFKBHh/BmdfWpvHajtZe++Xbuqc7DEM6/tIhQKD+p8xpj0YAEqRfozUWFdQ3Msg3lWZpBpE8f6WrahsHL2hNjy5oN1hAaQ3Tzq4vPLKRyTCzj3buOOw4MunzFm6FmASpH+DDH3Nu+VlxYP2wzmvZk7aUysma+pLczu4wNv5gsGA1xzyYLY+jOv2TQcxowEI++bL011ySLRR4DqmkFiZA2QiBIRVszqrEUNtpnvmksXxJ763ranimOnzw7qfMaYoWcBKkX608R3sEsOvpHV/+S10hOgthyqoz008Ga+kgkFVCwui63blPDGDH8WoFKkP018+0dgiqPuzCzOi01i2NoRYfvR+kGdr+tkhu9yrtny8xkznFmASpFER/HVN7ZQc7YJgMxggNJJ430v21ARkS61qME28y07rzR2v9raO3jqlWGTZMQY0w0LUCnQ3Baipd0Z+pwZEMZk9/x8tLf2VDatKDbv0UjlHc237chZWtsHPkRcRLjpmgti60+8tJ3Wto5Blc8YM3RG9rdfmoivPfWWtujAKOl/ipo6PpfSCc4Q8Y6wsvXI4AY3rLpgNiXjnVx/jU2t9uCuMcOYBagUqDnnHcHXew6+g94RfDNGbv+TV5dmvv2Da+YLBgOsver82Pq6F96xB3eNGaYsQKVAf0bw7ffk4BtpKY564h1uvvNYw4Cngo+66uLzYg/u1pxt4uVNewd1PmPM0LAAlQLeJr6igp4DVGNTK9V1zgOrwRE+QMKrpDCbWROdBK/hiPL2obpBnS8rM8iH3rc0tv6b57YQiQx8CLsxZmhYgEqBmsbEalDe/qeZUyYQDAZ8LVc6SeZDuwDXXraQ3BznXh+vrmfDtoODPqcxJrUsQKVAokPM9x/x9j+Njua9qBXl44mOHak80chZzz0biPzcbK7zPBf162e32Iy7xgwzFqBSoOtEhYnVoGaVjo4BElHj8rOY5860qwqbB9nMB3DD+5aS6dZCD1ad4Z3KqkGf0xiTOhagfBaOaGyiQoBxeb0EqBE8B1QikjmaD2BsQS5Xe5LIPv7M24M+pzEmdSxA+exsc3tsosLC3CBZwe5v+bnmNk7VONOfBwIZzJgyodv9RrILy8fFZtbdf7qJM42DT1V045Xnx7LB79p/gsqDJwd9TmNMaliA8lmiOfi8CWJnjLIBElEFOZksmFoQW0/GYImJEwq4fPmc2Pqvn90y6HMaY1LDApTPEn0Gav8ob96L8qY+2piEZj6Aj1x9QWwqjk07D3P4eE1SzmuM8ZcFKJ91nQeq5ywSFqAcF8wcTzDghJOjtS0cr2sZ9DmnTx7PyqXlsfVfP2e1KGOGA18DlIisEZFKEdknIl/v5v1sEfml+/6bIlLmbi8TkRYR2eq+fuA5ZrmIbHeP+Y70ltguDQykiW+0jeDzys0KsHR65xT3G5PQzAdw09WdSWRf27yPk2caknJeY4x/fAtQIhIAvgdcBywEbhWRhXG7fQ6oU9U5wLeAb3re26+qy9zXXZ7t/w+4E5jrvtb49RmSIZEmvqaWNk5UO3MhZWRkMGPq6Bsg4eVt5nvrQG1Snl+aM3MiS+eVAqDA757fOuhzGmP85WcNaiWwT1UPqGo78AiwNm6ftcBD7vJjwFW91YhEZApQqKpvqPOt9RPgw8kvevIk8pBu/ACJrMyep+MYDZZOH0d2pvOreaq+jaO1g2/mA7pMxfH8m5XU1jcl5bzGGH/4GaCmAUc961Xutm73UdUQUA8Uue+Vi8gWEXlJRC737O992rK7cwIgIneKyCYR2VRdXd3dLilR2yWTefcBypsgdjQ370VlBTO4YOa42Pqb+5MzqGHx3KnMmTERgFAozB9e3JaU8xpj/OFngOquJhTfVtPTPieAGap6AfAV4BciUpjgOZ2NqvepaoWqVpSUDM2gg+a2EM3uBHzBgFCQ033N6ECVDZCI583Nt/FAXVKa+eInNHz61V02LbwxaczPAFUFTPeslwLHe9pHRILAWKBWVdtUtQZAVTcD+4F57v6lfZwzbdQ1dWaQKMrveaLCA54cfLOmWw0KYNG0QvKynWfBas+1s/90cprjVi4ps2nhjRkm/AxQG4G5IlIuIlnALcC6uH3WAbe7yzcDz6uqikiJO8gCEZmFMxjigKqeABpF5GK3r+rTwO98/AyDUtPU+dd5T/1PLa3tnQMkRCibVtTtfqNNMJBBRVnndCPJSH0ENi28McOJbwHK7VP6IrAe2A08qqo7ReQeEbnR3e2HQJGI7MNpyosORV8NbBORd3AGT9ylqtFvqLuBB4B9ODWrp/z6DIOVyBDzg8dqYm2UpZPHj/oBEl4rvA/tHqglHElONvL4aeGffd2mhTcmHfn6baiqTwJPxm37hme5FfhYN8c9Djzewzk3AYuTW1J/JDLEfH+X5j3rf/KaP7mAsXmZ1Dd30NgaovJEIwunFQ76vMFggA9ftYz7H3sFgEef3sTqirmMLcgd9LmNMcljmSR8lMgQ864DJKz/ySsjQ6goT34zHzjTwk8udoJdc2s7P/v9m0k7tzEmOSxA+ajmXN9pjg54hpjbCL73usjTzPf2oTqa20JJOW9mZoDPffSy2Przb77LnkOnknJuY0xyWIDyUV9NfK1tHRw75UzMJ2ADJLpRXpLPpEInuDe3h3l+1+mknfvChTNYsbgstn7/Y68SiUSSdn5jzOBYgPJJJKLU9dHEdyhugER2VmaKSjd8iAjXL5sSW39mxyla3GfLkuEzN10am9rkwNFqnnvj3aSd2xgzOBagfHK2uYPooLOCnO4nKvRmMLcBEj27eE4RE6O1qLbk1qImFRXykauXxdZ//oc3aWxqTdr5jTEDZwHKJ4kNkLAUR4kIZMTVorafpDWJtaibrr4gNuz8XHMbv3jiraSd2xgzcBagfJLIM1A2B1TiLp49gZICpxbVlORaVFZmkM9+dFVs/dnXdnUZ/m+MGRoWoHxSc64zi0RRNwMk2to7qDrhDJsWoNxqUL0KBjK4ftnk2Pr6JNeiViyeyQULnMxcCtz/2CtJyf9njBk4C1A+6auJ7/Dx2tgAiakTx5GTbQMk+nLJnCKKC5x72dQW5vndyatFiQifvWkVgYDzX2Lv4dO88GZl0s5vjOk/C1A+6auJzwZI9F8wkMEHPX1R67edpK0jebWoqRPH8eErOwdM/PT3b9LUYtnOjRkqFqB80tczUNb/NDCXzCmKzavV1Bbmj0nsiwJnUsOicfkANJxr4ZEnNyb1/MaYxFmA8kmNp4mvuOC9WSS8GSRsio3EBQMZ3HCBf7WonOxM7vjIpbH1p17eweHjyZkw0RjTPxagfNDaHqa5reeJCts7Qhw9WRdbL59mAao/4mtRyRzRB3DJ+bNYOs+ZdkyB+35lAyaMGQoWoHzQZYBENxMVHj5eE0upM7VkLHm53Q9DN917T1/U9lNJrUWJCJ/96CoyMpz/Hu8eOMkrm/cm7fzGmMRYgPJBTR8DJLzNe+XW/zQgl84tit3bc60hXtid3OeWpk8ez4fetyS2/tBvN9Dc0t7LEcaYZLMA5QMbIOG/YCCDD57vX18UwMeuXc74wjwAzjY286v1m5N6fmNM7yxA+aCvZ6AqD3ZO62ApjgZu1bzOWlRja4gXk1yLys3J4o4Pdw6Y+MNL2zl07EwvRxhjkskClA+6ZJGIC1AnquupcqfYyAwGmFc2MaVlG0mCgQyuP78zu8TTPtSiVl04m4WznZpaJBLh33/4jCWTNSZFLED5oLaXiQrf3nUktrx0XqlNsTFIq+YVMz7fuYd+1KJEhDs/vjr2czpV08B/P/gc4bDNG2WM3yxA+SB+FJ/Xph2HY8sVi2emrEwjVWYgg+s9fVFPbztJWyi5tajpk8fzp7e9P7a+bU+VTRFvTAr4GqBEZI2IVIrIPhH5ejfvZ4vIL9333xSRMnf7NSKyWUS2u/9e6TnmRfecW91XWrWRRSLK2aaO2Pr4MZ01pOaWdnbuPx5bX75oRkrLNlJdNr9rLeqlJNeiAC5ZNoubr10eW1/3wju8tHFP0q9jjOnkW4ASkQDwPeA6YCFwq4gsjNvtc0Cdqs4BvgV8091+BviQqi4Bbgd+Gnfcbaq6zH0l9ynNQWpo6SDkzlQ4JidItjtbK8DWyqOxpqGyacUUjRszJGUcaVJRiwK45bqKLlPEf/+Rl9h3OK1+/YwZUfysQa0E9qnqAVVtBx4B1sbtsxZ4yF1+DLhKRERVt6hqtKqxE8gRkffmC0pDvTXvbd7Z2f9kzXvJ5a1FNbT4U4sSEb70ySspnTQegFAozL//aD1nG5uTfi1jjL8Bahpw1LNe5W7rdh9VDQH1QFHcPh8FtqiqN630j93mvb+V+DQNQ6ymywCJzgAViUS6DJBYscgCVDJlBjK4zlOL+u3m41TVJj9w5OVm8bXPX0tejvOzrTnbxL//8BlCPtTYjBnt/AxQ3QWO+IRmve4jIotwmv3+P8/7t7lNf5e7r091e3GRO0Vkk4hsqq5O3eyoPWWR2Hv4NA3nWgAYV5DH7Bn2gG6yXTavmEljnYp2eyjC95/bT3NbKOnXmTpxHF+545rYL2/lwZM88PirSb+OMaOdnwGqCpjuWS8Fjve0j4gEgbFArbteCvwG+LSq7o8eoKrH3H8bgV/gNCW+h6rep6oVqlpRUpK6YNBTFglv896FC2e8Jz+fGbysYAZfuHoO2ZnOr/Xphjbuf/GgL4leL1gwnU/eeHFs/dnXd7P+1Z1Jv44xo5mfAWojMFdEykUkC7gFWBe3zzqcQRAANwPPq6qKyDjgCeAvVfW16M4iEhSRYnc5E7gB2OHjZ+i3nrJIbNxxKLZs/U/+mTY+lzsuL4utbz9az7q34/8uSo61V57PqgvnxNYfePw1du0/4cu1jBmNfAtQbp/SF4H1wG7gUVXdKSL3iMiN7m4/BIpEZB/wFSA6FP2LwBzgb+OGk2cD60VkG7AVOAbc79dnGIjabvqgTtc2cuRELQCBQAbnzy8dkrKNFitmTWDN0s4ME7/fcoJ3jpxN+nVEhD+59QrK3XRVkUiE//jRM5ypO5f0axkzGvn6HJSqPqmq81R1tqr+s7vtG6q6zl1uVdWPqeocVV2pqgfc7f+kqvmeoeTLVPW0qjap6nJVXaqqi1T1z1Q1rXqnvWmOok18b+/0Zo+YRk62ZY/w20cqprFgakFs/YcvHuRUffJTFGVnZfK1z11L4ZhcwJmF998eeJq29o4+jjTG9MUySSRRW0eYpuhEhRnC2DwnEG3aeSi2z3IbvZcSgQzhf75/VqyZtbk9zPef25/0XH0AJRMK+OpnPxCbP+pg1Rm+/8hLNsmhMYNkASqJvCP4xrsTFba2dbBtz7HYdgtQqVOYm8ndV80mGHAGpByra+HBVw75EjgWzp7C525aFVt/dfM+/vNHz1hNyphBsACVRN0NkHinsiqWPWLGlAlMnFDQ7bHGH+Ul+XxyVecfBRsP1PHsjlO9HDFw1162kGsuXRBb37DtIH/znXXU1jf5cj1jRjoLUElU280zUF2Sw1rtaUhcNq+YKxZ0Pmrw2FtV7D7ekPTriAh3fuxybrhiaWzbgaPVfP2/f83BKptHypj+sgCVRF1G8OVnoaps3mXZy9PBLRdPZ9bEfAAiCvc9f6DLzytZMjIy+MxNl3Lnxy4nw33WreZsE3/97d91edTAGNM3C1BJVBPXxLf/SDX1jU72iIL8HObOTKvE66NKZiCDu6+aTWFuEHCynn//j/toD/kzr9O1ly3ir++6nlw3JVJbewffvP9pfv/CNhs8YUyCLEAlUXwT38adnbWn5YtmxkZ5maExPj+Lu66aTSDDqdkcqm7mkTeO9HHUwC07bzr/+uWPxPodFXjwt69z76MvW+4+YxJg35hJFJ/myNv/ZHM/pYd5kwv4+EWdD0q/XHmGxzdWEfJphtzpk8fzb1+5ifnlnQ8OP/v6bv753qdoamnr5UhjjAWoJFFV6jxNfIQ7OHTM6RgPBDJYNn96D0eaVLty4UQunjMhtv7UOyf5jycqOdPoT8AYW5DL3//JDVy2vDMt0rY9VfzVt37LyTPJH6xhzEhhASpJ6ps7JyrMzw6wo7Iq9t6i2VPJy83q6VCTYiLCpy6byaLSwti2/aebuOc3u9h8sM6Xa2ZlBvnzT13FJ66riG2rOlXH1//71+zc50+uQGOGOwtQSeJ9BqpoTLY176W57GCAP792LjetmIbbJUVze5j/98f9/Oy1w74MnhARPr6mgi9/+mqC7kzLjU2tfOP/ruM/f/wsJ6rrk35NY4YzC1BJ4s0iUZgTYNuezhqUZY9ITyLC9edP4Ws3nNdlcskXd1fzL+t2c7yuxZfrXrZ8Dvd88UOx/H0Ab2zdz5f+5Zfc/6tXYiM/jRntLEAliXeAREtzMx3uKK3SSeOZUjJ2qIplEjB70hj+7iMLubBsXGxbVW0L//S73bxSWe3LsPD55ZP597+4iYvPnxXbFolEePrVndx9zy949OlNtLZZmiQzulmASoLW9jCbD3X2XVSfro0t28O5w0NedpC7r5rNp1bNJNPN3dceivDQK4e5/4WDtLT7l2T2X7/8ERbO7pyuvq29g18+tYk/+ceHeea1XbFUWcaMNjIaHhqsqKjQTZs2+XLuc60hvr1+Lwero/nWlJZDe2htbATgH7+0tsuXj0l/VbXN3Pv8AU6c7Zyeo6QgmzuvnEV5Sb4v11RVNu08zM/WvUnVqa4DNaaWjOWTN17MyiVlNhOzGRFEZLOqVvS5nwWogWto6eC/n9pDVW1nn8Hq2YX87rfPA5Cfm82P//l2AgGrqA43bR1hHtlwlFcqO3PoicDCaYWsnl/C+TPGEvTh5xoOR3hxYyUPP7GRuobmLu/NL5/Mh963lAsXTic7y+YUM8NXogEqmIrCjES159r5r6cqOVXf+ezMbZfOoPpIZ2aCCxfOsOA0TGVnBrj98jIWTC3kJ68eorUjgirsrGpgZ1UDBTlBLp1XxOXzSpg8Lidp1w0EMrjq4gVcduEc/vDSdn7z3FZaWp3+zcqDJ6k8eJKszCAXLpjOJctms3zRjFg6JWNGGqtBDcDphlb+68k9sZF7IvCZ1WVcOreY//2fj7P/aDUAX/701V0ezjTD0+mGVn7x+hF2Hmugu/8u8yaP4fL5JSwvH09WMLl/kDSca+HxZ7bw1Ks7uu2LCgYDLJtfyiXLZlGxuIwxedlJvb4xfrAmPo9kBqjjdS3891N7ONvsjLAKujO3Li8fT11DM5//258AkCHCg/96B/m59oUxUlQ3tvFa5Rle3XMm9vP3yssKcPGcIlafV0zphLykXvt0bSN/3PAuG7YeeE8fVVRGRgZL503jkmWzWLG4jLEFud3uZ8xQswDlkawAdfhME996ei/nWkMAZAaEL1w9hyXTnWHkf9ywm+8//BIAi+ZM5Z4/vXHQ1zTpJxxRdlTV80rlGbYdOUukm/9Ck8ZmU1acz4yiPGYU5zGjKI/87OS0qB89WceGdw7wxtYDHD5e0+N+k4sLmTm1yH1NoGxaMZOKCmyghRlyFqA8khGg9p06x3fW76XZHW6cnZnBn14zh/OmdqbL+eYDT/PW9kMAfHrtJay98vxBXdOkv7NN7by2t4ZXKqs509j7/FLFBVnMKMpjZjRwFeUxNm9wgx1OVNfHglW0abk3OdmZzJxaRFksaBVROnk8eTlZFrhMyqRFgBKRNcC3gQDwgKr+W9z72cBPgOVADfAJVT3kvveXwOeAMPAlVV2fyDm7M9gAtft4A999dh9tHU4fQF5WgD9fM5dZE8fE9mnvCHHHXz1EW7vT9POdv76FaRPHdXs+M/KoKu8eb+SVyjO8fagulpexL+PyMikak0VhbiYFuUEKczPdl7uc4yznZgX6DCDVtY1seOcgb7xzgL2HThHpx//tzGCA8YV5jCvMc/4tyGNcYW7ntgLn3zF52WRnBS2YmUEZ8gAlIgFgD3ANUAVsBG5V1V2efb4ALFXVu0TkFuAjqvoJEVkIPAysBKYCzwHz3MN6PWd3BhqgauubeHX3SR556wTt4QiRiJIdgOvmFzImCG0dIdrbQ7S1h6itb4rNmDqlZCzf/Ztb+309MzK0hcIcq23h8JlmjtQ0c/hMM8frWhIOWt0JBoSCnCBZwQyygwGyghmxV2ZA3H87t2UAZxuaqK49x5naRqrrGjhV00BraweCUw4nxCg9hZro+91tz8wMkp0VIDszSGZmkJysTDKDAXKyg2RlBsjKDBIMZBAMBMgQCAQDZIgQCGQQyMggEBD33wwyRMjIEDJEnMAngggIQkaGk5JKRBCIBcZYgBS6BEtnn/jPEFd+C65JsaS8mOkTB5YlJx2Gma8E9qnqAbdAjwBrAW8wWQv8vbv8GPBdcX571gKPqGobcFBE9rnnI4FzJs0Pn97OLzZUxf4LBzXMFK3ndwd7zypQYbn3RrXsYIBZE8d0qWGHwhGOn23l8JkmjtQ0c7SmhSM1zQknpQ2FlbqmAaY+knyYkM+kCZPpCIVpbumgubWdlrZ2Wlo6aOvo6LYfrVdhoMV9OSV0XzbH1Whx11WzueuGZb5ew88ANQ046lmvAi7qaR9VDYlIPVDkbt8Qd+w0d7mvcwIgIncCdwLMmNH/bOKqyoG6UCw4ZWqYKXqWTHr/QgkEMnj/Ref1+3pmZAsGMmL9TlGRiFLd2MbZ5g4aWpxXY0vIWW4NxbY1tISSlF1dyAwGGVsQjBvhp4TDSkco3O0rFArTEYrQ0REiHIn0P5gZM0B+Bqju6tHxv9o97dPT9u4eMun2v4uq3gfcB04TX8/F7J6IcEvFZA4er4WMDFZOzKAwdzzZWUGyM4NOU4a7nJ3lNHFkZQaZWzaRSUWFfV/AjHoZGcKksTlMGtv3g75tHWEaW51A1R6K0B6O0OEud4SVtlC4871QhFBECUcUVWfUYUSVSEQJqxJRJzjGtrv/O1Sd/dVdhuiy8360T0sjSigccV9hQuEI4VCEjuiy+14kokQizgPOEXX/da+psWs75SLumqh2liO+THH/mzX6FRC/fRQMABtKE8cl91GK7vgZoKoA7zSypUD8zGzRfapEJAiMBWr7OLavcybNFRVzqFhSRjiiFOZaahkzdLIzA2RnBoa6GMaklJ95eDYCc0WkXESygFuAdXH7rANud5dvBp5X58+edcAtIpItIuXAXOCtBM+ZVPnZQQtOxhgzBHyrQbl9Sl8E1uMMCf+Rqu4UkXuATaq6Dvgh8FN3EEQtTsDB3e9RnMEPIeBPVDUM0N05/foMxhhjho49qGuMMSalEh1mbqm2jTHGpCULUMYYY9KSBShjjDFpyQKUMcaYtGQByhhjTFoaFaP4RKQaODyIUxQDZ5JUnOHK7oHdA7B7AHYPYPD3YKaqlvS106gIUIMlIpsSGRI5ktk9sHsAdg/A7gGk7h5YE58xxpi0ZAHKGGNMWrIAlZj7hroAacDugd0DsHsAdg8gRffA+qCMMcakJatBGWOMSUsWoIwxxqQlC1AeIrJGRCpFZJ+IfL2b97NF5Jfu+2+KSFnqS+mvBO7BV0Rkl4hsE5E/isjMoSinn/q6B579bhYRFZERN+Q4kXsgIh93fxd2isgvUl1GvyXwf2GGiLwgIlvc/w/XD0U5/SIiPxKR0yKyo4f3RUS+496fbSJyYdIL4UzzbC+c+aX2A7OALOAdYGHcPl8AfuAu3wL8cqjLPQT34P1Anrt892i8B+5+BcDLwAagYqjLPQS/B3OBLcB4d33iUJd7CO7BfcDd7vJC4NBQlzvJ92A1cCGwo4f3rweeAgS4GHgz2WWwGlSnlcA+VT2gqu3AI8DauH3WAg+5y48BV4mIpLCMfuvzHqjqC6ra7K5uAEpTXEa/JfJ7APCPwL8DraksXOEecM0AAAelSURBVIokcg/+J/A9Va0DUNXTKS6j3xK5BwoUustjgeMpLJ/vVPVlnIlke7IW+Ik6NgDjRGRKMstgAarTNOCoZ73K3dbtPqoaAuqBopSULjUSuQden8P5C2ok6fMeiMgFwHRV/UMqC5ZCifwezAPmichrIrJBRNakrHSpkcg9+HvgkyJSBTwJ/GlqipY2+vt90W++Tfk+DHVXE4ofg5/IPsNZwp9PRD4JVABX+Fqi1Ov1HohIBvAt4I5UFWgIJPJ7EMRp5nsfTi36FRFZrKpnfS5bqiRyD24FHlTV/xKRS4Cfuvcg4n/x0oLv34dWg+pUBUz3rJfy3ip7bB8RCeJU63urAg83idwDRORq4K+BG1W1LUVlS5W+7kEBsBh4UUQO4bS9rxthAyUS/b/wO1XtUNWDQCVOwBopErkHnwMeBVDVN4AcnCSqo0VC3xeDYQGq00ZgroiUi0gWziCIdXH7rANud5dvBp5Xt7dwhOjzHrjNW/fiBKeR1u8AfdwDVa1X1WJVLVPVMpx+uBtVddPQFNcXifxf+C3OgBlEpBinye9ASkvpr0TuwRHgKgARWYAToKpTWsqhtQ74tDua72KgXlVPJPMC1sTnUtWQiHwRWI8zgudHqrpTRO4BNqnqOuCHONX4fTg1p1uGrsTJl+A9+A9gDPArd3zIEVW9ccgKnWQJ3oMRLcF7sB74gIjsAsLAV1W1ZuhKnVwJ3oO/AO4XkS/jNG3dMZL+YBWRh3GacIvdfra/AzIBVPUHOP1u1wP7gGbgM0kvwwi6n8YYY0YQa+IzxhiTlixAGWOMSUsWoIwxxqQlC1DGGGPSkgUoY4wxackClEl7IhIWka0iskNEfiUief08/lw/939QRG7uZnuFiHzHXb5DRL7rLt8lIp/2bJ/an+v1Uo7L3UzhW0UkN+69L4nIbhH5uYjc2FvW9R7OHfuMIvKAiCxMRpkTvPYh99kpY3plz0GZ4aBFVZcBiMjPgbuA/46+6SbsFb9TzLgP477ngVz3mZCoO4AdJOeJ+tuA/1TVH3fz3heA69wsDvDeh0gTpqqfH+ixxvjJalBmuHkFmCMiZW4N4vvA28B0EblVRLa7Na1veg8Skf8SkbfFmcOqxN32P0Vko4i8IyKPx9XMrhaRV0Rkj4jc4O7/PhF5T4JYEfl7Eflfbo2kAvi5W+v5oIj8xrPfNSLy626Ov0qcOYW2izMHT7aIfB74OPANNyh79/8BzjQQ60Tky3G1uQfFmaPndRE54KkliYh8V5z5m54AJnrO92I0VZOInBORf3bvyQYRmeRun+2ubxSRe3qqlYrIb0Vks1vzu7Pbn6DjqyLylvua4yl7rOYavYaIZIjI991z/kFEnuyuhmtGHgtQZtgQJ//hdcB2d9N8nHT/FwAdwDeBK4FlwAoR+bC7Xz7wtqpeCLyE80Q8wK9VdYWqng/sxsmtFlWGkwj3g8APRCSnr/Kp6mM4Nazb3Brfk8CCaEDEedK+S23IPe+DwCdUdQlOq8bdqvoATq3oq6p6W9x17sKpob1fVb/VTVGmAJcBNwD/5m77CM79WoIzVcalPXyMfGCDe09edvcF+DbwbVVdQe+1w8+q6nKcQP0lEekp23+Dqq4Evgv8n17OB3ATzs9jCfB54JI+9jcjhAUoMxzkishWnC//IzgppwAOu/PQAKwAXlTVancqlJ/jTLgGEAF+6S7/DOfLG2CxW0vajtOctshzzUdVNaKqe3FyzJ3X30K7aW9+ijMlwzicL9b46UnmAwdVdY+7/pCn3AP1W7fsu4BJ7rbVwMOqGlbV48DzPRzbDkRriZtxAgNu2X/lLvc2e+6XROQdnByF0+k5gezDnn/7CjiXAb9yP9NJ4IU+9jcjhPVBmeEg1gcV5XQ70eTd1I/zRfN7PQh8WFXfEZE7cPKOxe/T03qifgz8Hmdiw1+5wdPLjwkvvRnmvedP5DN0ePLJhenHd4SIvA+4GrhEVZtF5EWcBKrd0W6WQ7h/NLv9ilnRUydaBjOyWA3KjBRvAleISLGIBHDm6nnJfS8DJ/s8wP8AXnWXC4ATIpKJU4Py+pjb9zEbp7+nMsFyNLrnBcCtrRwH/gYnIMZ7FyiL9sMAn/KUO5leBm4RkYA4s56+v5/HbwA+6i73lCR5LFDnBqfzcKYi6cknPP++4S4fApa7y2txE5Pi/Lw+6v48JtH1DwkzglkNyowIqnpCRP4Sp/lHgCdV9Xfu203AIhHZjDMLcvTL8W9xAtthnH6tAs8pK3ECxSTgLlVtdWttfXkQp8+qBacm0YLT3FjiNrnFl7tVRD6Dkx0+iDPNww/i90uC3+D0z20H9tD/IPjnwM9E5C+AJ3DuY7yngbtEZBvO/dvQzT5R2SLyJs4fD7e62+4HficibwF/pLOG/DjOtBY73LK/2cP1zf/f3h0aMQhEQRjeMyi6iUwRWLBpJS1g0wVVxNAKAiZyETyPgnmQ/9M3c3Lnzdu5uxleMwcOFg270fZn93BS0XD82XYppZXU2W5OvL+2PUfp4ivpGfso3BgTFHCgmNoWbX8HXdlDUh+7oUnS6+T7hyiaVJLehNN/YIICAKRESQIAkBIBBQBIiYACAKREQAEAUiKgAAAprXcn+xVU+ezpAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p1_marginal = suite.Marginal(1, label='p1')\n", "p2_marginal = suite.Marginal(2, label='p2')\n", "\n", "print('post mean p1', p1_marginal.Mean())\n", "print('MAP p1', p1_marginal.MaximumLikelihood())\n", "\n", "print('post mean p2', p2_marginal.Mean())\n", "print('MAP p2', p2_marginal.MaximumLikelihood())\n", "\n", "thinkplot.Pdf(p1_marginal)\n", "thinkplot.Pdf(p2_marginal)\n", "\n", "thinkplot.decorate(xlabel='Probability of finding a bug',\n", " ylabel='PMF')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }