Skip to content

Instantly share code, notes, and snippets.

@thomasgilgenast
Created September 20, 2016 01:55
Show Gist options
  • Save thomasgilgenast/935e1a36e910bd820467a6fae021745d to your computer and use it in GitHub Desktop.
Save thomasgilgenast/935e1a36e910bd820467a6fae021745d to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.stats as stats\n",
"from scipy.special import binom, gamma, factorial"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"### define some variables\n",
"# wikipedia parametrization\n",
"wiki_p = 0.5\n",
"wiki_r = 10.0\n",
"\n",
"# actual indisputable mean and variance\n",
"mu = wiki_p*wiki_r / (1-wiki_p)\n",
"sigma_2 = wiki_p*wiki_r / (1-wiki_p)**2\n",
"\n",
"# scipy parametrization\n",
"scipy_p = 1 - (sigma_2 - mu) / sigma_2\n",
"scipy_n = mu**2 / (sigma_2 - mu)\n",
"\n",
"# mixture parametrization\n",
"mixture_m = mu\n",
"mixture_phi = 1 / wiki_r\n",
"\n",
"### define the pmfs\n",
"# wikipedia pmf\n",
"def wiki_pmf(k, r, p):\n",
" return binom(k+r-1, k) * (1-p)**r * p**k\n",
"\n",
"# scipy pmf\n",
"def scipy_pmf(k, n, p):\n",
" return binom(k+n-1, n-1) * p**n * (1-p)**k\n",
"\n",
"# mixture pmf\n",
"def mixture_pmf(k, m, phi):\n",
" r = 1 / phi\n",
" return gamma(r + k) / (factorial(k) * gamma(r)) * (m / (r+m))**r * (r/(r+m))**k"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fcc0ccf5a10>]"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdUVEcbx/Hv0ARUFHvvHXvvYjdGY41i773XGBO7MbH3\nFnuMCZbYYuxGLLH3CnYFewMVpe68fyxBQPOKClzK8zmHk925s3d/EHn2MnfuXKW1RgghRMJgYXQA\nIYQQMUeKvhBCJCBS9IUQIgGRoi+EEAmIFH0hhEhApOgLIUQCEqmir5Sqo5RyV0pdUUp9857tlZRS\nJ5VSgUqpxhG2tQt5nYdSqm1UBRdCCPHx1Ifm6SulLIArQHXgHnAccNFau4fpkwVwAAYDm7XW60Pa\nHYETQHFAASeB4lprn6j/VoQQQnxIZI70SwNXtda3tdaBgCvQIGwHrfUdrfUFIOInSG1gp9baR2vt\nDewE6kRBbiGEEJ8gMkU/I+AZ5rlXSFtkRHzt3Y94rRBCiCgWmaKv3tMW2bUbPue1QgghophVJPp4\nAVnCPM+EeWw/MrwA5wiv3Ruxk1JKPgiEEOITaK3fd3D9nyJzpH8cyKWUyqqUsgFcgM3/p3/YADuA\nmkqpZCEndWuGtL1Day1fWjNq1CjDM8SWL/lZyM9Cfhb//+tTfLDoa62Dgd6YT8JeBFy11peVUmOU\nUvUAlFIllVKeQFNggVLqfMhrnwPjMM/gOQqM0eYTukIIIQwQmeEdtNbbgbwR2kaFeXwCyPwfr10O\nLP/khEIIIaKMXJEbyzg7OxsdIdaQn8Vb8rN4S34Wn+eDF2fFSAildGzIIYQQcYlSCh0NJ3KFEELE\nE1L0hRAiAZGiL4QQCYgUfSGESECk6AshRAIiRV8IIRKQSF2cJUREfkF+eDzxACCRVSLypcr3zvZb\n3rfeaRdCGEuKvogUf3+YNw8OHYKgIHiZ6DZ78hcFwPF1Vtod3xSu/6PE99hcsD813M0fDNmyQf/+\nkDVrTCcXQoQlF2eJD/L2hgYNYP/+MI0pPaCP+Sg+zxPwmBP+NdtSZqNuCxuY4xHaliwZrFhh3pcQ\n4vPJxVkiSgUGBzL575/pU/IgN/Z7fvgFH+DjAw0bQufhF1l0YmkUJBRCfCwp+uI/Xb4SyNQ1M0iR\n05lt1MaRZ5+9zyQ2d9n+sgxd/+pEy9WdeBP4JgqSCiEiS4q+eK8TJ2Bp6d+5sMKdPbmC2VrhMluo\nhx2v3+nrjw1nKRzu6ybZ37NXEwUblOFuKl8A1p5fxfq9t6L3GxFChCMncsU7du/SHKz3EzMChgOw\nYyVU6ARpfA/z9eU/qTCrOTpZIn68VRiA1NmzcmN9+PvqWAbcouTdvgxfD0uXwpYtUKlsYw443Q3t\ns3yzPz4T5jJ70lR6D0qE+qiRSSHEp5ATuSLUxUcX2bjVl386e7M1uHb4bakUpTskZUq1jfSoU/Wj\n9msyQfcJB1gSWAWThfn/c69jMGerefsRyjCn0d/MX2FP0qRR8q0IkSDIiVzxyQ55HqL0gkp8f7ku\n2xwzM5veodv8SMRMtR7XpgdpU63UR+/bwgIGdEtFpsTm2T4lvRRTw9w08yxFWLXBnhLlXtFsZTcZ\n5xciGsmRvsA/KIA043PzQt0xN/hkhiUHcX0xhNrsYECOzfxwoDIZMnze+7wKeEWPDYN5M60+k4/1\nJju3OEEJKnIQfxJBk1ZYYcuKxkto2VLGeoT4kE850pein8AFBUHPnrDor+PQvirYmE+y8qgA1st2\n0zT/C+btyUvy5FH3ngEBMKrPM5x+7sv3jOc22aD0HCi+GBYfhiA7+vSBKVPAxibq3leI+EaKvvgo\nfn7Qt8Ft/t4ZyHVyQc6d0LIeWAYC4PiqPHdG7ydJYstoef/Vq6FTJ/B1PAItvjIX/Oc5ASjIeUoU\nDGD8thJkyhQtby9EnCdj+iJSTNqE56OX9Kh4nlE7y7OTWqTlAVyvBRtWAKC0BRObd4q2gg/QvDns\nOvgUqxbNYPPi0ILvgA/raczCC+X5scAKpu1cTbApONpyCJGQSNFPgBYeXEOJseUYd64SGblHDm6y\nnTo44AMXWlAraA4bXDbQpWTHaM9SunBy1rZZTPOiX4W0aJbTntxcw9cuAM8v2zPosAsj94yP9ixC\nJAQyvJPAvHwJfQrsJKjkFzy3N7HRFaxN5m3D+YF0M4fTt2/M59Ia5syBP/vvYaepBkczQvOv4fa/\n5xK04p9OBymfuXzMhxMilpLhHfFB08a+4kevdizbbK70nb8Ck4I5Fn0ovGqYIQUfQCno0wdG7a9G\nxWotqdgxTMEHGh/JhPuuYsaEEyIekaKfgHh5wYZZnniTHGsTrFkLV1JCpVoVybttJi4tjf/nULac\nicAm9wkKOZWQ/A1scIVvdqRjdG9fLl0yNp8QcZ3xv+UixowcCWcD8lOI83RnPq8C01DftQfJmmSj\nYlU/o+MBYGlhycY2q3CwSEM+r6ScXgg33ftTkYN4vknF11+Dr6/RKYWIu2RMPwHQWjNv1zZ6f1Eb\nTG9n4yThJf36KcbPSGJguve79Ogyc3rZc3fdGTYTfgH+Dh1g7sI32FnbGZROiNhB5umL99p6dStf\n/vYlPCgCuyaZp2YCjo5w/br5v7HRq1dQqhS4u4dt1VB2Blm+WsHNYaewUPLHqki45ESueEewKZie\nG4aan6Q7C21qQ+vakPYsI0bE3oIPkCQJrF0Ldv8e0Ns9A5eGUOh3Hs3egIe7/PMV4mPJb0085x8U\niM3xqhBg/7Yx107oXozstf4yLlgkFSxonspJxqPQrTg8z0HypZtZdX8Qk+ru5fW7y/sLIf4PGd6J\n535ZbiJ/hzJYJb1NS+e8uBc7BBYmUlplxnOoR5wYF/fyuUvuKSXxWzef0u7pWU1zsnGb+6RjUosz\nTP8trdERhTCEjOmLcN68gSGZfmfOs5ahbZtS56Bjk9xM796KtkXaGJju49x7+gKX0i/ZeSMntviH\ntu+iBteX/EaT5pA6cWoDEwoR82RMX4Qzd5o/g54ND9cW+LgYG1psp03h1gal+jQZUjowd2NGplsO\nCdfukHE3357PQ4OVLWV9HiEiQYp+PLX+9B4eTZhJdm6FtgVixe5qP1K5svkIIa4pVAjSzhuFG1XQ\nwLRyULEjeCf35vDD3Yzd+6PREYWI9WR4Jx4yaRP5hjck4M1f/LTHRPMLoIDZqg81L80iXz6jE346\nraFLs4t4WRVnR76AtxveJKeG7zJ2zW5oXDghYpgM7wgAbly34OaUzSTfNJmx5e0p0wW2ZrXnTtsR\ncbrgg3mNniw9toQv+HdLwcJT7J7TkF9/NS6bEHFBpIq+UqqOUspdKXVFKfXNe7bbKKVclVJXlVKH\nlVJZQtqtlFLLlVLnlFIXlVLDovobEO/69lvzHbHO3hrI5UU+WB3piUsjO5K3djU6WpQY7jyIEikr\nm58c7g9LD4J3dgC6d494MZcQIqwPFn2llAUwB6gNOAEtlFIRjxc7Ac+01rmBGcCkkPavARutdWGg\nJNDt3w8EET2OHIF168I0aCsOn59Lv+Re9KwQt07e/hcrCys2t/udno4bYcd0CH57T0VfX2jWzDxz\nSQjxrsgc6ZcGrmqtb2utAwFXiLAYivn5ipDH64BqIY81kFgpZQnYA/7Ai89OLd7r2evnDB78bnv6\n9DBssC2OdrH48tuPlCFpBub0aUDLlu9uu33tBoVH1eOo19GYDyZELBeZop8R8Azz3Cuk7b19tNbB\ngI9SKgXmD4DXwH3gFjBFa+39mZnFe/gH+VNgZkn+yfIVpA6//vDYsZA4sUHBopFSsGAB5Mnztq1A\nxiUk7ZYH9ewo7hfs//vFQiRQVpHo874zwxGn2kTso0L6lAaCgHRASuCAUmq31vpWxB2OHj069LGz\nszPOzs6RiCb+NefofB4G3IC8NyD3X3CqM7iNoUCWdLRvb3S66JM0qXl9ntKlTZQp3ojr5TYz/m9o\ne/YJnfc+o5UHWEXmX7kQcYCbmxtubm6ftY8PTtlUSpUFRmut64Q8HwZorfXEMH22hfQ5GjKUc19r\nnUYpNQc4rLVeFdJvCbBNa70uwnvIlM3PVP2nfuz1m0W4yVsBiemVeypz2nczLFdMqfBjJzIf2MT8\nvU9xDLk1wAWc2Df9NL36WxsbTohoEl1TNo8DuZRSWZVSNoALsDlCnz+BdiGPvwb+Dnl8h5DxfaVU\nYqAsIHMrotiLF1Dkh4KcXAg1rofZYONLtTJpDMsVk/7qNwP765tCCz5AQS7iNWwODx8al0uI2CZS\nF2cppeoAMzF/SCzRWv+klBoDHNdab1FKJQJWAsWAp4CL1vpWSKFfBhQI2dVSrfW09+xfjvQ/w9ih\nr+gyOTfpeYAGduSCzjWTkyqfE6f7HIiTV99+iqNH4XLZ9rQPnVMAL0hK17ZrGDYxHUXTFTUwnRBR\nTxZcS4Du3oVl2cbwfdDo0LY32DK4oTsjV9qSNknCWoFyYKuHjPwtL8nx4ZpdMr6qXojLJf4hr0Nx\nLvY/iqWF5Yd3IkQcIVfkJjA+fj6MGKlJHXQvXPssiwEMmZ41wRV8gG9npGW83TjaFq9Inj4WXC55\nEJTG4+VJFhxfZHQ8IQwnRT8O67LmG5bpynTP1IGSHGcvzjwmFa96fUO2bEanM0bq1LB9xHZWfnUQ\nbf/87YYrdXl0pKZxwYSIJWR4Jw5r0jSY9dd/gaojwass7JlAwWBb9t/MHKtvgxjdfj6xmG5/dTE/\neZ4Nts8Ej/okS6bw8IC0Ce8PIBFPyZh+AvLgAWTObF5jB+vXUHYGlJtGzeQ92PntOKPjGcqkTZSY\nXZUzG5zh4DAIent3sPbtYdkyw6IJEaVkTD8BWbEipOADBNrDgeGkWXuJ/o2rGJorNrBQFpzss5eO\nOcaEFnyFic4soujyfhw6ZHBAIQwkR/pxkNbmpQeuXQvf/t13MH68MZlio8ePzT+nXN7HmUsvSnMc\ngDoVp9J+XmZcCn1tcEIhPo8c6ScQi7acZti1Dozje7JxM7S9UycDQ8VCqVPDD+M1P9OV0hzH0wGa\nfQ3nSwxj/w5bo+MJYQgp+nHMmQdn6HaqONN6Lidx2R84bp+DndSknvMrsmc3Ol3s0627YkbeaUys\nAMW6Q/7HcG1uIA4jbvDokdHphIh5MrwTx3TZ0IfF5+aEPrcJgirujlT6cgsj2pc3MFnsVXthKx7u\n288f273IGTKL0wcHRjbzYObqdMaGE+IzyPBOAnDD3Qb7gLf/jwOsYFfB55iy7jUwVey2tv18nPUJ\nUj13CG27Qh52rPHm8GEDgwlhACn6cYjWkGdhDR5M0SzaDGW8zO1KW9KlZAdjw8ViDokcGD4zLRPs\nxvGUFHRlIWU5ggf56DzIkx1XdxsdUYgYIyuNxyHHjkGG24dICnQ+Zf4am/Yr9OIvyZA0g9HxYrU0\naSDrTz3J068Vz0gJlgFQbjKXKo+j8e+23Bl8hZT2KY2OKUS0kzH9OKRzZ1iyBPLgQSeW0J7ljCy0\nkQXnZCw/MoKDoVQpOP1iJ3zRB1JdCd3WOn8XVjb72cB0Qnw8GdOPx16+BFdX8+Mr5OUbJpEJL8oN\nLGdssDjE0hKKDhgJbWqHK/g8LMTTvW2MCyZEDJKiHwcEBAdQZk5NfAvMBdu3txi2c7Ch6dcJY638\nqNK9Wj1Cby/m5wDbZpBi4U7SLrgmJ3VFgiBFPw6wsrAieP83kHU/9M8GjdpClgO0bKXj5Q3Po1Pp\njKVp49QF60ttUHMu0fWoLVdMTiyjI4vbHSA42OiEQkQvGdOPA86ehaL/3vTJ/jEUWQnFF/NV8bJs\n6rjU0GxxkUmbmDfXgqR92tGOX0Lbz1KYY/NP0aW73GhFxA2yymY89U2XZ2RZPIIldOI0xQEoVlyz\n56A3jnYJeA3lzxAUBN3y72fJtfAL1HXJMIUqm7PRukQTg5IJEXlyIjceevMGWPUrvZjHKUpwghI0\nYzVdOisp+J/Bygo6rajMKloC8NIGRlaFda2HsHT9GYPTCRF9pOjHYsGmYOb85kHrN29v81eCU+Sz\nvk7LlgYGiyfKl4eDX45lXnEr8vaBW8nh7EJNhsmF8PIyOp0Q0UOKfiy24/oOhnrlo1eHC6wsDG+s\nIAhLXjRqT7JkRqeLHwLazmVq4fT8+Rv02FAWF59/WBXYjO+/NzqZENFDxvRjsZpLGrHba2Po82R+\nUOFcNlr220qrWvkNTBZ/+Ab48sOgIK7O2ck6mgLm4VGl4ORJKFbM2HxC/D8yph+PmLQJrxuWWJre\ntvnYwtbStwhIfdS4YPFMYpvEDBmbjL9TfM2/BR/M6xwNHBLItafXjQsnRDSQoh9LBQVa8GzxOhyn\nnaXK7tpkeWZeJskWB5o5yR2fopKjI4wcGbZFQ76NuBVwwnlRXQKDA42KJkSUk+GdWOqPP6Bp0zAN\nKogMOdbyzTI/+laSFTWjWkAAODnBtdfHodZgyLY/dNvM2nPoW7aXgemEeL9PGd6RVTZjqcWLIzRo\nKyqWaEHfSobEifdsbKDggGFcezwxXHsSP0uO/mMBZQ0KJkQUk+GdWOj2bdix4932Ll1iPktC4lK5\neOhjq2DodwRuzQwmzRgbXrwwMJgQUUiGd2KZZ2+eUXpiK66vaw/uDSE4EQDZs8O1a2AhH9PRRmtN\nkVkVSHTIi992e5L7mbn9Pun4efBVRk1OYmxAISKQ2TvxgK1FYlLuLINjiYkwMDPUHgipL9GpkxT8\n6KaU4kjP3ZT1O0imZ7ah7el5gOX0KXh6GhhOiCgiZSSWcduTiImH3Xj2y2n+XpyI8kEnsWhbDZ8i\n442OliDYW9szaGYWZlsOCNc+IHgyEwZe47b3bYOSCRE1pOjHMn9Nv4Iz+wCo+tyLf/bsp8upXQyr\n3tPgZAlHtmzwqvcwHpEagEAFvYsV4+cslZi8Y5Wx4YT4TDJ7JxZ58ACy7VkSru0AFflyQCFS2BkU\nKoEaONqBHxaPIUW6RYyt9Qo/PwWrN3HxVGl0M/MVu0LERVL0YwnfAF/mLX5KL708XPsfyTsz5Qtj\nMiVkyZOD23cHOXH/JeyaZD6pjsLtLmzZAvXrG51QiE8jwzuxxKrzv/FDQA56NH/EtlwQrMCbZDh2\n/Ror+Wg2xLw2Q8i96yK4NyLsEg1DhkCgXKQr4iiZshlL5JtWGo+Xx0OfZ/aBgqdLMmL2JsoVzGBg\nsoRt0yZo2PDd9jlzNF27B2FtaR3zoYQIEW1TNpVSdZRS7kqpK0qpb96z3UYp5aqUuqqUOqyUyhJm\nW2Gl1CGl1AWl1FmllM3HBEwIXgW84tnDROHaPJPBNucT+NpfMiiVAPjqK6hcOXxb3vSrGHu6CIO2\nDTcmlBCf4YNH+kopC+AKUB24BxwHXLTW7mH69AAKaa17KqWaA4201i5KKUvgFNBKa31BKeUIeEc8\nrE/oR/rPn0P69OCf1B2KL4YiKyDxE9JYZ+f+t9ewUDIKZ6QTJ6BUKUjncJQc1VtyqMgNACy1DVf7\nuZPdMbvBCUVCFV1H+qWBq1rr21rrQMAVaBChTwNgRcjjdUC1kMe1gLNa6wsAWuvnCbq6/4cNG8Df\nH3iSD3ZOgWl3SbJtDVPr/igFPxYoUUJTuPc3ePcpG1rwAbQOZP2pfQYmE+LjRaaiZATCXovoFdL2\n3j5a62DARymVAsgDoJTarpQ6oZQa8vmR4x9X1wgNwTZ0LPM1rYs2NySPCE8pRbZSd/ALM3z/lTuc\nnGvD7bk1jAsmxCeITNF/358OEY/WI/ZRIX2sgApAC6AS0EgpVfVjQ8Znjx7B9d03sSA4XHuLFgYF\nEu81o/4ELLUNhe9ZsXc5bHKFok/9Kb7he06cMDqdEJEXmcmAXkCWMM8zYR7bD8sTyAzcCxnHd9Ba\nP1dKeQH7tNbPAZRSW4HiwN6IbzJ69OjQx87Ozjg7O0f+u4ijAoMDabd8NLPSLaHofcVamrOa5tzP\nUpYyZeTqn9gku2N29rY+wsZiR3B+9fbq6HpsoW1fb/76J7lcsCWinZubG25ubp+1j8icyLUEPDCf\nyL0PHANaaK0vh+nTEygYciLXBWgYciI3ObAbqAgEAduAaVrrbRHeI0EO9b/wf0HV9r3xTrMSSw0u\nF6DBhURsav2EsdNkRcfYaP7sICr3LUIerjCPnoxhFM9Jweq1QTRsZMLGUianiZgTLSdyQ8boewM7\ngYuAq9b6slJqjFKqXki3JUAqpdRVoD8wLOS13sA04ATmWTwnIhb8hOzFYwfauybj2ixY9Qf4WkOt\nNha4FxlqdDTxH7r0sGJM1mUU5AL9mclzUkCOXbQ+WIxfz6w2Op4QHyQXZxlo+uQgWg7NSFoehbb1\nzPAHY65UJnXiVAYmE//PX39BvXpAKnfzrRVTucPOKUzs0IChQ2WMR8QcWU8/jrm2ZF+4gu+DAxk6\n1JWCH8vVqB1A5q59oUMluFkV5l4E94b88IPi0aMPv14II0nRN8j163DBw4q9OGMKmfy0gUY0aWX7\ngVcKo9lYWtOyTk7U/EtweFDo3c1evjAxfvhrg9MJ8f/J8I4BXgW8IsePpXj8T3244EK6+2lpxjqe\n5yrNL1fLGR1PRFL37rBwoflxRQ4wnQEcowyVz8+lYEFjs4mEQYZ34og/Pf7kMe5QYTJ0K8GDPlWZ\nVfUpqTokNzqa+Ahjx0K2JE9YS1MOUJmSnOQLhwU0nF2fPz22GB1PiPeSom+Alcc3h29IeRWqjONq\nxrHGBBKfJE0a6P1tUopyhlc2MLIqOPUxcT3DFnpuHExgsKy/LGIfKfoGKHZrBazaAmdbg//b+fid\ny7gYmEp8ip4DrelSqQ55+sC4KvAmZKkGLz8PtnhsNzacEO8ht+eIYVrDWlcbuPql+cvqDeTeSonW\nG6mTq47R8cRHUlYBXKzzF4/DHNQXuw/9tmfhXvIvoIBx2YR4HzmRG8NOnYKjJXqg0LjiwgEqoZUl\nnp6QMeIydiJOWH1hDS5/NCfdS5iwB2qfTc0o/QMbHTty5boljo5GJxTx1aecyJUj/Ri24ZeXDGc5\ndvjRnYXcIz39Sx8hY8YsH36xiJWaOX3NGY9npG55kgcBqcjHt7zEAZ7DuHEwbZrRCYV4S470Y9Cp\ne2eYU34fS2/3D227RVa2z7tJ9x5yJWdc164d/PJL+DZra7h4EXLnNiaTiN9kymYst3zf36xrPoSK\nHWFOaXiQBNYoF5o0lYIfH0yYAPb24dsCk12m4vy6uN1yMySTEBFJ0Y9B9nvbc2+q5tsDcDQj5OsN\ns3rvRdvLtfvxQcaM8M2/d5C2fwJ1e0OHyjw6XJOA6+UNzSbEv6Tox5DgYPBcewTbYM2XV2HlBtg3\nJReN8n5DKntZaye+6DsggGR1pkGv/KAVFnPP0/2ILb91O0Zw8IdfL0R0k6IfQ/btg9+865KBe/Rg\nHvuozCbdmnGtGst9cOMR60SBFKpxFpbtp/a2Lzn3ujrz6Unf631ZsVSqvjCenMiNIV27wqJF4dsa\nfKXZuEnG8+MbkwlaFLnE6gtO4dr7OyxlnFcHkiY1KJiId+REbiy17ORKfr0xBZLdCdfu0kIKfnxk\nYQF9FxTgd8JfYd3bbxj1fvqGm89vGpRMCDnSjxE5JhXh5ptz5id3KsAFF+xufM3jW2lJnNjYbCL6\n9Kp3myl/5cNG+bGiKHxXDR4khbpZv+av9muMjifigU850peiH80uPb6E0zynd9qVtuBq3yvkTJHT\ngFQiJty6BROqt+Z4zVWcSR9+27HOxyiVsZQhuUT8IcM7sVAq68wU/2sgKa+VAJNlaHsG2zzkcMxh\nYDIR3exSPWRpm7XhC/6LDLBhBf63ShiWSyRsUvSj2cG/k7Lg+H6e/HqSy1McaLGlEvaeZelUqiVK\nyZh+fJY2SVp6lwq5+jrQDtxGwuwrcLYtgwZaYDIZm08kTDK8E8361L3G7G1vr8E3oRjY1JPpazJI\n0U8AXvi/4KtZ37Bv/HB4kTnctsWLoVMng4KJeEGGd2KZly8h5a7V4dr2U5k6nTJKwU8gHBI5sGfg\nfApny/zOtqFD4ckTA0KJBE2KfjTRWrN64wsaB4Uv+lsSu1C9ukGhhCEsLWHGjHfbbYP/ocxPX3D8\n7vGYDyUSLCn60cTjqQe9rmRkdPPzrHGC19YQhCU0aYK1tdHpREyrWhVatzY/trR4TeUKdfHrVpES\nx+7hc7WwseFEgiJFP5qktcyH3YyLPPLowMRiKckwCGo2SU2WRu5GRxMGmTIFsub/g+zdU2CfbRvH\nFsGaA+dwbbefgACj04mEQop+NNmwAXxeZeGfM0s59esTrGZfwPNpVzLkfmx0NGGQpCleY2o5ji5u\n6dm6CnI+N7d/c6cnMyf6GRtOJBhS9KOJq2v45099nahXZAxNnRobE0gYzt7anpvfnsZdrUbz9kR+\nbq4RMPYnbt0yLptIOKToR4NHj2DPnnfbXVzebRMJi6Wlou+vpVlI93DtXfVU6o0bzrar2w1KJhIK\nKfpRLDA4kCG/rsRk7ROuPWtWKFPGoFAiVilaFO50n8BD0gAwOUcRMvZMycUsP9JhXU9eB742OKGI\nz6ToR7HdN3bzy8u2MCQtNG8MTmvA2hcXF5Cp+eJfwyclZ0CmH8nWuDxD254lMOVtAB4G3GTe4SUG\npxPxmZXRAeKbpcdCBvOt/CH/BvNXQGL8i44ChhiaTcQeSZJoDveezW2/M28b/Rxgz4/ce9kNKhuX\nTcRvcqQfxUy3S1DwfoTPUhtfiuROY0wgESsppZjVcNzbhnMtYY4HHO/JrBmWnDtnXDYRv0nRj2Lp\nluXj/MIg3GfDmL2Q57EFltqGRvkaGh1NxDL189ajY/4BWP++C9avglfpAPP9lAd2eSkLsoloIQuu\nRaHLl+FYgXa045fQtmW0JcuJsVQvkdXAZCI2mzABvvvO/DgTnsykH1m5zekFx+jU1ULWaRL/SRZc\nM9ia5a9pzPpwbcdzt5GCL/6vwYMhX17NAKZxmfw0ZgO5E51iwZbq9Nw02Oh4Ip6Roh9FTCbN9t+e\nsZNa+GMDwD3Sk7d7VYOTidjOxgbmL1CU4SiJ8WVdASjQC/JbHMK0uoPR8UQ8E6nhHaVUHWAG5g+J\nJVrriRF3fajZAAAgAElEQVS22wC/ACWAJ0BzrfWdMNuzABeBUVrrae/Zf5wf3qkwtz6HdqWCc61J\ndqsIX+sNWCkTo+51I106o9OJuKB9y8PctajE/WTBLNgCFe+AK81Jt9cVZ2ej04nYKFqGd5RSFsAc\noDbgBLRQSuWL0K0T8ExrnRvzh8OkCNunAVs/Jlhck/PSQnhUCGoOxWdAERbX8uBMw9KkTRu3P8xE\nzEncbD2BXl9yeoG54AO4sJpf2u6WBdlElInM8E5p4KrW+rbWOhBwBRpE6NMAWBHyeB0QumK8UqoB\ncB3zkX68FBAAf7lmgMMD4eeT8MtuCLLlWem+aKToi8iZ23AyrVzWc85UHIBArJjAt7h6lmfKFIPD\niXgjMkU/I+AZ5rlXSNt7+2itgwFvpVQKpZQ9MBQYA8TbKQjbt8OzZ2EanuTH7vB4TvQ6gIWS0yYi\n8jp1tWRuwQXsxZminOE7JvAGO0avX0mL37oZHU/EA5GpSO8r1hEPXyP2USF9xgDTtdav/6NfvLBq\n1bttDRpA0qQxn0XEbRYW0O/XUtS03MslnCCVO7SrRmD9trhe/Zk/PbYYHVHEcZFZhsELyBLmeSbg\nXoQ+nkBm4J5SyhJw0Fo/V0qVAZoopSYBjkCwUuqN1npexDcZPXp06GNnZ2ec48iZq8M3zpN105+s\n4TSraMU2viCARKF3SRLiYxUpAj37vWH2mR+gwiSwDAzdNnnrGurnrWdgOmEkNzc33NzcPmsfH5y9\nE1LEPTCP098HjgEttNaXw/TpCRTUWvdUSrkADbXWLhH2Mwp4GZ9m7wSbgkk1ITPBvg9pftFEq/NQ\n8HZyOjps4o/HleW2iOKTjdnzE6MPfvu2Idga9kwg9aXenL5gS8aIA6wiQYqW2TshY/S9gZ2YT8a6\naq0vK6XGKKX+PeRYAqRSSl0F+gPDPi563LT31l68g+/z0tbE4hJQtT0U6+/Nw05/YGkl19CLTzek\ncl/SWGczP3lUABYdo/Gh7JzyzsXIppdkiQbxyWQZhs+w5J8NfLOmA09ThF87v5hjFU71dTMmlIg3\ntl3dTp9Z23i4YDizgobRgeUAnKYoByYfpe9gG2MDCsPJMgwxzOef+pyflYhDi6HnMUgZcrq6e4VW\nxgYT8cIXuetwZNRMWibZF1rwAYpxBv9hI/lt3zGCTcHGBRRxkhT9z3B4mTtJ8aWcF8zdCpenJqVl\nwDqaOX1tdDQRT6RKBU1cv8aV5qFt/pbwoPpEumz9Eo+Ht4wLJ+IkKfqf6NIlWOdekLQ8pAW/sYUv\nWR/cgpFtm5DcNrnR8UQ8Uqu24my3+XiSiQtpoHQXuOkIbnMSs2xilg/vQIgwZEz/E333nXlJ3LBK\nltAcPxEvL0UQBvN9baJa0z7cKDiPibuh4umctGEVxyjD7t1QvfqH9yHin08Z05fbJX6CVed+Z47H\naUjbGh4WDm1v3UYKvogetraanC0sqTGgE0FPLSjONHxJAkC7dnDuHKRIYXBIESfI8M4nsHhQkhfe\nVtCiPvQoBBUmopLfoXnzD79WiE9haWHJb21m4Tj0Z7rxc2jBB7h7F9r1vsdhzyMGJhRxhRT9T7Bv\nfW7YMwFm3oS/5oHjDax6FeOePmV0NBHPDRxsQdWIt2jIv54tmQvzxYpGPPZ9bEguEXfImP5HCgiA\ntim2YOnrw0Ya8prEACxZ4U/7NtaywJqIdp6eULgweL9+CXX6Q/GlodtqZ2nM9g5/GJhOxCSZpx8D\ntm2Dfr4/sIrWPCQtK2lNQdtrNGucSAq+iBGZM8Pg6cege7FwBT+TD7CqDsEydV/8H1KlPoJvgC87\nF1ylHOax0yT40ppVVKttTZIkH3ixEFGofBVfcLwR+rzFeTg/H6btmc7U8W8MTCZiOyn6H2Ho9hGs\nyleKb6vDxdTmtv1UonZXufG5iFlVs1elT8lBWPvbs+oP+O0PSO4HBbiM/ZhvOHHC6IQitpIx/UgK\nNgWT8odM+JgehLYVeQDJ3JuwZtPPpHWQ+XIiZvkH+bP574e8qj2aDiwLt61zxm3M9KhD4sQGhRMx\nQsb0o9G1Z9fw83sdru1sOthfZT0B2tegVCIhS2SViK9rZeH2wJlcJ0e4bfVfDSPvDzU4++CsQelE\nbCVFP5KS+OfFNPEupX4fTrkLmbANua9F8ZRVyJwss7HhRII2/MekjMvzK8Ehv869C5WlUde73N1f\nk1vHCxqcTsQ2MrwTSVOmwJAhb58nSeRFhgquTJ/vRN08XxgXTAjA3R0WlRzDH3V2cDvdc1i/Cu4X\nJ3VqOH8e0qY1OqGIDp8yvCNFP5KKFoWzEf5SHjUKwtzlUQjDaK3JOqEEnofKw65JEGgfuq1WLfjr\nL7CSRVfiHRnTjyYXLrxb8AFaybL5IpZQSnGy/07qMidcwQfYuRMGDA7E/Ym7QelEbCJF/wOevH5C\nX9fJ4OAVrr10acid26BQQrxH6sSpWLoUUqcO3+5ge5Ntz/JQcl5Fbjy/8f4XiwRDiv4HvA7w4/K1\nE6juhaFdNSi2BGy95ShfxEpp08KKFaBC/uDPkmIXKTvn43rOW/jqp1RbXJ9XAa+MDSkMJUX/A26d\ny8Ti1a94Pu0VQ4/5kiHPz9A/K34FFxgdTYj3+uIL88SDfNnm8aJLbW6mCgjdlnFnVu7elMn7CZmc\nyP2AwW0f8dPKDFjxdkGTTpX3M25LTjIkzWBgMiH+270X98kyNQfBFn4A2AbC8o3Q9KIFXTJuY+r5\nWjg6GhxSfDY5kRvF/P3BYt2acAX/PAWp1qWiFHwRq2VwSM/k2pMASPdSsX8ZNL8IlpiYfvdr+ta7\nQWCgwSGFIWQS1/+xbRs0fvNruLa11q0Y2lDukCViv/7leuPj68eVHqkpca8jYP5rejnt+f1QFpL0\nhbFTnpA6cSpjg4oYJUf6/+HEvRN0+ceJQ+WOci/p23afL1vKipoiTlBKMbrmEEb82Z4RtlMIwpLu\nzKc/MwnGkgXuoyk2vSYmbTI6qohBMqb/H5r91p61V1cAYGGCajeh5Lk8lBx/nCb1HAxOJ8TH2blD\nM7juJc6bnMDqDTRsD8nuoNZsYOuadNSpY3RC8SlkTD+KBJmC2O9xPvS5yQJ254SfGl3hWrKFBiYT\n4tPUqq3oOtMJktyHDlXAZAUr9qJfpqN5c7h0yeiEIqbIkf57vHoFmTJrfFLuhPJTIecuACx1Iu4N\nuUOaxGkMTijEx/ML8iPtWCde7OsA+78DzAeIqXhMjfSXmH2uCqlkeD9OkSP9KLJ8Ofh4K7heG1bu\nhHnnsDjbgbYFO0vBF3GWrZUtJ/rsoYbN9/xb8J24wDFK8/P9evRouINfTv9mbEgR7eRIP4LgYMib\nF65fD9/esSMsXqxRSmbuiLjt+XMoWxZyXtmKKy448JIjmaC+iwVPE8OWlltk5dg4Qo70o8Cff4LD\n9VP8O73tXwMHIgVfxAuOjrBlCzjZ3cSBl/xWCJzbw5MkJrQy0ez3NrJUQzwmRT+MrVe38t2GJrim\nLMExStMcVywJonZtcHIyOp0QUSd3bqizpQd1q5akVRPwD7liJ+Vr6LesOHu3y1IN8ZUU/TDUgxLk\n83SnYkcY1+IEPbK2YCntGTTI6GRCRL38pR5wwPnm2+eP4egiSHy7JC1cNG5uxmUT0UfG9MPoU/8W\nM7bkJMDKxC9FYHo5eEoujozZQc4UOT68AyHimAO3D+C8rDrlr1vzx9pABvkv4VfaAGBnB5s2Qc2a\nBocU/0nG9D/DnTuQ46/ZWGLCLgi6nYQ1cwvQOv8MMjlkNDqeENGiUtZK7OuwjzxeJ2jkvze04AO8\neQP12lyn8LSqvAl8Y2BKEZWk6If4ebIPnfSicG3Lkg7ip05fksgqkUGphIh+FbOWY/qm/LwuWiH8\nhmS3CWhRnYurXdi+xc6YcCLKSdEHDlw7zeKVmlGM4RZZAXhAWtIMaEUiqfciAXBwgD17oGTJfxu8\nzDcNOjyQsscKcqdJf1b/Lmv0xAeRGtNXStUBZmD+kFiitZ4YYbsN8AtQAngCNNda31FK1QB+AqyB\nAGCo1nrve/Zv2Jj+S/+XpPkpE34vksLRPlic7EhjPzeSWb/hp3tt5QpFkaD4+ED1Bg84WagKnOpM\nqUNV2E0NHHjJCtphsWwJbdpbGh1ThIiWMX2llAUwB6gNOAEtlFL5InTrBDzTWufG/OEwKaT9MVBP\na10EaA+s/JhwMeHnk0vw4wU43IWawzANzM66Lw7wpmtFKfgiwUmWDFzGriPbizYUPVSDHdTGgZcA\ntGMFlh1a8dWUkZx5cMbgpOJTRWZ4pzRwVWt9W2sdCLgCDSL0aQCsCHm8DqgOoLU+q7V+EPL4IpBI\nKWUdJcmjiPt5O3iV9m2DjS+Umc3DogONCyWEgQZX7s3F+d9TvLwdb3g7lq+Bi9VW86fvOCr8XI0T\n904YF1J8ssgU/YyAZ5jnXiFt7+2jtQ4GvJVSKcJ2UEo1BU6HfHDEClrDuaXdYPpt2LgUHhYM3Tay\nphR9kXDZ28PcPfkYWWU/d8gMwPjKMKGyeftr/Zwey+YYmFB8qsjcOet940URB+Aj9lFh+yilnIAf\ngf+c8Tt69OjQx87Ozjg7O0ci2uc5fBisjv2DBWUxnekAZ9pDjt00+vZPKmWpFO3vL0RsZmsL83bm\none9/Ti+Ls6kas/fbvSox4k1PzNBwfDhxmVMaNzc3HD7zKvmPngiVylVFhitta4T8nwYoMOezFVK\nbQvpc1QpZQnc11qnCdmWCdgDtNNaH/mP9zDkRG6v2teYvTMPt8jGTPqxlI7kKZ6UEydAltkRwuz0\n3QsUX1QYVMjv6LVa4LoJgmwBGDECxoyR3xkjRNfFWceBXEqprCGzdFyAzRH6/Am0C3n8NfB3SKDk\nwBZg2H8VfKNcuRpM/p0zsUCTg5vMpD9/Uj9kYTWj0wkRexTLWJDF9ZeCVnCzKqzeEFrwAcaNg2HD\n4Nnr5/9nLyK2+GDRDxmj7w3sBC4Crlrry0qpMUqpeiHdlgCplFJXgf7AsJD2XkBOYIRS6rRS6pRS\nyvA5MQ9ePaDEysw8rLyQx/Zv29c7dqZZM+NyCRFbdSrRni0tt9LOZjME2r+z/c+1s8n6kxPeb14Y\nkE58jAS59s7z5/Ct0wCCS85gXQFodhFaHknN8QFeDB5uE2M5hIhrTCbo2xfmzn3bVijLHO4378Ok\nNTk4nGsXE9fmwNHRuIwJiay9E0lL5gcw4v4aFv0J7nMg/Sv4sv0bLCsuMzqaELGahQXMnm2+vwRA\ngUyLud+8L7//AR1u32DSnuIMzbeZU6eMzSn+W4I70g8IgNzZgyh/by2DmEpJTuKLPSO7XuX7GXY4\n2skhihAfojV0GnGS9QHl+W1jAHWvht8+xXIo96ZWp/aXUDtXLWNCJgBypB8Ja9fCnXtWuNKCUhyn\nMvsYoGbS65sMUvCFiCSloG3nlzSwW8Gbq43f2V4l8Qpm3GtBnV/rMHznaIJNwQakFO+ToI709950\no2u3YK7tqkbYSwuaNIF166L97YWIlzZu0BxtMYOx/kOxJghfS0vSt8/Py8wXALB6k54DLc9RtrDh\nczjiHTnS/z+01nRbP4hrFWpAt+JQeCVYBgBvxyeFEB+vYSNF54sD6JLLDS8yUqlGhdCCj8mSoN9X\nU7NCKlavNjanMEswR/r7bu3DeYVz+MYXGch8vw8eK/pjZ2373tcJISLnzRtoOmg5W9N2eNu4czIc\nGgyAJUEM6vKScXMcsZFJclFCjvT/j6CnWclwrAEEhrkZhMM9/AouxNoyMqtRCCH+Hzs7GDo0Ow4W\nIQsYXm4Ih97eYHo839N9UXFald9KG9cexIYDzoQowRT97bMTc23rDq5Os6HmnipYvEwDwDDn/lhZ\nSNEXIipUyVYF9/6nqZ+lLTnOL+Pfc2f12cwwJpLJ4hZ3i9Tj5qK77NhubNaEKkEU/fv3wWHVfOzw\nI9cbH3Ye2Me5GUlpabucLiU6Gh1PiHglfdL0bO6wgtOHk9O0KWThNitCVmn5tgYk99fs3/snT+q2\nYfywVwTLxJ4YFe/H9E0mTdcqV5h1sBj2vL258/BEUxn2aCAODtHytkIIzPP5504PQA8eQsZ8sxhY\nG04uhJQhv4rz6MGykvPoO+EMeQsEUDpjaWMDxzEyph+Bj58POSaUptzthuEKvjfJUF06S8EXIpop\nBb0H2pB+Vz9a1Xdg2Vq70ILvSSZGMYYTJ020XduFMovL0HJ1Jx77PjY2dDwXr4v+vZvJsFjdhm+b\nezCyKgSGfLejkk6n/0ip+ELElHpVMvB7i42synaKCzgRiBXNWMMTUkOxpZDRfBeu3y+sYsb8lwTG\nmlstxT/xdngnIADKloXTpyFVkrNkbVADC/sntF1fg+zLd/JlPVk/WYiYFhwMk0b5su+Hg+ygNtg+\nhz55IPETc4d9I6i9tzyZMlvgsrQWNWoYmze2+5ThnXhb9IcPhx9/DNtionypNpytvRvP4e6y5IIQ\nBjp4EHr1gnMZ+0KZ2eZG7yzYzT2BR2BxMuPFBhqys/Y0vhj/lKK505AlWRZjQ8dCUvSBIFMQhw5a\n4uysiLjLvHnh70PeZEiRPEreSwjx6Z77viTr1Fy81I/MDWvWMu7SGb7nh9A+PpaJyNkzMXUdFrHw\nu8bY2f3HzhIoOZELDNo6nDorGqLtwp8MsrKCX39FCr4QsYRj4qTcHHyJDgV7kMmvFikuOTOQaeH6\nLCzrT9mnz2j5w88UKAAbNvDOwZz4OPHqSH/vzb1UW1HdfC/PV2lh01K4WheA8ePhu+8++y2EENEg\nyBTE2dNWzO1wgq7ne1OWo9xLCoV7wJHF0PPZDnZhXqK5Zk2oM/APCudJRvXs1VEJ+P6mCX54p/mc\nsax5Oip84/EelH0xhYN77bG0/Oy3EEJEI5MJfv3FxMl+v3C3Rnfy+PhTak9DGrPhbSeLQOifHRzu\nks2+EDs7rid3ylzGhTZQgh7e8fKCbEMd2fYrpHsZZkOp+eTtNUwKvhBxgIUFtG1vQbdL1dnjlJ6k\nB/q8M+SD01pwuAuAzyMPJpa8xJLFGj8/AwLHQfHiSN9kgi7lLzLnaEns8OOJPXSpDxvzQzKLdFwb\ndJ5U9rKWtxBxiX+QPzeuJqJfP9i1699WDV1LQYaTAIz9G0bsh4sUYEmS/qTo25qvu73BZPeQ/Knz\nG5Y9piTY4Z1Zk/2pNLQsxTgT2uZNUhq2GcGwEYWpk7t2VMQUQhhAa9i40Xzfi1umA9CxMgCJgsBz\nGqR+/bbvY1JRwbkeqqQl69svxsnJoNAxJMEN75i0ianb1jLxex9eEP4K2++Sz+eP6UOk4AsRxykF\njRrBpUtQsvvC0Pa2Z8MXfABbS2+ul9jKlRUDKFgQateGHTtkxk9Ycbro3/d+zohtU7nXrC1Vk/7K\nUCYSgDW/40ID15akTGl0QiFEVLGzA7dBC5nzxTwyJMqDt/fv/I4LQbw9YTegUBlMD4vBY/Mh/s6d\nUKcOFCxoosCEWrRb252/b/5t1LcQK8Tp4Z3Bg2Hq9CCoPB5Kzoe/5lP0cnbqdMvGjwvkilsh4iut\nNUopTp6EFeM9ybJpNp30QjL0SIvfztlwPfxf+EWyTeNse/MNXdI/LcaIfCdp1FiRLp0R6aNOghrT\n//tvqF49TEOmI9C4Ncl9KnFt1lxSOthHbUghRKzl5QXfzPoLV7+hmGZf4N+bt/yrVIt0HM/7EIBZ\nW+GrY1lYTxOOl+pBuTa5yVX5JMlTBFAuczkD0n+6BFP0nz+HwoXN/6PDsrJ/Ra/lc5naZDCWFjJH\nU4iE5MGrB3g88OTCjlLMmAHXrpnbs6bcwe0+dUL7eU2FjCHTuquxh71Ug0ZtyG5Tjj5letKkCWSJ\nI8v8JIii73phNWt+sKO56yr6M4MHpA/dNmkSDBkSXSmFEHFFcDBs2QLTp8ODJNXwKLUXgHKecGiJ\nuc9jUpGe+wRbBsPgdDD3Erwy15O2BU5QqE1RXhSaR+Z0dlTKWom8KfPGuqt/P6Xox6mbwy49so5O\nO1ywzAOFneH0/t30MC1mI41wdjZP6RJCCEtLaNAAvvpKU+1nB67ct0ArE7kvFeI117DnDRtoRDBW\nkGMHPCoUWvDT8oBll0rjNSIFeYb44J8oCJvgJExP95wqVazIn998EVlcFSeKvtawZrXmdM8HZG2r\nuO2oGeMM23M9Y+X6xii/7cxYUVuuuhVChKOUYm+3jTzyfcRG901krvEFEzc48mzVNg48CFm6ocA6\nuNQ09DWN2IAFmtM5n+KfyNxW9eYrUo1rSUHWkDIlVKgA+Stc47zjODa2X4G1tQHf3CeK9UX/3j3o\n2RPyb/qJ2Qyn4WaoYb7HMkczQeHuFnTK9CTOjMEJIWJemsRp6FqiCwBfVAA9uSnnz8PqdQH8FLQZ\n095xoX2bsg6AdQXevr7SbbiM+Qrfp09h82bY7LkPsgdRuq87FZ2ek+aLEpR3tiFdvptc9D6OU2on\n8qTMg7Vl7PpEiLVFf8+Nv/l1yy02jOyIjw/8Q2cGMo3qN5/Q6RQsKW7uF2BtQaPacXzelRAiRill\nngySt4Cm3I0V5GidiXXrYN1azc0L2XlgmYzNeX1C+1e6A2OpFH4nWQ/AnUq095tPv5OzeHPSlqPj\nyzCoghU7au4BoJaaRLeCQyhYEHLmNA87+fj5kMQmiWGTTWLlidybN6FlXw+OFKgAP58A72wANMcV\nV1rgZWtP/l7WWCSzZbXLCrnqVggRZTw84EfXjaygEQDWQQqvnxKTI+g+viR527FvTvh9Mycft6E4\np0Obv2oBf+YNefLrVrj2BQCJEkH+/PDKuSlFU9albZGOFCwIWbOazxHce3mPFHYpsLWyjXTWOD97\nJzgY5s0IYMb3T7jhlwEq/gTZ98DKnZjn3Wq+ZzyPa7eh4eQnlMmVU257KISIcq8CXuF2y42/rx3g\nqucLSt2dzb5/rDhyBF6/BpLegx6FSDL5Ot46JZaYADApSDUUnv97h6/pt8En/Nhznu52LN6ssb1X\niAsU5KqNE6+yFWRFk368SHQFB4u0TC90kHJ5c5ElCyRObH7d1adXyZY8W7jhojhd9P/45zRbugfR\n73wn/ElEeQ4RbAF0LgPHe8HpjqRNC3PnQpMmRicWQiREgYFw6hTM2buG3Q9XkWr5XIZ7D6ESB8jE\nXS6kgUI9Qzr7J4UffQh7oZiFhR+JvrXj8SRIHPh2vxqw+M4OrN+YG370Bv9kAPRPvAidLj2Lmren\nvfdBcuXKR5Ys5r8QSpWKw0U/de8U3FzwnMRB5jyDmcxUBkPas9C2Jk2fnGHhlAykSGFwWCFEghcY\nHMhzv+ekskvD5ctw8IDGffstznqPZm/VX8ydPMvCksPhXpct5TYsWtXl+qzw+7tsn4ICQ5+Zn/g5\nwE/m8wn2+OJLEvysINkwePMDvNQO3CAHJTiJxjJ65ukrpeoAMzAv0LZEaz0xwnYb4BegBPAEaK61\nvhOy7VugIxAE9NNa73zfezxO9YwfqsAE8/kPxjGCTTTA36YI9QuMptsXT0iRIsPHfG9CCBEtrC2t\nSZM4DQBOTuDkpKB7dvbebE/200k45XmRlHYlydsTLlwwfz17BvZpjpHlkTUQGG5/J5KnBkKKvk/W\n0PYs3AHA0wEyvQALDcl4QUqeoj9xvcwPFn2llAUwB6gO3AOOK6U2aa3dw3TrBDzTWudWSjUHJgEu\nSqkCQDMgP5AJ2K2Uyv1fl99OqgBNL0Hx+7COpjTu6Mj3MyBp0p7v6x4vubm54ezsbHSMWEF+Fm/J\nz+Kt2PyzqJq9KlWzV32nXWt4+BAm78mEh+coBlt0JPD0BexuXCSX/wX+TvYc8DB39n5b9LNyG4Db\nySHr28lE3OHT56hH5ki/NHBVa30bQCnlCjQAwhb9BsC/N6ddB8wOefwV4Kq1DgJuKaWuhuzv6Pve\nKNiC/7V3v6FVV3Ecx98fmrDKWla2RUuLLOhJCCuScmAPCulBRiCtSS2LCLI/D5UgpEetIMMe9KRW\nmGikD0oTUomCioxGtbI2GyKuLd1cZtYebLn27cHvd7e7u3u3e0f8ztz5vuDA7x5+d/fs7Mv399u5\n55wfj9y3gBX7dtPy3hoebix21vw2lwM6a94XE7wvJpyPfSFBXR28uu7xvNqrMbub3l44+sVOLj32\nFX+N9XNT7RJqG6GnB/7oW0zb2GN01hzisrN9DPMP1YzQw9KSnzWTcpL+NUBv3us+ksRd9Bwz+1fS\nWUmXp/X5g1q/pXVTVJ+D4QVwpE68/tmVNN5Y7q/gnHPnJynZ3K11XTOtNDM8OszI6Ag16azN0dEG\nTpxo4+ihd7jk9yG2NGzgTPcgJ3tHuXkAuroq/8xykn6xLwkKh2dKnVPOewGo//kFLlz1Ce83t0Xx\nbEvnnCtUXVU9aZ5+VVVyUXhpyfq8s2rHj2a1/5uZTVuAFcD+vNebgI0F53wM3J4eXwCcKnYusD93\nXsH7zYsXL168VF5myuGFpZw7/XZgmaSlwEmgCXio4JyPgBaSsfq1QO55ZHuBHZJeIxnWWQZ8U/gB\nlU45cs45NzszJv10jP5p4CATUza7JL0ItJvZPqAN2J5+UXua5MKAmXVK2gV0ksxRempWz0V0zjn3\nv5gTi7Occ85lI/ijACStlnREUrekjaHbE5Kk45J+kPS9pCnDYPOZpDZJA5J+zKtbJOmgpF8kHZBU\nE7KNWSnRF5sl9Un6Li2rp/sZ84WkekmfSuqUdFjSs2l9dLFRpC+eSesrio2gd/rpwq9u8hZ+AU0F\nC7+iIekY0GBmZ0K3JWuSVgJDwLtmdkta9zJw2sxeSW8IFpnZppDtzEKJvtgM/G1mW4I2LmOS6oA6\nM+uQtBD4lmRd0Hoii41p+uJBKoiN0Hf64wu/zOwckFv4FSsR/m8ShJl9CRRe7NYA29LjbcD9mTYq\nkBJ9AcWnQM9rZtZvZh3p8RDQRbK6P7rYKNEXuXVPZcdG6ARTbOFX0cVbkTDggKR2SU+EbswccJWZ\nDW/aNz8AAAGDSURBVEAS8MDiwO0JbYOkDklvxTCcUUjSdcBy4GugNubYyOuL3O4GZcdG6KRf9uKt\nSNxhZrcC95L8EVeGbpCbM94AbjCz5UA/ENswz0KSLV6eS+9yo80TRfqiotgInfT7YNLOQfUkY/tR\nSu9YMLNB4AOmbncRmwFJtTA+nnkqcHuCMbPBvOnObwK3hWxPliRVkSS57Wa2J62OMjaK9UWlsRE6\n6Y8v/Eq3Z24iWdAVHUkXpVdwJF0M3AP8FLZVmROT//vbCzyaHrcAewrfMI9N6os0seU8QFyx8TbQ\naWZb8+pijY0pfVFpbASfp59OL9rKxMKv1qANCkTS9SR390ayaG5HTH0haSewCrgCGCDZtfVDYDdw\nLfArsNbM/gzVxqyU6Iu7SMZwx4DjwJO5Me35TNKdwOfAYSa2HnieZGX/LiKKjWn6opkKYiN40nfO\nOZed0MM7zjnnMuRJ3znnIuJJ3znnIuJJ3znnIuJJ3znnIuJJ3znnIuJJ3znnIuJJ3znnIvIfTpht\ntxU/43UAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fcc0cc91110>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# plot pmfs to see that they are indeed just different parametrizations of the same distribution\n",
"ks = np.arange(26)\n",
"plt.plot(ks, wiki_pmf(ks, wiki_r, wiki_p), c='b', lw=5)\n",
"plt.plot(ks, scipy_pmf(ks, scipy_n, scipy_p), c='r', ls='--', lw=4)\n",
"plt.plot(ks, mixture_pmf(ks, mixture_m, mixture_phi), c='g', ls='-.', lw=6)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(array([ 8., 47., 152., 244., 451., 602., 761., 876., 938.,\n",
" 904., 887., 802., 692., 558., 477., 404., 305., 253.,\n",
" 185., 144., 88., 61., 49., 35., 37.]),\n",
" array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.,\n",
" 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.,\n",
" 22., 23., 24., 25.]),\n",
" <a list of 1 Patch objects>)"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGAhJREFUeJzt3XmUVPWZxvHv2zStyCa7yiIuiHGJxoii4lhRA24Jjh5N\nNAmiHsd9magJcaKAJjPqmGT0mIlmXAJGRyXRQBKNjIEyxxVRRFEEVLSBlkaw2QSkm37nj7pg2XZL\n3aq+VdX9ez7n1LH6V/eter1e67l7mbsjIiLhqSh1AyIiUhoKABGRQCkAREQCpQAQEQmUAkBEJFAK\nABGRQG03AMzsXjOrNbPXs8Z6mNl0M1tgZk+ZWfes1+4ws0Vm9pqZHZw1fo6ZLYxqxrT+v4qIiMSR\nyxbA/cCoJmPjgKfdfSgwA/gJgJmdCOzl7kOAC4G7ovEewA3AMOBwYHx2aIiISPFtNwDc/Vmgrsnw\naGBS9HxS9PfW8clR3UtAdzPrRyZAprv7GndfDUwHTii8fRERyVe+xwD6unstgLsvB/pG4/2BJVnT\nLY3Gmo4vi8ZERKREWvsgsDXztzczTjQuIiIlUplnXa2Z9XP3WjPbBVgRjS8FBmZNNwCoicZTTcZn\nNvfGZqZgEBHJg7s3t7Ldoly3AIzPr8VPA8ZGz8cCU7PGxwCY2XBgdbSr6Cngm2bWPTog/M1orFnu\nroc748ePL3kP5fLQvNC80Lz48kc+trsFYGYPkVl772Vm1cB44GZgipmdB1QDZ0Rf3E+Y2Ulm9g7w\nCXBuNF5nZjcBs8ns+pnomYPBIiJSItsNAHc/u4WXjm9h+staGP8d8LtcGxMRkWTpSuAylkqlSt1C\n2dC8+IzmxWc0Lwpj+e47SoqZebn1JCJS7swMT+ggsIiItDMKABGRQCkAREQCpQAQEQmUAkBEJFAK\nABGRQOV7LyARNmyAzZvj1ZhBd/0ShEhZUABIXjZvhkGDoKEhXt2GDfDkk3Dcccn0JSK5UwBIXhoa\nMl/mGzbEqzv7bKitTaYnEYlHxwBERAKlABARCZR2AUlxLV8O9/4Fnnk5Xt3IkXD66cn0JBIoBYAU\n13vvwuBGOOSQ3GsWLIA771QAiLQyBYDkr3ELTH4wXs36frD//nDhUbnXzJwJc+bE+xwR2S4FgORn\n40b41OHpp+PVdbuYqdUH8cF/xKh5bxD/tOYAYkSGiORAvwcgedmwcgO9+8AG3ylW3YsvwtSp258u\n2wezavnw1Rpm1n0tXqFIQPL5PQBtAUhRDR+eecQx85c13DirAf7xj3iF3bvDQQfFqxEJiAJAMqvl\nl18er6a+IxBz90++Bg6Ejkvhp1fHq5s9O3MAeeDAZPoSaeMUAAKzZsGee8I11+Res6kCRu6YXE9Z\nKnfpzRsVvTmtd7wtgJ07TOHXa+vplFBfIm2dAkAy+vWDYcNyn34DEGtvY/5GjIBJk+DTT+PVXTzt\nG/xb7SfstX8yfYm0dQoAYUujMevDwfBC7jWbNiXWzheYwcknx6+71j5p/WZE2hEFgHDfi/tx/RPD\n2GNpvLozz0ymHxEpDgWAsLG+kjP3nsMdLxxT6lZEpIh0MzgRkUApAEREAqUAEBEJlAJARCRQCgAR\nkUApAEREAqUAEBEJlAJARCRQCgARkUApAEREAqUAEBEJVEEBYGb/ambzzOx1M3vQzKrMbLCZvWhm\nC8zsf82sMpq2ysweNrNFZvaCmQ1qnX8FERHJR943gzOz3YDLgX3dfbOZPQKcBZwE/MLdp5jZb4Dz\ngbujf37s7kPM7DvArcB3C/43EPkyr7wCWxbGqzniCOjSJZl+RMpIoXcD7QB0NrNGoBNQA3yDTBAA\nTALGkwmA0dFzgD8Adxb42SJfrmtX+NOfYPqHudcsXgzf+x5MnJhcXyJlIu8AcPcaM/sFUE3m96Gm\nA68Cq929MZpsKdA/et4fWBLVbjGz1WbW090/zrt7kS/RsW9Pzu84mc5VMYoq3+HaRTNIJdWUSBkp\nZBfQzmTW6ncH1gBTgBObmdS3ljR9i6zXPmfChAnbnqdSKVKpVL5tSsCmToV33olX8+j1G3lq8T4K\nACl76XSadDpd0HsUsgvoeOC9rWvwZvY4cCSws5lVRFsBA8jsFoLM1sBAoMbMOgDd3L2uuTfODgCR\nfA0dmnnE8fodq1i7Npl+RFpT05XjiXnstizkLKBqYLiZ7WhmBhwHvAnMBM6IpjkHmBo9nxb9TfT6\njAI+W0RECpR3ALj7LDIHc+cAc8ns0vktMA74oZktBHoC90Yl9wK9zWwRcFU0nYiIlEhBZwG5+0Sg\n6XbHYuDwZqb9FNDPiIuIlAn9KHx74w6NjdufrmnNF47Ri0h7p1tBtDennAIdO0JVVe6PPz0O3bqV\nunMRKTJtAbQ31dUwdy4ceGDuNXcAMU+XFJG2TwHQzjy6ZhQv/ueu0Dv3mjlz4uWFiLQPCoB25qYV\nF3JChbPrgNxrBgyAU09NricRKU8KgHZozLfqOPD0PqVuQ0TKnA4Ci4gESgEgIhIoBYCISKAUACIi\ngVIAiIgESgEgIhIoBYCISKAUACIigVIAiIgESgEgIhIoBYCISKAUACIigdLN4ESaWLVsE2//+2Ox\nanp/dTd6nzI8oY5EkqEAEMlyyGmDmfSjHfjHLbnXNDQavmEj725Jri+RJCgARLKMunAwb18Yr2bl\nglXs+5WOyTQkkiAdAxARCZQCQEQkUAoAEZFAKQBERAKlABARCZQCQEQkUAoAEZFAKQBERAKlABAR\nCZQCQEQkUAoAEZFAKQBERAKlABARCZQCQEQkUAXdDtrMugP3AAcAjcB5wELgEWB34H3gTHdfE01/\nB3Ai8Akw1t1fK+TzRcpCVRX4Zth773h1lZXw2GOw337J9CWyHYX+HsDtwBPufoaZVQKdgeuAp939\nVjP7MfATYJyZnQjs5e5DzOxw4C5AP6EkbV/XrtCjEf72t3h1F10ECxcqAKRk8g4AM+sKHO3uYwHc\nvQFYY2ajgWOiySYBM4FxwGhgcjTtS2bW3cz6uXttAf2LlIeKivhbAF26JNOLSI4K2QLYE1hpZvcD\nBwGzgauAbV/q7r7czPpG0/cHlmTVL4vGFAAtWLcOGhri1WyhAtBvE4rI9hUSAJXAIcCl7j7bzH5F\nZk3fW5jemhlradrgLVwIBx0EO+4Yr65r40706bk6maZEpF0pJACWAkvcfXb09x/JBEDt1l07ZrYL\nsCJr+oFZ9QOAmubeeMKECduep1IpUqlUAW22TXXvr+GrO9by0rDL4hU+/zzs9koyTYlI2Uin06TT\n6YLeI+8AiL7gl5jZPu6+EDgOeDN6jAVuif45NSqZBlwKPGJmw4HVLe3/zw6AYC1ZApu2wDXXxKvr\n2hWGDk2mJxEpG01XjidOnBj7PQo9C+gK4EEz6wi8B5wLdAAeNbPzgGrgDAB3f8LMTjKzd8icBnpu\ngZ/d/nXoACNHlroLEWmnCgoAd58LDGvmpeNbmD7m/gwREUmKrgQWEQmUAkBEJFAKABGRQCkAREQC\nVehZQCIC1NfDM8/Eq6lctS9HevNXSIoUgwJApEDdumXO1r3hhnh18166jgfnzOOEf06mL5HtUQCI\nFKiqCqZMiV936q7z2VTfofUbEsmRjgGIiARKASAiEigFgIhIoBQAIiKBUgCIiARKASAiEigFgIhI\noBQAIiKBUgCIiARKASAiEigFgIhIoBQAIiKBUgCIiARKASAiEigFgIhIoPR7ACIlVP1RJ958M17N\noEHQtWsy/UhYFAAiJXJ0r7e4e9rJ3P1/63OuWbexAyNGwEOPdUqwMwmFAkCkRK6+xrj6oe/HqvnL\nO/ty1/NjgEOTaUqCogAQKZWxYzOPOG6YBXc3JtGNBEgHgUVEAqUAEBEJlHYBFcGyZTBjRryad17Y\nGVidSD8iIqAAKIobb4S5c2HIkBhFK3fk4r6PAQck1ZaIBE4BUASNjXD++XDBBTGKXnoXrvgrcENS\nbYlI4HQMQEQkUAoAEZFAKQBERAKlABARCZQCQEQkUAWfBWRmFcBsYKm7f9vMBgMPAz2AV4EfuHuD\nmVUBk4GvAyuB77h7daGf3ya8+y7MmAL3T8u9Zt066NUruZ5EJHitcRrolcBbQLfo71uAX7j7FDP7\nDXA+cHf0z4/dfYiZfQe4FfhuK3x++atZlrkI4Ibb4tXtuWcy/YiIUGAAmNkA4CTg58APo+FjgbOi\n55OA8WQCYHT0HOAPwJ2FfHab07cPHHlkqbsQEdmm0C2AXwHXAt0BzKwXUOfuW29XuBToHz3vDywB\ncPctZrbazHq6+8cF9iASjJ071zNj5dfZY494db17QzoNnTsn0pa0UXkHgJmdDNS6+2tmlto6HD2y\nedZrn3uLrNdEJAcj9q9jQeoiGu75Xay6ww6DtWsVAPJ5hWwBHAV828xOAjoBXYH/ArqbWUW0FTAA\nqImmXwoMBGrMrAPQzd3rmnvjCRMmbHueSqVIpVIFtCnSvgzstBJibgFU6qYv7U46nSadThf0Hnkv\nFu5+HXAdgJkdA1zt7t83s0eAM4BHgHOAqVHJtOjvl6LXW7w/ZnYAiIjIFzVdOZ44cWLs90hivWAc\n8LCZ3QTMAe6Nxu8FHjCzRcAqQjkDSKQ1VVTAnDkwZky8urV3Qn0noGMibUnb1CoB4O7PAM9EzxcD\nhzczzafAma3xeSLBGjkSbrsN6uvj1f1+A3y8EQb1S6YvaZO0Z1CkLamshLPO2v50TZ1b2/q9SJun\nW0GIiARKASAiEigFgIhIoBQAIiKBUgCIiARKASAiEigFgIhIoBQAIiKBUgCIiARKASAiEigFgIhI\noBQAIiKBUgCIiARKASAiEijdDlokBBUVMPJ4qFoVr+7OO+HUU5PpSUpOASASgE4DenFK15ep6ui5\nF9XU8POpz3Ksvv/bLQWASADSz1SwbFlVrJrfXlTHc9UDOTahnqT0FAAiARg0KPOI44lua5NpRsqG\nDgKLiARKASAiEigFgIhIoBQAIiKBUgCIiARKASAiEigFgIhIoHQdQFxz50JNTbya9QCdk+hGRCRv\nCoC4Dj0Ujj8+Xo1dCXvtnUw/IgmqXd+ZN9+MV9O3L/Tpk0w/0roUAHE1NMCTT8aruQDYLZFuRBJz\n2G5LGTdzFDPPzL1m82bYYQeYNy+5vqT1KABiWk4/6ubHq6mrS6YXkSR9a58FfOuAxXD99TnXVFfD\niBEJNiWtSgEQ0z4spP9p8WrM4Nprk+lHRCRfCoCY1tGN+TG3AEREypFOAxURCZS2AESkZatWwaJF\nuU9fUwmNu6N1y7ZBASAizTvkEPjRj+Cvf829ZmMfWPFHYNfE2pLWY+4xfiKuCMzMy62nbGZQxu2J\nlFT1C8sYcTRUN/QvdSvBMTPc3eLU5L2dZmYDzGyGmb1lZm+Y2RXReA8zm25mC8zsKTPrnlVzh5kt\nMrPXzOzgfD9bREQKV8iOugbgh+6+H3AEcKmZ7QuMA55296HADOAnAGZ2IrCXuw8BLgTuKqhzEREp\nSN7HANx9ObA8er7ezOYDA4DRwDHRZJOAmWRCYTQwOZr+JTPrbmb93L22gP5FpNy4w7PPxqupqoJh\nwzL7WKVoWuUgsJkNBg4GXgS2fam7+3Iz6xtN1h9YklW2LBpTAIi0Fz17QtV6GDcuXt28efDnP8PR\nRyfTlzSr4AAwsy7AH4Aroy2Blg6RNhftzU47YcKEbc9TqRSpVKrALkWkKDp1gj6d4m8BjBwJGzcm\n01M7lU6nSafTBb1HQQFgZpVkvvwfcPep0XDt1l07ZrYLsCIaXwoMzCofADR7X+XsABARkS9qunI8\nceLE2O9R6NUa9wFvufvtWWPTgLHR87HA1KzxMQBmNhxYrf3/IiKlk/cWgJkdBXwPeMPM5pDZnXMd\ncAvwqJmdB1QDZwC4+xNmdpKZvQN8ApxbaPMiUl66doX166Ei5qpllf2Fl0+fzYHJtCUtKOQsoOeA\nDi283Owvprj7Zfl+noiUvx49YOXK+HWj+szjw493UAAUmW4FISKtKu7aP4A1fz6IJEx3bBIRCZQC\nQEQkUAoAEZFAKQBERAIV7EHgxYvhxhvj39rZaES5KSLtQbAB8Pe/w7vvwnnnxasbM+l4Mjc5FRFp\n24INAIB99oGxY2MWnTsziVZEgtatcgPn/3IYXSbHq7voIrjyymR6CkHQAUD9Zli7qdRdiATv/v1v\nY+nZP4Kjjsq55skn4bnnFACFCDcAXn8dHpgNj18Vr+5AXaso0tq69uzIV246G7p0ybnmjTWjoM/F\nwJDkGmvnwg2ANWtgyBBYsLbUnYjIfffBkiXbny7bze9BeiUKgPyFGwAiUj66dYP9949X03tNMr0E\nROcziogESgEgIhIoBYCISKAUACIigVIAiIgESgEgIhIonQYqIm3Xp5/Cs8/Gq+nVC77ylWT6aWMU\nACLSNu0+mOpNdTxy/iOxyvotfpHUx4/Fuuq4vVIAiEibdMTpu7HnrN14bMuNseoeX7iZlXUb6abv\nfwWAiLRNAwfCQw/Fr+v+6MbYvwPSXukgsIhIoBQAIiKBUgCIiARKxwBEJDh//ltHOveJVzNqFOy0\nUzL9lIoCQESCctkO9/DHiUdDheVc88bqgVRf6lx5864JdlZ8CgARCcrPZx4J7y2KVXP1T5fQMKMe\npnWK92FDh2YeZapdBMCqVVBdHa+m+uPOwCeJ9CMiZeyIIzKPGCqmfsBD0+t57ZIYv1rWUM+RFY9y\ncc31MRssHvMyOyHWzDxuT8cP+YDqZRV07hDjB97r67nqqNmc8/cxMTsUkdCsWAFPPRWvZtmcFUz6\n9Trmf7pXMk01YWa4e+77tWgnWwDrP1zH5IvnM/zkXvEKD/l2Mg2JSLvSty/84Afxaub3WMekX+fx\nYXV18X8fOU/tIgCAzM2djj2g1F2IiGzz0ZaeXB9zD9Dk2zZTvemr7GQbY1Tltyen/QSAiEgZ2WeP\neq6vvJm1fz08Vt0lHd7kgidGssM/5V73+0ue56LJcTtUAIiIJKLDfkO58sHDYO2aeIWdhsCoYbEu\n092hY2O8z4goAEREkmAGp59e6i6+VNFvBWFmJ5jZ22a20Mx+XOzPFxGRjKIGgJlVAHcCo4D9gbPM\nbN9i9tCWpNPpUrdQNjQvPqN58RnNi8IUexfQYcAid/8AwMweBkYDb2dP9LOfxXvTpZv7Aitap8My\nkk6nSaVSpW6jLGhefEbz4jOaF4UpdgD0B7JPcF1KJhQ+Z9NjT8R600s6zuFrQ79ZWGciIoEpdgA0\nd5XaF05g/dm/xLyvw04D4ahD82xJRKRt67VrVV51Rb0VhJkNBya4+wnR3+MAd/dbsqYpr3tTiIi0\nEXFvBVHsAOgALACOAz4EZgFnufv8ojUhIiJAkXcBufsWM7sMmE7mDKR79eUvIlIaZXc3UBERKY6y\n+k1gXST2GTN738zmmtkcM5tV6n6KyczuNbNaM3s9a6yHmU03swVm9pSZdS9lj8XSwrwYb2ZLzezV\n6HFCKXssFjMbYGYzzOwtM3vDzK6IxoNbNpqZF5dH47GWjbLZAoguEltI5vhADfAy8F13f/tLC9sp\nM3sP+Lq715W6l2IzsxHAemCyu381GrsFWOXut0YrBz3cfVwp+yyGFubFeGCdu/+ypM0VmZntAuzi\n7q+ZWRfgFTLXEZ1LYMvGl8yL7xBj2SinLYBtF4m5ez2w9SKxUBnl9d+naNz9WaBp8I0GJkXPJwGn\nFrWpEmlhXkDzp1S3a+6+3N1fi56vB+YDAwhw2WhhXvSPXs552SinL5jmLhLr38K0IXDgKTN72cwu\nKHUzZaCvu9dCZuEH+pS4n1K71MxeM7N7Qtjl0ZSZDQYOBl4E+oW8bGTNi5eioZyXjXIKgJwuEgvI\nke5+KHASmf+gI0rdkJSN/wb2cveDgeVAaLuCugB/AK6M1n6D/Z5oZl7EWjbKKQCWAoOy/h5A5lhA\nkKI1Gdz9I+BxmrllRmBqzawfbNv/2f5u/pQjd/8o64ez/wcYVsp+isnMKsl84T3g7lOj4SCXjebm\nRdxlo5wC4GVgbzPb3cyqgO8C00rcU0mY2U5RsmNmnYGRwLzSdlV0xue3CqcBY6Pn5wBTmxa0Y5+b\nF9GX3FanEdaycR/wlrvfnjUW6rLxhXkRd9kom7OAIHMaKHA7n10kdnOJWyoJM9uDzFq/k7lY78GQ\n5oWZPQSkgF5ALTAe+BMwBRgIVANnuPvqUvVYLC3Mi2+Q2efbCLwPXLh1H3h7ZmZHAf8A3iDz/4YD\n15G5o8CjBLRsfMm8OJsYy0ZZBYCIiBRPOe0CEhGRIlIAiIgESgEgIhIoBYCISKAUACIigVIAiIgE\nSgEgIhIoBYCISKD+HwhJUfUNR9SYAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fcc0cb72690>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"### simulate samples directly vs from a hierarchy\n",
"# directly\n",
"direct_contacts = stats.nbinom.rvs(n=scipy_n, p=scipy_p, size=10000)\n",
"\n",
"# hierarchy\n",
"gamma_shape = scipy_n\n",
"gamma_scale = (1-scipy_p) / scipy_p\n",
"hierarchy_lambdas = stats.gamma.rvs(a=gamma_shape, scale=gamma_scale, size=10000)\n",
"hierarchy_contacts = stats.poisson.rvs(mu=hierarchy_lambdas, size=10000)\n",
"\n",
"# histogram\n",
"bins = np.linspace(0,25,26)\n",
"plt.hist(direct_contacts, color='r', histtype='step', bins=bins)\n",
"plt.hist(hierarchy_contacts, color='b', histtype='step', bins=bins)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2 (SageMath)",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.10"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment