{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Think Bayes\n", "\n", "This notebook presents code and exercises from Think Bayes, second edition.\n", "\n", "Copyright 2016 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 math\n", "import numpy as np\n", "\n", "from thinkbayes2 import Hist, Pmf, Suite, Beta\n", "import thinkplot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The Alien Blaster problem\n", "\n", "In preparation for an alien invasion, the Earth Defense League (EDL) has been working on new missiles to shoot down space invaders. Of course, some missile designs are better than others; let's assume that each design has some probability of hitting an alien ship, x.\n", "\n", "Based on previous tests, the distribution of x in the population of designs is well described by a Beta distribution with parameters 5, 10.\n", "\n", "Now suppose the new ultra-secret Alien Blaster 9000 is being tested. In a press conference, an EDL general reports that the new design has been tested twice, taking two shots during each test. The results of the test are confidential, so the general won't say how many targets were hit, but they report: \"The same number of targets were hit in the two tests, so we have reason to think this new design is consistent.\"\n", "\n", "Is this data good or bad; that is, does it increase or decrease your estimate of x for the Alien Blaster 9000?" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.3333333333333333" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8XdWZ7//Po+reZLkXufcuN4xpphgCGAhJgBQywwyTTJi0O/nd5Hcn+SXczMxlZjJ5TW4qCQmEFCAQBgNOTLEx2Bj33uUuFyw3uclWe35/nO2tI0XV1tE+kr7v10sv9tp77a3nbGQ9WmuvvZa5OyIiIskmJeoAREREqqMEJSIiSUkJSkREkpISlIiIJCUlKBERSUpKUCIikpSUoEREJCkpQYmISFJSghIRkaSUFnUATaF79+6ek5MTdRgiIgKsWbPmuLtn11WvVSSonJwcVq9eHXUYIiICmNn++tRTF5+IiCQlJSgREUlKSlAiIpKUlKBERCQpKUGJiEhSUoISEZGkpAQlzU5JSRlaCVqk5WsV70FJ81VcUsrbH2xn1aZ9nDh9nuOnz3HxUglZXdozadQAJo3qz4QR/WjbJiPqUEWkkSlBSVIqLS1j0YodvPjGGk6cPv8Xx0+cPs9by7fx1vJttGuTwcP3zGTOjJGYWQTRikgiKEFJ0jlSUMi//GwBhwsKqz1uQHwH34WLxfzkuSUsW7ubzz94PT26dWySOEUksZSgJKnsO3Scx3/yOoVni8J9nTu25d45kxg3vA9ZXTrQNjOdnfuPsW7rAd5bk0fBqbMAbNyZz5f/9QW+8vAcpo7NiegTiEhjsdbwsDk3N9c1F1/y27r7CP/65J+4cLEYgPS0VB64Yypzrx1Dm8z0as+5VFzCcwtW8+riDWGrKjU1hf/nkdvIHTOwiSIXkYYwszXunltXPY3ik6SwdfcRHv/xa2Fyatcmg29/4S7umTOxxuQEkJmRzsP3zORfvnIvPbM6AVBWVs6/PbWQNVvqNR+liCQpJSiJ3OmzF/jer96kpLQMiHXpffdL8xg5uFe9rzE8pyeP/8PdlZLUE08tZN22gwmJWUQSTwlKIuXu/NevF3H67AUAOnVoyz9/6R4G9slq8LW6d+3Adx67KxwkUVZWzveefpOjx880aswi0jSUoCRSL725jo0784HY6Lwvffomemd3vuLrZXfryHf+4W66d+0AQNHFYr739JuUlJQ1Rrgi0oSUoCQyW/IO89zrK8PyfbdMZuLI/ld93R7dOvK1v7qV1NTYj/eegwX8ev7yq76uiDQtJSiJRElJGT/+/TvhyLtRg3vzidvrHNRTb0MH9uDheTPD8oJ3N7N8/Z5Gu76IJJ4SlETitSUbw2dD7dtm8pWH54QtnsZyx3VjmT5+UFj+8e/f4WThX85KISLJSQlKmtzpsxd48Y21YfkTt+eS1aVDo38fM+PvH7yB7K6xQRMXLhbz9H+rq0+kuVCCkib321dXcvFSCQD9enbltlmjE/a9OrTL5AsP3RCWl63NY8OO/IR9PxFpPEpQ0qT25h9n8YrtYfmz915DWlpqQr/nuOF9uXbK0LD8iz+8p1F9Is2AEpQ0GXfnl39cFg6MmDJ6IJNGXf2ovfr47D3XhEtyHC4o5OW31zXJ9xWRK5fQBGVmc81sh5nlmdnXqzmeaWbPB8dXmFlOsH+ama0PvjaY2b31vaYkr827DrN19xEAUlJSePjemXWc0Xi6dmrHQx+ZGpZfenOdXuAVSXIJS1Bmlgr8CLgdGA08aGZVHzY8Apxy96HA94Engv2bgVx3nwjMBX5mZmn1vKYkqZffqmi1zJkxgr49ujTp95977RgG988GYutN/X7ByjrOEJEoJbIFNQ3Ic/c97l4MPAfMq1JnHvBMsP0iMMfMzN0vuHtpsL8NFcv/1OeakoT25h8PBycYMO+miU0eQ0pKCo/cNyssL1uTx/7DJ5o8DhGpn0QmqL5A/Eyd+cG+ausECakQyAIws+lmtgXYBHwuOF6faxKc/6iZrTaz1QUFBY3wceRq/DGu9TRz0pCrms7oaowc3CtchsOJjSgUkeSUyARV3drbVRefqrGOu69w9zHAVOAbZtamntckOP9Jd89199zs7OwGhC2N7UhBIcvX7Q7L9908KcJo4KE7p4U/SGu27mf7nqORxiMi1UtkgsoH4odo9QMO11THzNKAzsDJ+Aruvg04D4yt5zUlybyyaH34V8TEkf0Z1K97pPEM7JPFtVOGheXfvraC1rBwp0hzk8gEtQoYZmaDzCwDeACYX6XOfODhYPt+YJG7e3BOGoCZDQRGAPvqeU1JIqfOXGDxyp1h+d6bm/7ZU3U+cXsuKSmxH/+tu49o3SiRJJSwBBU8M3oMWAhsA15w9y1m9riZ3R1UewrIMrM84KvA5WHj1wIbzGw98DLw9+5+vKZrJuozyNV7a/k2SoOFCIcO6MGYoX0ijiimd3Znbpk5Kiz/YeEataJEkkxaIi/u7guABVX2fStu+yLwsWrOexZ4tr7XlORUXl7O28srZo248/pxmFX3GDEa9982mbc+2EZZWTk7933I1t1HkiaBiohmkpAEWr89n4JTZ4HYnHjTJwyq44ym1a1ze26cNiIsx7+nJSLRU4KShHlr+bZw+8ZpI8hIT2iD/YrcM2diOKJv3baD7M0/Hmk8IlJBCUoS4mTheVZt2heWb75mVM2VI9Q7uzMzJg4Jyy+/vT7CaEQknhKUJMSiFTsoDwYdjB7Sm349u0YcUc3uixtZ+P7aPI4UFEYYjYhcpgQljc7deev9iu69W5K09XTZ4P7ZTBjRD4i99f3KIrWiRJKBEpQ0uqqDI2ZMGBxxRHW775aK2S3eWbmTM+eKIoxGREAJShJg8cod4fYNU5NzcERVY4b2CWc6Lykt46244fEiEg0lKGlUFy+VsHLj3rB804wRtdROHmbGR64bG5b/vHQzZWXlEUYkIkpQ0qhWbtpLSTBzRP9eXRnYJyviiOpv1qShdOrQFoATp8+zYtPeOs4QkURSgpJGtXRNxazl8ROyNgfp6ancOqti/cvXl2yKMBoRUYKSRnPmXBHrtldMujp7ytAIo7kyt80aHU4iu33PUfYc1FpiIlFRgpJGs3z9HsrLY89thuf0pGdWp4gjarhundszc2LFqMPX390cYTQirZsSlDSa99bkhdvNsfV02Z3Xjwu331uzi8KzGnIuEgUlKGkUBSfPsm3PEQBSzLhm0pA6zkhewwb2YEgw5LysrJx3Vu2s4wwRSQQlKGkUS9dWtJ7GDe9Ll47tIozm6pgZt11bMVjizfe3aq0okQgoQUmjWLq2YvTe7GY2eq86syYNpW2bDACOFBSyedfhiCMSaX2UoOSqfXjiDPsOxZapSE1NYdr4nGgDagRtMtO5Prci0S5ctjXCaERaJyUouWorN+4Lt8cP70v7tpnRBdOI4t+JWrlprwZLiDQxJSi5aivipjaaPj65Vs29GgP7ZDE8pycQGyzx9gean0+kKSlByVUpPFvE9mD0ngFTx+VEGk9ju/WailbUW8u3abCESBNSgpKrsmrzPi7/yh4xuFezHr1XnVmTh9AuGCzx4YkzbNx5KOKIRFoPJSi5KvHdezPGJ/+6Tw2VkZ7GDdOGh2V184k0HSUouWIXiorZsCM/LE+f0HKeP8WbM2NkuL1i417Onr8YYTQirUdCE5SZzTWzHWaWZ2Zfr+Z4ppk9HxxfYWY5wf5bzGyNmW0K/ntT3DnvBNdcH3z1SORnkJqt3XYgXDMpp293enTrGHFEiZHTt3u4mGFpaVmll5JFJHESlqDMLBX4EXA7MBp40MxGV6n2CHDK3YcC3weeCPYfB+5y93HAw8CzVc77pLtPDL6OJeozSO0qj97LiS6QJjBnekUratGKHbXUFJHGksgW1DQgz933uHsx8Bwwr0qdecAzwfaLwBwzM3df5+6XX93fArQxs5bxck0LUVJSxtqtB8JySxpeXp1rpwwlLS0VgD0HC8IXk0UkcRKZoPoCB+PK+cG+auu4eylQCFRdgvWjwDp3vxS371dB9943zcyq++Zm9qiZrTaz1QUFWtOnsW3dc4SLl0oA6JnViQG9u0UcUWJ1aJdZKQmrFSWSeIlMUNUljqovkdRax8zGEOv2+7u4458Muv5mB1+fru6bu/uT7p7r7rnZ2dkNClzqtnZLRetpypgB1PB3QosSP1hiyaqdlJSURRiNSMuXyASVD/SPK/cDqs64GdYxszSgM3AyKPcDXgY+4+7hTKTufij471ngd8S6EqWJrdm6P9yeMmZghJE0nfHD+9K9awcAzl24xKot+6INSKSFS2SCWgUMM7NBZpYBPADMr1JnPrFBEAD3A4vc3c2sC/A68A13X3a5spmlmVn3YDsduBPQkqdN7PCx0xwpKAQgMyOdMUP6RBxR0zAzbpg2IiwvVjefSEIlLEEFz5QeAxYC24AX3H2LmT1uZncH1Z4CsswsD/gqcHko+mPAUOCbVYaTZwILzWwjsB44BPw8UZ9Bqrcmrntv/PC+pKenRhhN07ppekWCWr/tIKfOXIgwGpGWLS2RF3f3BcCCKvu+Fbd9EfhYNed9F/huDZed0pgxSsPFj96bPHpAhJE0vZ5ZnRg9pDdbdx+h3J331uzi7hsnRB2WSIukmSSkQYouFrNld8WjxNaWoIBKUx+9s1LLwYskihKUNMiGHfnh7BED+2SFgwZak5kThpAevBO1//AJvRMlkiBKUNIg8c+fclvJ6L2q2rXNYFrcO1FqRYkkhhKU1Ju7s25b633+FO/GuNF8767ZFbYqRaTxKEFJve07dCIctdahXSbDc1rvPL0TRvSla6fY2leFZ4tYv/1gHWeISEMpQUm9rY1rPU0Y2Z+UlNb745OSksJ1ucPC8mJ184k0utb7G0YabMP2irWfJo/qX0vN1iH+pd1Vm/dxvuhSLbVFpKGUoKReLl4qYfveo2F5/Ih+EUaTHAb07kZO3+5AbJ2oDzbsiTgikZZFCUrqZXPe4XAgwIDe3ejWuX3EESWHG6ZWvBO1ZNWuCCMRaXmUoKReNsQNApg4Ut17l82aPCSckn9L3mEKTp6NNB6RlkQJSupl/ba4BKXnT6FundtX6u58b42WgxdpLEpQUqdjJ89yOJi9PD0tlVGDe0UcUXK5Pq6b793VO3GvuuyZiFwJJSipU3z33thhfchIT+gcw83O9PGDwnty8Ogp9h06EXFEIi2DEpTUKb57b8IIde9V1SYzvdJy8O+u1mAJkcagBCW1KisrZ+POQ2FZz5+qF9/N996aXZSXa+ojkaulBCW1yjtwjAsXi4HYgIB+PbtEHFFyGj+8L507tgXg1JkLbNp1uI4zRKQuSlBSq/Vxs0dMHNkfM6ulduuVmprC7MkVUx+pm0/k6ilBSa027qxIUBM0e0StZk8ZGm5/sGEPl4pLIoxGpPlTgpIaXbxUws59x8LyuOF9I4wm+Q0ZkE2f7M5A7N6t2rw/4ohEmjclKKnRlrzD4cP+gX2ywmcsUj0zY3bcDOfvqZtP5KooQUmNNu6oGL03Xq2nepk9pSJBrd12kDPniiKMRqR5U4KSGsU/f9Ls5fXTO7szwwbGFnIsLy/n/XWa4VzkSilBSbUKzxZx4MhJILY43+ghvSOOqPmIX8jw3TXq5hO5UglNUGY218x2mFmemX29muOZZvZ8cHyFmeUE+28xszVmtin4701x50wJ9ueZ2Q9M454TYlPcy7kjcnrSJjM9wmial1mThpIS/Fju2HuUD0+ciTgikeYpYQnKzFKBHwG3A6OBB81sdJVqjwCn3H0o8H3giWD/ceAudx8HPAw8G3fOT4BHgWHB19xEfYbWbMOOiu49jd5rmM4d21aacUMznItcmUS2oKYBee6+x92LgeeAeVXqzAOeCbZfBOaYmbn7One//Cr+FqBN0NrqDXRy9+UemzL618A9CfwMrZK7V2pB6f2nhrtuSuXRfJrhXKThEpmg+gIH48r5wb5q67h7KVAIZFWp81FgnbtfCurnxx2r7ppylY4eP0PBqdjCe5kZ6QwdkB1xRM3P1HE5ZGbEukXzP9QM5yJXIpEJqrpnQ1X/jKy1jpmNIdbt93cNuOblcx81s9VmtrqgoKAe4cpl8a2nsUP7kJaWGmE0zVNshvOcsLxk1c7oghFpphKZoPKB+Kmv+wFVZ9AM65hZGtAZOBmU+wEvA59x991x9eP7m6q7JgDu/qS757p7bna2WgANET97uZ4/Xbn4d6KWrs3TDOciDZTIBLUKGGZmg8wsA3gAmF+lznxigyAA7gcWububWRfgdeAb7r7scmV3PwKcNbMZwei9zwCvJPAztDruzuZdcS/ojlCCulITRvSjU4eKGc43a4ZzkQZJWIIKnik9BiwEtgEvuPsWM3vczO4Oqj0FZJlZHvBV4PJQ9MeAocA3zWx98NUjOPZ54BdAHrAb+FOiPkNrdODISc6evwhApw5tGdC7W8QRNV+pqSlcO3lIWF6iqY9EGiSha3e7+wJgQZV934rbvgh8rJrzvgt8t4ZrrgbGNm6kcln89EZjhvbR8hpX6brcYSx4dzMQm+H87z4+O1weXkRqp5kkpJIteRXdUOOG9YkwkpZh6IAe9OreCdAM5yINpQQlobKycjbHJaixGiBx1cyM63LjloNXN59IvSlBSWhv/nGK4pZ3v7y2kVyd+Ln51m47ED7jE5HaKUFJKH54+dhhev7UWOJnOC8rK2fZ2t11nCEioAQlcSo/f1L3XmPSDOciDacEJQCUlpaxdfeRsKznT42r6gznR49rhnORuihBCQC79h+juKQUgJ5ZnejRrWPEEbUsnTu2ZdKoAWH53dWa+kikLkpQAsCmXZXff5LGF9/NpxnOReqmBCUAlabhGa/uvYSYOm5gOMP54YJCdh/QJMYitVGCEopLStm+92hYHqMXdBMiMyOdmRMHh+Ul6uYTqZUSlLBj74eUlcVm2u7bowvdOrePOKKWK76bb+na3ZSWlkUYjUhyqzVBmdnTcdsP11JVmrFK6z9peHlCjRvWh66d2gFw5lwR63fk13GGSOtVVwtqQtz2lxIZiESn8vRG6t5LpJSUlEqtqHdWqptPpCZ1JSgNM2rhLl4qYdf+Y2F5rEbwJdz1U0eE26s27+N80aUIoxFJXnXN+9/PzH5AbKn1y9shd/9iwiKTJrFtz9FwpdcBvbuFC+xJ4gzs042cvt3Zd+g4paVlLF+/h5tnjoo6LJGkU1eC+lrc9upEBiLR2LKr8vx70jRumDqcpw8dB2DJqp1KUCLVqDVBufszTRWIRGOjBkhE4topQ/n1K8spd2fr7iN8eOIMPbM6RR2WSFKpNUGZ2fzajrv73bUdl+R2vugSew7GXhY1NINEU+raqR0TR/Vn7dYDQKwV9fG5uRFHJZJc6urimwkcBH4PrCD2e0xaiG17joajYHL6dadDu8xI42ltrs8dHiaod1fv4mO3TdESJyJx6hrF1wv4f4GxwH8BtwDH3X2Juy9JdHCSWJvjuve0vEbTmzY+h7ZtMgA4UlDIjr0fRhyRSHKpNUG5e5m7/9ndHwZmAHnAO2b2D00SnSTUprj59zS9UdPLSE/jmripj95ZtSPCaESST51THZlZppndB/wG+ALwA+CPiQ5MEuvs+YvsD0aRpZgxenDviCNqnW6aPjLcXrp2d7jkiYjUPdXRM8D7wGTgO+4+1d3/t7sfqu08SX5b8g6Hz58G98+mXduMSONprUYM6kmv7rHRe0UXi1m5cV+0AYkkkbpaUJ8GhhOb5mi5mZ0Jvs6amZYEbcYqL++u7r2omBk3xrWi3v5ge4TRiCSXup5Bpbh7x7ivTsFXR3ev86UNM5trZjvMLM/Mvl7N8Uwzez44vsLMcoL9WWa22MzOmdkPq5zzTnDN9cFXj4Z9ZIHKE8SOG9Evwkjk+txh4fDYTTvzOX7qXKTxiCSLurr42pjZl83sh2b2qJnVNSw9/txU4EfA7cBo4EEzG12l2iPAKXcfCnwfeCLYfxH4JvCPNVz+k+4+Mfg6VkMdqcHpsxc4ePQUAKmpKYwc1DPiiFq37G4dGTc89keCA++s0gSyIlB3F98zQC6wCbgD+F4Drj0NyHP3Pe5eDDwHzKtSZ17wPQBeBOaYmbn7eXdfSixRSSOLXz13+MCe4SqvEp0bpw8Pt99ZuUPLwYtQd4Ia7e6fcvefAfcDsxtw7b7EXvK9LD/YV20ddy8FCoGselz7V0H33jethjcbgxbfajNbXVCgpbXjbY6bf0/Dy5PD9PGDaJMZ+0NB70SJxNSVoEoubwQJpCGqSxxV/yysT52qPunu44gly9nEBnL85UXcn3T3XHfPzc7OrjPY1iS+BTV+uF7QTQaZGenMmjQkLGuwhEg9FiyMH7kHjG/AKL58oH9cuR9wuKY6wfOtzsDJ2i56eYi7u58FfkesK1Hq6fipcxwpKAQgPS2V4QP1/ClZzJlRMZpv2brdFF0sjjAakejVNYovtcrIvbQGjOJbBQwzs0FmlgE8AFSdfHY+cHkp+fuBRV5L57uZpZlZ92A7HbgT2FxHHBInfnj5yMG9SE9PjTAaiTc8pyf9enYF4FJxCcvW7Y44IpFo1TmTxJUKugQfAxYC24AX3H2LmT1uZpdnQX8KyDKzPOCrQDgU3cz2Af8JfNbM8oMRgJnAQjPbCKwHDgE/T9RnaIm0vEbyMjPmzKxoRb21fFuE0YhEr97Dxq+Euy8AFlTZ96247YvAx2o4N6eGy05prPhaG3evNEBCz5+Sz/W5w/nNqysoKytn1/5jHDx6iv69ukYdlkgkEtaCkuTz4Ymz4UugmRnpDO7XPeKIpKrOHdsydWxOWH5brShpxZSgWpFKw8uH9iYtTc+fklH8YIl3Vu2ktLQswmhEoqME1Yro+VPzMHFkP7K6tAdis86v3Lwv2oBEIqIE1Uq4e6X59/T8KXmlpKRUmkD2rffVzSetkxJUK3HgyEnOnCsCoEO7THL61mfCDonKnBkjw7fYN+zI5+hxLR4grY8SVCuxcUfc7OXD+1HDDFGSJHp068ik0QPC8lvvb40wGpFoKEG1Eurea35uuaZi8v+3V+zQYAlpdZSgWoHS0jI2xy9QqATVLEwZPSAcLHHmXBErNu2LNiCRJqYE1QrkHSjgUnFs3t/srh3DJcYluaWmpjBnxqiw/MayLRFGI9L0lKBagY0788PtccP76vlTM3LzzJGkBP+/Nu86zKFjpyOOSKTpKEG1Anr+1HxldenAlDEDw7KGnEtrogTVwl28VMKOfRWL340boQTV3Nw6q2KwxKIV2ykuaejSbCLNkxJUC7dtz1HKysoB6N+7G106tos4ImmoSaP60zMr9tzw3IVLvLdmV8QRiTQNJagWblPc8yd17zVPZsZt144Jywve3UIty6aJtBhKUC3chrgXdMeP6BdhJHI1bpo+gvRgct99h46zY++HdZwh0vwpQbVghWeL2HfoOBCb32304N4RRyRXqmP7Nlw/dXhYfv3dTRFGI9I0lKBasI07Krr3RuT0pF3bjAijkat1x3Vjw+0PNuzlZOH5CKMRSTwlqBZs3faD4faEkerea+4G9sli9JBYK7i8vJyFyzQ/n7RsSlAtlLuzYXtFC2rSyP4RRiONZe7silbUm+9vpaRE8/NJy6UE1UIdOHKS02cvALHlNQb31/LuLcH0cTl06xybn6/wbBFL1+ZFHJFI4ihBtVDrt8dPb9SPlBT9r24J0tJSmTu7Ysj5/MUbNORcWiz91mqh1m+reP40aZSeP7Ukt14zmoz0NCDWUt4YN5WVSEuiBNUCFZeUsnXPkbA8YYSeP7UkHdu3Yc6MiiXhX128IcJoRBInoQnKzOaa2Q4zyzOzr1dzPNPMng+OrzCznGB/lpktNrNzZvbDKudMMbNNwTk/ME3N/Re27j4SLm7Xr2dXunftEHFE0tg+cv24cEn4ddsOcuDIyUjjEUmEhCUoM0sFfgTcDowGHjSz0VWqPQKccvehwPeBJ4L9F4FvAv9YzaV/AjwKDAu+5jZ+9M1b/Og9DS9vmXpnd2ba+EFh+dXFGyOMRiQxEtmCmgbkufsedy8GngPmVakzD3gm2H4RmGNm5u7n3X0psUQVMrPeQCd3X+6xJ8O/Bu5J4Gdoliq9/6TpjVqsu2+cEG4vWb2TU2cuRBiNSONLZILqCxyMK+cH+6qt4+6lQCGQVcc18+PK1V0TADN71MxWm9nqgoKCBobefJ04fY6DQXdPamoKY4b2iTgiSZQRg3oydEAPAMrKylmwRNMfScuSyARV3bOhquNh61Pniuq7+5PunuvuudnZ2bVcsmVZu/VAuD16SG/aZKZHGI0kkpkxb05FK+pPS7dwvuhShBGJNK5EJqh8IH74WD/gcE11zCwN6AzU9rQ3P7hObdds1eIT1JTRA2upKS3BjPGD6JPdGYCii8X8eemWiCMSaTyJTFCrgGFmNsjMMoAHgPlV6swHHg627wcWeS1vHbr7EeCsmc0IRu99Bnil8UNvnkpKyiotrzF5zIAIo5GmkJKSwn23TA7Lr72ziUvFJRFGJNJ4EpaggmdKjwELgW3AC+6+xcweN7O7g2pPAVlmlgd8FQiHopvZPuA/gc+aWX7cCMDPA78A8oDdwJ8S9Rmam617joS/nHp17xT+ZS0t2+wpQ8NXCc6cK+Kt5dsjjkikcaQl8uLuvgBYUGXft+K2LwIfq+HcnBr2rwbGVnestVu7paJ7b/LoAegVsdYhLS2VeTdN4KmXlgHwyqL13DZrNGnBAocizZVmkmhB1mzdH25PGaPnT63JzTNH0alDWwBOnD7PktU7I45I5OopQbUQh4+d5khBIQAZ6WnhukHSOmSkp3HXDePD8ktvrAtnExFprpSgWoj40XsTRvQLJxOV1mPutWPo0C4TgA9PnOGdVWpFSfOmBNVCrKny/Elan3ZtM5h308Sw/IeFa7SgoTRrSlAtwMVLJWzZXfE6mBJU63XHdWPDZ1HHT53j7Q80ok+aLyWoFmD99oOUlZUDMKB3N81e3oq1yUznnjkVragX31hDcUlphBGJXDklqBZgxca94fbUsTnRBSJJYe61o+nSsR0Ap85c4I1lWyOOSORnwjB9AAAVcklEQVTKKEE1c6WlZazeXDG8fHrcEgzSOmVmpHPfLZPC8ktvrqPoYnGEEYlcGSWoZm7L7iNcCH75dO/agcH9u0cckSSDW64ZRVaX9kBsdon/fnt9xBGJNJwSVDO3YkNF9960cTmaPUKA2HtRD31kWlh+ZdEGThaejzAikYZTgmrG3J2VmyoSlLr3JN71U4eT0zfWoi4pLeP3r6+KOCKRhlGCasZ27T8WrqLaoV0mowZr9gipYGZ8Zt6MsLx4xXb2H65tNRuR5KIE1YytjB+9Ny6H1FT975TKJozox6RRsWXZHHh2/vJoAxJpAP1Ga6bcnQ82qntP6vbpu2eES1Gv23aw0rRYIslMCaqZOnj0VDg5bGZGOhNG9KvjDGmtBvbJ4sbpI8PyUy8t1RRI0iwoQTVTH2zYE25PHj1Ak8NKrT5113TatckA4OjxM7yyeEPEEYnUTQmqGXJ3lq3dHZZnqHtP6tC5Y1seurNi2PmLC9dw7OTZCCMSqZsSVDO0//AJ8j88BcTed8kdq8UJpW63zRpdadj5My+/H3FEIrVTgmqG3l29K9yeNj6HNpnpEUYjzUVKSgp/e/+1YfmDjXs1YEKSmhJUM+PuLF2bF5avmzIswmikuRk5uBfXTx0eln/6/BIuFGmePklOSlDNzNbdRzhxOjZlTYd2mRq9Jw328LyZdGzfBoATp8/zm1dXRByRSPWUoJqZ99ZUdO/NmjSUtLTUCKOR5qhzx7b8zUcruvoWLtvC5l2HIoxIpHpKUM1IaWkZy9dXDC+fPWVohNFIczZr8pBKa4f95LklXCouiS4gkWokNEGZ2Vwz22FmeWb29WqOZ5rZ88HxFWaWE3fsG8H+HWZ2W9z+fWa2yczWm9nqRMafbNbvyOfchUtAbGmNkYN7RRyRNFdmxqMfn13p3Sh19UmySViCMrNU4EfA7cBo4EEzG12l2iPAKXcfCnwfeCI4dzTwADAGmAv8OLjeZTe6+0R3z01U/MkofvTe7MlDtbSGXJVundvz2XtnhuUF727WqD5JKolsQU0D8tx9j7sXA88B86rUmQc8E2y/CMyx2G/decBz7n7J3fcCecH1Wq2z5y9WWtp9dq5G78nVu2n6SHLHVLxH939/u5jTZy9EGJFIhUQmqL7AwbhyfrCv2jruXgoUAll1nOvAG2a2xswerembm9mjZrbazFYXFBRc1QdJBotX7qC0NDZ/2pD+2QzskxVxRNISmBl//+ANdOnYDoitvvvD3y7G3SOOTCSxCaq6/qeqP/U11ant3FnuPplY1+EXzOy66r65uz/p7rnunpudnV3fmJOSu/Pmsq1h+dZZVXtKRa5c545t+eKnbwrL67Yd5PUlmyKMSCQmkQkqH+gfV+4HHK6pjpmlAZ2Bk7Wd6+6X/3sMeJlW0PW3Je8wh4OZy9tkpnPtZI3ek8Y1YUQ/7r5xQlh+5pUP2Lb7SIQRiSQ2Qa0ChpnZIDPLIDboYX6VOvOBh4Pt+4FFHutbmA88EIzyGwQMA1aaWXsz6whgZu2BW4HNCfwMSeGN97eF29fnDtfURpIQD31kGoP7x3obysvL+Y9fvcnJwvMRRyWtWcISVPBM6TFgIbANeMHdt5jZ42Z2d1DtKSDLzPKArwJfD87dArwAbAX+DHzB3cuAnsBSM9sArARed/c/J+ozJIPCs0WVlta4ddaoCKORliw9PZWv/fWt4SwTp89e4N9/+YbWjpLIWGt4GJqbm+urVzfPV6ZeWbSBX78SW6Z72MAe/J+v3hdxRNLSbdp5iO/86NXwoe8t14zi7z5+nV5rkEZjZmvq85qQZpJIYu7Om+/HDY64RoMjJPHGDe/Lp+dVvB/15vvbeGWRFjiUpqcElcRWbtoXLuvetk0G10waEnFE0lrcfeN4ZsUNxnl2/gcsXZNXyxkijU8JKkm5Oy+9sTYs3zZrtAZHSJMxMx576AZGD+kd7vvBbxexaacmlZWmowSVpDbtPMTug7EXjNPSUrnzhvERRyStTUZ6Gv/zb+bSr2dXAMrKynniqYXsPtD8X3yX5kEJKkn98a114fZN00fQtVO7CKOR1qpDu0z+6XN3hD9/RReL+c6PX2PPQSUpSTwlqCSUt/9Y2JWSYsa8myZGHJG0ZtndOvLNz3+E9m0zAThfdIlv/+g19uYfjzgyaemUoJJQfOtp1uSh9OreKcJoRGBgnyy+/YU7qySpV8nbfyziyKQlU4JKMvsPn2Rl3Kzl9948KcJoRCoM7p/Nt79wZ7iG1LkLl/jWD1/VEh2SMEpQScTd+dXLy8IXJHPHDGRgn26RxiQSL5ak7gpbUpeKS/jXn/+ZxSt2RByZtERKUElk1eb94bMnAx66c3q0AYlUY8iAbP75y/fQvWsHIDZv3w9/t5jfL1ilZTqkUSlBJYmSkjKefvn9sHzbtWPUepKk1b9XV/71K/cyoHfFz+iLC9fwzz9bwNnzFyOMTFoSJagk8eo7G/nwxBkA2rfN5IE7pkYckUjtunVuzz9/6R7GD+8X7lu37SBf+/eXNAxdGoUSVBI4WXieF+NmjXjgjtxwRmmRZNaubQb/9Lk7uHdOxasQBafO8vXvv8xLb66lrKw8wuikuVOCipi789Pn3uVScQkQ6zq5bdaYiKMSqb/U1BQ+dfcMvvbXt4bTcZWVlfO711byv/7rvzl07HTEEUpzpQQVsdfe2cSarfvD8l/fN4vUVP1vkeZnxoTB/Ns/fpShA3qE+3btP8ZXn/gDv1+wKvwjTKS+9JswQnn7j/Hsqx+E5btuGM/4Ef1qOUMkufXt0YV/+fI9PPiRaeEfWqWlZby4cA1f/JfneX/9bo30k3pTgorIhaJi/vOZt8I++iH9s/nUXRpWLs1famoK9986mX/7H/cxJFhCHuD4qXN871dv8j+/90fWbNmvRCV10oq6ESgpKePff/lG2LXXJjOd//ja/fTO7hxxZCKNy91ZtGI7v3l1JWfOFVU6NmxgD+6ZM5Fp43JISdHfyq1JfVfUVYJqYiUlZfzHr95g9ZaK505fefhmro1bHE6kpTlfdIk//HkNf166hZLSskrHsrt2ZO7sMcyZMVKjV1sJJag4yZKgSktjLaf45HTfzZP4pLr2pJU4WXiel99axxvvb6O0SqJKTU1h8qgBXD91OLljBpKenhpRlJJoSlBxkiFBnTpzgf/7m0Vs2JEf7rt3zkQ+edd0zCzCyESa3onT51i4dCsLl23h3IVLf3G8TWY6k0YNYOrYgUwePUAtqxZGCSpO1Alq+fo9/PT5JZX+Id4zZyKfUnKSVq64pJSla/JYuGwreQeqX7rDgIF9uzN2aB/GDOvDiJyedO7YtmkDlUalBBUnqgS1N/84f3xrHe+v211p//23TuaBO6YqOYnEOXTsNO+u3sV7q3eF037VJKtLe4YN6MHAvln069WV/r260bt7J9LS1C3YHCRFgjKzucB/AanAL9z9/1Q5ngn8GpgCnAA+4e77gmPfAB4ByoAvuvvC+lyzOk2ZoE6ducDmnYdYuGwr2/YcqXQsq0t7vvDQjUzQu04iNXJ3Dh49xarN+1i1aR95+49Rn99SBnTv2pFe2Z3I7tqRrK7tyercnm6d29O5Q1s6dWxL5w5tyMxIT/RHkDpEnqDMLBXYCdwC5AOrgAfdfWtcnb8Hxrv758zsAeBed/+EmY0Gfg9MA/oAbwHDg9NqvWZ1rjRBHSko5GTh+Ur7ysudktIySsvKKS4u5dSZCxSevcCxU+fI23+sxr/8Zk8Zxt9+7NpwHR0RqZ8LRcVs23OELXmH2bbnKHvzj//FSMCGSEtLpUPbTNq3zaBd2wzaZKbTNjOdzIx0MjPSyEhPJTM9jdS0VNKDr7TUlNhXWgoplkJqqpFiKViKkWJGamoKZmAWK0Ns+/K+y30l1fWaxO9qTr0qfXt2oUvHdld0bn0TVNoVXb1+pgF57r4nCOg5YB4Qn0zmAd8Otl8Efmix/0PzgOfc/RKw18zygutRj2s2mj+/t4XXlmy84vNTUlKYOXEwd14/juE5PRsxMpHWo13bDKaMGciUMQOB2Dx/B4+eZPfBAg4eOUX+h6c4cOQkJ06fr+NKMaWlZZw+e4HTZy8kMuwW70ufvonrcofXXfEqJDJB9QUOxpXzgarjqcM67l5qZoVAVrD/gyrn9g2267omAGb2KPAowIABA67sE1yB9LRUhg3swbjhfbl55ii6dW7fZN9bpDVITU0hp293cvp2r7S/uKSUD0+c5cMTZzh+8hwnC89zovA8J0+f58z5i5w5V0ThuSLNsN6MJDJBVddWrdqfWFOdmvZX97p5tX2U7v4k8CTEuvhqDrNmPbt3ZNTg3rFA4yIKm/1pqXTp2JYundrRrXM7+vfqxuB+3fWgViQCGelp9O/Vlf69utZYx90pLinlfFEx5y5couhiMUWXSii6WMKl4hKKS8q4VFxKcWlprCu/pIyS0nJKy8ooKy+ntKycsrJyyt0pK3O8PLZdXu64g+N4XBmg3MurxFA5nurjvOrbkXCdr7B7ryESmaDygf5x5X7A4Rrq5JtZGtAZOFnHuXVds9Hccd047rhuXKIuLyJNzMyCZ03p6t1oBhI5AdYqYJiZDTKzDOABYH6VOvOBh4Pt+4FFHvuTYj7wgJllmtkgYBiwsp7XFBGRFiBhLajgmdJjwEJiQ8J/6e5bzOxxYLW7zweeAp4NBkGcJJZwCOq9QGzwQynwBXcvA6jumon6DCIiEh29qCsiIk2qvsPMNce9iIgkJSUoERFJSkpQIiKSlJSgREQkKSlBiYhIUmoVo/jMrADYX2fFmnUHjjdSOM2V7oHuAegegO4BXP09GOju2XVVahUJ6mqZ2er6DIlsyXQPdA9A9wB0D6Dp7oG6+EREJCkpQYmISFJSgqqfJ6MOIAnoHugegO4B6B5AE90DPYMSEZGkpBaUiIgkJSUoERFJSkpQccxsrpntMLM8M/t6Ncczzez54PgKM8tp+igTqx734KtmttXMNprZ22Y2MIo4E6muexBX734zczNrcUOO63MPzOzjwc/CFjP7XVPHmGj1+LcwwMwWm9m64N/DHVHEmShm9kszO2Zmm2s4bmb2g+D+bDSzyY0ehLvrK/YcLhXYDQwGMoANwOgqdf4e+Gmw/QDwfNRxR3APbgTaBdufb433IKjXEXgX+ADIjTruCH4OhgHrgK5BuUfUcUdwD54EPh9sjwb2RR13I9+D64DJwOYajt8B/AkwYAaworFjUAuqwjQgz933uHsx8Bwwr0qdecAzwfaLwBwzsyaMMdHqvAfuvtjdLwTFD4B+TRxjotXn5wDgfwP/BlxsyuCaSH3uwd8CP3L3UwDufqyJY0y0+twDBzoF252Bw00YX8K5+7vEFpKtyTzg1x7zAdDFzHo3ZgxKUBX6AgfjyvnBvmrruHspUAhkNUl0TaM+9yDeI8T+gmpJ6rwHZjYJ6O/urzVlYE2oPj8Hw4HhZrbMzD4ws7lNFl3TqM89+DbwKTPLBxYA/9A0oSWNhv6+aLCELfneDFXXEqo6Br8+dZqzen8+M/sUkAtcn9CIml6t98DMUoDvA59tqoAiUJ+fgzRi3Xw3EGtFv2dmY939dIJjayr1uQcPAk+7+/fMbCbwbHAPyhMfXlJI+O9DtaAq5AP948r9+Msme1jHzNKINetrawI3N/W5B5jZzcD/Au5290tNFFtTqesedATGAu+Y2T5ife/zW9hAifr+W3jF3UvcfS+wg1jCainqcw8eAV4AcPflQBtik6i2FvX6fXE1lKAqrAKGmdkgM8sgNghifpU684GHg+37gUUePC1sIeq8B0H31s+IJaeW9twB6rgH7l7o7t3dPcfdc4g9h7vb3VdHE25C1Offwn8TGzCDmXUn1uW3p0mjTKz63IMDwBwAMxtFLEEVNGmU0ZoPfCYYzTcDKHT3I435DdTFF3D3UjN7DFhIbATPL919i5k9Dqx29/nAU8Sa8XnEWk4PRBdx46vnPfh3oAPwh2B8yAF3vzuyoBtZPe9Bi1bPe7AQuNXMtgJlwNfc/UR0UTeuet6D/wH83My+Qqxr67Mt6Q9WM/s9sS7c7sFztv8PSAdw958Se+52B5AHXAD+qtFjaEH3U0REWhB18YmISFJSghIRkaSkBCUiIklJCUpERJKSEpSIiCQlJSiRgJmVmdl6M9tsZn8ws3YNPP9cA+s/bWb3V7M/18x+EGx/1sx+GGx/zsw+E7e/T0O+Xy1xzA5mJF9vZm3j9ufUMpP148EL25jZlxt6r0TqQwlKpEKRu09097FAMfC5+IPBC4kJ/zfj7qvd/YvV7P+pu/86KH4WaJQEBXwS+I/gsxfVM8ZvuftbQfHLgBKUNDolKJHqvQcMDVoR28zsx8BaoL+ZPWhmm4KW1hPxJ5nZ98xsrcXWysoO9v2tma0ysw1m9lKV1sbNZvaeme00szuD+jeY2V9MRGtm3zazfwxaXbnAb4NWz0fM7OW4ereY2R+rOX+OxdYu2mSxtX4yzexvgI8D3zKz31ZzH1LN7OdBC+uNyy2sy60/M/sisUS52MwWN+gOi9RBCUqkimCexduBTcGuEcSWFZgElABPADcBE4GpZnZPUK89sNbdJwNLiL15D/BHd5/q7hOAbcTmcLssh9iEux8BfmpmbeqKz91fBFYDn3T3icTe6B91OSESe6P/V1U+UxvgaeAT7j6O2Cwyn3f3XxCbsuZr7v7Jar7dMGLLaowBTgMfrRLLD4jNv3aju99YV+wiDaEEJVKhrZmtJ/bL/wCxqa0A9gfr3QBMBd5x94JgyZXfElvYDaAceD7Y/g1wbbA9NmglbSLWnTYm7nu+4O7l7r6L2Fx2IxsadDC9zrPEln7oAszkL5dBGQHsdfedQfmZuLhrs9fd1wfba4glVJEmobn4RCoUBS2SUDDf4Pn4XQ243uV5xJ4G7nH3DWb2WWLzm1WtU1O5vn4FvEpsAcU/BMkz3pUurBk/W30Z0LamiiKNTS0okYZZAVxvZt3NLJXYmkBLgmMpxGa5B3gIWBpsdwSOmFk6sRZUvI+ZWYqZDSG2vPiOesZxNrguAO5+mFhX2z8RS4hVbQdyzGxoUP50XNxXq1IsIo1FLSiRBnD3I2b2DWAxsVbJAnd/JTh8HhhjZmuIrbb8iWD/N4kltv3EnmvF/zLfQSxR9AQ+5+4Xg1ZbXZ4m9syqCJgZjL77LZDt7lurifuimf0VsVno04gtJ/HT+n/yWj0J/MnMjug5lDQmzWYu0kIE70utc/en6qws0gwoQYm0AEGr7TxwSwtc5VhaKSUoERFJShokISIiSUkJSkREkpISlIiIJCUlKBERSUpKUCIikpT+f7Yyfgr1p2ICAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Solution\n", "\n", "# Here's the prior\n", "\n", "prior = Beta(5, 10)\n", "thinkplot.Pdf(prior.MakePmf())\n", "thinkplot.decorate(xlabel='Probability of hit',\n", " ylabel='PMF')\n", "prior.Mean()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Solution\n", "\n", "# And here's the likelhood function\n", "\n", "from scipy.stats import binom\n", "\n", "class AlienBlaster(Suite):\n", " \n", " def Likelihood(self, data, hypo):\n", " \"\"\"Computes the likeliood of data under hypo.\n", " \n", " data: number of shots they took\n", " hypo: probability of a hit, p\n", " \"\"\"\n", " n = data\n", " x = hypo\n", " \n", " # specific version for n=2 shots\n", " likes = [x**4, (1-x)**4, (2*x*(1-x))**2]\n", "\n", " # general version for any n shots\n", " likes = [binom.pmf(k, n, x)**2 for k in range(n+1)]\n", " \n", " return np.sum(likes)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VOeZ///3rY6ERJEEoosOovdmG9vYDiYOJC4JTtzWTrxpm02c8nM2m3yz2e/mG+8m9m7W3iTeOLZT3Cu2cacYsCmiiQ4CBBJVNNEk1J7fHzOMRrIqaOaMRp/XdXFd55x5ZuaeIzG3nm7OOURERCJNjNcBiIiI1EcJSkREIpISlIiIRCQlKBERiUhKUCIiEpGUoEREJCIpQYmISERSghIRkYikBCUiIhEpzusAwiEjI8NlZ2d7HYaIiABr16495pzLbKpcu0hQ2dnZ5Obmeh2GiIgAZravOeXUxCciIhFJCUpERCKSEpSIiEQkJSgREYlISlAiIhKRlKBERCQiKUE1obyiki35BymvqPQ6FBERT50oOceewmKqq6vD8n7tYh7UpXr0mcUsW5tPZWUVP//W5xg1pJfXIYmIeGbpmp389Y1VJCclcNvsCcy9ZkxI3081qEbEx8VSWVkFwOb8gx5HIyLira27DwFwvqyclA4JIX8/JahGjBjUM3C8VQlKRNqxqqrqQIICyBnYs5HSrUMJqhHBCWrnvqPqhxKRdmtv0THKLlQA0LVTClkZaSF/TyWoRnRJS6ZnZicAKiur2LXvqMcRiYh4Y0tw7WlQD8ws5O+pBNWEnKBa1BY184lIOxXczTEiDM17EOIEZWazzWyHmeWb2YP1PJ5oZs/7H19lZtn+6+lmttjMzprZo3Wec7uZbTKzPDN7x8wyQvkZRgzqEThWghKR9qi6unb/04jBbTxBmVks8BhwI5AD3G5mOXWK3QecdM4NAh4BHvJfLwN+CvygzmvGAf8FXOOcGw3kAd8O1WeA2h2BO/YeoaKiKpRvJyIScfYdPMH5snIAOqfWdH2EWihrUJOBfOfcHudcOfAcMK9OmXnA0/7jl4BZZmbOuXPOueX4ElUw8/9LMV8DaBoQ0mpNRpeOgc7Aisoq8verH0pE2pfNu2q+ZsPV/wShTVC9gMKg8yL/tXrLOOcqgRIgvaEXdM5VAN8ANuFLTDnAE/WVNbP7zSzXzHKLi4sv9TMAtWtRwR2FIiLtwdbd4e9/gtAmqPpSrLuEMjWFzeLxJahxQE98TXw/rq+sc+5x59xE59zEzMwmdxZuVK1+qF3qhxKR9sM550n/E4Q2QRUBfYLOe/Pp5rhAGX//UifgRCOvORbAObfbOeeAF4DprRVwQ4JH8m3feziwuoSISLTbf+gEZ89fACCtYwd6d+8ctvcOZYJaAww2s/5mlgDMBxbUKbMAuNt/fCuwyJ94GnIAyDGzi1Wi64FtrRhzvbp1TSWzSyrgWzx2d+HlNRmKiLQVtfqfBoav/wlCmKD8fUrfBt7Fl0RecM5tMbNfmNlcf7EngHQzywceAAJD0c2sAHgYuMfMiswsxzl3EPgX4CMzy8NXo/plqD5DsOBq7ZZ89UOJSPtQa/5TUHdHOIR0NXPn3EJgYZ1rPws6LgNua+C52Q1c/z3w+9aLsnlGDurJktU7AN98qJuvHxfuEEREwso5V2tgWPDyb+GglSSaKbgGtXX3IfVDiUjU23/oBGfO+Wb7pHXsQN8eXcP6/kpQzdStayrd033zocorKtmpdflEJMrl7TgQOB4xqGdY+59ACapFgjcs3LTzQCMlRUTavuDvudEebNiqBNUCowYrQYlI+1BVVc2WoAm6XuworgTVAiOHBO8PdSSwN4qISLTZXVgc+I4LXvItnJSgWqBzajJ9/J2EVVXVbNtz2OOIRERCIy+olWjk4F5h738CJagWC26H3bxLzXwiEp027SwKHHvR/wRKUC02MqgfKk/9UCIShcorKtm+90jg3Iv+J1CCarERg3oEVrjdW1gcWKNKRCRa7Nh7JDDXs1e3znTtlOJJHEpQLZTSIZGBfbsBvmXX1cwnItFmU53+J68oQV2C2v1Q2n5DRKJLXlD/k1fNe6AEdUlGasKuiESp86Xl7N7v27HBgJFh3P+pLiWoSzB8QBZxcbEAFB05ybGTZz2OSESkdWzOP0i1f9ej7N4ZpKYkeRaLEtQlSIiPI2dAzbLzeTuKGiktItJ2bNhWGDgeO7S3h5EoQV2yMcNqfnAblKBEJEps3FGToMYM69NIydBTgrpEY4MS1MbthTS+EbCISOQ7fOw0h4+dBnwtRcP6Z3kajxLUJerXM51OqR0AOHv+AnsKj3kckYjI5Qnurhg5uCfx8bEeRqMEdcnMjDFDg5v5ChspLSIS+TZuD2re87j/CZSgLsvYoPbZjdvVDyUibVdVVXWt5du87n8CJajLEjyBbfvew9p+Q0TarPz9RzlfVg5AeucUenfv7HFESlCXpWunFPoGbb+xOV+rSohI27QhqBVozNA+nmyvUZcS1GWq3cynfigRaZs2Bg2QCJ5G4yUlqMs0drj6oUSkbTtXeoFdBb7tNQzv9n+qSwnqMg0fkEW8f9mjA0dPUXzijMcRiYi0zKadBwLLGw3ok0laxw4eR+SjBHWZEuLjyBlYs+zR+m1q5hORtiX4e2tsBIzeu0gJqhWMz+kbOF63db+HkYiItIxzrtb3VvD3mdeUoFpB8A80b+cBKiqqPIxGRKT59h86wYmSc4BvQ9bB/bp5HFGNkCYoM5ttZjvMLN/MHqzn8UQze97/+Cozy/ZfTzezxWZ21swerfOcBDN73Mx2mtl2M7sllJ+hOXp260xWRhoAF8or2LbnkMcRiYg0z9otNbWnMcN6ExsbOfWWkEViZrHAY8CNQA5wu5nl1Cl2H3DSOTcIeAR4yH+9DPgp8IN6XvonwFHn3BD/6y4NQfgtpmY+EWmL1m+r+b6aEEHNexDaGtRkIN85t8c5Vw48B8yrU2Ye8LT/+CVglpmZc+6cc245vkRV173A/wNwzlU75yJildZxw5WgRKRtOVd6ge17DgfOg7/HIkEoE1QvIHhIW5H/Wr1lnHOVQAmQ3tALmtnFtTf+1czWmdmLZta9gbL3m1mumeUWFxdf6mdotpGDe9Yabn7k+OmQv6eIyOXYsL0oMLx8YJ/MwA4NkSKUCaq+dTLqbprUnDLB4oDewArn3HjgE+DX9RV0zj3unJvonJuYmZnZnHgvS0J8XK21+dZv1XBzEYlswa094yKseQ9Cm6CKgOAB9b2BuovVBcqYWRzQCTjRyGseB84Dr/rPXwTGt0awrUH9UCLSVjjnIrr/CUKboNYAg82sv5klAPOBBXXKLADu9h/fCixyjWxN63/sDeBq/6VZwNbWDPpyBLff5u0soryi0sNoREQatqfwGCVnSgFITUliUN/QtzS1VMgSlL9P6dvAu8A24AXn3BYz+4WZzfUXewJIN7N84AEgMBTdzAqAh4F7zKwoaATg/wf83MzygDuB74fqM7RUVkYavbr5uskqKqvYvEurm4tIZFq7dV/geNzwPsTERM7w8oviQvnizrmFwMI6134WdFwG3NbAc7MbuL4PuKr1omxd44b35cDRUwCs3bIvomZli4hcFDz/aXyEjd67KPJSZhs3cWS/wPGazQU00mIpIuKJEyXnyN9/FIAYs4gcIAFKUK1u+IAskpMSADh+6hz7Dh73OCIRkdqCB3ENH9iDjsmJHkbTMCWoVhYXF8v4ETV/jazeVOBdMCIi9VgT9L00aWS2Z3E0RQkqBCaNyA4c527e13BBEZEwu1BeUWv33OBuiUijBBUC43JqRsTsLiwOrBQsIuK1jTsOUFHp23Ghd/cu9Mjs5HFEDVOCCoGUDomMGFSziaFqUSISKXI3FwSOJ0Vw7QmUoEJm4oiaH7wSlIhEAuccuVtqvo8mRnD/EyhBhUzwD37jziLKLlR4F4yICLBr39HA6hFpHTswJDtyNiesjxJUiGRlpNGnR1cAKiuranVKioh4Ibg1Z8KIvhG5ekSwyI6ujZscVItaE9TuKyLihdW1+p+yPYujuZSgQqjWqhKbCqiqqvYwGhFpzw4ePUXhId9mEfFxsYwZ2tvjiJqmBBVCg/t1o2unFADOnr/A1t2HPI5IRNqrVXl7A8djh/UhKTHew2iaRwkqhMyMKaP7B86Df0FERMIp+Ptn6pj+jZSMHEpQIRacoFZu3KPFY0Uk7I6dPMuuff7FYWNimDAisuc/XaQEFWI5QQsxnjx9PvBLIiISLsG1p5GDepKakuRhNM2nBBVisbExTB5VuxYlIhJObbF5D5SgwmLKmNr9UGrmE5FwKTlTytZ83+7eBkwale1pPC2hBBUGY4b0DoyYOXzsNPv9Qz1FREItd0sBF/8kHtI/KzCyuC1QggqD+PjYWp2Sn6iZT0TCZOXGttm8B0pQYRP8ixH8CyMiEirnSi/UWmYteFRxW6AEFSbjh/clPi4WgMJDJyg8fNLjiEQk2uVu3hdYwSa7Vwbd09M8jqhllKDCJCkxngk5NVvBf7x+t4fRiEh7sGJdzffMjHEDPYzk0ihBhdH08YMCxx+v363RfCISMmfPX2DDjsLA+XQlKGnMhJy+JMTHAVB05CT7D6mZT0RCY3Xe3kDz3sA+mWRltK3mPVCCCqukxPhaK5x/vD7fw2hEJJp9vCGoeS+o9aYtUYIKs+lja6rZK9TMJyIhcOZcGRt3HAicTxs7wMNoLp0SVJhNGNGXxATfpN1DxSUUHDjucUQiEm1W5e2lutrXvDe4Xze6dU31OKJLowQVZgnxcUwOWmpkxTo184lI66o9eq9tNu9BiBOUmc02sx1mlm9mD9bzeKKZPe9/fJWZZfuvp5vZYjM7a2aPNvDaC8xscyjjD5UZ49XMJyKhUXKmlE07aybnTh/XNpv3IIQJysxigceAG4Ec4HYzy6lT7D7gpHNuEPAI8JD/ehnwU+AHDbz2zcDZUMQdDmOH9iE5KQGAoyfOkL9fW3CISOv4ZMOewNp7wwZkkd65o6fxXI5Q1qAmA/nOuT3OuXLgOWBenTLzgKf9xy8Bs8zMnHPnnHPL8SWqWsysI/AA8H9DF3poxcfH1lrh/KPcXR5GIyLR5KO1Nd8nV7TR0XsXhTJB9QIKg86L/NfqLeOcqwRKgPQmXvdfgd8A5xsrZGb3m1mumeUWFxe3JO6wuGrC4MDx8nW7qays8jAaEYkGh4+dZsfew4Bv59y2uHpEsFAmKKvnWt3OluaUqSlsNhYY5Jx7tak3d8497pyb6JybmJmZ2VTxsBs5uCdd0pIBOH22lLydB5p4hohI4z7K3Rk4HjesD2kdO3gYzeULZYIqAvoEnfcGDjZUxszigE5AY5slTQMmmFkBsBwYYmZLWinesIqJieHKoFqUmvlE5HI451gW9D1y1cTBjZRuG0KZoNYAg82sv5klAPOBBXXKLADu9h/fCixyjQxpc879zjnX0zmXDVwB7HTOXd3qkYdJ8C/Qqry9lF2o8DAaEWnLdu8v5mBxCQCJCfFMGtWviWdEvpAlKH+f0reBd4FtwAvOuS1m9gszm+sv9gSQbmb5+AY+BIai+2tJDwP3mFlRPSMA27zsXun07t4FgPKKSlblaZ8oEbk0wYMjpo7pH1gQoC2LC+WLO+cWAgvrXPtZ0HEZcFsDz81u4rULgJGXHaSHzIyrJg3mmTdXA7Bs7S5mThricVQi0tZUVVWzPGjSf7R8j2glCY8F90Nt2FbIqTONDk4UEfmUvJ0HKDlTCkDn1GRGDe7pcUStQwnKY926pjJ8QA/AN3xxWa6WPhKRllm6pmb03hXjBxETEx1f7dHxKdq4mZNqalGLVm3X0kci0mznSi+wcuOewHnw90lbpwQVAaaPG0h8XCwA+w+dYG/RMY8jEpG2YsW63VT4J/r37dGV/r0zPI6o9ShBRYCUDolMHVOzoOPi1Ts8jEZE2pLg74trpwzDrL71D9omJagIce2UoYHjj3J3UVGhpY9EpHFFR06ys+AI4Jv8Hw2Tc4MpQUWIUUN6kdnFt6nY2fMXWLOlwNuARCTiLVlVU3uaNLIfnVLb9tJGdSlBRQgzY+bkmrkLi1epmU9EGlZVVc2SoNF71wS1wkQLJagIcs3kml+w9Vv3c6LknIfRiEgk27C9kJOnffMm0zp2YNywPk08o+1RgoogWRlp5AysmRO1ZPXOxp8gIu3WoqBWlqsnDSHOPxI4mihBRZhrpwwLHH+4cpvmRInIp5ScKWXN5oLA+dWTo695D5SgIs60sQPo4N8O/vCx02zeVXeHEhFp7xat2k5VVTUAg/t1o1/Prh5HFBqNJigzeyro+O5GikorSUqMZ2bQUNH3Pt7qYTQiEmmcc3zwybbA+WdmjPAwmtBqqgY1Juj4H0MZiNS4YUbNziKr8vYGFoEUEdm08wCHj50GIDkpgenjBjTxjLarqQSlDhAP9OuZzuB+3QDfUFKtLCEiF733cU3taeakIVGx71NDmtoPqreZ/RawoOMA59x3QhZZO/eZGSPYte8oAO9/vJV5146JqiVMRKTlSs6UsnpTzcam108f7mE0oddUgvph0HFuKAOR2qaPG8CfXlnB+bJyDh87zaadBxg9tLfXYYmIhz5cWTM4Ykh2d/r1TPc4otBqNEE5554OVyBSW2JCPDMnDeHtZZsBX7VeCUqk/XLO8eHKmua9G6bnNFI6OjSaoMxsQWOPO+fmtm44Euz66cMDCWpV3l5OlJyja6cUj6MSES9s2F5Ua3DEjPEDPY4o9Jpq4psGFALPAqvw9UVJmPTrmc6wAVls33OY6upq3vt4K/NvnOR1WCLigbc/2hw4vnbKMBLim/r6bvuaGsWXBfwTMBL4L+B64JhzbqlzbmmogxOYc9WowPH7K7ZRWaltOETam0PFJazbui9wPvvK6J37FKzRBOWcq3LOveOcuxuYCuQDS8zsH8ISnTBlVDZd0pIBOHXmPCs37m3iGSISbd5ZtiUw52fc8D70yOzkaTzh0uRSR2aWaGY3A38FvgX8Fngl1IGJT1xcbK2Ju299tMnDaEQk3MouVLBo1fbAeXCrSrRraqmjp4GPgfHAvzjnJjnn/tU5dyAs0QngW1kiNtb3o9pZcITd+4s9jkhEwuWj3F2cLysHfDsejBsefdtqNKSpGtSdwBB8yxx9Yman/f/OmNnp0IcnAJ1Tk5k+tmbEzsJlmxspLSLRwjnHwqBWk9lXjGxXE/ab6oOKcc6lBv1L8/9Ldc6lhStIgTlXjQwcL1+Xr/X5RNqBzbsOUnj4JOCbG3nt1OjcVqMhTTXxJZnZd83sUTO738yif1xjhBqS3Z2BfTIBqKys4p3lWzyOSERC7Y3FeYHjqycNIaVDoofRhF9TTXxPAxOBTcAc4Dchj0gaNPeamsXl31m+hfKKSg+jEZFQKjx8krVBQ8s/e3X7GRxxUVMJKsc5d4dz7g/ArcCVLXlxM5ttZjvMLN/MHqzn8UQze97/+Cozy/ZfTzezxWZ21sweDSqfbGZvmdl2M9tiZr9qSTxt3dQx/Unv7FtJ4vTZUpau0ZbwItHqzSU1taeJI/rRq1tnD6PxRlMJquLigXOuRX+um1ks8BhwI5AD3G5mdRePug846ZwbBDwCPOS/Xgb8FPhBPS/9a+fcMGAcMMPMbmxJXG1ZXFwsN109OnC+YNFGbQkvEoVKzpSyJOgP0LnXjmmkdPRqcsPC4JF7wOgWjOKbDOQ75/Y458qB54B5dcrMw9eMCPASMMvMzDl3zjm3HF+iCnDOnXfOLfYflwPrgHa1gup1U4cHtoQ/WFzC2q37PY5IRFrb28s3B1aNGdgnk5yBPTyOyBtNjeKLrTNyL64Fo/h64VvH76Ii/7V6y/hraCVAs9aPN7POwOeADxt4/H4zyzWz3OLi6Jk3lNwhgeun1ewBs2DRRg+jEZHWVl5RyTvLagZBzb2m/e4F1+RKEpehvjtatz2qOWU+/cK+0YTPAr91zu2pr4xz7nHn3ETn3MTMzMwmg21L5lw1khj/L+yW/IPk+zc2FJG2b8nqnZw552s8yujSkWljo3dL96aEMkEVAcFTnnsDBxsq4086nYATzXjtx4Fdzrn/bIU425zMrqlMG1czcffVDzd4GI2ItJaqqmpeX1Tz//mzM0cFVpFpj0L5ydcAg82sv5klAPOBuvtLLQDu9h/fCixyTfT6m9n/xZfIvtvK8bYpX5g1NnC8auMeio6c9DAaEWkNH6/fHdjzKaVDYq3m/PYoZAnK36f0beBdYBvwgnNui5n9wswubnT4BJBuZvnAA0BgKLqZFQAPA/eYWZGZ5ZhZb+An+EYFrjOzDWb21VB9hkjWv3cGE3L6Ab420Vc/UC1KpC1zzvHyB+sD53NmjgwMiGqvQroyhHNuIbCwzrWfBR2XAbc18NzsBl62ffYW1uPm68cFJvJ9lLuLL904kW5dUz2OSkQuxZrN+yg85OvhSEyI57PtaNXyhrTfxs0oMGxAFiMG9QSgurqa19UXJdImOed4+b11gfPPzMghNSXJw4gigxJUG3fLDeMDxx+s3M7J0+c9jEZELsWmnQfI3+8bjRsXF8vnrhndxDPaByWoNm70kF61FpF9Y7HmRYm0NS+/X1N7mjVlGF07pXgYTeRQgmrjzKxWLertZVu0FYdIG7J51wE27/LNwIkxY96s9rmsUX2UoKLA5FHZZPfKAHyz0F9TX5RIm+Cc4/m3cwPnV08eSvd0bbV3kRJUFDAzvnTjxMD528s2c6LknIcRiUhzbNp5gK27DwEQExPDbbMneBxRZFGCihKTRvZjgL8vqqKySrUokQjnnOO5oNrTrKlDNU2kDiWoKGFmzA+qRb27YivHT531MCIRacyG7UXs2HsYgNjYGG65fnwTz2h/lKCiyPicvgzu1w3wjeh7+b31TTxDRLzg63taEzi/ftpwMlV7+hQlqChiZsyfMylw/sHKbYF1vUQkcqzeVMCufTXznm6+fpzHEUUmJagoM2Zob4YP8G1uVlVVzTNvrfY4IhEJVlVVzd/eWBU4nz1jBOmdO3oYUeRSgooyZsadc6cEzlesy2dPYfRs2CjS1i1evYMDR08B0CEpgVtuUO2pIUpQUWho/yymjO4fOP/LglWNlBaRcLlQXsFzC2v6nj4/ayxpHTt4GFFkU4KKUl++aXJg2fe8nUVs3FHkaTwiAm8t3RxYL7NLWjI3zdSK5Y1RgopSvbt3YVbQZmd/fn0lTewFKSIhdOZcGa8G7ff0xdkTSUqM9zCiyKcEFcW+OHsC8XGxABQcOMaS1Ts9jkik/Xr+7VzOl5UD0DOzE7OmDvM4osinBBXF0jt3ZO41NQtP/u3NVZT6/4OISPjsP3SCd5dvCZzfMXcqsbH6+m2K7lCUu/n6cXRJSwbg5OnzvPK+Ju+KhJNzjidf+ZhqfxP7yME9mTwq29ug2gglqCiXlBjPHZ+rGXa+YEkeR45r8q5IuORu2UfeTt8gJQPuvXkGZtb4kwRQgmoXZk4awqC+NUsg/fm1TzyOSKR9qKys4qlXPw6cXz8jh3490z2MqG1RgmoHzIx7b54eOF+Zt5dNOw94GJFI+/Dm0k2B5caSkxKYf+OkJp4hwZSg2omh/bO4YsKgwPn/vriMysoqDyMSiW7HTp7lhXfWBs6/OHsinVI1KbcllKDakbvmTiUxwTfv4sDRU7y2aKPHEYlErz+9soIL5RUA9Mnqwo1XjvA4orZHCaodSe/ckduDVjt/6d21GjAhEgK5W/axKm9v4Pz+L15FnH9OojSfElQ7M+eqkWT3ygB8O+8+8dIKrTAh0ooulFfwxEvLA+dXTx5KzsAeHkbUdilBtTOxsTH8/RevDKzTt3Zr7b/0ROTyvPzeeo6eOANASodE7po71eOI2i4lqHZoSHZ3rptes07fH19azrnSCx5GJBIdCg4c49UPNwTO75w7RQMjLkNIE5SZzTazHWaWb2YP1vN4opk97398lZll+6+nm9liMztrZo/Wec4EM9vkf85vTTPeLslXbqr5j3Py9HmeelVzo0QuR1VVNY89u5Tq6mrAN3L2uqAFm6XlQpagzCwWeAy4EcgBbjeznDrF7gNOOucGAY8AD/mvlwE/BX5Qz0v/DrgfGOz/N7v1o49+qSlJfO3WKwPni1Zt15YcIpfh9UUbA5uDxsXF8q0vX60VIy5TKGtQk4F859we51w58Bwwr06ZecDT/uOXgFlmZs65c8655fgSVYCZ9QDSnHOfOF/P/p+Bz4fwM0S1aWMHMHXMgMD5755dStmFCg8jEmmbDhw9xfPv5AbOvzR7Ir26dfYwougQygTVCygMOi/yX6u3jHOuEigBGlsHpJf/dRp7TQDM7H4zyzWz3OJibXnekK/ddgUdkxMBKD55hr++od13RVqiurqax55ZEpj4PqBPJvOuHdP4k6RZQpmg6qvb1h3P3Jwyl1TeOfe4c26ic25iZmZmIy/ZvnVOTebem2cEzt9etllNfSIt8NqHG9mx9zAAMTExfOv2mdpKo5WE8i4WAX2CznsDBxsqY2ZxQCfgRBOv2buJ15QWumriYCbk9Auc//dfF3HmXFkjzxARgL1Fx3ju7TWB81tvGB+YZyiXL5QJag0w2Mz6m1kCMB9YUKfMAuBu//GtwCLXyKxR59wh4IyZTfWP3rsLeL31Q29fzIxvfnkmqSlJgG9U3+MvLtMEXpFGlFdU8l9/+ZCqKt+ovUF9u3HL9eM8jiq6hCxB+fuUvg28C2wDXnDObTGzX5jZXH+xJ4B0M8sHHgACQ9HNrAB4GLjHzIqCRgB+A/gjkA/sBt4O1WdoTzqnJvPN268OnH+8fjfL1u7yLiCRCPe3N1ZTePgkAAnxcXznzmu1nFEriwvlizvnFgIL61z7WdBxGXBbA8/NbuB6LjCy9aKUiyaPyubaKcNYtGo7AI+/uJyh/bPonp7mcWQikWXD9kLeXJoXOL/n89M0ai8E1JMntdx78/RAQiotK+c3T76vbTlEgpwoOcd//WVR4Hx8Tl9umFF3iqe0BiUoqaVDUgLfvWsWMTG+X43dhcX8ZYGGnouAb7WIR57+gNNnSwFf07gm5IaOEpR8ypDs7tzxuSkF97RvAAAWS0lEQVSB8zeX5rF6U4F3AYlEiBfeXcvW3YcA35yX7941i86pyd4GFcWUoKRec68ZXWvo+aN/WxxYoVmkPdq4o4iX363ZIfe22RMZNaTedQKklShBSb3MjH+44xrSO6cAcK70Ag/98d3ADqEi7cnRE2d4+Kn3A6sCjBzck9s+M97TmNoDJShpUGpKEt+/5/rArPiCA8f43XMfaX6UtCsXyit46I/vcva8b0uazqnJ/OOdNf20Ejq6w9Koof2zuC9oKaRla3fx5pJNHkYkEj7OOX733EcUHDgG+Db8/OG9N9C1U4rHkbUPSlDSpBtm5DBr6rDA+Z9f/4RNOw94GJFIeLy5ZFOtCetfveUKhg3I8jCi9kUJSppkZnzt1isZ3K8bANXO8R9/eo8DR095HJlI6KzZXMDTr30cOL9u2nDNdwozJShplvj4WH547w10SfMNqT1XeoF/+/1CSs6UehyZSOvbW3SMR57+MDAoYkh2d756yxWextQeKUFJs6V37siPv3YjCfG+FbKOHD/NQ0+8S3lFpceRibSeYyfP8m9/WBgYsdqtayoPfnU28fFaZy/clKCkRQb2zeR7d18X2Jhrx97D/PffFmtkn0SFc6UX+OXjb3Py9HkAkpMS+Ke/n0On1A4eR9Y+KUFJi00elc1dn58WOP94/W7++NJyJSlp08orKvnV/77DvoPHAd/mgz+89wb6ZHXxOLL2SwlKLsnnrh7NjVfWLCr/zvItvPDO2kaeIRK5qqqqefipDwLLGAF840tXMXpo70aeJaGmBCWXxMy475YZTB83MHDthXdyeXvZZg+jEmk55xy/f/4j1mwuCFy743NTuDZoaoV4QwlKLpmZ8Y93XMuYoL8y//jScj5cuc3DqESazznHEy+vCOyBBjD3mjF84TrtjBsJlKDkssTFxfKj+z7DoL7dAtd+9+xSFq/a4WFUIk1zzvHkqx/XqvVfPXkod82b6mFUEkwJSi5bUmI8//z1OWT3ygDAAY89s5ila3Z6G5hIA5xzPP3aJ7y1tGbZrhnjB/HN+TO1t1MEUYKSVpGaksTPv3UT/XqmA74k9d9/XaSalEQc5xxPvfoJbyyp2bJ92tiB/OMd1wYWRpbIoJ+GtJqLSapvj66AL0k9+sziWn+linipurqa/3l2KW8urUlOU0f357t3KjlFIv1EpFWldezAz7/1uUBzH8CfXlnBS++t0zwp8VRlZRUPP/1hrQERU0f353t3X0dcnFaJiERKUNLqOqV24Bf/8DmGZHcPXHv2rdU8+erHVFdXexiZtFfnS8v55eNv88mG3YFrV08eygP3XK/kFMGUoCQkUjok8n++eROjh9QMQX9r6SZ+8+T7WrtPwur4qbP8829fZ+OOosC1OVeN5NtfvlrNehFOPx0JmaTEeH58/2ymju4fuLYyby8/++8FWgVdwmLfweP8+JFXA8sXAdw2ewL33jxDo/XaACUoCamE+Dh+cO8N3DRzdODarn1HefDhV2p9aYi0ttWbCvin/3yN46fOAb619b51+9XMv3GSklMboQQlIWdm/N3N031/tfqvHT1xhgcffpVPNuzxNDaJPs45Xnx3LQ/98R3KLvi2zLg4V0/LF7UtSlASNp+dOYoffXU2iQnxgG/16F8/+R7PvrVagyekVZwvLefXT77PcwvXBK5ldknll9/9fK0luaRtUIKSsJo8KptfPfAFsjLSAtdeem8d//I/bwb24BG5FHuLjvGj37zMyo01tfKRg3vy7z+4OTCBXNqWkCYoM5ttZjvMLN/MHqzn8UQze97/+Cozyw567Mf+6zvM7DNB179nZlvMbLOZPWtmSaH8DNL6+vboyq8euLnWX7Sbdx3k+//+InlBI61EmsM5x3srtvLgI69yqLgkcH3OVSP56dc/S1pHbTbYVlmoJk+aWSywE7geKALWALc757YGlfkmMNo593Uzmw98wTn3JTPLAZ4FJgM9gQ+AIUAWsBzIcc6VmtkLwELn3FONxTJx4kSXm5vb6p9RLk91dTUvvLuWl95Zy8XfQgPmzRrL/BsnaYttaVLJmVL+8MJHrMrbG7iWmBDPN+fP5IoJgzyMTBpjZmudcxObKhfKGtRkIN85t8c5Vw48B8yrU2Ye8LT/+CVglvmG18wDnnPOXXDO7QXy/a8HEAd0MLM4IBk4GMLPICEUExPD/Bsn8bNv3hT4K9cBr324gR/95mWN8pNG5W7Zx/ceeqFWcurboyv/8cNblJyiRCgTVC+gMOi8yH+t3jLOuUqgBEhv6LnOuQPAr4H9wCGgxDn3Xn1vbmb3m1mumeUWFxe3wseRUBk9tDe/+dGttSb17j90gh/++mVeem8dlZVVHkYnkebs+Qs89swS/t/jb9eaT3fDjBwe+v7N9OrW2bvgpFWFMkHVN9GgbntiQ2XqvW5mXfDVrvrja/pLMbM76ntz59zjzrmJzrmJmZmZLQhbvNC1Uwo/++Znue+WGcT7l56pqqrm2bdW84Nfv8zOgiMeRyhec86xYv1uvvPL52qtp9c5NZmf/P0c/v6LV5EQH+dhhNLaQvnTLAL6BJ335tPNcRfLFPmb7DoBJxp57nXAXudcMYCZvQJMB/4aig8g4WVmzLlqFKOH9ua3f1nE7kJfzbfw0An+6ZFX+cwVI5g/ZxKpKRoX094cKi7hyVc+Zu3WfbWuTxs7kPtvu0IDIaJUKBPUGmCwmfUHDgDzgS/XKbMAuBv4BLgVWOScc2a2AHjGzB7GV1MaDKwGqoGpZpYMlAKzAI1+iDK9u3fhVw98gbeWbuaZt1ZTXlGJA95ZvoXl6/L5yk1TuG7aMGJiNEsi2pVdqOCV99fz2qINVFXVzJXrkpbM1267kilBy2hJ9AnZKD4AM5sD/CcQC/zJOfdvZvYLINc5t8A/RPwvwDh8Naf5zrk9/uf+BLgXqAS+65x723/9X4Av+a+vB77qnLvQWBwaxdd2HT1xhj88/xEbthfWup7dK4O75k3V5MsoVV1dzZLVO3nu7TWBpYrA1/Z/3fTh3Dl3KikdEr0LUC5Lc0fxhTRBRQolqLbNOceqvL089eonFJ88U+ux0UN6c+fcKQzoo37GaOCcY+3W/fx1wUoKD5+s9digvt342q1XMKhfN4+ik9aiBBVECSo6lFdU8uoHG3j1g/VU1BnZN3lUNl+cPZH+vTMaeLZEMucc67bu54V31pK//2itxzqlduCOm6ZwzZShWuQ1SihBBVGCii7HT53l+bdzWbRy+6eGhU4amc0XrhvL0P5ZnsQmLVNdXc3qTQW88v76wKCYixIT4vn8rDHMvWYMSYnxHkUooaAEFUQJKjoVHj7Jc2+tZmXQRM2LhvbPYt61Y5g0sp8GU0Sg8opKFq/awRtL8motTwQQFxfLDdOHc8sN4+mcmuxRhBJKSlBBlKCi276Dx3nhnbW1Fgm9qFvXVG6YkcOsqcM0FDkCHCou4b0VW/lw5XbOldYe23QxMX1+1ljSO3f0KEIJByWoIEpQ7cP+Qyd4fdFGlq3dVWtIMvi+/KaM7s+sqcMYNbinalVhVF5Ryeq8Ahat2l5r2/WLkpMS+MyMHObMHEXXTikeRCjhpgQVRAmqfTlRco6FSzfx/ifbOHv+0zMQ0junMHPiEK6YMIi+Pbqq4z0EnHNs3X2IZWt3sWLdbs6XlX+qTLeuqcy5ahTXTRtGh6QED6IUryhBBVGCap/KKypZsW437yzf8qmRYRf17t6FGeMHMmV0fyWry1RdXc22PYdZnVfAivX59e7vZcD4nH7MvnIE44b30f1up5SggihByd6iYyxevYOPcndx5lxZvWW6dU1l0qhsxg3vy4hBPbSuWzOcPX+BjTuKWLd1P7mbC+qtsQJkZaRx9eShzJw0hG5dU8McpUQaJaggSlByUWVlFWu37mfF+t2s2VRAeUVlveXi42LJGdiDUUN6MXJwTwb0ziQ2Vv1WF8or2LH3CFt2H2Lj9kLy9x391FD/i9I6dmDamAFcMWEQwwdkqbYkAUpQQZSgpD5lFypYu3U/q/L2sm7rfkrr6Se5KDEhnqHZ3RnSvztDs7szuF+3qF+01jnH8VPn2FFwhF0FR9hRcITdhcWfGoASrEtaMpNGZTN1zABGDuqppC71UoIKogQlTamsrGLL7kOs3bKPDdsKOXD0VJPPyejSkf69MujfO4O+PbrSO6sLPTLSiItrezsBl1dUcuDIKfYfOsH+QyfYU3iMvQeONdgcepEBA/t2Y+zwPkwa0Y+BfTNVU5ImKUEFUYKSlio+cYa8nUVs3nWQrbsPcezk2WY9LzY2hqz0NLIyOtEjsxPd0lPJ7JpKZpeOpHdOITUlyZMv8Orqak6dKeX4qbMUnzxL8YmzHDl2mkPFJRwqLuHYyTMNNtXV1SerCzkDezJicE9GD+kV9TVJaX1KUEGUoORyHT1xhp17j7Cj4DA7C46y98CxRpu6GhIbG0Pn1A50Tk0mrWMSHZOTSE1JJLlDIslJCXRIjCcxIY74+FgS4uOIi40hxgwzIybGqK52VDtHdbWjorLK96+ikrILlZwvK6fsQgXnSi9w5vwFzpwt4/S5Mk6dPk/JmfPNTkDBkhLjGdQ3k6HZWQzO7saQft3plKoJz3J5mpugNExJpBm6dU2lW9dUrpgwCPA1CR44eoq9RccoOHCcwsMnKTx8otbWEPWpqqrm+KlzTZYLNwOyMjvRJ6sLfbK60q9XOgN6Z5CVkaYmO/GMEpTIJYiLi6Vfz3T69Uyvdb20rJzDx05zMKjp7Ji/Se14yblGB2KEWsfkRNI7d6Rb11QyunQks2sqPTJ9TZHd01M1rF4ijn4jRVpRh6QE+vfOaHDbjwvlFZw8XUrJmfOBZrgz58soLaugtKyc82XllFf4mu3KK6uorKym2lVTXe1wDmJiDDOIsRji42L9TYGxJCbEkZyUQFJiPCkdEklNSaRjchJpKUl0Tkumc2qHNjl4Q9o3JSiRMEpMiCcrI56sjDSvQxGJeJqkICIiEUkJSkREIpISlIiIRCQlKBERiUhKUCIiEpGUoEREJCIpQYmISERSghIRkYjULhaLNbNiYN9lvEQGcKyVwmmrdA90D0D3AHQP4PLvQT/nXGZThdpFgrpcZpbbnJV3o5nuge4B6B6A7gGE7x6oiU9ERCKSEpSIiEQkJajmedzrACKA7oHuAegegO4BhOkeqA9KREQikmpQIiISkZSgREQkIilBBTGz2Wa2w8zyzezBeh5PNLPn/Y+vMrPs8EcZWs24Bw+Y2VYzyzOzD82snxdxhlJT9yCo3K1m5sws6oYcN+cemNkX/b8LW8zsmXDHGGrN+L/Q18wWm9l6//+HOV7EGSpm9iczO2pmmxt43Mzst/77k2dm41s9COec/vn64WKB3cAAIAHYCOTUKfNN4Pf+4/nA817H7cE9uAZI9h9/oz3eA3+5VOAjYCUw0eu4Pfg9GAysB7r4z7t5HbcH9+Bx4Bv+4xygwOu4W/keXAWMBzY38Pgc4G3AgKnAqtaOQTWoGpOBfOfcHudcOfAcMK9OmXnA0/7jl4BZZmZhjDHUmrwHzrnFzrnz/tOVQO8wxxhqzfk9APhX4N+BsnAGFybNuQdfAx5zzp0EcM4dDXOModace+CANP9xJ+BgGOMLOefcR8CJRorMA/7sfFYCnc2sR2vGoARVoxdQGHRe5L9WbxnnXCVQAqSHJbrwaM49CHYfvr+gokmT98DMxgF9nHNvhjOwMGrO78EQYIiZrTCzlWY2O2zRhUdz7sHPgTvMrAhYCPxDeEKLGC39vmixuNZ8sTauvppQ3TH4zSnTljX785nZHcBEYGZIIwq/Ru+BmcUAjwD3hCsgDzTn9yAOXzPf1fhq0cvMbKRz7lSIYwuX5tyD24GnnHO/MbNpwF/896A69OFFhJB/H6oGVaMI6BN03ptPV9kDZcwsDl+1vrEqcFvTnHuAmV0H/ASY65y7EKbYwqWpe5AKjASWmFkBvrb3BVE2UKK5/xded85VOOf2AjvwJaxo0Zx7cB/wAoBz7hMgCd8iqu1Fs74vLocSVI01wGAz629mCfgGQSyoU2YBcLf/+FZgkfP3FkaJJu+Bv3nrD/iSU7T1O0AT98A5V+Kcy3DOZTvnsvH1w811zuV6E25INOf/wmv4BsxgZhn4mvz2hDXK0GrOPdgPzAIws+H4ElRxWKP01gLgLv9ovqlAiXPuUGu+gZr4/JxzlWb2beBdfCN4/uSc22JmvwBynXMLgCfwVePz8dWc5nsXcetr5j34D6Aj8KJ/fMh+59xcz4JuZc28B1GtmffgXeAGM9sKVAE/dM4d9y7q1tXMe/B94H/N7Hv4mrbuiaY/WM3sWXxNuBn+frb/A8QDOOd+j6/fbQ6QD5wH/q7VY4ii+ykiIlFETXwiIhKRlKBERCQiKUGJiEhEUoISEZGIpAQlIiIRSQlKxM/Mqsxsg5ltNrMXzSy5hc8/28LyT5nZrfVcn2hmv/Uf32Nmj/qPv25mdwVd79mS92skjiv9K5JvMLMOQdezG1nJ+hf+CduY2Xdbeq9EmkMJSqRGqXNurHNuJFAOfD34Qf+ExJD/n3HO5TrnvlPP9d875/7sP70HaJUEBXwF+LX/s5c2M8afOec+8J9+F1CCklanBCVSv2XAIH8tYpuZ/Q+wDuhjZreb2SZ/Teuh4CeZ2W/MbJ359srK9F/7mpmtMbONZvZyndrGdWa2zMx2mtlN/vJXm9mnFqI1s5+b2Q/8ta6JwN/8tZ7PmtmrQeWuN7NX6nn+LPPtXbTJfHv9JJrZV4EvAj8zs7/Vcx9izex//TWs9y7WsC7W/szsO/gS5WIzW9yiOyzSBCUokTr86yzeCGzyXxqKb1uBcUAF8BBwLTAWmGRmn/eXSwHWOefGA0vxzbwHeMU5N8k5NwbYhm8Nt4uy8S24+1ng92aW1FR8zrmXgFzgK865sfhm9A+/mBDxzeh/ss5nSgKeAr7knBuFbxWZbzjn/ohvyZofOue+Us/bDca3rcYI4BRwS51Yfotv/bVrnHPXNBW7SEsoQYnU6GBmG/B9+e/Ht7QVwD7/fjcAk4Alzrli/5Yrf8O3sRtANfC8//ivwBX+45H+WtImfM1pI4Le8wXnXLVzbhe+teyGtTRo//I6f8G39UNnYBqf3gZlKLDXObfTf/50UNyN2euc2+A/XosvoYqEhdbiE6lR6q+RBPjXGzwXfKkFr3dxHbGngM875zaa2T341jerW6ah8+Z6EngD3waKL/qTZ7BL3VgzeLX6KqBDQwVFWptqUCItswqYaWYZZhaLb0+gpf7HYvCtcg/wZWC5/zgVOGRm8fhqUMFuM7MYMxuIb3vxHc2M44z/dQFwzh3E19T2z/gSYl3bgWwzG+Q/vzMo7stVKxaR1qIalEgLOOcOmdmPgcX4aiULnXOv+x8+B4wws7X4dlv+kv/6T/Eltn34+rWCv8x34EsU3YGvO+fK/LW2pjyFr8+qFJjmH333NyDTObe1nrjLzOzv8K1CH4dvO4nfN/+TN+px4G0zO6R+KGlNWs1cJEr450utd8490WRhkTZACUokCvhrbeeA66Nwl2Npp5SgREQkImmQhIiIRCQlKBERiUhKUCIiEpGUoEREJCIpQYmISET6/wH5AeHNfBN/fgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Solution\n", "\n", "# If we start with a uniform prior, \n", "# we can see what the likelihood function looks like:\n", "\n", "pmf = Beta(1, 1).MakePmf()\n", "blaster = AlienBlaster(pmf)\n", "blaster.Update(2)\n", "thinkplot.Pdf(blaster)\n", "thinkplot.decorate(xlabel='Probability of hit',\n", " ylabel='PMF')" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XlYXdd56P/ve85hBjFPAiEkgdA8WFiDB3mQHcuOY2VwUjuje924TZO2ado+N7n9Jbc3bZ7nl3vb5ja/pk3TZk4TO0njRElsy7Elz9aALCSBRkCAQEgMAsQM57B+f+yjffZBIEDiDMD7eR4e7+ls1jlGvKy13v0uMcaglFJKRRtXpBuglFJKjUcDlFJKqaikAUoppVRU0gCllFIqKmmAUkopFZU0QCmllIpKGqCUUkpFJQ1QSimlopIGKKWUUlHJE+kGhENWVpYpLi6OdDOUUkoBhw8fbjfGZE923bwIUMXFxVRUVES6GUoppQARaZjKdTrEp5RSKippgFJKKRWVNEAppZSKShqglFJKRSUNUEoppaKSBiillFJRSQOUmnV8Ph+6ErRSc9+8eA5KzV4+n4/a2lqam5vp7++nv78fr9dLYmIi+fn55Ofnk5eXR0xMTKSbqpSaYRqgVFQaHR2lrq6O6upq+vv7rznf399PbW0ttbW1xMTEsHHjRpYuXYqIRKC1SqlQ0AClok5PTw+vvvoqPT09U7p+ZGSEgwcP0tDQwJYtW0hKSgpxC5VS4aABSkWVzs5OXnnlFQYHB+1j8fHxrFy5ktzcXBITE4mJiaG9vZ2Wlhbq6+vtHtalS5d47rnnuO222ygoKIjUW1BKzRCZD5PN5eXlRmvxRa/R0VHO1LfyyoEqjh0/gcdliI1xk5ocx3vu28qa1SvxeMb/W8rr9XL8+HFOnTplH3O5XNxxxx0apJSKUiJy2BhTPtl12oNSEdPS1s2v9lZy8Hg9re1ddHR02Nl5Ii4yMzM53XqIu8q7eOCO1RTlZ1xzD4/Hw8aNG1m0aBFvvfUWfX19jI6O8sYbb3DnnXeycOHCcL8tpdQM0R6UCjtjDL976yTfffYthke8+Hw+WlvbGB31AeByucnKyrwmM297eSlPffBOEuJjx71vf38/L730En19ff77uNi+fTv5+fmhfUNKqWmZag9KA5QKq84r/fzrT17l8Amr2r4x0NHRjgsfyxamkJORRNmKVfiMiyMnGrnQ1h30+vzsVP7iiftZUpg17v37+vp4+eWX7SAVExPDzp07SU5ODu0bU0pNmQYoBw1Q0aG9s5e//qdf0t7Zax9LjDGsXRRDQXYiLhHuvvtuu8djjKHq7AWee+04B4/X269xu108+f7beeCO1eN+n76+Pl566SU7eSI9PZ37778ft9sdujenlJqyqQYorSShwqKnb5C//dffBgWn7bcU864NC1iUk4RLhFWrVgUNx4kIa5cX8N//YCef/dgO4mKtIT+fb5Rv/ex1fvPKsXG/V1JSEnfccQcul/Xj3dnZSWVlZQjfnVIqFDRAqZAbHBrhK//2HE2XOgGrB/Tfn3wXSzNG8LitH8Hs7GzWrl074T3uLC/l7//qAxQXBIb2vvvsW7z09slxr8/MzGTDhg32/pkzZzh//vxMvB2lVJhogFIh5fX6+D/feZGzDa0ACPBnH9tBsmeQ3l6rNxUbG8ttt91m93gmsjAnja/82S7KluTZx7759Ku8XnF23OuXL19OYWGhvX/gwAEGBgZu8h0ppcJFA5QKqWdfrqTyVKDn8geP3smmlQVUV1fbx9asWUNiYuKU7hcfF8Nf/+GDdpKEAb7+o70cP9N8zbUiwpYtW+x7j4yM8M4779zEu1FKhZMGKBUyDRc6+Nmew/b+o++6hZ13rubo0aN4vV4AFixYQGlp6bTum5QQx5c+9W4W5aUDMGoMX/vBS3ReubZmX2xsLFu3brX3GxsbuXjx4o28HaVUmGmAUiHh9fr4+o/24fONAlC6OIffe7Cczs5O6urq7OtuueWWSYf2xrMgOYEv/fHDLEhOAKC7Z4D/+4OXGB0dveba3NxcFi9ebO9XVFTg8/mm/T2VUuGlAUqFxH/97gj1ze0AxHjc/MlH70VEOHw40KNauHDhTT1Em5GaxJ9/fAdX65dXnb3AMy8cHvfajRs32g/+9vT0cPLk+MkVSqnoEdIAJSI7ReS0iNSIyOfHOR8nIs/4zx8QkWL/8c0iUun/Oioi75vqPVXknWtq5+cvBuZ6PvzwZgpy0mhtbaWtrQ2w5oc2btx4099rXVkhj+7cZO//157DHD3ddM11CQkJrFu3zt6vrq62kzSUUtEpZAFKRNzAN4AHgVXA4yKyasxlTwKdxpgS4GvAV/3Hq4ByY8wGYCfwbyLimeI9VQQZY/jeL9+yh9rKluTx8F1W+viJEyfs65YuXcqCBQtm5Ht+6IFNrF1uFYY1wL/+5FUGh0auua60tJT0dP+81egox46N/xyVUio6hLIHtRmoMcbUGWOGgaeBXWOu2QV837/9c2CHiIgxpt8Y4/Ufj8f6vTPVe6oIOnammaqzFwBwifDpD9+Ny+Wis7MzKDlh5cqVM/Y9XS4Xn/34DpIT4wBo6+zh6ecOXXOdiLBpU6C31dDQQFdX14y1Qyk1s0IZoAoA55ORTf5j417jD0jdQCaAiGwRkWrgOPBH/vNTuSf+1z8lIhUiUnF1WEmFljGGH/36gL2/Y9sKCnLSgODeU1FRESkpKTP6vdNSEnnivbfZ+7955Ri1jdf+f8/Ozg6qcH706NEZbYdSauaEMkCNt/b22MJ/E15jjDlgjFkN3Ap8QUTip3hP/K//ljGm3BhTnp2dPY1mqxv19tE66s5bQSHG4+aDD1i9lZ6eHhobG+3rVq0Kzajs3ZuXBw31/cvTr9pZhE7r16+3ty9cuID+AaNUdAplgGoCFjn2C4ELE10jIh4gFbjsvMAYcxLoA9ZM8Z4qAny+UX7ym4P2/kPb15CZZlUQd2bM5eXl2fNAM01E+MMPbSfGYxWFrW9u59fj1OtLS0sLSjs/evQo86FoslKzTSgD1CGgVESWiEgs8Biwe8w1u4FP+LcfBfYaY4z/NR4AEVkMlAH1U7ynioBXDp22l8ZIiI/lffdZGXoDAwOcO3fOvi5Uvaer8rNT+dDOQJHkn75weNwHeNeuXYuI1SFva2ujpaUlpO1SSk1fyAKUf87oM8Ae4CTwU2NMtYh8WUQe8V/2bSBTRGqAzwFX08bvAI6KSCXwLPDHxpj2ie4ZqvegpmZ0dJSf7wmklb93xwZSkuIBqK2ttTP6MjIyyMnJCXl7HrlnnV1lYmh4hGeevzZhIiUlhWXLltn7VVVV2otSKsqE9DkoY8xzxpjlxphlxpiv+I99yRiz2789aIz5oDGmxBiz2RhT5z/+Q2PMamPMBmPMLcaYX17vniqyDh6vp/VyDwDJiXF2WrkxhtraWvu6srIyu9cSSh6Pm4/v2mbvv/TWSc5f7LzmujVr1thVLDo6OnQuSqkoo5Uk1E377avH7e133baK+DirYkNLS4u9aGBsbCyLFi0a9/WhsHHlItYttyqZG+CHv9p/zTUJCQksWbLE3ndmGiqlIk8DlLopdefbOFFrzd+4XK6gVW6dvaclS5aEdUVbEeHju7baaZ+HTzRwbJwKE87nsVpaWujsvLanpZSKDA1Q6qY4s+S2bVhKVrqVuTcwMEBzc2AJjJKSkrC3bUlhFndtLrP3v/+r/dfMM6WkpAT17LQXpVT00AClbtjl7j7ePBLoJb3n7sCKuHV1dXYwyM7OnrGyRtP1+EO3BqWdv3207pprnJmFjY2N9PT0hK19SqmJaYBSN2zPG9X2g7BlS/IoXZwLWMkRNTU19nWR6D1dlZWebCdtgJV2PrYXlZGRQV5eYJVerXSuVHTQAKVuyMiIjz1vBobDHnb0ni5evBix5IjxPHLveuJircSN8y2XJ+1F1dfXMzQ0FLb2KaXGpwFK3ZCDVfX09A0CkJmWxJa1gWw454KE4U6OGM+C5AQeujOQvDFeLyonJ8eucOHz+YISPJRSkaEBSt2QVw6etrfv2bICt9v6UfJ6vUHJEUuXLg1728Yzthf1VmVwL0pEKCsLJFScPXt23NV5lVLhowFKTdvl7j6OnAgUf73HkSnX1NRkL6e+YMEC0tLSwt6+8SxITuDd29fY+z97oeKaXlRRURFxcdaSHf39/TQ1XZuWrpQKHw1QatpePXTGLiG/alk+eVmBDL2GhgZ7u7i4OLwNm8R77lkX6EVd7LxmLsrtdgcldJw+fRqlVORogFLTYoxh34HAL+57t6ywt4eGhoKKrjorhkeDsXNRv3yp8ppeVGlpqV2Oqb29ncuXg4rrK6XCSAOUmpazDa00t1qr0MbFxrBtQ2COqbGx0f6Fn5mZSXJyckTaeD0P37MOj/+5qNrzbVTXBK/WkpCQEJR1eObMmbC2TykVoAFKTcveA6fs7ds2LrXr7kF0D+9dlZaSyL1bAnNmv3y58pprnMkSDQ0NDA4OhqVtSqlgGqDUlA2PeHnznUD6tXN4r6+vz64GLiIUFRWFvX1T9cg96+0afUdOnqfhQkfQ+czMTDIyMgBrKRHnelZKqfDRAKWm7FBVA/2DwwDkZS1g5dJA9QVn7yk3N5f4+Piwt2+q8rNT2bI+MDQ5thclIkHJEjU1NbpWlFIRoAFKTdlbjrp7d5aXBq3t5AxQ0ZYcMZ737dhgb79xuMZez+qqxYsXExNjDV/29vbS2toa1vYppTRAqSkaHBrhcHUgCN2+MdDD6O3tpavLSpxwuVwUFhaGvX3TVbI4hzWlCwEYNYbfOKqyA3g8nqB5tLNnz4azeUopNECpKTp8opERr/UAbmFuur2kOhD0QGtubi6xsbFhb9+N2HVvoBf18v5T9A8MB513DvM1NzdrsoRSYaYBSk2Jc3jvto3Lgs45A9Rs6D1dtXHlIgpyrEoXg0Mj7DsY/GBuWloamZmZgJUs4awxqJQKPQ1QalJjh/ecAWpwcNDO3oPZFaBEhIe2B6qwP/fa8Wvq72myhFKRowFKTep6w3vOwrBZWVlRnb03nrs3Lycx3hqSvNh+hcOOGoMQnCzR19fHpUuXwt5GpeYrDVBqUtcb3jt//ry9Hel1n25EfFwM99+20t7/7avHg8673W6WLAksJaLLcCgVPhqg1HVdb3hvZGQkqEcxGwMUwIN3rsHlT5k/fqb5mgd3nUuGNDU16WKGSoVJSAOUiOwUkdMiUiMinx/nfJyIPOM/f0BEiv3H7xeRwyJy3P/fex2vecV/z0r/V04o38N8d73hvQsXLthzNmlpaSQlJUWkjTcrOyOFzesCvaSxvaj09HR7McPR0dGgZ76UUqETsgAlIm7gG8CDwCrgcRFZNeayJ4FOY0wJ8DXgq/7j7cB7jDFrgU8APxzzuo8YYzb4v/QJyhA6cCxQ5mfbxuDFB53Ze7O193TVw3cFkiVeqzhLb39wL2nZskDPUbP5lAqPUPagNgM1xpg6Y8ww8DSwa8w1u4Dv+7d/DuwQETHGHDHGXC0zXQ3Ei0hcCNuqxuH1+oIWJtzq6GX4fD4uXAhUAp9N2XvjWbE0j+KCLABGvL6gJUXASpZwuax/Lp2dnXR2doa9jUrNN6EMUAXAecd+k//YuNcYY7xAN5A55poPAEeMMc4/ab/rH977ojjr7TiIyFMiUiEiFc40aDV1J+su2rX3MtOSWLww8L+mra0Nr9cLQFJSEqmpqRFp40wRER50rBW1583qoJTy2NjYoCCsvSilQi+UAWq8wDH2IZLrXiMiq7GG/f7Qcf4j/qG/O/1fHxvvmxtjvmWMKTfGlGdnZ0+r4cpSURWYa7l1TXFQ7T1n76mgoIAJ/k6YVe64pcROOW9p6+bo6eAl353DfPX19fbS9kqp0AhlgGoCnBMThcCFia4REQ+QClz27xcCzwIfN8bYub3GmGb/f3uAH2MNJaoZZoyhorre3i9fUxx03vn808KFC8PUqtCKj4vhHsdaUS+8Xh10Pjc3l8TERACGh4eDPgOl1MwLZYA6BJSKyBIRiQUeA3aPuWY3VhIEwKPAXmOMEZE04LfAF4wxb169WEQ8IpLl344BHgaqQvge5q3m1i4utl8BrJVz15QEglBPTw+9vb2AVVQ1J2fuJFI+cEdgmK+iqp42R5VzEQl6JkqH+ZQKrZAFKP+c0meAPcBJ4KfGmGoR+bKIPOK/7NtApojUAJ8DrqaifwYoAb44Jp08DtgjIseASqAZ+PdQvYf5zDm8t3FFITExbnvf2XPIzc3F7XYzVxTkpLFuuTXXZIDfvXUy6LzzmaiLFy8yMDAQzuYpNa94QnlzY8xzwHNjjn3JsT0IfHCc1/0d8HcT3HbTTLZRje9QVb29PXZ4zzn/NFeG95weuGMVx85Y80+/e/skH9q5CY/HCsLJyclkZ2fT1taGMYaGhgZWrFhxvdsppW6QVpJQ1+jpG+R03UXAymK5ZVVg+faRkZGg4rBzMUDduqaYjFTroeMrvQMcOF4fdN45zKfLwSsVOhqg1DXeOdFop1KWFueSmpJgn7t48WJQ9YirSQNzidvt4r5tgfp8L74ZnCxRVFRkD2t2dXXpM1FKhYgGKHWNQ475p/I1wcu3z/Xhvat2bF1hPwNRdfYCF1q77HMxMTEUFAQe6dNelFKhoQFKBfH5Rjl6KvB8dfnqYnvbGENLS4u9P5cDVFZ6MptWB4Lz9ZIlGhoarllHSil18zRAqSBnG1qDqkcU5QeKw3Z1ddlZa7GxsWRlZUWkjeFy/+2B0pH7Dp5mZCTwYG5eXh4JCdbQ5+DgIBcvXgx7+5Sa6zRAqSBHHL2n9WWLJqwekZeXNyeqR1zPLSsXkZWeDFiJI87CuSLC4sWBHpY+E6XUzNMApYJUngwEqA0rgyuUO3sJc3l47yqXy8WOrYEU8j1jkiWcw3zNzc0MDw+HrW1KzQcaoJStp2+Q2kZr9RIB1pcFiqN6vV7a29vt/by8vHA3LyLu27bSXszwRG0LTZcCGXupqamkpaUB1jpRztWFlVI3TwOUsh093WSnl5csziE5MbDCyaVLl+xEgNTUVHv+Za7LSE0KymR8+e1TQeedz0TV19eHq1lKzQsaoJSt8tTUhvfy8/PD1qZocP9twckSXm8gWcI5D9Xa2kpfX19Y26bUXKYBSgFWCrlz/umWlUVB553p5fMtQG1YUUhmmlVZoqdvMKiyREJCQtBwpy4Hr9TM0QClAGhsuUznlX4AkhLiKCkKrKHV19dHT49V1dvtdjPf1tdyuVzc60iWePnt4GeiiouL7e1z584FLXSolLpxGqAUAJWnAovzrSsrtJc3h+DeU05OzpyqXj5V924JVJY4drqJVscyHIsWLbI/kytXrtDV1TXOHZRS06UBSgHB6eUbVxYGnXMGqPmSvTdWTkaKPS9ngJf3B5IlPB5P0HLwmiyh1MzQAKUYGh6hujbwEO76skCCxOjoKJcuXbL359v8k9OOrYECsnv3n8LnC5Q3cmbzNTQ06DCfUjNAA5TiZN1F+5dtYW66XT0B4PLly4yMjABWQsCCBQsi0sZocOuaxSxIttLrL3f3BVXdyM3NJT4+HoCBgYGgoK6UujEaoBTHzwRWyF1XVhB0bmz23lwvb3Q9Ho+be7eU2fsvOQrIulyuoJRzHeZT6uZpgFIccwSotcuD55+czz/N1/knJ2c23+HqBjvzEYKfiTp//jxerzesbVNqrtEANc/19A1y7ry1Qq4Aq0sCc0xer5eOjg57Pzc3N9zNizoFOWmsXGp9RqPG8MrB0/a5jIwMUlJSAOuza25uHvceSqmp0QA1z1WdvWCXN1pWlENSQqC8UWtrqz3Zn5aWZs+xzHf3bQv0ovbuP2V/RiIS9EyUDvMpdXM0QM1zQfNPy4Pnn5zDe9p7Cti2YSkJ8bEAXGjr5mRd4HNyDvO1tLQwNDQU9vYpNVdogJrnjp8JPKC79joBSuefAuJiY7jjlmX2/kuOyhIpKSlkZmYCVvmoxsbGsLdPqblCA9Q81tHVy4W2bsDKUFuxNBCEBgcH6e62zokIOTk5EWljtLrP8UzUW0dq6RsI9JR0mE+pmRHSACUiO0XktIjUiMjnxzkfJyLP+M8fEJFi//H7ReSwiBz3//dex2s2+Y/XiMjXZT7nPd8k5/DeiiW5xMZ47H3nczxZWVl4PB5UwLKibIryMwAY8fp443CNfa6oqMhOx29vb6e3tzcibVRqtgtZgBIRN/AN4EFgFfC4iKwac9mTQKcxpgT4GvBV//F24D3GmLXAJ4AfOl7zr8BTQKn/a2eo3sNcN9X0cp1/upaIcN+2QC/KWfooPj5eK5wrNQNC2YPaDNQYY+qMMcPA08CuMdfsAr7v3/45sENExBhzxBhztfZONRDv723lAwuMMW8bK3XqB8B7Q/ge5ixjzIQJEsYYnX+agu3lpXg8VpHY2vNt1DcHVhweO8ynpY+Umr5QBqgCwLkGdpP/2LjXGGO8QDeQOeaaDwBHjDFD/uubHOfGu6eaggtt3VzuthbXS4iPZdmiwBIavb299PdbD6B6PB570l8FS0mKZ8u6QA0+Zy+qsLDQHhbVCudK3ZhQBqjx5obG/hl53WtEZDXWsN8fTuOeV1/7lIhUiEhFW1vbFJo7v1Q5ek9rShbidgd+FJzzTzk5OUFLb6hgOxyVJV49dJbhEat6xNgK5+fOnQt725Sa7UL5m6cJcK4bXghcmOgaEfEAqcBl/34h8CzwcWNMreN652TJePcEwBjzLWNMuTGmfL4tsDcVVTWBj211ycKgczr/NHXrlheQnW5Vj+gbGOLgsXr7nPOZqMbGRh3mU2qaQhmgDgGlIrJERGKBx4DdY67ZjZUEAfAosNcYY0QkDfgt8AVjzJtXLzbGtAA9IrLVn733ceBXIXwPc5IxhmpHgFpTujDoXGtrq72v80/XJyLcuzVQQNY5zJeXl0dcnFWZQyucKzV9IQtQ/jmlzwB7gJPAT40x1SLyZRF5xH/Zt4FMEakBPgdcTUX/DFACfFFEKv1fVx/E+RTwH0ANUAs8H6r3MFc1t3bR3TMAWMu7FxcE5pi6u7vt6gdxcXGkpqZGpI2zSdBqu2cCq+1qhXOlbk5IJxeMMc8ZY5YbY5YZY77iP/YlY8xu//agMeaDxpgSY8xmY0yd//jfGWOSjDEbHF+t/nMVxpg1/nt+xui4ybRVnXEO7wUvoeEc3svJyZnXy2tMVVZ6sr3aLgT3opzZfOfPn8fn84WzaUrNajr7PQ9db/7JObyn809T51yGY9+BU4yOWgtAZmRkkJxsLQCpFc6Vmh4NUPPM9eafRkdHNUDdoM1riklJsqq9d3T1UXnKehpCK5wrdeM0QM0zTZe6uNJrzT8lJ8axeGFg/qmzszNoeferaxupyXk8bu6+dbm9P9Ewn1Y4V2rqNEDNM9Vng4f3nHNMziyz3NxcnX+aph2O0keHqurtRBRnhfPR0VGtcK7UFGmAmmeuN/80NkCp6VmUl87yYutz8/lGeeXQGfucDvMpNX0aoOaRyeafnBU3NEDdmIlW29UK50pNnwaoeWTs/NPV5SIAOjo67BTopKQkkpKSItLG2e72jSXExcYA0HSpk9PnrF5pfHw8+fn59nXai1Jqchqg5pHpzD+pGxMfN2a13f2B1Xa1wrlS06MBah7R+afwcK4T9eY7tfQPDANQUFBgVzjv6enh8uXLEWmfUrOFBqh5whjDidrgChJX+Xw+2tsDaxlpgLo5pYtzWOQfPh0e8fLGO9Zqux6Ph0WLAhUndJhPqevTADVPXGjrDqq/53z+qb293a58kJKSQkJCQkTaOFeICPc5Kku89Pb4w3wNDQ32566UutZ1A5SIfM+x/YnrXKqinHP+adWyievvae9pZtx16/Kg1XbPNVk91NzcXPsPgKGhIVpaWiLWRqWi3WQ9qPWO7T8LZUNUaFU7hvdWLssPOqfljWZeSlI8W9cHVtv93VtWL0pEgiqc60KGSk1ssgClaUZzgDGGEzWBv9RXOwKU1+ulo6PD3tcANXPudyRLvH74LEPDVhmpJUsCgau5uZnh4eGwt02p2WCyAFUoIl8Xkf/PsW1/haOB6uZd6ujhcncfYKVBLynMss+1tbXZ6c6pqan2Anvq5q0uWUh+trWeVv/gMG8dqQMgLS2NtLQ0wHpA+vz58xFro1LRbLIA9VfAYaDCse38UrPACUd6+cqlebjdgf/tml4eOiLCDkeyxO8cyRLOXpQO8yk1Ps/1Thpjvh+uhqjQqa4NDO+tWqbPP4XTPVvK+PFvDzE6OsrpcxdpbLlMUX4GixcvprKyEmMMbW1t9Pb22utGKaUsk2Xx7b7eV7gaqW6Oswe1yjH/NDw8HPSwaE5OTljbNR+kpSSyeU0gKeLlt61lOBISEsjLy7OP6zNRSl1rsiG+bUAh8Drw98A/jPlSUa69s5fWyz0AxHjclBRl2+ecxWHT09OJjY0Ne/vmg/tvX2Vv7zt4muERL3DtMJ+WPlIq2GQBKg/4H8Aa4J+A+4F2Y8yrxphXQ904dfOc1ctXLM2zn80BHd4Ll/VlheRkWIs/9g0M8XallSxRWFhITIxVWLa3tzeomodSapIAZYzxGWNeMMZ8AtgK1ACviMifhKV16qadCJp/Cn7+yRmgdHgvdESE+28L9KL2vHkCALfbHVT6SJMllAo2aakjEYkTkfcDPwI+DXwd+EWoG6ZmxkTzT0NDQ3R1dQHWL1ANUKF179YyXC7rn9vpcxdpuGDN/S1dutS+prGx0V7yRCk1eZLE94G3gFuA/2WMudUY87fGmOawtE7dlM4r/Vxo6wbA7XbZq71CcPWI9PR0e6hJhUZaSiJb1jkrS1i9qKysLDt7b2RkhKampoi0T6loNFkP6mPAcqwyR2+LyBX/V4+IXAl989TNcA7vLV+cS2xM4KkCLW8Ufu+6LVBZ4tVDZxgcGkFEgnpRtbW1kWiaUlFpsjkolzEmxfG1wP+VYoxZMNnNRWSniJwWkRoR+fw45+NE5Bn/+QMiUuw/niki+0SkV0T+ecxrXvHfs9L/pWNTEzh5nfknZ4FYZ7qzCp21ywuCKku8ecRahsNZ4fzSpUv09/f9REOMAAAgAElEQVRHonlKRZ3JhvjiReSzIvLPIvKUiFz3wd4xr3UD3wAeBFYBj4vIqjGXPQl0GmNKgK8BX/UfHwS+CPzlBLf/iDFmg/+rdYJr5j1nBt8qx/pPg4ODXLlidYBdLhdZWVnXvFbNvLHJEi++aVWWSEpKCurFarKEUpbJhvi+D5QDx4GHmN6zT5uBGmNMnTFmGHga2DXmml3+7wHwc2CHiIgxps8Y8wZWoFI3oKdvkMYWayLeJcKKJYFekjN7LzMz017lVYXePZuX26WmahpbqW20nkVzDvPV1dXpM1FKMXmAWmWM+agx5t+AR4E7p3HvAsBZBbPJf2zca4wxXqAbyGRy3/UP731RnAsbOfh7fBUiUuF8IHW+OFkXGMJbuiib+LhAEoSml0fOguQEbt+4zN5/4Y1qwHom6uofCvpMlFKWyQLUyNUNfwCZjvECx9g/C6dyzVgfMcasxQqWd2Ilclx7E2O+ZYwpN8aUZ2dnj3fJnHa9+SdngNL5p/Dbecdqe/v1w2fp7R/C4/FQVFRkH6+rq4tE05SKKpMuWOjM3APWTSOLrwlY5NgvBC5MdI1/fisVuMx1XE1xN8b0AD/GGkpUYzjnn1aXBgrE9vf309vbC1gPimZmTqXDqmbS8uJciguseb8Rr499B04DsGxZoGfV2NjIyMjIuK9Xar6YLIvPPSZzzzONLL5DQKmILBGRWOAxYGyB2d3A1aXkHwX2musMvouIR0Sy/NsxwMNA1STtmHcGBoepO28NawpMOP+UlZWF2+0e+3IVYiLCzjuclSWqMcaQmZnJggXWPyuv10tjY2OkmqhUVJi0ksSN8g8JfgbYA5wEfmqMqRaRL4vII/7Lvg1kikgN8DnATkUXkXrgH4EnRKTJnwEYB+wRkWNAJdAM/Huo3sNsdercJXuctGhhJsmJgUUItf5edLhzUykJ8VZx3pa2bo6daUZEgnpR+kyUmu9Cmr5ljHkOeG7MsS85tgeBD07w2uIJbrtppto3Vznnn1Y70suNMRqgokR8XAz3bF7Oc69ZAwAvvF7F+rJCiouLOXr0KKOjo3R0dNDd3U1qamqEW6tUZISsB6Uipyqo/l5g/qmvr89+CNTj8ZCRkRH2tqmAnXeusbcPHa+n7XIP8fHxFBQEkl21F6XmMw1Qc8zwiJeaxsCzy84MPmfvKTs72y5eqiKjICeNtcutYGSAPf6Uc+czUfX19YyOjkaieUpFnP6GmmNOn7uEz2f9QivISSM1JcE+5yxvpMN70eGh7Wvt7d+9fZLhES/5+fkkJiYCVtV5LSCr5isNUHNMde346eVj55/0+afoUL66iOx0azHD3v4hXj98FhEJWm1Xh/nUfKUBao45UeNIkHDMP3V3dzM0NARAbGwsaWlpYW+bupbL5eLB7YG5qN+8chxjTFA238WLF+1n15SaTzRAzSHDI15O1wd6Sc4CsWOH9yaoEKUiYMfWFfZSKI0tlzlR20JSUhL5+YH/fzU1NZFqnlIRowFqDjnb0IrXa63Imp+dSkZqkn1Oh/eiV3JiHHdvXm7vP/fqcQBKSkrsY3V1dZosoeYdDVBzSPUEy7uPjo7qAoVRzpksceDYOVov97Bw4UJNllDzmgaoOcS5gu4aR4JER0cHXq9V6zcxMdFeYlxFj0V56axbXghYKefPv1aFy+UKSjk/e/ZshFqnVGRogJojvF4fpxxLbDgf0B07vKfzT9Hp3XcHp5z3DwyzbNky+/9Xa2urvdCkUvOBBqg5oqaxjRH//FNu5gKy0gO9JH3+aXbYtKqIhf4l4QcGh3l5/ykSExNZuDDwx4amnKv5RAPUHOF8/smZvef1euno6LD3NUEieokI77lnvb3/21eP4/ONXpMs4fP5ItE8pcJOA9QcMdHzT21tbXb2V2pqKvHx8WFvm5q6u24ttavPt3X2sP/YOfLz80lKsjIyh4eHqa+vj2ALlQofDVBzgNfrC1ri3VlBQof3Zpe42JigIrK79x4FoLS01D525swZrrNsmlJzhgaoOaCuqZ2hYWv11ez0FHIyUuxzzgClw3uzw847VuN2W/80axpbOX3uEkuXLrUXl+zq6qK9vT2STVQqLDRAzQFVZ8effxocHKSrqwuw5jdycnLC3jY1fekLEtleHugx/fLlSuLi4iguLraPnTlzJgItUyq8NEDNAVVnm+3ttaWBtYScvaesrCxiYmLC2i514x5xJEscqqrn/MVOli8PVJs4f/48AwMDkWiaUmGjAWqWu978U0tLIHFCh/dml6L8DMpXL7b3f7W3krS0NLKzswGrOr0+uKvmOg1Qs1xNYxvDI1aViNzMBfb8kzEmqAflLDyqZof33bfR3n6t4iztnb1Bvaja2lpNOVdzmgaoWe7YmUB9ttUlwctrDA4OAtbyGrq8++yzYmkeK5daf1j4fKP8et8xCgsLSUiwFqEcHBykoaEhkk1UKqQ0QM1yzgKxa5ePP7yny2vMXu+7b4O9/bu3T9I3MByUcn7q1ClNOVdzlgaoWWx4xMupc4E6e2scCRLOAKXDe7PXLauKKMq3er9DwyM8/3oVJSUldsp5d3d3UK1FpeYSDVCz2Olzl+z1nxY61n/y+Xy0tbXZ12mAmr1EhPc75qJ+++pxRo0EVTk/depUJJqmVMiFNECJyE4ROS0iNSLy+XHOx4nIM/7zB0Sk2H88U0T2iUiviPzzmNdsEpHj/td8Xebx2FWVY3hvzfJA76m1tdUub7RgwQJ7TSE1O922cRm5mQsA6O0f4vnXqykrK7PPt7S00N3dHanmKRUyIQtQIuIGvgE8CKwCHheRVWMuexLoNMaUAF8Dvuo/Pgh8EfjLcW79r8BTQKn/a+fMt352OH4m8PzTmgmef9L08tnP7XbxgXcFelG79x0lJjaewsJC+5j2otRcFMoe1GagxhhTZ4wZBp4Gdo25Zhfwff/2z4EdIiLGmD5jzBtYgcomIvnAAmPM28aaGf4B8N4QvoeoNTg0Qk1jYJXcNSXjJ0jo8N7ccFf5crLTrUcIevoGeeGNalasWGGfr6+v1wd31ZwTygBVAJx37Df5j417jTHGC3QDmZPc07nu9Xj3BEBEnhKRChGpcM7HzBWnzl3E57OG8RblZ5CaYqUe9/f328M9LpdLyxvNER6Pm/ffH+hF/WrvUVIWpNqPD4yOjmr5IzXnhDJAjTc3NDYfdirX3ND1xphvGWPKjTHlV5++n0uqzjjLGwV6TxcuBOalsrOz8Xg8YW2XCp17t5TZC1Fe6R3gxbdOsnLlSvv82bNnGR4ejlTzlJpxoQxQTcAix34hcGGia0TEA6QClye5Z6Fjf7x7zguVpwMdybXLAx+JM0A5V2JVs5/H4w7K6Pvly5Xk5OaRkmIN/Y2MjGj5IzWnhDJAHQJKRWSJiMQCjwG7x1yzG/iEf/tRYK+5zlOHxpgWoEdEtvqz9z4O/Grmmx7dunsGONdkLbfgEmF1ydVqA76gZ2I0QM09925ZQWaa9ThBd88Az79ezapVgdyj06dP4/V6I9U8pWZUyAKUf07pM8Ae4CTwU2NMtYh8WUQe8V/2bSBTRGqAzwF2KrqI1AP/CDwhIk2ODMBPAf8B1AC1wPOheg/Rypm9V1qcS1KCfwXWtjb7l1NycrL9l7WaO2Ji3Dz6rk32/rMvVZKdm28/SjA0NERtbW2kmqfUjArpBIUx5jnguTHHvuTYHgQ+OMFriyc4XgGsGe/cfHHUMby3vmzi4b15/IjYnHbvljJ+tbeSi+1X6BsY4jevHGfTypUcPnwYsFLOS0tLcbn0OXw1u+lP8CxjjOHo6UBypDNANTcHelY6vDd3eTxuHnvwVnv/168cJyt3IXFxVk+6v7+fc+fORap5Ss0YDVCzTHNrFx1dfQAkxMdSuthKI+/p6aG3txcAt9ut6eVz3B2bSoJq9P3y5aNBz0VVV1fb1USUmq00QM0yR085svdKF+J2W/8LncN7eXl5djFRNTeJCI+/e7O9/8Ib1aRl5hEbGwtAX1+f9qLUrKcBapY5FjT/FMji1+G9+efWNYvtHrTPN8rPXjwS9FxUVVWVLmioZjUNULOI1+vj+NlAT2ldWYH/uDeoerkGqPlBRPjYI1vt/dcqzuKKTw+ai6qrq4tU85S6aRqgZpEzDa0MDY8AkJ2eQn52KmDV3rs635CamqrVy+eR1SUL2by22N7/0a8PXjMXpb0oNVtpgJpFgtLLVxTaaeRNTYHjBQXjliZUc9hHH9lqp5SfrGuhczCG+Ph4AAYGBqipqYlk85S6YRqgZpGjpwLp5ev86eWjo6NB80+LFi265nVqbivISePBO1fb+z/+7SGWlwX3okZGRiLRNKVuigaoWaK7Z4CaBmt5DQHW+RcobG1ttX/5JCYmkp6eHqkmqgj64AObSIy3Mvgutl+hpmUoqLrEyZMnI9k8pW6IBqhZ4sjJRrtse9nSPFKSrCGc8+cDvarCwkKtHjFPpSTF86Gd5fb+z158h6Ilpfb+qVOndL0oNetogJolKqob7e1NqxYDVlUJ5/Cec4VVNf88eOdqCnLSABgYHOa1oy2kpVn7Pp+PY8eORbJ5Sk2bBqhZwOv1UemYf9q02gpQHR0d9l/FsbGxzMV1r9TUeTxunnz0Dnv/tYqzJKYHkmbq6uro6uqKRNOUuiEaoGaBU+cuMjBoLUSXlZ5MUb41zzQ2e0+Lg6r1ZYVs27DM3n/2lVPk5OTa+5WVlZFollI3RH+jzQKHxwzviQjGmKD5J83eU1c98d5txMXGAHC+5TJtAwn2uZaWlqCyWEpFMw1Qs8Dh6gZ7e9PqIgCuXLliF4f1eDzk5eVFpG0q+mSlJ/PBB26x93e/eoL0rHx7//Dhw/rwrpoVNEBFuZa2bppbrXmDGI+btf70cmfvKT8/X4vDqiDvuXsdhbnWUPDQ8Ahvn+rG47GWf+vt7eXUqVORbJ5SU6IBKsq9cyIwvLdueSGxMR6MMTQ0BHpVOrynxvJ43Hz6w3dz9aGDqpoWBt0Z9vnq6mr6+voi0zilpkgDVJQbb3ivq6uLK1euANbaT1reSI1neXEu775rnb3/u4MNeOKSACvt/J133olU05SaEg1QUWxgcJiqmsCE9i2rrABVX19vHyssLLSHbpQa6/F330p2egoAfQPDVDV57XNNTU2aMKGimgaoKFZR1YDPZ1UpX7wwk+yMFIwxNDYGhv2Ki4sj1Do1G8THxfCpx++y96tqW+kaDmT1HTp0SOv0qailASqKvX00sJbPbRutZ1va2tro7+8HrIdzNXtPTWZ9WSE7tgaKx75+vI1BrzU71d/fr89GqailASpKDQ6NBCVIbNuwFAge3isqKtKHc9WU/P77biMvawEAQyM+jtQPM2qs6o41NTVcunQpks1Talz62y1KVVQ3MOK1nlUpys+gICeN0dHRoPRyHd5TU5UQH8uffWwHrqtriLX1U99u7PMHDx7E6/VO9HKlIiKkAUpEdorIaRGpEZHPj3M+TkSe8Z8/ICLFjnNf8B8/LSIPOI7Xi8hxEakUkYpQtj+S3j5Sa29f7T21tLQwPGyVPEpMTCQrKysibVOz0/LiXB59YBMAIvBOzRXar1jzT729vRw9ejSSzVPqGiELUCLiBr4BPAisAh4XkVVjLnsS6DTGlABfA77qf+0q4DFgNbAT+Bf//a66xxizwRhTzhw0ODTC4aDhPWv+yTm8t3jxYl1aQ03bo++6hdLFOdaOCG+d7mNgyOo5nTlzRrP6VFQJZQ9qM1BjjKkzxgwDTwO7xlyzC/i+f/vnwA6xfuvuAp42xgwZY84BNf77zQuHTzTaw3uFueksyktnaGgoqDisDu+pG+F2u/jsx++zFzcc8gpvneqx56P279/P4OBgJJuolC2UAaoAOO/Yb/IfG/caY4wX6AYyJ3mtAV4UkcMi8tRE31xEnhKRChGpaGtru6k3Em5vVway97b6h/fOnTvH6KiVcp6RkWGv86PUdOVlLeBPP3YvYA31tfcaKmu7AWv13f3792OMud4tlAqLUAao8cafxv7UT3TN9V57uzHmFqyhw0+LyPbxvrkx5lvGmHJjTPlsWidpaHgkqHrEbRuWYYyhpqbGPlZSUhKJpqk55NY1xbxvxwbAqkZysnmY+otW8eGWlhbOnDkTyeYpBYQ2QDUBziJxhcDYAW77GhHxAKnA5eu91hhz9b+twLPMsaG/iupGhkesOYGF2akU5afT2tpKT08PYFUuX7x4cSSbqOaIx9+9mTWlCwGIi4vjzZNX6LgyBMCRI0eYbSMPau4JZYA6BJSKyBIRicVKetg95prdwCf8248Ce401trAbeMyf5bcEKAUOikiSiKQAiEgS8C6gKoTvIez27g9Umb59UwkiEtR7Ki4u1tJGaka43S7+/BP3kZlm1eeLjU/khYpL9A95Mcbwxhtv2Cs2KxUJIQtQ/jmlzwB7gJPAT40x1SLyZRF5xH/Zt4FMEakBPgd83v/aauCnwAngBeDTxhgfkAu8ISJHgYPAb40xL4TqPYRbe2cvR/1Luwtw75YVDA4OBiVH6PCemklpKYn8j6ceJC42BhHBHZvM8wdb8PpGGRwc5PXXX9e1o1TEhPQ5KGPMc8aY5caYZcaYr/iPfckYs9u/PWiM+aAxpsQYs9kYU+d47Vf8ryszxjzvP1ZnjFnv/1p99Z5zxb6Dp+2JtrXLC8nJSAlKjsjMzCQ9PT1yDVRzUnFBFp974j4Ea5mOwdE4Xqy4wKgxdHR0cPjwYU2aUBGhlSSihDEmaHhvx9YVmhyhwqZ89WI+8d7bAGs+6tIVYd+RixhjqK2t1QUOVURogIoSVWcv0HrZSoRIjI9l87pimpub7WXdY2JiKCoqimQT1Rz38N1reWj7GgCSk5OpvTjIG1WtGGOorKwMWiRTqXDQABUl9h4I/IW6vbyUGI+b6upq+1hJSYkmR6iQEhH+2/tvZ3t5KSKQnp5OVX0Ph053ANZDvFpUVoWTBqgo0DcwFPRw7o6tK7h06RKXL18GwOVyUVZWFqnmqXlERPjMh+/h1jXFiAiZmZkcOn2ZAyfb8Pl8vP766/bPpVKhpgEqCrxeUWOXNiouyGLpomxOnDhhn1+6dCkJCQkTvVypGeV2u/jcE/exdnkBLpeLzMwsDp2+zFvVbQwPD7Nv3z4NUiosNEBF2OjoKL9+JVBF+t4tZXR0dNhDKSLCypUrI9U8NU/Fxnj4wid3sr6sEI/HTWZmJu+cvcyrRy8xNDTEvn376OzsjHQz1RynASrC9h87x8X2K4CVHHHvlrKg3lNRURHJycmRap6ax+JiY/jCJx/k1jXFxMTEkJmZyfFzXTx/sJn+gUH27t1LR0dHpJup5jANUBFkjOHZlwLLbT945xqGBvuDHsxdtWrsCiVKhU9MjJu//P372bZhGbGxsWRmZlLT3Mt/vdZAV08/e/fu1SU6VMhogIqg42eaqTtv1TuL8bh5aPsa3nnnHfv8woULtWq5ijiPx83nPrGDh7avITY2lqysTFq7hnh67zkudvTy2muvUVdXN/mNlJomDVAR5Ow93btlBb1XLgel8a5fvz4SzVLqGi6Xiyc/cAf/7f23ExcbS1ZWFv1Dozyzr56qc50cOHCAY8eOacUJNaP0wZoIqW1s49gZayjPJcLDd6+hYv/r9vnS0lLtPamo8+671pKbtYB//N5LuFwu2ts7eOnwBS609zPiHeXy5cts27aNuLi4SDdVzQHag4qQn+05bG/fdssyutpb6OvrAyA2Npa1a9dGqmlKXVf56sV89S/ez+KFWWRnZxEXF8eJhi5+8nId71TV8sILL2gaupoRGqAioPLUeQ5V1dv7D96+MqhqxNq1a/UvUBXVFuWl89W/eD93by4jMzOT5ORkOnuH+dkr9bx44BzPv7CH6upqu9CxUjdCh/jCbGTEx7d//oa9f1d5KZeaavB6rUUKFyxYoEVh1awQHxfDn370XtaULuQ7v3iL2NhYOjs7OXymg9oLPTRe6qW5uZmtW7eyYMGCSDdXzUIaoMJs9ytHudDWDUBCfCzbVmVRcybw3NOmTZtwubRjq2YHEWHH1pWsXV7Iv/zkFd6pbqCzs5Ou3mF+9WYjx+o6OXf+EttuXc+qVau0nqSaFv1pCaO2yz387IXA3NPDd66griZQJLasrIy8vLxINE2pm5KTkcL//OOH2fPGCX706/20tl/mypUezrX00HCpl8On29i+8TTbb9/CokWLEJFIN1nNAhqgwsQYw3d+8aZdc68wN414bxsD/jH6jIwMNmzYEMkmKnVTRISdd65m64Yl/OjXB9jzehVdXV2MjAxzrK6T6vou3jjaxL23LuOObeXk5+droFLXJfPhuYXy8nJTUVER0TY8+9IRfvTrA4AVrN57eyExo9ZaTx6Ph507d5KSkhLJJio1o842XOK7v3iLd6rruHLlip0wISKUFqRw963LuH/7rRQWFmqgmmdE5LAxpnzS6zRAhd7B4/X87/94AYMVnJbmxrJpaaJ9/rbbbmPx4sURa59SoWKMofJUE//56/1Unqijr68v6GHe7LR4Ni7P5ZH7bmXt6hWavTpPTDVA6RBfiDVc6OD//uBlOzilxPrYUBxYOmPVqlUanNScJSJsXLmIDSsKOXyikWdfPMz+yrP+Z/4MbV2DvHiwgZcqGlmSl8KW9Ut4YPtGSpYW43a7I918FWHagwqh8xc7+co3n6Otswefz4cZGWDXtlwS4qy/C1auXMn69et1eEPNK+ea2vnFixX87q0qrlzpueZZKZdLWJSTzMaVRdxRvoLy9SuIj4+PUGtVKOgQn0MkAtSbR2r5xo9fYWh4hIGBAfp6r/Do9iKyUq1/aBqc1HzX2z/Ea4dO86uXKjhR08zIyPC41yXGeVhamMGa0kJuWbOUDauXsSBFl6CZzTRAOYQzQA0MDvPTFw6ze99RhodH6O3tYWR4iHdvLaQ4z/pHtXr1atauXavBSSm/5tYu9r1dzctvHedM/UV8Pt+E14oIOemJFBdksnRRDksW5VJavJBli/OJjY0JY6vVjYqKACUiO4F/AtzAfxhj/t8x5+OAHwCbgA7g94wx9f5zXwCeBHzAnxpj9kzlnuMJR4BqvdzD869VseeNKjq7e+jr62N4eJjUpFge3lZIVmo8iYmJbNmyRZ91Uuo6Wi9f4e13zvBGxSmOnT5Pd0//lF4nQGpKAtnpyeRkppCVnkJ2xgJyslLJTEshOyOVzIwFpKUk4Xbrw/CRFPEAJSJu4AxwP9AEHAIeN8accFzzx8A6Y8wfichjwPuMMb8nIquAnwCbgYXAS8By/8uue8/x3GiAqqm/QEurVfTSGIMx4PX66BscYmBwhK6efs6db6eh5TKX2nsYGh4K+stvSV4yD9xaQFysm8WLF1NeXk5sbOy026HUfDU6Okrd+TYOHj1D5YlznKm/xKWOnptc1kOIi/WQGB9DQnwM8bExJMR5iI+LJT7OQ1ysh7jYGGI8bmJjPMTEuIn1eHB73MR43HjcLjweNy6XC4/bhcvlwu0SRASXS3C7XCCCyz9C4na5sAdLRBABITB64hxImU1VZJYvKSA3O/2GXhsNWXybgRpjTJ2/QU8DuwBnMNkF/I1/++fAP4s17rULeNoYMwScE5Ea//2Ywj1nzI9++SovvDH9W6cnx7KxNJO1S9MpKiqirKyMrKysELRQqbnN5XJRsjiXksW5fPiROwHoHxii6nQ91WfPU9/UyvmWy1xou0J3zwBTC1uGoeERhoZH6LwSytbPbZ974n4+9PD2kH6PUAaoAuC8Y78J2DLRNcYYr4h0A5n+4/vHvLbAvz3ZPQEQkaeApwCKiopu7B1Mg4hQmJ3IpuXZbFpdRH5+PsuWLSMhIWHyFyulpiwxIY7NG8rYvKEs6Pjg0DCNzZdoaG6jpfUyrR1X6OjspaO7l57eQXr6hugdGGZo2DvFQKYiLZQBarwMgLE/FxNdM9Hx8fq/4/6sGWO+BXwLrCG+iZs5sYLcDEoWZYEEGuR2CfFxMfYQQWFuGsUFWSxblENWVgYZGRmzqpuu1FwRHxfL8qWLWL500YTXGGPwer109/TTdaWXnt5+evsG6e0fpH9gyMq6HRxmcHiEkREfwyNehke8eL0+vL5RRrw+RkdHGR01eH2j9vaosb6uTgWMjvq3/d/z6m8p49+4OkI50VDlbAigWemhr1AfygDVBDh/UgqBCxNc0yQiHiAVuDzJaye754z55OMP8MnHHwjV7ZVSYSYixMTEkJWRSlZGaqSboyYRyj/1DwGlIrJERGKBx4DdY67ZDXzCv/0osNdYf1LsBh4TkTgRWQKUAgeneE+llFJzQMh6UP45pc8Ae7BSwr9jjKkWkS8DFcaY3cC3gR/6kyAuYwUc/Nf9FCv5wQt82hjjAxjvnqF6D0oppSJHH9RVSikVVlNNM9fZfKWUUlFJA5RSSqmopAFKKaVUVNIApZRSKippgFJKKRWV5kUWn4i0AQ03cYssoH2GmjNb6WegnwHoZwD6GcDNfwaLjTHZk100LwLUzRKRiqmkRM5l+hnoZwD6GYB+BhC+z0CH+JRSSkUlDVBKKaWikgaoqflWpBsQBfQz0M8A9DMA/QwgTJ+BzkEppZSKStqDUkopFZU0QCmllIpKGqAcRGSniJwWkRoR+fw45+NE5Bn/+QMiUhz+VobWFD6Dz4nICRE5JiIvi8jiSLQzlCb7DBzXPSoiRkTmXMrxVD4DEfmQ/2ehWkR+HO42htoU/i0Uicg+ETni//fwUCTaGSoi8h0RaRWRqgnOi4h83f/5HBORW2a8EcZepnh+f2GtL1ULLAVigaPAqjHX/DHwTf/2Y8AzkW53BD6De4BE//an5uNn4L8uBXgN2A+UR7rdEfg5KAWOAOn+/ZxItzsCn8G3gE/5t1cB9ZFu9wx/BtuBW4CqCc4/BDwPCLAVODDTbdAeVMBmoMYYU2eMGQaeBtdFPBIAAAZXSURBVHaNuWYX8H3/9s+BHSIiYWxjqE36GRhj9hlj+v27+4HCMLcx1KbycwDwt8D/BgbD2bgwmcpn8EngG8aYTgBjTGuY2xhqU/kMDLDAv50KXAhj+0LOGPMa1kKyE9kF/MBY9gNpIpI/k23QABVQAJx37Df5j417jTHGC3QDmWFpXXhM5TNwehLrL6i5ZNLPQEQ2AouMMb8JZ8PCaCo/B8uB5SLypojsF5GdYWtdeEzlM/gb4KMi0gQ8B/xJeJoWNab7+2LaQrbk+yw0Xk9obA7+VK6Zzab8/kTko0A5cFdIWxR+1/0MRMQFfA14IlwNioCp/Bx4sIb57sbqRb8uImuMMV0hblu4TOUzeBz4njHmH0RkG/BD/2cwGvrmRYWQ/z7UHlRAE7DIsV/ItV12+xoR8WB166/XBZ5tpvIZICL3AX8NPGKMGQpT28Jlss8gBVgDvCIi9Vhj77vnWKLEVP8t/MoYM2KMOQecxgpYc8VUPoMngZ8CGGPeBuKxiqjOF1P6fXEzNEAFHAJKRWSJiMRiJUHsHnPNbuAT/u1Hgb3GP1s4R0z6GfiHt/4NKzjNtXkHmOQzMMZ0G2OyjDHFxphirHm4R4wxFZFpbkhM5d/CL7ESZhCRLKwhv7qwtjK0pvIZNAI7AERkJVaAagtrKyNrN/BxfzbfVqDbGNMyk99Ah/j8jDFeEfkMsAcrg+c7xphqEfkyUGGM2Q18G6sbX4PVc3osci2eeVP8DP4PkAz8zJ8f0miMeSRijZ5hU/wM5rQpfgZ7gHeJyAnAB/yVMaYjcq2eWVP8DP4C+HcR+XOsoa0n5tIfrCLyE6wh3Cz/PNv/BGIAjDHfxJp3ewioAfqB35/xNsyhz1MppdQcokN8SimlopIGKKWUUlFJA5RSSqmopAFKKaVUVNIApZRSKippgFLKT0R8IlIpIlUi8jMRSZzm63unef33ROTRcY6Xi8jX/dtPiMg/+7f/SEQ+7ji+cDrf7zrtuNNfkbxSRBIcx4uvU8n6y/4HthGRz073s1JqKjRAKRUwYIzZYIxZAwwDf+Q86X8gMeT/ZowxFcaYPx3n+DeNMT/w7z4BzEiAAj4C/L3/vQ9MsY1fMsa85N/9LKABSs04DVBKje91oMTfizgpIv8CvAMsEpHHReS4v6f1VeeLROQfROQdsdbKyvYf+6SIHBKRoyLyX2N6G/eJyOsickZEHvZff7eIXFOIVkT+RkT+0t/rKgf+09/rebeIPOu47n4R+cU4r98h1tpFx8Va6ydORP4A+BDwJRH5z3E+B7eI/Lu/h/Xi1R7W1d6fiPwpVqDcJyL7pvUJKzUJDVBKjeGvs/ggcNx/qAxrWYGNwAjwVeBeYANwq4i8139dEvCOMeYW4FWsJ+8BfmGMudUYsx44iVXD7apirIK77wa+KSLxk7XPGPNzoAL4iDFmA9YT/SuvBkSsJ/q/O+Y9xQPfA37PGLMWq4rMp4wx/4FVsuavjDEfGefblWItq7Ea6AI+MKYtX8eqv3aPMeaeydqu1HRogFIqIEFEKrF++TdilbYCaPCvdwNwK/CKMabNv+TKf2It7AYwCjzj3/4RcId/e42/l3QcazhtteN7/tQYM2qMOYtVy27FdBvtL6/zQ6ylH9KAbVy7DEoZcM4Yc8a//31Hu6/nnDGm0r99GCugKhUWWotPqYABf4/E5q832Oc8NI37Xa0j9j3gvcaYoyLyBFZ9s7HXTLQ/Vd8Ffo21gOLP/MHT6UYX1nRWq/cBCRNdqNRM0x6UUtNzALhLRLJExI21JtCr/nMurCr3AB8G3vBvpwAtIhKD1YNy+qCIuERkGdby4qen2I4e/30BMMZcwBpq+3+wAuJYp4BiESnx73/M0e6bFdQWpWaK9qCUmgZjzP/f3h3aIAwFYRz/LmEFwhBMwADMgMaxAoI9CAkJhgEwIJmA4EiQqE7QUHeIO4EgFIF4lP9Ptk+0pl8u9/quMrO5pKOiKtm7+y5v15KGZnZSTFue5PWFIthuir7W88f8qgiKgaSZuzdZtbXZKHpWd0mj3H23ldR398uL527MbKo4hb6nGCex/PzN31pJOphZRR8K38Rp5kBH5P9SZ3dfty4GfgABBXRAVm21pHEHpxzjTxFQAIAisUkCAFAkAgoAUCQCCgBQJAIKAFAkAgoAUKQHBjBFKfA6J3EAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Solution\n", "\n", "# Now let's run it with the specified prior and \n", "# see what happens when we multiply the convex prior and \n", "# the concave posterior.\n", "\n", "pmf = Beta(5, 10).MakePmf()\n", "blaster = AlienBlaster(pmf)\n", "thinkplot.Pdf(blaster, color='gray')\n", "blaster.Update(2)\n", "thinkplot.Pdf(blaster)\n", "thinkplot.decorate(xlabel='Probability of hit',\n", " ylabel='PMF')" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.3333333333333333, 0.3175635713858314)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "# The posterior mean is lower\n", "\n", "prior.Mean(), blaster.Mean()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.3076923076923077, 0.28)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solution\n", "\n", "# So is the MAP\n", "\n", "prior.MAP(), blaster.MAP()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# So if we learn that the new design is \"consistent\",\n", "# it is more likely to be consistently bad (in this case)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part Two" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose we\n", "have we have a stockpile of 3 Alien Blaster 9000s and 7 Alien\n", "Blaster 10Ks. After extensive testing, we have concluded that\n", "the AB9000 hits the target 30% of the time, precisely, and the\n", "AB10K hits the target 40% of the time.\n", "\n", "If I grab a random weapon from the stockpile and shoot at 10 targets,\n", "what is the probability of hitting exactly 3? Again, you can write a\n", "number, mathematical expression, or Python code." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.23054197320000014" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "k = 3\n", "n = 10\n", "x1 = 0.3\n", "x2 = 0.4\n", "\n", "0.3 * binom.pmf(k, n, x1) + 0.7 * binom.pmf(k, n, x2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The answer is a value drawn from the mixture of the two distributions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Continuing the previous problem, let's estimate the distribution\n", "of `k`, the number of successful shots out of 10. \n", "\n", "1. Write a few lines of Python code to simulate choosing a random weapon and firing it.\n", "\n", "2. Write a loop that simulates the scenario and generates random values of `k` 1000 times. \n", "\n", "3. Store the values of `k` you generate and plot their distribution." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def flip(p):\n", " return np.random.random() < p\n", "\n", "def simulate_shots(n, p):\n", " return np.random.binomial(n, p)\n", "\n", "ks = []\n", "for i in range(1000):\n", " if flip(0.3):\n", " k = simulate_shots(n, x1)\n", " else:\n", " k = simulate_shots(n, x2)\n", " ks.append(k)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what the distribution looks like." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1000, 3.784)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE5JJREFUeJzt3X+sX3ddx/Hny9YNhYDTNRrWQQsUZYJseikoCgrbKGI2oltoFZxmZoIWUUQzdBlal8iYGDRbdJVVp0DnHCTeaHEibMQf2da7H250s9IV3K7FcLUTFXCj69s/vqfuu8u3vbdrT89nvc9HctNzPudzznn3pO2rn/M9389JVSFJUmu+ZugCJEmaxICSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNWn50AUcLSeffHKtWrVq6DIkSQu4/fbb/72qVizU77gJqFWrVjEzMzN0GZKkBST5l8X08xafJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUm9BlSSdUl2JtmV5OIJ29+e5N4kdyf5eJJnj217NMld3c90n3VKktrT22PmSZYBVwFnAbPA9iTTVXXvWLc7gamq+lKStwDvAd7QbftyVZ3eV32SpLb1OYJaC+yqqt1V9QhwHXDueIequqmqvtSt3gKs7LEeSdKTSJ8BdQrw4Nj6bNd2MBcCHx1bf0qSmSS3JHn9pB2SXNT1mZmbmzvyiiVJzehzJolMaKuJHZM3AlPAK8ean1VVe5I8B/hEknuq6v7HHaxqM7AZYGpqauKxpSO18bKtg5z3yks2DHJeqRV9jqBmgVPH1lcCe+Z3SnIm8KvAOVX18IH2qtrT/bobuBk4o8daJUmN6TOgtgNrkqxOcgKwHnjc03hJzgCuZhROnx9rPynJid3yycDLgfGHKyRJx7nebvFV1b4kG4EbgWXAlqrakWQTMFNV08AVwNOAP0sC8EBVnQO8ALg6yX5GIfrueU//SZKOc73OZl5V24Bt89ouHVs+8yD7/QPwoj5rkyS1zZkkJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNWj50AdLBbLxs6yDnvfKSDYOcV9LjOYKSJDXJgJIkNcmAkiQ1qdeASrIuyc4ku5JcPGH725Pcm+TuJB9P8uyxbRck+XT3c0GfdUqS2tNbQCVZBlwFvBY4DdiQ5LR53e4EpqrqO4AbgPd0+34j8C7gpcBa4F1JTuqrVklSe/ocQa0FdlXV7qp6BLgOOHe8Q1XdVFVf6lZvAVZ2y68BPlZVe6vqIeBjwLoea5UkNabPgDoFeHBsfbZrO5gLgY8ezr5JLkoyk2Rmbm7uCMuVJLWkz4DKhLaa2DF5IzAFXHE4+1bV5qqaqqqpFStWPOFCJUnt6TOgZoFTx9ZXAnvmd0pyJvCrwDlV9fDh7CtJOn71GVDbgTVJVic5AVgPTI93SHIGcDWjcPr82KYbgbOTnNQ9HHF21yZJWiJ6m+qoqvYl2cgoWJYBW6pqR5JNwExVTTO6pfc04M+SADxQVedU1d4kv8Eo5AA2VdXevmqVJLWn17n4qmobsG1e26Vjy2ceYt8twJb+qpMktcyZJCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTVo+dAGSFrbxsq2DnPfKSzYMcl4JHEFJkhplQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkprk96A0kd+7kTQ0R1CSpCYZUJKkJvUaUEnWJdmZZFeSiydsf0WSO5LsS3LevG2PJrmr+5nus05JUnt6+wwqyTLgKuAsYBbYnmS6qu4d6/YA8BPAOyYc4stVdXpf9UmS2tbnQxJrgV1VtRsgyXXAucD/B1RVfbbbtr/HOiRJT0J93uI7BXhwbH22a1uspySZSXJLktdP6pDkoq7PzNzc3JHUKklqTJ8BlQltdRj7P6uqpoAfBd6X5LlfdbCqzVU1VVVTK1aseKJ1SpIa1GdAzQKnjq2vBPYsdueq2tP9uhu4GTjjaBYnSWpbnwG1HViTZHWSE4D1wKKexktyUpITu+WTgZcz9tmVJOn411tAVdU+YCNwI3AfcH1V7UiyKck5AElekmQWOB+4OsmObvcXADNJ/hG4CXj3vKf/JEnHuV6nOqqqbcC2eW2Xji1vZ3Trb/5+/wC8qM/aJEltcyYJSVKTDChJUpMMKElSkwwoSVKTDhlQSf5obPmC3quRJKmz0AjqxWPLb+uzEEmSxi0UUIczNZEkSUfNQt+DWpnkdxnNq3dg+f9V1c/1VpkkaUlbKKB+aWx5ps9CJEkad8iAqqprj1UhkiSNO2RALfSq9ao65+iWI0nSyEK3+L6b0UsHtwK3MvkdT5IkHXULBdS3AGcBGxi9OPAvga1VteOQe0mSdIQO+Zh5VT1aVX9VVRcALwN2ATcneesxqU6StGQt+LqN7sWBr2M0iloF/C7wkX7LkiQtdQs9JHEt8ELgo8CvV9WnjklVkqQlb6ER1JuALwLPB96W5MDMEgGqqp7eZ3GSpKVroe9BOdu5JGkQC93iewrwZuB5wN3AlqradywKkyQtbQuNkK4FpoB7gB8E3tt7RZIksfBnUKdV1YsAklwD3NZ/SZIkLTyC+sqBBW/tSZKOpYVGUC9O8l/dcoCv69Z9ik+S1KuFnuJbdqwKkSRpnI+RS5KaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkpq04Bt1JQlg42VbBznvlZdsGOS8Gp4jKElSkwwoSVKTeg2oJOuS7EyyK8nFE7a/IskdSfYlOW/etguSfLr7uaDPOiVJ7ektoJIsA64CXgucBmxIctq8bg8APwF8aN6+3wi8C3gpsBZ4V5KT+qpVktSePkdQa4FdVbW7qh4BrgPOHe9QVZ+tqruB/fP2fQ3wsaraW1UPAR8D1vVYqySpMX0G1CnAg2Prs13bUds3yUVJZpLMzM3NPeFCJUnt6TOgMqGtjua+VbW5qqaqamrFihWHVZwkqW19BtQscOrY+kpgzzHYV5J0HOgzoLYDa5KsTnICsB6YXuS+NwJnJzmpezji7K5NkrRE9BZQVbUP2MgoWO4Drq+qHUk2JTkHIMlLkswC5wNXJ9nR7bsX+A1GIbcd2NS1SZKWiF6nOqqqbcC2eW2Xji1vZ3T7btK+W4AtfdYnSWqXM0lIkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKatHzoAvR4Gy/bOsh5r7xkwyDnlaSDcQQlSWqSASVJapIBJUlqkgElSWpSrwGVZF2SnUl2Jbl4wvYTk/xpt/3WJKu69lVJvpzkru7n9/usU5LUnt6e4kuyDLgKOAuYBbYnma6qe8e6XQg8VFXPS7IeuBx4Q7ft/qo6va/6JElt63MEtRbYVVW7q+oR4Drg3Hl9zgWu7ZZvAF6dJD3WJEl6kugzoE4BHhxbn+3aJvapqn3AF4Bv6ratTnJnkk8m+b5JJ0hyUZKZJDNzc3NHt3pJ0qD6DKhJI6FaZJ/PAc+qqjOAtwMfSvL0r+pYtbmqpqpqasWKFUdcsCSpHX0G1Cxw6tj6SmDPwfokWQ48A9hbVQ9X1X8AVNXtwP3A83usVZLUmD4DajuwJsnqJCcA64HpeX2mgQu65fOAT1RVJVnRPWRBkucAa4DdPdYqSWpMb0/xVdW+JBuBG4FlwJaq2pFkEzBTVdPANcCfJNkF7GUUYgCvADYl2Qc8Cry5qvb2VaskqT29ThZbVduAbfPaLh1b/l/g/An7fRj4cJ+1SZLa5kwSkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQm9fo9KEk62jZetnWQ8155yYZBzruUOYKSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNWn50AW0ZONlWwc575WXbBjkvJLUMkdQkqQmOYKSpMPk3ZZjwxGUJKlJBpQkqUm9BlSSdUl2JtmV5OIJ209M8qfd9luTrBrb9s6ufWeS1/RZpySpPb19BpVkGXAVcBYwC2xPMl1V9451uxB4qKqel2Q9cDnwhiSnAeuBbweeCfxNkudX1aN91StJTzbH+2dhfY6g1gK7qmp3VT0CXAecO6/PucC13fINwKuTpGu/rqoerqrPALu640mSlohUVT8HTs4D1lXVT3XrbwJeWlUbx/p8qusz263fD7wU+DXglqr6QNd+DfDRqrph3jkuAi7qVr8V2NnLb2ZxTgb+fcDzt8Rr8RivxYjX4TFeC3h2Va1YqFOfj5lnQtv8NDxYn8XsS1VtBjYffmlHX5KZqpoauo4WeC0e47UY8To8xmuxeH3e4psFTh1bXwnsOVifJMuBZwB7F7mvJOk41mdAbQfWJFmd5ARGDz1Mz+szDVzQLZ8HfKJG9xyngfXdU36rgTXAbT3WKklqTG+3+KpqX5KNwI3AMmBLVe1IsgmYqapp4BrgT5LsYjRyWt/tuyPJ9cC9wD7gZ58ET/A1cauxEV6Lx3gtRrwOj/FaLFJvD0lIknQknElCktQkA0qS1CQD6ihYaEqnpSDJqUluSnJfkh1J3jZ0TUNLsizJnUn+YuhahpTkG5LckOSfuj8f3z10TUNJ8gvd349PJdma5ClD19QyA+oIjU3p9FrgNGBDN1XTUrMP+MWqegHwMuBnl+h1GPc24L6hi2jA7wB/VVXfBryYJXpNkpwC/BwwVVUvZPTw2Pphq2qbAXXkFjOl03Gvqj5XVXd0y//N6B+hU4atajhJVgKvA94/dC1DSvJ04BWMntilqh6pqv8ctqpBLQe+rvve59fj9zsPyYA6cqcAD46tz7KE/2EG6GalPwO4ddhKBvU+4JeB/UMXMrDnAHPAH3a3O9+f5KlDFzWEqvpX4LeAB4DPAV+oqr8etqq2GVBHblHTMi0VSZ4GfBj4+ar6r6HrGUKSHwI+X1W3D11LA5YD3wn8XlWdAXwRWKqf057E6O7KakZvaXhqkjcOW1XbDKgj57RMnSRfyyicPlhVHxm6ngG9HDgnyWcZ3fJ9VZIPDFvSYGaB2ao6MJq+gVFgLUVnAp+pqrmq+grwEeB7Bq6paQbUkVvMlE7Hve41KdcA91XVbw9dz5Cq6p1VtbKqVjH68/CJqlqS/1Ouqn8DHkzyrV3TqxnNELMUPQC8LMnXd39fXs0SfWBksfqczXxJONiUTgOXNYSXA28C7klyV9f2K1W1bcCa1Ia3Ah/s/gO3G/jJgesZRFXdmuQG4A5GT73eidMeHZJTHUmSmuQtPklSkwwoSVKTDChJUpMMKElSkwwoSVKTDChpTJJK8t6x9Xck+bWjdOw/SnLe0TjWAuc5v5s1/KZ57d9/sJnVuymITuuWf6XvGqXFMKCkx3sY+OEkJw9dyLhu1vzFuhD4mar6gcXuUFU/VVUHvkBrQKkJBpT0ePsYfXnyF+ZvmD8CSvI/3a/fn+STSa5P8s9J3p3kx5LcluSeJM8dO8yZSf626/dD3f7LklyRZHuSu5P89Nhxb0ryIeCeCfVs6I7/qSSXd22XAt8L/H6SKyb8/p429m6mD3YzGpDk5iRTSd7NaLbtu7rtT03yl0n+sTvPG57YZZUOnzNJSF/tKuDuJO85jH1eDLwA2MtotoT3V9Xa7sWNbwV+vuu3Cngl8FzgpiTPA36c0czWL0lyIvD3SQ7Mcr0WeGFVfWb8ZEmeCVwOfBfwEPDXSV5fVZuSvAp4R1XNTKjzDODbGc0X+feMZgD5uwMbq+riJBur6vTuPD8C7Kmq13XrzziMayIdEUdQ0jzdLOx/zOjlcou1vXsn1sPA/cCBgLmHUSgdcH1V7a+qTzMKsm8DzgZ+vJsi6lbgm4A1Xf/b5odT5yXAzd3Eo/uADzJ679JCbquq2araD9w1r7ZJ7mE06rs8yfdV1RcWcQ7pqDCgpMnex+iznPF3F+2j+zvT3Ro7YWzbw2PL+8fW9/P4OxXz5xYrRq9seWtVnd79rB57T9AXD1LfpNe8LMZ4nY+ywF2UqvpnRqO0e4Df7G4hSseEASVNUFV7gesZhdQBn2X0jzWM3uvztU/g0Ocn+Zruc6nnADsZTTT8lu51JSR5/iJe6ncr8MokJ3cPUGwAPvkE6pnkK2O1PBP4UlV9gNHL9pbqqzI0AD+Dkg7uvcDGsfU/AP48yW3Axzn46OZQdjIKkm8G3lxV/5vk/Yxutd3RjczmgNcf6iBV9bkk7wRuYjSa2lZVf/4E6plkM6PP4O5gdKvziiT7ga8AbzlK55AW5GzmkqQmeYtPktQkA0qS1CQDSpLUJANKktQkA0qS1CQDSpLUJANKktSk/wPq9TqXkUIK1QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pmf = Pmf(ks)\n", "thinkplot.Hist(pmf)\n", "thinkplot.decorate(xlabel='Number of hits',\n", " ylabel='PMF')\n", "len(ks), np.mean(ks)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mean should be near 3.7. We can run this simulation more efficiently using NumPy. First we generate a sample of `xs`:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Hist({0.4: 706, 0.3: 294})" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xs = np.random.choice(a=[x1, x2], p=[0.3, 0.7], size=1000)\n", "Hist(xs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then for each `x` we generate a `k`:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "ks = np.random.binomial(n, xs);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the results look similar." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.629" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE3RJREFUeJzt3XuwXWV9xvHvYyJ4Gy1Kph0JmqDxErVCe4xa66WKGKsDTgtj0qrYwaFa47W2g8qgjcxUpLbWgVZSSaVeQERnzNQoOgpOqwPkcBEMNBqihdPY8dhQbb2AIb/+sVfs5riTcwJnZb9Jvp+ZPVnrXe9a+5c1OefJu/ba70pVIUlSa+437gIkSRrFgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1aeG4C5gvRx55ZC1ZsmTcZUiSZnHttdf+oKoWzdbvoAmoJUuWMDk5Oe4yJEmzSPLvc+nnJT5JUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKTDChJUpN6DagkK5NsSbI1yRkjtr81yc1Jbkzy5SSPHtp2d5IbuteGPuuUJLWnty/qJlkAnA+8EJgCNiXZUFU3D3W7Hpioqp8keR3wPuDl3bafVtWxfdUnSWpbnzNJrAC2VtU2gCSXACcBvwioqrpiqP9VwCt6rEcHsTVnXzzuEgA478zV4y5BOmj0eYnvKOD2ofWprm1PTgM+P7T+gCSTSa5K8rJROyQ5veszOT09fd8rliQ1o88RVEa01ciOySuACeC5Q82PqqrtSY4BvpLkpqq69R4Hq1oHrAOYmJgYeWxJ0oGpzxHUFHD00PpiYPvMTkmOB94JnFhVd+5ur6rt3Z/bgCuB43qsVZLUmD4DahOwLMnSJIcBq4B73I2X5DjgAgbh9P2h9iOSHN4tHwk8i6HPriRJB7/eLvFV1c4ka4DLgQXA+qranGQtMFlVG4BzgYcAn0oCcFtVnQg8EbggyS4GIfreGXf/SZIOcr0+D6qqNgIbZ7SdNbR8/B72+zrwlD5rkyS1zZkkJElNOmieqCu1qIXvZ/ndLB2oHEFJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKa1GtAJVmZZEuSrUnOGLH9rUluTnJjki8nefTQtlOTfLt7ndpnnZKk9vQWUEkWAOcDLwaWA6uTLJ/R7Xpgoqp+HbgMeF+378OBdwFPB1YA70pyRF+1SpLa0+cIagWwtaq2VdVdwCXAScMdquqKqvpJt3oVsLhbfhHwparaUVV3AF8CVvZYqySpMX0G1FHA7UPrU13bnpwGfH5f9k1yepLJJJPT09P3sVxJUkv6DKiMaKuRHZNXABPAufuyb1Wtq6qJqppYtGjRvS5UktSePgNqCjh6aH0xsH1mpyTHA+8ETqyqO/dlX0nSwavPgNoELEuyNMlhwCpgw3CHJMcBFzAIp+8PbbocOCHJEd3NESd0bZKkQ8TCvg5cVTuTrGEQLAuA9VW1OclaYLKqNjC4pPcQ4FNJAG6rqhOrakeS9zAIOYC1VbWjr1olSe3pLaAAqmojsHFG21lDy8fvZd/1wPr+qpMktcyZJCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTVo47gIk7V9rzr543CVw3pmrx12CDgCOoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElN8ou62mctfNET/LKndLDrdQSVZGWSLUm2JjljxPbnJLkuyc4kJ8/YdneSG7rXhj7rlCS1p7cRVJIFwPnAC4EpYFOSDVV181C324BXA28bcYifVtWxfdUnSWpbn5f4VgBbq2obQJJLgJOAXwRUVX2327arxzokSQegPi/xHQXcPrQ+1bXN1QOSTCa5KsnLRnVIcnrXZ3J6evq+1CpJakyfAZURbbUP+z+qqiaAPwA+kOQxv3SwqnVVNVFVE4sWLbq3dUqSGtRnQE0BRw+tLwa2z3Xnqtre/bkNuBI4bj6LkyS1rc+A2gQsS7I0yWHAKmBOd+MlOSLJ4d3ykcCzGPrsSpJ08OstoKpqJ7AGuBy4Bbi0qjYnWZvkRIAkT0syBZwCXJBkc7f7E4HJJN8ArgDeO+PuP0nSQa7XL+pW1UZg44y2s4aWNzG49Ddzv68DT+mzNklS25zqSJLUJANKktQkA0qS1KS9BlSSjwwtn9p7NZIkdWYbQT11aPlNfRYiSdKw2QJqX2Z+kCRp3sx2m/niJB9kMG3R7uVfqKo39laZJOmQNltA/dnQ8mSfhUiSNGyvAVVVF+2vQiRJGrbXgJrtSbZVdeL8liNJ0sBsl/ieyeCZThcDVzP6ERqSJM272QLq1xg8sn01g+cyfQ64uKo273UvSZLuo73eZl5Vd1fVF6rqVOAZwFbgyiRv2C/VSZIOWbPOZt49l+klDEZRS4APAp/ptyxJ0qFutpskLgKeDHwe+Iuq+uZ+qUqSdMibbQT1SuDHwOOANyXZPbNEgKqqh/ZZnCTp0DXb96Cc7VySNBazXeJ7APBa4LHAjcD67lHukiT1arYR0kXABHAT8LvA+3uvSJIkZv8ManlVPQUgyYXANf2XJEnS7COon+9e8NKeJGl/mm0E9dQkP+qWAzywW/cuPklSr2a7i2/B/ipEkqRh3kYuSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqUq8BlWRlki1JtiY5Y8T25yS5LsnOJCfP2HZqkm93r1P7rFOS1J7eAirJAuB84MXAcmB1kuUzut0GvBr4xIx9Hw68C3g6sAJ4V5Ij+qpVktSePkdQK4CtVbWtqu4CLgFOGu5QVd+tqhuBXTP2fRHwparaUVV3AF8CVvZYqySpMX0G1FHA7UPrU13bvO2b5PQkk0kmp6en73WhkqT29BlQGdFW87lvVa2rqomqmli0aNE+FSdJalufATUFHD20vhjYvh/2lSQdBPoMqE3AsiRLkxwGrAI2zHHfy4ETkhzR3RxxQtcmSTpE9BZQVbUTWMMgWG4BLq2qzUnWJjkRIMnTkkwBpwAXJNnc7bsDeA+DkNsErO3aJEmHiIV9HryqNgIbZ7SdNbS8icHlu1H7rgfW91mfJKldziQhSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWrSwnEXoL1bc/bF4y4BgPPOXD3uEiQdYhxBSZKa5AhK0ti1cKXAqwTtcQQlSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWpSrwGVZGWSLUm2JjljxPbDk3yy2351kiVd+5IkP01yQ/f6UJ91SpLa09vjNpIsAM4HXghMAZuSbKiqm4e6nQbcUVWPTbIKOAd4ebft1qo6tq/6JElt63MEtQLYWlXbquou4BLgpBl9TgIu6pYvA16QJD3WJEk6QPQZUEcBtw+tT3VtI/tU1U7gh8Ajum1Lk1yf5KtJnj3qDZKcnmQyyeT09PT8Vi9JGqs+A2rUSKjm2Od7wKOq6jjgrcAnkjz0lzpWrauqiaqaWLRo0X0uWJLUjj4Dago4emh9MbB9T32SLAQeBuyoqjur6r8Aqupa4FbgcT3WKklqTJ8BtQlYlmRpksOAVcCGGX02AKd2yycDX6mqSrKou8mCJMcAy4BtPdYqSWpMb3fxVdXOJGuAy4EFwPqq2pxkLTBZVRuAC4GPJtkK7GAQYgDPAdYm2QncDby2qnb0VaskqT29BRRAVW0ENs5oO2to+WfAKSP2+zTw6T5rkyS1zZkkJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTTKgJElNMqAkSU0yoCRJTep1JglJOlCtOfvicZfAeWeuHncJY+UISpLUJANKktQkA0qS1CQDSpLUJANKktQkA0qS1CQDSpLUJANKktQkA0qS1CQDSpLUJANKktQkA0qS1CQDSpLUJANKktQkA0qS1CQDSpLUJANKktQkA0qS1CQDSpLUJANKktQkA0qS1KSF4y6gJWvOvnjcJXDemavHXYKkRrXwOwr23+8pR1CSpCYZUJKkJhlQkqQm9RpQSVYm2ZJka5IzRmw/PMknu+1XJ1kytO3tXfuWJC/qs05JUnt6C6gkC4DzgRcDy4HVSZbP6HYacEdVPRb4G+Ccbt/lwCrgScBK4O+640mSDhF9jqBWAFuraltV3QVcApw0o89JwEXd8mXAC5Kka7+kqu6squ8AW7vjSZIOEamqfg6cnAysrKrXdOuvBJ5eVWuG+nyz6zPVrd8KPB14N3BVVX2sa78Q+HxVXTbjPU4HTu9WHw9s6eUvs2+OBH4w7iIOAJ6nufE8zZ3nam5aOE+PrqpFs3Xq83tQGdE2Mw331Gcu+1JV64B1+15af5JMVtXEuOtonedpbjxPc+e5mpsD6Tz1eYlvCjh6aH0xsH1PfZIsBB4G7JjjvpKkg1ifAbUJWJZkaZLDGNz0sGFGnw3Aqd3yycBXanDNcQOwqrvLbymwDLimx1olSY3p7RJfVe1Msga4HFgArK+qzUnWApNVtQG4EPhokq0MRk6run03J7kUuBnYCby+qu7uq9Z51tQlx4Z5nubG8zR3nqu5OWDOU283SUiSdF84k4QkqUkGlCSpSQbUPJltWicNJDk6yRVJbkmyOcmbxl1Ty5IsSHJ9kn8edy2tSvIrSS5L8m/dv6tnjrumFiV5S/cz980kFyd5wLhrmo0BNQ/mOK2TBnYCf1pVTwSeAbzec7VXbwJuGXcRjftb4AtV9QTgqXi+fkmSo4A3AhNV9WQGN66tGm9VszOg5sdcpnUSUFXfq6rruuX/YfDL5KjxVtWmJIuBlwAfHnctrUryUOA5DO4Ipqruqqr/Hm9VzVoIPLD7zumDOAC+W2pAzY+jgNuH1qfwl+6sutnrjwOuHm8lzfoA8OfArnEX0rBjgGngH7tLoR9O8uBxF9WaqvoP4K+A24DvAT+sqi+Ot6rZGVDzY05TM+n/JXkI8GngzVX1o3HX05okLwW+X1XXjruWxi0EfgP4+6o6Dvgx4GfAMyQ5gsFVnaXAI4EHJ3nFeKuanQE1P5yaaR8kuT+DcPp4VX1m3PU06lnAiUm+y+CS8fOTfGy8JTVpCpiqqt2j8MsYBJbu6XjgO1U1XVU/Bz4D/NaYa5qVATU/5jKtk4DucSoXArdU1V+Pu55WVdXbq2pxVS1h8O/pK1XV/P9497eq+k/g9iSP75pewGAGGt3TbcAzkjyo+xl8AQfAzSR9zmZ+yNjTtE5jLqtVzwJeCdyU5Iau7R1VtXGMNenA9gbg491/DrcBfzTmeppTVVcnuQy4jsGdtNdzAEx55FRHkqQmeYlPktQkA0qS1CQDSpLUJANKktQkA0qS1CQDShqSpJK8f2j9bUnePU/H/kiSk+fjWLO8zyndrN5XzGh/3p5mRe+mCFreLb+j7xqluTCgpHu6E/i9JEeOu5Bh3Yz5c3Ua8CdV9Ttz3aGqXlNVu7/gakCpCQaUdE87GXyB8S0zN8wcASX53+7P5yX5apJLk3wryXuT/GGSa5LclOQxQ4c5Psm/dP1e2u2/IMm5STYluTHJHw8d94oknwBuGlHP6u7430xyTtd2FvDbwIeSnDvi7/eQoWcnfbybVYAkVyaZSPJeBjNe39Btf3CSzyX5Rvc+L793p1Xad84kIf2y84Ebk7xvH/Z5KvBEYAeD2Qw+XFUrugcyvgF4c9dvCfBc4DHAFUkeC7yKwezST0tyOPC1JLtnml4BPLmqvjP8ZkkeCZwD/CZwB/DFJC+rqrVJng+8raomR9R5HPAkBnNFfo3BzB7/untjVZ2RZE1VHdu9z+8D26vqJd36w/bhnEj3iSMoaYZudvV/YvCAt7na1D3r6k7gVmB3wNzEIJR2u7SqdlXVtxkE2ROAE4BXdVM/XQ08AljW9b9mZjh1ngZc2U3+uRP4OIPnIs3mmqqaqqpdwA0zahvlJgajvnOSPLuqfjiH95DmhQEljfYBBp/lDD9baCfdz0x3aeywoW13Di3vGlrfxT2vVMycW6wYPK7lDVV1bPdaOvSsnh/vob5Rj3iZi+E672aWqyhV9S0Go7SbgL/sLiFK+4UBJY1QVTuASxmE1G7fZfDLGgbP1rn/vTj0KUnu130udQywhcEkw6/rHkNCksfN4aF7VwPPTXJkdwPFauCr96KeUX4+VMsjgZ9U1ccYPPDOR1lov/EzKGnP3g+sGVr/B+CzSa4BvsyeRzd7s4VBkPwq8Nqq+lmSDzO41HZdNzKbBl62t4NU1feSvB24gsFoamNVffZe1DPKOgafwV3H4FLnuUl2AT8HXjdP7yHNytnMJUlN8hKfJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSQaUJKlJ/wfUjicJZ8h1owAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pmf = Pmf(ks)\n", "thinkplot.Hist(pmf)\n", "thinkplot.decorate(xlabel='Number of hits',\n", " ylabel='PMF')\n", "np.mean(ks)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One more way to do the same thing is to make a meta-Pmf, which contains the two binomial `Pmf` objects:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Pmf({0: 0.028247524900000005, 1: 0.12106082100000018, 2: 0.2334744405, 3: 0.26682793200000016, 4: 0.20012094900000013, 5: 0.10291934520000007, 6: 0.03675690899999999, 7: 0.009001692000000002, 8: 0.0014467004999999982, 9: 0.00013778100000000015, 10: 5.904899999999995e-06}) 0.3\n", "Pmf({0: 0.0060466176, 1: 0.04031078400000004, 2: 0.12093235199999994, 3: 0.21499084800000012, 4: 0.2508226560000002, 5: 0.20065812480000034, 6: 0.11147673600000013, 7: 0.04246732800000004, 8: 0.010616832, 9: 0.0015728640000000028, 10: 0.00010485760000000014}) 0.7\n" ] } ], "source": [ "from thinkbayes2 import MakeBinomialPmf\n", "\n", "pmf1 = MakeBinomialPmf(n, x1)\n", "pmf2 = MakeBinomialPmf(n, x2)\n", "\n", "metapmf = Pmf({pmf1:0.3, pmf2:0.7})\n", "metapmf.Print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's how we can draw samples from the meta-Pmf:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "ks = [metapmf.Random().Random() for _ in range(1000)];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here are the results, one more time:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.708" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE3BJREFUeJzt3X2wXVV9xvHvYxBQGS1Kph0JmiBRiVqlvUatFd8QY3XAaWFMrIodHKo1vtZ2UBmwkZmC1NY60EoKqdQXENEZMzWKjoLT2gFyQQQDjYZo4TZ2vDZU6xsY8usfZ8ceLje5N+HunJXk+5m5k73XXuucX/YkebL22WftVBWSJLXmIaMuQJKk6RhQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQm9RpQSZYl2ZhkU5Izpzn+ziS3JbklyVeSPH7o2H1Jbu5+1vZZpySpPelrJYkk84BvAy8BJoD1wIqqum2ozwuB66vqZ0neBLygql7VHftJVR022/c74ogjauHChXP5W5Ak9eDGG2/8YVXNn6nfQT3WsBTYVFWbAZJcAZwM/Cqgquqaof7XAa/Z0zdbuHAh4+PjezpckrSXJPmP2fTr8xLfkcBdQ/sTXdvOnA58YWj/0CTjSa5L8srpBiQ5o+szPjk5+eArliQ1o88ZVKZpm/Z6YpLXAGPA84eaH1dVW5IcDXw1ya1Vdcf9XqxqNbAaYGxszFVvJWk/0ucMagI4amh/AbBlaqckJwDvBU6qqnt2tFfVlu7XzcC1wHE91ipJakyfAbUeWJxkUZKDgeXA/e7GS3IccDGDcPrBUPvhSQ7pto8AnsvQZ1eSpP1fb5f4qmpbkpXA1cA8YE1VbUiyChivqrXABcBhwKeTANxZVScBxwIXJ9nOIETPG777T5K0/+vtNvO9bWxsrLyLT5Lal+TGqhqbqZ8rSUiSmmRASZKaZEBJkprU5/egpAPeynMvH3UJXHjWilGXIO0RZ1CSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCb5wELtF1p4MCD4cEBpLjmDkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1qdeASrIsycYkm5KcOc3xdya5LcktSb6S5PFDx05L8p3u57Q+65Qktae3gEoyD7gIeBmwBFiRZMmUbt8AxqrqN4GrgA90Yx8NnAM8C1gKnJPk8L5qlSS1p88Z1FJgU1Vtrqp7gSuAk4c7VNU1VfWzbvc6YEG3/VLgy1W1taruBr4MLOuxVklSY/oMqCOBu4b2J7q2nTkd+MLujE1yRpLxJOOTk5MPslxJUkv6DKhM01bTdkxeA4wBF+zO2KpaXVVjVTU2f/78PS5UktSePgNqAjhqaH8BsGVqpyQnAO8FTqqqe3ZnrCRp/9VnQK0HFidZlORgYDmwdrhDkuOAixmE0w+GDl0NnJjk8O7miBO7NknSAeKgvl64qrYlWckgWOYBa6pqQ5JVwHhVrWVwSe8w4NNJAO6sqpOqamuS9zMIOYBVVbW1r1olSe3pLaAAqmodsG5K29lD2yfsYuwaYE1/1UmSWuZKEpKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCb1+jwo7Z9Wnnv5qEsA4MKzVoy6BEk9cgYlSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSDyyUDjAtPHDSh01qNpxBSZKaZEBJkppkQEmSmmRASZKa1GtAJVmWZGOSTUnOnOb48UluSrItySlTjt2X5ObuZ22fdUqS2tPbXXxJ5gEXAS8BJoD1SdZW1W1D3e4EXg+8a5qX+HlVPaOv+iRJbevzNvOlwKaq2gyQ5ArgZOBXAVVV3+uObe+xDknSPqjPS3xHAncN7U90bbN1aJLxJNcleeV0HZKc0fUZn5ycfDC1SpIa02dAZZq22o3xj6uqMeDVwIeSPOEBL1a1uqrGqmps/vz5e1qnJKlBfQbUBHDU0P4CYMtsB1fVlu7XzcC1wHFzWZwkqW19BtR6YHGSRUkOBpYDs7obL8nhSQ7pto8AnsvQZ1eSpP1fbwFVVduAlcDVwO3AlVW1IcmqJCcBJHlmkgngVODiJBu64ccC40m+CVwDnDfl7j9J0n6u18Viq2odsG5K29lD2+sZXPqbOu7fgKf1WZskqW2uJCFJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWrSLgMqyUeHtk/rvRpJkjozzaCePrT9tj4LkSRp2EwBtTvPb5Ikac7MtFjsgiQfZvDwwR3bv1JVb+2tMknSAW2mgPqzoe3xPguRJGnYLgOqqi7bW4VIkjRslwGVZJdPwK2qk+a2HEmSBma6xPcc4C7gcuB6Bp9FSZLUu5kC6jeAlwArgFcDnwcur6oNuxwlSdKDtMvbzKvqvqr6YlWdBjwb2ARcm+Qte6U6SdIBa6YZFEkOAV7OYBa1EPgw8Nl+y5IkHehmukniMuCpwBeAv6iqb+2VqiRJB7yZZlCvBX4KPBF4W5IdK0sEqKp6ZJ/FSZIOXDN9D8rVziVJIzHTJb5DgTcCxwC3AGuqatveKEySdGCbaYZ0GTAG3Ar8HvDB3iuSJImZP4NaUlVPA0hyKXBD/yVJkjTzDOqXOza8tCdJ2ptmmkE9PcmPu+0AD+v2vYtPktSrme7im7e3CpEkaZi3kUuSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKa1GtAJVmWZGOSTUnOnOb48UluSrItySlTjp2W5Dvdz2l91ilJak9vAZVkHnAR8DJgCbAiyZIp3e4EXg98csrYRwPnAM8ClgLnJDm8r1olSe3pcwa1FNhUVZur6l7gCuDk4Q5V9b2qugXYPmXsS4EvV9XWqrob+DKwrMdaJUmN6TOgjgTuGtqf6NrmbGySM5KMJxmfnJzc40IlSe3pM6AyTVvN5diqWl1VY1U1Nn/+/N0qTpLUtj4DagI4amh/AbBlL4yVJO0H+gyo9cDiJIuSHAwsB9bOcuzVwIlJDu9ujjixa5MkHSB6C6juCbwrGQTL7cCVVbUhyaokJwEkeWaSCeBU4OIkG7qxW4H3Mwi59cCqrk2SdICY6Ym6D0pVrQPWTWk7e2h7PYPLd9ONXQOs6bM+SVK7XElCktQkA0qS1CQDSpLUJANKktQkA0qS1CQDSpLUJANKktSkXr8HJUmzsfLcy0ddAheetWLUJWgKZ1CSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJh006gK0ayvPvXzUJQBw4VkrRl2CpAOMMyhJUpMMKElSk3oNqCTLkmxMsinJmdMcPyTJp7rj1ydZ2LUvTPLzJDd3Px/ps05JUnt6+wwqyTzgIuAlwASwPsnaqrptqNvpwN1VdUyS5cD5wKu6Y3dU1TP6qk+S1LY+Z1BLgU1Vtbmq7gWuAE6e0udk4LJu+yrgxUnSY02SpH1EnwF1JHDX0P5E1zZtn6raBvwIeEx3bFGSbyT5WpLnTfcGSc5IMp5kfHJycm6rlySNVJ8BNd1MqGbZ5/vA46rqOOCdwCeTPPIBHatWV9VYVY3Nnz//QRcsSWpHnwE1ARw1tL8A2LKzPkkOAh4FbK2qe6rqvwGq6kbgDuCJPdYqSWpMnwG1HlicZFGSg4HlwNopfdYCp3XbpwBfrapKMr+7yYIkRwOLgc091ipJakxvd/FV1bYkK4GrgXnAmqrakGQVMF5Va4FLgY8l2QRsZRBiAMcDq5JsA+4D3lhVW/uqVZLUnl6XOqqqdcC6KW1nD23/Ajh1mnGfAT7TZ22SpLa5koQkqUkuFitJ02hhoeYDfZFmZ1CSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJhlQkqQmGVCSpCYZUJKkJh006gJasvLcy0ddAheetWLUJUhSEwwoSdpHtPCfaNh7/5H2Ep8kqUkGlCSpSQaUJKlJBpQkqUm9BlSSZUk2JtmU5Mxpjh+S5FPd8euTLBw69u6ufWOSl/ZZpySpPb0FVJJ5wEXAy4AlwIokS6Z0Ox24u6qOAf4GOL8buwRYDjwFWAb8Xfd6kqQDRJ8zqKXApqraXFX3AlcAJ0/pczJwWbd9FfDiJOnar6iqe6rqu8Cm7vUkSQeIVFU/L5ycAiyrqjd0+68FnlVVK4f6fKvrM9Ht3wE8C3gfcF1VfbxrvxT4QlVdNeU9zgDO6HafBGzs5Teze44AfjjqIvYBnqfZ8TzNnudqdlo4T4+vqvkzderzi7qZpm1qGu6sz2zGUlWrgdW7X1p/koxX1dio62id52l2PE+z57manX3pPPV5iW8COGpofwGwZWd9khwEPArYOsuxkqT9WJ8BtR5YnGRRkoMZ3PSwdkqftcBp3fYpwFdrcM1xLbC8u8tvEbAYuKHHWiVJjentEl9VbUuyErgamAesqaoNSVYB41W1FrgU+FiSTQxmTsu7sRuSXAncBmwD3lxV9/VV6xxr6pJjwzxPs+N5mj3P1ezsM+ept5skJEl6MFxJQpLUJANKktQkA2qOzLSskwaSHJXkmiS3J9mQ5G2jrqllSeYl+UaSfx51La1K8mtJrkry792fq+eMuqYWJXlH93fuW0kuT3LoqGuaiQE1B2a5rJMGtgF/WlXHAs8G3uy52qW3AbePuojG/S3wxap6MvB0PF8PkORI4K3AWFU9lcGNa8tHW9XMDKi5MZtlnQRU1fer6qZu+38Z/GNy5GiralOSBcDLgUtGXUurkjwSOJ7BHcFU1b1V9T+jrapZBwEP675z+nD2ge+WGlBz40jgrqH9CfxHd0bd6vXHAdePtpJmfQj4c2D7qAtp2NHAJPCP3aXQS5I8YtRFtaaq/hP4K+BO4PvAj6rqS6OtamYG1NyY1dJM+n9JDgM+A7y9qn486npak+QVwA+q6sZR19K4g4DfAv6+qo4Dfgr4GfAUSQ5ncFVnEfBY4BFJXjPaqmZmQM0Nl2baDUkeyiCcPlFVnx11PY16LnBSku8xuGT8oiQfH21JTZoAJqpqxyz8KgaBpfs7AfhuVU1W1S+BzwK/M+KaZmRAzY3ZLOskoHucyqXA7VX116Oup1VV9e6qWlBVCxn8efpqVTX/P969rar+C7gryZO6phczWIFG93cn8OwkD+/+Dr6YfeBmkj5XMz9g7GxZpxGX1arnAq8Fbk1yc9f2nqpaN8KatG97C/CJ7j+Hm4E/GnE9zamq65NcBdzE4E7ab7APLHnkUkeSpCZ5iU+S1CQDSpLUJANKktQkA0qS1CQDSpLUJANKGpKkknxwaP9dSd43R6/90SSnzMVrzfA+p3arel8zpf0FO1sVvVsiaEm3/Z6+a5Rmw4CS7u8e4PeTHDHqQoZ1K+bP1unAn1TVC2c7oKreUFU7vuBqQKkJBpR0f9sYfIHxHVMPTJ0BJflJ9+sLknwtyZVJvp3kvCR/mOSGJLcmecLQy5yQ5F+6fq/oxs9LckGS9UluSfLHQ697TZJPArdOU8+K7vW/leT8ru1s4HeBjyS5YJrf32FDz076RLeqAEmuTTKW5DwGK17f3B1/RJLPJ/lm9z6v2rPTKu0+V5KQHugi4JYkH9iNMU8HjgW2MljN4JKqWto9kPEtwNu7fguB5wNPAK5JcgzwOgarSz8zySHA15PsWGl6KfDUqvru8JsleSxwPvDbwN3Al5K8sqpWJXkR8K6qGp+mzuOApzBYK/LrDFb2+NcdB6vqzCQrq+oZ3fv8AbClql7e7T9qN86J9KA4g5Km6FZX/ycGD3ibrfXds67uAe4AdgTMrQxCaYcrq2p7VX2HQZA9GTgReF239NP1wGOAxV3/G6aGU+eZwLXd4p/bgE8weC7STG6oqomq2g7cPKW26dzKYNZ3fpLnVdWPZvEe0pwwoKTpfYjBZznDzxbaRvd3prs0dvDQsXuGtrcP7W/n/lcqpq4tVgwe1/KWqnpG97No6Fk9P91JfdM94mU2huu8jxmuolTVtxnM0m4F/rK7hCjtFQaUNI2q2gpcySCkdvgeg3+sYfBsnYfuwUufmuQh3edSRwMbGSwy/KbuMSQkeeIsHrp3PfD8JEd0N1CsAL62B/VM55dDtTwW+FlVfZzBA+98lIX2Gj+Dknbug8DKof1/AD6X5AbgK+x8drMrGxkEya8Db6yqXyS5hMGltpu6mdkk8MpdvUhVfT/Ju4FrGMym1lXV5/agnumsZvAZ3E0MLnVekGQ78EvgTXP0HtKMXM1cktQkL/FJkppkQEmSmmRASZKaZEBJkppkQEmSmmRASZKaZEBJkpr0fwb3KGpLCf2CAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pmf = Pmf(ks)\n", "thinkplot.Hist(pmf)\n", "thinkplot.decorate(xlabel='Number of hits',\n", " ylabel='PMF')\n", "np.mean(ks)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This result, which we have estimated three ways, is a predictive distribution, based on our uncertainty about `x`.\n", "\n", "We can compute the mixture analtically using `thinkbayes2.MakeMixture`:\n", "\n", "\n", " def MakeMixture(metapmf, label='mix'):\n", " \"\"\"Make a mixture distribution.\n", "\n", " Args:\n", " metapmf: Pmf that maps from Pmfs to probs.\n", " label: string label for the new Pmf.\n", "\n", " Returns: Pmf object.\n", " \"\"\"\n", " mix = Pmf(label=label)\n", " for pmf, p1 in metapmf.Items():\n", " for k, p2 in pmf.Items():\n", " mix[k] += p1 * p2\n", " return mix\n", " \n", "The outer loop iterates through the Pmfs; the inner loop iterates through the items.\n", "\n", "So `p1` is the probability of choosing a particular Pmf; `p2` is the probability of choosing a value from the Pmf.\n", "\n", "In the example, each Pmf is associated with a value of `x` (probability of hitting a target). The inner loop enumerates the values of `k` (number of targets hit after 10 shots)." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.700000000000003" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADk1JREFUeJzt3X+s3fVdx/Hny1aY2+Is42q0P9aS1blO59CuTIlohEGJhu4PiG2ypTOYxmSd81cMUwJJxxKmRrcEVAh0km22w7LojakiAdQ/FGyBOVawoXQI16KwFZnRDSy8/eN8t5zd3O5+b++5Pb3383wkNz3fX+f7+abN8377Ped8T6oKSVIbvmPcA5AknT5GX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSHLxz2A6c4999xau3btuIchSYvKQw899OWqmphtvTMu+mvXruXgwYPjHoYkLSpJ/q3Pel7ekaSGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGnHGfyFW7dt6wZ8H3cdO12xZ8H9KZzDN9SWqI0Zekhhh9SWqI1/R1UqfjGjt4nV06nTzTl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1Jakiv6CfZnORwkiNJrplh+a8neSzJF5Lcm+RNQ8u2J3mi+9k+ysFLkuZm1ugnWQbcDFwObAC2JdkwbbVHgI1V9XZgH/C73bbnANcDFwCbgOuTrBjd8CVJc9HnTH8TcKSqjlbVy8BeYMvwClV1f1X9bzf5ALCqe3wZcE9VHa+qF4B7gM2jGbokaa76RH8l8MzQ9FQ372SuBv76FLeVJC2g5T3WyQzzasYVk/cCG4Gfnsu2SXYAOwDWrFnTY0iSpFPR50x/Clg9NL0KODZ9pSSXAL8DXFFVL81l26q6tao2VtXGiYmJvmOXJM1Rn+gfANYnWZfkLGArMDm8QpLzgVsYBP+5oUV3A5cmWdG9gHtpN0+SNAazXt6pqhNJdjKI9TJgd1UdSrILOFhVk8DvAa8H/jwJwNNVdUVVHU/yEQa/OAB2VdXxBTkSSdKs+lzTp6r2A/unzbtu6PEl32bb3cDuUx2gJGl0/ESuJDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ3p9Mbq0FO28Yc+C7+Oma7ct+D6kufBMX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSG9bq2cZDPwCWAZcFtV3Tht+UXAx4G3A1urat/QsleAR7vJp6vqilEMvFXeDljSfMwa/STLgJuBdwNTwIEkk1X12NBqTwPvB35zhqf4WlW9YwRjlSTNU58z/U3Akao6CpBkL7AF+Gb0q+qpbtmrCzBGSdKI9LmmvxJ4Zmh6qpvX12uSHEzyQJL3zGl0kqSR6nOmnxnm1Rz2saaqjiU5D7gvyaNV9eS37CDZAewAWLNmzRyeWpI0F33O9KeA1UPTq4BjfXdQVce6P48CfwecP8M6t1bVxqraODEx0fepJUlz1Cf6B4D1SdYlOQvYCkz2efIkK5Kc3T0+F7iQodcCJEmn16zRr6oTwE7gbuBx4M6qOpRkV5IrAJK8M8kUcBVwS5JD3eZvBQ4m+RfgfuDGae/6kSSdRr3ep19V+4H90+ZdN/T4AIPLPtO3+0fgR+Y5RknSiPiJXElqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqSK/oJ9mc5HCSI0mumWH5RUkeTnIiyZXTlm1P8kT3s31UA5ckzd2s0U+yDLgZuBzYAGxLsmHaak8D7wf+bNq25wDXAxcAm4Drk6yY/7AlSaeiz5n+JuBIVR2tqpeBvcCW4RWq6qmq+gLw6rRtLwPuqarjVfUCcA+weQTjliSdgj7RXwk8MzQ91c3rYz7bSpJGrE/0M8O86vn8vbZNsiPJwSQHn3/++Z5PLUmaqz7RnwJWD02vAo71fP5e21bVrVW1sao2TkxM9HxqSdJc9Yn+AWB9knVJzgK2ApM9n/9u4NIkK7oXcC/t5kmSxmDW6FfVCWAng1g/DtxZVYeS7EpyBUCSdyaZAq4CbklyqNv2OPARBr84DgC7unmSpDFY3melqtoP7J8277qhxwcYXLqZadvdwO55jFGSNCJ+IleSGmL0JakhRl+SGmL0JakhRl+SGtLr3TuS5m/nDXsWfB83Xbttwfehxc0zfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIb4xejz5JddS1pMPNOXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqSK/oJ9mc5HCSI0mumWH52Uk+2y1/MMnabv7aJF9L8vnu509GO3xJ0lzM+oncJMuAm4F3A1PAgSSTVfXY0GpXAy9U1ZuTbAU+BvxCt+zJqnrHiMctSToFfc70NwFHqupoVb0M7AW2TFtnC3BH93gfcHGSjG6YkqRR6BP9lcAzQ9NT3bwZ16mqE8CLwBu7ZeuSPJLk75P81DzHK0mahz43XJvpjL16rvMssKaqvpLkx4G/SPK2qvrqt2yc7AB2AKxZs6bHkCRJp6LPmf4UsHpoehVw7GTrJFkOvAE4XlUvVdVXAKrqIeBJ4Aen76Cqbq2qjVW1cWJiYu5HIUnqpU/0DwDrk6xLchawFZicts4ksL17fCVwX1VVkonuhWCSnAesB46OZuiSpLma9fJOVZ1IshO4G1gG7K6qQ0l2AQerahK4HfhUkiPAcQa/GAAuAnYlOQG8AvxyVR1fiAORJM2u15eoVNV+YP+0edcNPf46cNUM290F3DXPMUqSRsRP5EpSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDWk120YJC0+O2/Ys+D7uOnabQu+D42WZ/qS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNMfqS1JAl981ZfluQJJ2cZ/qS1BCjL0kNMfqS1JAld01f0nj4etri4Jm+JDXE6EtSQ3pFP8nmJIeTHElyzQzLz07y2W75g0nWDi37cDf/cJLLRjd0SdJczRr9JMuAm4HLgQ3AtiQbpq12NfBCVb0Z+EPgY922G4CtwNuAzcAfdc8nSRqDPmf6m4AjVXW0ql4G9gJbpq2zBbije7wPuDhJuvl7q+qlqvoScKR7PknSGPR5985K4Jmh6SnggpOtU1UnkrwIvLGb/8C0bVee8mglqeO7hU5Nqurbr5BcBVxWVb/UTb8P2FRVHxxa51C3zlQ3/SSDM/pdwD9V1ae7+bcD+6vqrmn72AHs6CbfAhwewbH1dS7w5dO4v9PN41vcPL7F63Qf25uqamK2lfqc6U8Bq4emVwHHTrLOVJLlwBuA4z23papuBW7tMZaRS3KwqjaOY9+ng8e3uHl8i9eZemx9rukfANYnWZfkLAYvzE5OW2cS2N49vhK4rwb/hZgEtnbv7lkHrAf+eTRDlyTN1axn+t01+p3A3cAyYHdVHUqyCzhYVZPA7cCnkhxhcIa/tdv2UJI7gceAE8AHquqVBToWSdIset2Goar2A/unzbtu6PHXgatOsu1HgY/OY4wLbSyXlU4jj29x8/gWrzPy2GZ9IVeStHR4GwZJakjT0Z/t9hKLWZLVSe5P8niSQ0k+NO4xjVqSZUkeSfJX4x7LqCX5niT7kvxr93f4E+Me0ygl+bXu3+UXk+xJ8ppxj2k+kuxO8lySLw7NOyfJPUme6P5cMc4xfkOz0e95e4nF7ATwG1X1VuBdwAeW2PEBfAh4fNyDWCCfAP6mqn4I+FGW0HEmWQn8CrCxqn6YwRtEto53VPP2pwxuNTPsGuDeqloP3NtNj12z0aff7SUWrap6tqoe7h7/N4NoLJlPQydZBfwccNu4xzJqSb4buIjBu+Koqper6r/GO6qRWw58V/e5ntcyw+d3FpOq+gcG71wcNnx7mjuA95zWQZ1Ey9Gf6fYSSyaKw7q7np4PPDjekYzUx4HfAl4d90AWwHnA88Anu8tXtyV53bgHNSpV9e/A7wNPA88CL1bV3453VAvi+6rqWRichAHfO+bxAG1HPzPMW3JvZUryeuAu4Fer6qvjHs8oJPl54LmqemjcY1kgy4EfA/64qs4H/ocz5NLAKHTXtrcA64AfAF6X5L3jHVU7Wo5+r1tELGZJvpNB8D9TVZ8b93hG6ELgiiRPMbgs97NJPj3eIY3UFDBVVd/4n9k+Br8ElopLgC9V1fNV9X/A54CfHPOYFsJ/Jvl+gO7P58Y8HqDt6Pe5vcSi1d3a+nbg8ar6g3GPZ5Sq6sNVtaqq1jL4e7uvqpbMmWJV/QfwTJK3dLMuZvCp9qXiaeBdSV7b/Tu9mCX0QvWQ4dvTbAf+coxj+aZmvxj9ZLeXGPOwRulC4H3Ao0k+38377e7T1TrzfRD4THdCchT4xTGPZ2Sq6sEk+4CHGbzL7BHO0E+v9pVkD/AzwLlJpoDrgRuBO5NczeAX3Yx3LTjd/ESuJDWk5cs7ktQcoy9JDTH6ktQQoy9JDTH6ktQQoy9JDTH6ktQQoy9JDfl/THSwO2mSHxwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from thinkbayes2 import MakeMixture\n", "\n", "mix = MakeMixture(metapmf)\n", "thinkplot.Hist(mix)\n", "mix.Mean()" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.23054197320000014" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mix[3]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: Assuming again that the distribution of `x` in the population of designs is well-modeled by a beta distribution with parameters α=2 and β=3, what the distribution if `k` if I choose a random Alien Blaster and fire 10 shots?" ] }, { "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": 1 }