gibbs_sampler.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Franklin/anaconda/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n"
     ]
    }
   ],
   "source": [
    "# Gibbs sampler for the change-point model described in a Cognition cheat sheet titled \"Gibbs sampling.\"\n",
    "# This is a Python implementation of the procedure at http://www.cmpe.boun.edu.tr/courses/cmpe58n/fall2009/\n",
    "# Written by Ilker Yildirim, September 2012.\n",
    "\n",
    "from scipy.stats import uniform, gamma, poisson\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "from numpy import log,exp\n",
    "from numpy.random import multinomial\n",
    "\n",
    "# fix the random seed for replicability.\n",
    "numpy.random.seed(123456789)\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50\n"
     ]
    }
   ],
   "source": [
    "# Hyperparameters\n",
    "N=100\n",
    "a=2\n",
    "b=1\n",
    "\n",
    "# Change-point: where the intensity parameter changes.\n",
    "#n=int(round(uniform.rvs()*N))\n",
    "n = 50\n",
    "print str(n)\n",
    "\n",
    "# Intensity values\n",
    "#lambda1=gamma.rvs(a,scale=1./b) # We use 1/b instead of b because of the way Gamma distribution is parametrized in the package random.\n",
    "#lambda2=gamma.rvs(a,scale=1./b)\n",
    "\n",
    "lambda1 = 3.479722 # We use 1/b instead of b because of the way Gamma distribution is parametrized in the package random.\n",
    "lambda2 = 0.3134539\n",
    "\n",
    "lambdas=[lambda1]*n\n",
    "lambdas[n:N-1]=[lambda2]*(N-n)\n",
    "\n",
    "# Observations, x_1 ... x_N\n",
    "x=poisson.rvs(lambdas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At iteration 0\n",
      "At iteration 1\n",
      "At iteration 2\n",
      "At iteration 3\n",
      "At iteration 4\n",
      "At iteration 5\n",
      "At iteration 6\n",
      "At iteration 7\n",
      "At iteration 8\n",
      "At iteration 9\n",
      "At iteration 10\n",
      "At iteration 11\n",
      "At iteration 12\n",
      "At iteration 13\n",
      "At iteration 14\n",
      "At iteration 15\n",
      "At iteration 16\n",
      "At iteration 17\n",
      "At iteration 18\n",
      "At iteration 19\n",
      "At iteration 20\n",
      "At iteration 21\n",
      "At iteration 22\n",
      "At iteration 23\n",
      "At iteration 24\n",
      "At iteration 25\n",
      "At iteration 26\n",
      "At iteration 27\n",
      "At iteration 28\n",
      "At iteration 29\n",
      "At iteration 30\n",
      "At iteration 31\n",
      "At iteration 32\n",
      "At iteration 33\n",
      "At iteration 34\n",
      "At iteration 35\n",
      "At iteration 36\n",
      "At iteration 37\n",
      "At iteration 38\n",
      "At iteration 39\n",
      "At iteration 40\n",
      "At iteration 41\n",
      "At iteration 42\n",
      "At iteration 43\n",
      "At iteration 44\n",
      "At iteration 45\n",
      "At iteration 46\n",
      "At iteration 47\n",
      "At iteration 48\n",
      "At iteration 49\n",
      "At iteration 50\n",
      "At iteration 51\n",
      "At iteration 52\n",
      "At iteration 53\n",
      "At iteration 54\n",
      "At iteration 55\n",
      "At iteration 56\n",
      "At iteration 57\n",
      "At iteration 58\n",
      "At iteration 59\n",
      "At iteration 60\n",
      "At iteration 61\n",
      "At iteration 62\n",
      "At iteration 63\n",
      "At iteration 64\n",
      "At iteration 65\n",
      "At iteration 66\n",
      "At iteration 67\n",
      "At iteration 68\n",
      "At iteration 69\n",
      "At iteration 70\n",
      "At iteration 71\n",
      "At iteration 72\n",
      "At iteration 73\n",
      "At iteration 74\n",
      "At iteration 75\n",
      "At iteration 76\n",
      "At iteration 77\n",
      "At iteration 78\n",
      "At iteration 79\n",
      "At iteration 80\n",
      "At iteration 81\n",
      "At iteration 82\n",
      "At iteration 83\n",
      "At iteration 84\n",
      "At iteration 85\n",
      "At iteration 86\n",
      "At iteration 87\n",
      "At iteration 88\n",
      "At iteration 89\n",
      "At iteration 90\n",
      "At iteration 91\n",
      "At iteration 92\n",
      "At iteration 93\n",
      "At iteration 94\n",
      "At iteration 95\n",
      "At iteration 96\n",
      "At iteration 97\n",
      "At iteration 98\n",
      "At iteration 99\n",
      "At iteration 100\n",
      "At iteration 101\n",
      "At iteration 102\n",
      "At iteration 103\n",
      "At iteration 104\n",
      "At iteration 105\n",
      "At iteration 106\n",
      "At iteration 107\n",
      "At iteration 108\n",
      "At iteration 109\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEPCAYAAACzwehFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8VNX1wL+HQAIECDsBwhpWQUBsQeqGu+Lexapgq1Zt\nq1XrVrfaQG2t9WertmpbqrVapVp3QHEHW1xAZd9kDXsCBEIgIfv5/XFnkpnJzGQmmWQmyfl+PvNJ\n3nv33nfufe/dc+85dxFVxTAMwzDC0SreAhiGYRiJjykLwzAMo1ZMWRiGYRi1YsrCMAzDqBVTFoZh\nGEatmLIwDMMwaiWhlIWI3CwiKz2/m+Itj2EYhuFIGGUhIqOAHwHfAMYB54nI4PhKZRiGYUACKQtg\nJLBIVUtUtQL4L/DtOMtkGIZhkFjKYhVwooh0EZH2wBSgX5xlMgzDMIDW8RbAi6quE5GlwC6gAtgJ\n2FokhmEYCUDCKAsR6QNMBNJUtVRE1gBpQcKZAjEMw6gDqip1jZtIZiiANkCqiAzCmaBeDxZo6tTp\nwGFcx8P7O8zUqdNR1Rbzy8rKirsMifKzsrCysLII/6svCaMsVHUXUALsA9YDn6nq7GBh77//SjIz\ns4BCz5lCMjOzuP/+KxtBUsMwjJZHwigLEekMbAV6AG2BQhG5PFjYQYMG8P77NzJ16sMATJ36MO+/\nfyODBg1oNHkNwzBaEg3qsxCRLkA/VV0RQfDTcb2KD3F2pa7AFBHprqp/8g04ffp0AIYMAVjA889n\nxVDqpsPkyZPjLULCYGVRjZVFNS25LBYsWMCCBQtilp7Ewpbll6DIAuACnCL6CtgDfKKqt9YSbwLw\nNPBNnDnqGeAi4GhV3e4TTn1lFoEYZ8HwYcuWrdx33z/ZubOSvn1bcf/9V1oPzjCaICKC1sPB3RA9\nizRVLRCRa4DnVDVLRGrtWajqYhF5BVgKlAE5wCpfRWE0Llu2bOWMM/7Mpk0zgFSgkM8/zzKTn2G0\nQBrCZ9FaRHoDlwBzo4moqjNUdaSqjgG2A7MaQD4jQu67758+igIglU2bZnDfff+Mo1SGYcSDhuhZ\nzADeBRaq6hee9Z02RBpZRNJw5qiLgONF5CtVXeQbxuuzcEz2/IxYs3NnJdWKwksqu3ZVxkMcwzCi\noCn4LI5X1U9qOxcm/j+Bg8Bw4DygvaoW+FyPic+irrb4UPEa07bfWPeaNm0GL7xwO/4Ko5CpUx9u\nsYMKvJgvx2hq1Ndn0RATP5ZEci5E3E7AJuDfwA9DhFFfAg4jYvPmbM3MvE3hsDpVc1gzM2/TzZuz\n6xTv448X1im9ulBX2RP9Xk0JKxejKeKpO+tet9cnsl9CMAm4DedruNXnNx1YHmEaY4EvgGJgGTAT\naBcQJqAAoi+0qVOn+3zoWvXBT506vU7xBg78dp3Sqwt1lb2ubN6c7bmnu7dViI3/DAwjFtRXWcTS\nZ5EMdMD5QTr6nC8AvhthGq1x+1lswM21uADYDfjZPOrrs6irLT5UvPz81DqlVxca248waNAAnn8+\nixdeoMWbnryYL8doCsTaZxEzZaGqHwMfi8g/VXVrHZPZAZQDE1X1gIicANwZGMhXWcyYEf1N+vZt\nhVsqxN8W36dP+MFhoeJ17lxIfn706dWFuspuxA57BkZTYPLkyX6TEmfUpbL0oSEc3MOA24GB+Cgj\nVT01wvhHgJPVzbvIwjm47/S5rr4yi8DmzdXOxrS0AlRbU1DQPqTjMdj8gczM4PMHfB2ZaWkFLF1a\nybZtD/jF+8c/Lubqq1/3S69fv1s45pi0sHIEEnivYPmIRvZI7xWJfMEGEtTXyRtrJ3Es0oskjVg9\nA8NoTBLRwb0c+CkwATjW+4si/g5cs60IWIKb5BfWZ1HtbMxWuCUix2Mktvhgjsx+/a7VCy64vUY8\n3/QuvPDn2r//TVE5QP3vFT4f9fUj1MVBG+gbqq+TN9ZO4likF00a5ssxmhokioO7KkH4qp7xe3v+\n9sA5uU8IuB5QAOrzcUfveAznIA/nyAwVz1t5RCuHf5zI4tfFuV9bvkIReK/6Onlj7SSORXqxKBfD\nSFTqqywaYlLeHBG5HrcXRYlPD2Z/hPFzRWQJrofxOq6HstA3QGgHdzAHbN0dj7F2hIeL5x8ntvmI\nhXyxTiPWTuJEyJNhJBIJ6+D24Yeev3f4nFNgcG0RPXtv3wisAboAZ+JmhPtR08Ht62yMneOxNkdm\nMPs2DIjKAepNY82alTHPRyj7e035tgJPsXr1VqZNmxGRrT9cHiOx+wePv5YtW1ZxyilZUfscYpFe\nXZ+3YSQisXZwx9wMVZ8fbq7GIdzQ2QLgriBhArpWGmBnjsxn4Rs/FOFs2DXve9hzHLntO7yfovZ8\n1FX2mteiv1d9JyjWjL9GW7f+YVTPLtbpRVdm1c/bMJoCJKDP4gfBfhHGfRkYB5wMzA4RJqAAajqX\ngzmgQxdg+AIO5cj095X4H4eL50tNG3m2wi8jzkdd/S2B8tXVPxIsj9HY/X3jx2JiYyzSC/XcQuXL\nlIXRVEhEZfFnn9/fgc3AKxHEOxd4Eljk07PIChJOs7Kyqn4w36cwNOj/4QuwbuH8Kw3/X6TpT578\nq1rTCBe/LmmfcsqvIspHsHCRyBHpfQPj1yVeQ6YXmN9InpVhJBLz588PqCvrpyxi7rNQ1Rt9jz3b\npb4YQdTjcQsHVuC2VW0F/ExE5qnqYt+AHRasZ9q000lP78mXHIa5c+GEE4DONRLd9dLLzPrL6+Tl\nKampRagmUVSUQt7IAdz9+5/itTl77dFdV2+hf1ppVfpe0giWfiEnsIQ0DnqOi53oc3HydK4pD//7\nHxx04c/VDaTyCtCWhZzAQToTaDM/gf/B3IPk5Ozh+ec/IC9P6dZN+N5jdwFj/GTfubOSU9ts59pL\nTvBLG/Ck36aG3+ME/kcauVXhvPno072EoPjI7+VcYOvy/tz7f6/X8L14059cucE9p6ob+5eP119Q\nXZ57gPfovGQPD0++nO89dhcDxtbM76TytbQrLaeoKIVu3YRp007nXHrSp0cJC4L4H05pvR3mzvUr\nz9TUItZ07c/uI11r+CK870/fZVs4l6NrlGdNh7hhJAYJPymvxg1E2uA2MRoeRZyTgV8AvYCfquoX\nPtd0DmfTvv16jps4kI/mt+O8c4FHH0WGDsGbHe9kvY+PuZjuB3sBpcAW3GK2SfycByDzOTZtepjN\nm6snWT3KPWSyvir99u3bATDkrUfZqEN8ZITMzNu5cVMRmWwDKmjffj1FRSOr5GHIkJqT2W6+GTZt\nAqCo6AifL8qmqGgYP+fPbKI3mZlZbNr0cFWcx+Rmrj1lVVU4SAIqeKR/Hz7a9rSf7JDKo1zP0e3f\nZezYdJYvz6mK481v4MSxx+Rmjm4/1y/t9u3XM3TuP+h3ysl++VX1l9/L3Lfgkf69+Gjbn3A74z4G\n3O8nj29ZQnX5+D6rM874s095rql6Vt78PrXgVwA++d3Ho1xEJul+shcVjWTsR3dwyrVzakyc++rE\nPNps3eZTnu69+DlvsIkxVeG870Uk708Df0KGERMScVLeHGC25/cWzgz1YBTxW+HMUGXA74Jc97M/\nhzLZeG3Otc1dqBlO/dIPlrb3OJh9OxozlGrtaYSTrzbZI/GbhMpHsHChqOn3qPa9ROM3isSPEukz\nDVW2qvV5L2rmy8xQRlOBBPRZnOzzOx7IiCJuBvARsNrzWw8cFRBGIUshSwcOPEnD+SzC2Zl9f5HY\nt8MpgXA+hkgrk3DphctHXWSP5L7hwgW7Vh97frAyC5WnSJ9puPTr+15EWmaGEU+ags/iYxHpBXzT\ncyriXfJwiwjeqqrLRKQDkI0bTXWXf7DpQCHHHy9kZ08OmVjNcfM1bdiQmvALw4WSL7FkDy5fLNOr\nzlPtzzQc9X0vaqRfXg6tG2LKkmHUnVj7LGJeo4jIJcBi4Hu4fbgXiUikS5SX4wzD4BzdAAdrBnN2\nZeeIDM39919JZmYW7gMHuM/n/0LPtWDh1tKhw2Vs2lTEtGkz2LKlrovoOrZscZPdTjkli4suuoUL\nL7yDU07JijjtmvLVTXZfOWKRL19CyRd430jzHypPhw/n07//PdT2TAPxlSHSNMKVu5eRrIFvfCNc\n0RhGVDTkd1ov6tMtCfbDLSTY0+e4B5FvfnQ0bvHAZcBa4ADQISBMjfkOXoL9X9scjGDhOnS4ysdW\nHXzyVTRmqEgWCIwmH3WRPZJJZfUxQ4WSL5IFEoOVWbg8+S7mGO6ZBn8G0aVRm08piTLVtDTVvXtD\nF5BhREhD7sJIPc1QDaEsVgYctwo8Fybu00AusAr4ErgwSJiAAqj9/2jDRTL5qm7pRbYwYV3yUZvs\nkUwqq4+yCPV/NA7pUGnEoswaqtxBVadMUX311aBlYxjR0JC7MNZXWTSEofUdEXkXt482wPeBtyOM\n+wxuYt4C4Jeq+mawQPXdKa92ai4mF5v0soJei+1CdcFlD7VIXkMTyQKJtROLxf0acIHAk0+GBQvg\n29+OTXpGiyWWi1km7EKCIjIE6KWqd4jIt4ETPJc+A16IJA1VXSgirwLFqvpYqHD13SmvdhrSWdvQ\nzuh6OmtjTKTO5PDEoswasNxPPhmuuy42aTVBIt10KhbhIr3mu3lYJBui1SdfsSSaRT5rky9hFxLE\nzVs+Osj5o4E5EaZxPM6xfQRYivNfnB0QJqBrVfv/0YaLZMG4uqdXN59FfWVvaJ9FqP8jWbSwtjyG\nW8yxPuUSi3IHVS0tVR05UrWoKGj5NGfqtmhm3cJFfs3XNxb5hmh1kTfWRLrIZ902MKufGSqWyuKL\nMNci8ll4wg4AVoS5HlAAtf8fbbhoJ9tFk16oBQJjpSwidYSHmgAYiroqi8D7RuqQri1Pkdw31mlE\nEq6lEamNPRbhIr8W/UZidc1XQxD95NTI5EskZbEhzLWNEaZxNrARt2nSnSHC1GkhwYZQKvEPN79B\n7xUqTrhriVAW0SrYWIZLBObPn99o94p0wcZYhKvLNZgf0X3rmq+GJHh+ai8LL4k8Ke9LEblWVf/u\ne1JErgG+qi2yiLQCHsdtnvRX4DIReVNV1wWGbXifRVNhAbF37ofe5CceNtzIWUBDlEWkRGovv+66\n05k584MGDbd8+WLGjRsf0/uGurZmzSbCbToVy3BFRUooe3779oUE940twM0PDn/fwLLIzl4V8l4X\nXXRLvcqstnL3X8wy+rKYNm0G1113Ok899TE7d2pVegmzkKBn1vbruBXXvMrhG0AycLGq5tQS/zjg\nNc9hN5zf4j1VvSQgnPrK7LtQX6j/m2+46ahOj+m9Ahcm9F1YLzPz9hrn33//RgYPHpBQZRFpeg1T\nZvs4n5u4nX24OaY7cR99EnOYwKOt91Fe/gS+Cy6ez7+5nelV4aAIYTuzuYmHmUb1wowuzvkcx+08\nXhVOGQKUM4cUHmYecC9uxLpbzPF8nuAOfl0VzivTHE7mYYqrwsFaWrf+PeeUT+F2/uAnOxTxVit4\nqPJj/BeL9M1vkp9M/un75vdDbufXfrLDTubw04D8OpkuanUDt1SWBClPPPmtlr28/AlP+tcDA3iG\nHVxFhSe/F/Ew51WF8y/PcPlNBW6NMr/hytM3v/f6ye6e13/4ZdtbKS4eTNu2mykuru39qc5vzffn\nXqq/4T+g9VhIMOarzorIKcBoz+FqVf0ownjfAc5S1es8x9OACap6U0A4UxYNqCymTp3BCy/cTvAR\nSzVHFE2d+jAvvJCVUGXR2MrCv8xmkM40hrITNxL8UsCttpvDy2zgwapw4OKkcztDOaMqnDdeDoPY\nwL+rwnnjpFPAUDYEpP8MOdzMBsbhelhv4X1W/ulXx/GXh4D07/GT3T/9atld+l8wlD8D19A7/T52\n59wTIn3X+k3nek5Kv7QqnH9+h1WFg/sZOPA7FGc/5ClPV4puLM015NDJI4+XtQwc+EuKsx/nquNv\nQkli1iftmHr8YZQkXv3k15QNvJvs7OeClGdgfiGHQsoGPkV29qtUv/u+z7f6WUVXnoH5DSzP3Qxl\nBdCOl17M4amn5vH+B1fRuu9vmL/z9aDl6Q1f8/0Z5jlXCHRILGVRV6JRFvGQzzAMo6lTH2WRSKuf\n7QT6+xxneM75UZ/MGoZhGHUjMZZVdXwBDBGRASKSjOuzzY6zTIZhGAaNrCxEpJWILBGRYErgBCAd\n+Bq3/3auqq5tTPkMwzCM4DS2Gepm3H6ZnUJc/0hVL2hEeQzDMIwIaLSehYhkAFNwwxxCBqsljbNF\nZJ2IrBeRO2MqYIIjIhki8pGIrBaRlSJyk+d8FxF5T0S+FpF3RSQt3rI2BoG91BZcDmki8rKIrPW8\nGxNbcFncIiKrRGSFiLwgIsktqSxE5GkRyRWRFT7nQuZfRO4WkQ2ed+fM2tJvTDPUI8AdQLjRTJNE\nZJmIvCUiR/le8Jm0dxYwCjdpb0SDSZt4eHcRHAVMAm7w5P8u4ANVHY7bkvbuOMrYmHh7qV5aajk8\nBrytqiOBscA6WmBZiEgf4EZgvKqOwVlNLqNllcUzuPrRl6D599SvlwAjgXOAJ0UkbGO9UZSFiJyL\n80Esw/Ueggn1FdBfVcfhlMIbAdcn4JYU2aqqZcCLwIUNKHZCoao5nvJDVQ/jNofKwJXBs55gzwIX\nxUfCxiNEL7UllkMn4ERVfQZAVctV9SAtsCw8JAGpItIaN+lgJy2oLFR1IW7DOF9C5f8C4EXPO5ON\n2/56Qrj0G8tncTxwgYhMAfoBIiLPqeoPvAFU9bCI/ElEzsEzA0ZEuqrqfk+QvsD2wHkWIvJgI+Uh\nEfEqyxzfRkELmotyOvjlt0WWQ4h8tsiyAPb7/H++52+LK4uAPPrm36tM+uK2j/Cy03MuJI3Ss1DV\ne1S1P/BnYD6wz1dRAIjIZUCmqg7Fda27+yiKwPQa/Ldhg7J1a8Pfpz6/rKysuMsQ6e+dd5TDh+se\n/6abFFC2b2/6ZWHvRfMui9JSJTc3/nkP/Hmrz7rW4/FwcL/lc+7HIuLdNeZ64GgRWQpcC+zyrDfl\nJXDSXkj27YMPP6y7rF9+CRMnwpVX1j2NWJCdDUOHwqmnwr33wuzZ8MwzcO21cNRR8Mor/ktSePnw\nQ9i8udHFDck778DZZ8Nvf1u3+M8/D2+95Taie/nl2MpmtCxU4auv4IEH3PfVEOn/6EfQrx9873vw\n8cfBv9FYcfvtMGEC3HUXvP8+HDlSa5Q9nr87cVYeL0EnQfsSDwf3MmAxgKr+TVVneq7nA5eq6jGq\n+i1gE/7doi+AIbXdZMUK+OY34fLL4YYboKQkOiE/+wymTIG//x22b4ePIlrZqmG46y64+GK44w5I\nSoInnoAPPoBx4+Bf/4KcHJg3zz/Ozp2uUr34YigurpnmokUN+/IGsmsXXHUV/Pvf8Le/waZN0cVf\nuhRuuQXeeMNtRvfiiw0jZ6x58UX4z39qD/fRR/DTn8Ktt8I998Cjj0JubsPLF469e6P/bhKd4mK4\n+27IzITvfx9WrYLTT4fdu8PHW74cDh2K/D5//KNLe9s2t4Hij38MY8fCk0/CwYP1y0Mgzz0Hc+fC\n734HKSluBe4hQ5wyDIN3q+rZwKWeEWODcHXr4rAxG6n7cyFuBbClwGbg6yBhPgEO4XbHW4Lb12J8\nQJizwX8/C9+1+19/XbVHD9VZs1Tz81Uvvlh1wgTV7GyNiI8/dvHnzXPHzz+v+q1vqVZWRhbfl8pK\n1SefVH3zTdWKiprXDh4MH//TT1X79lU9fDh0mAcemK/Dh7uN2rx897uq992n+p3vqN50k3/43/xG\nFVT//Ofo8vHss6p/+5vqkSORx1NVLS9XPeUU1Rkz3PFvf6t60UWRx9+3T3XgQNWXXnLHpaWq3bur\nbt5cM2xj7uFQG2+9pZqerpqZqfqzn6mWlAQPt3evaq9eqr/7neof/uDK58orVTt3Vr30UtX581U3\nbFBdvVp16VLVQ4dqplFervrtb6tOnqx6xx2urN58c36NcMXFqmefrTp7dnjZN21ysg8f7u4fa2bP\nVh0zRvXll8N/V0uW1C5rJHjfi1tvVT3rLFeO3vv+5jeqRx+tun9/aBnatXP1SCR1wLx5qr17q27d\nWn2uokL1/ffdd9m5s+qPfqT6z3+qvv226pdfqubl1UynstKV/X/+o/rFFy5M4P1XrnTfwsqV/udf\ne82df/vt4PtZAJ21uj6921PPrgXO1Nrq8doCxOIHPABs9yiK3bitU98KCPMm8KXP8Trcnt6BaQV9\nUH/4g2pGhurixf6F/vDDqj17qj76aOjKbtUq1csvd4X8wQfV58vLVY86yhV8tMyY4V7EY49VHTZM\n9S9/cWnfeKNq//6qbduqPvdc8LgVFaoTJ7qXKhyVlapnnqn62GPu+K23VIcMcfncv9/dZ+5cd+3h\nh1WHDlVduNApxK++qj0P+fmq3/ue6tixquee6yqR3/5W9cCByMvglFNcOao6uQYNch9Pbaxc6eT9\n5S/9z193neqDD4aOV17uKqJgFfTu3arvvlv77qcvvugqFl8lHClLlrj36NNPXTmdf77qpEmq27fX\nDDttmuott9Q8v3+/e1/HjVMdPNjt2Dp8uOpxx6mWlfmH/ctfXIPmnXdc5XfuuU5J7d7tH+6GG1S/\n8Q33jRQUBJd9zx5X5k8+6Rpe/fo55RWYVm289ZZL55NP/M/n5Djl+MgjLm8TJ7pvYvNmpxTXrFF9\n4gnV8eNVBwxw+bj7bv+KsrhY9bbbXL4j5aOPVPv0cY0PXyornRKZNKlmo2zvXifDv/7lyu2RR8Lf\nY906913973+hw+TkqD70kOoPfuAU9/jxqmlpTtm/9ZZ7Z19+2dUZI0Y4JTVunGqnTu5ZPPCAy8Oh\nQ+56qPrhk09cOT/1lP956rn5UTycLGfiTE7fDDj/C9zwWoDjgM9DxK/R4njySfdR7dgRvPCWL3cf\nbb9+qjNnqq5Y4Votjz2meuGF1a27YK1978OLpncxc6aTJyfHxfv4Y3efiRPdB71qlWstDhyo+utf\n10x71ix3z8AeSTBWrXIv6fbtriJ+993qa//7n6vgZ8xw17wV1osvOqXiW2mUlKguW+Yq6XXrVD/8\n0MW5/vpqJbtypeoVV7j7zZwZXr5PP3X33rnT//xrrzkFHFjpBea/e/fgyvTDD1WPOSZ4vMpKJ2+3\nbu5DXLu2+tpLL7lGw4QJqh07qp53nuo//lGz7HftcuG++U1XUflSWupa/JdeGrx3s22b6w2+8kr1\nuYoK95Gnp/u31OfNc88/XM/Rl4oK1TPOcO+Ll5wcV07Ll/uH/fWvXevdq9RfeMFVvAcOuMo/mII6\ndMjl+d57q88VFKjefLNqhw5OyZxzjlPee/cGl7G83MXv29cp9PR0pwRUXTmfd151mVZUuIp49Ohq\nxTB8uOoll6i+9567vnevq6h/8hN3vH69e/bnneee8Zo1tZdbfr5rNIVq8FVWql59tZPj7bfdcVmZ\n6qmnqt55pwuzebN7Jz77LHga+/e7BuHMmbXLE8jBg67XPmGCakqKqyPeeMP/26qsdI0Qb69z9GjX\nQwnHunWuPDdtqj7XZJQFzj+yFDcsdoPn3I+B6zz/n4zb8KgEtzbUd0Kkoz17qr76qiuA5593L6dv\noYTi00/dBzdypOqUKa5i+dvfVAsLQ8epqHDa/fXXa09f1T3o3r2rP5Jw7N7tKrWrrnIffHa2qwD6\n91f9738ju5+qy0efPq4SC2TGDKckt2zxP3/ddaqXXebuec897sMeOVJ11Cj34o8c6V/p+bJsmWvl\nTppUs6JSdcplxAinaAOprFQ97bRq05Qv5eWuIhs82N0jGOXlTrmvX1/z2kMPud5cfr5reXbr5lro\nl17qPpxFi1y4vDynkMaPd5WhV2FUVjqlfu+9qrm5rky9Pc3yctf7nDLFyd61q+rtt7sK66WXnNwD\nB7oeXDDefddVOH/8o6uEBwzwV+yRsGOHS8ObjyuucDIEUlmp+vOfux7H4sVOoXjLc88el8aSJdXh\n9+93Ld2rrgreKKqocN/XG2+o/vSnLr6voi0vV/38c1fBnnaaKztV1b/+1fUw9u5V/fvf3XcUyiQX\nioMHnYnttNNcPp54wt33iSfcO+jttYbiiiuczOGorHTf94gR7l5XXOF67L5pv/mm+y4DeyelpS7f\nwRRwtOTm1t4ozclx73a4OstXNl+ahLIAUoBFHmWx2mOOOiogTAfgSdzkkE3A1hBp6XXXZWlqapYe\nd1yWdukyX1etqr3g6sPcua5yOPNM18r54x9rvjSqrtLo0cPZGSPl0CHXLR092lXqaWmuBRENe/eq\nnnSSaxUHUlkZ3PxWVOQq1i5dnG8jklaaLxUVrjLo0aNmD+Cee1zXOhRbtzqFdM891R9HUZHrdp9y\nSmgbspcbblC9/37/c//+tys/X3PPunXuQ7755uCmpwMHXEv1rrucHLNmOWVZXOyuv/++Uxi5uU65\nTp5cnc6uXarXXOMaKhdc4Hqm//1v+I9982ZXYQ4c6J55XXj5ZVcBz5nj8hvMj6Hqns8Pf6jaqlVN\nc8VTT7mWbGmpaw336uUq1EjNbl995Vr8J51Ubb4dNcr1JgIr7zvvdD2W7t21zt9pUZHzxyxd6p+/\nyZOd+dnL7t2u0r78cve76CJXVpH23srKnFI755zgvoRf/MJ9p++9544rK917ce65tSuteBDrPbgb\ns2fR3vM3yeO/eCTg+jlePwYwESgGugZJR1WdKWn8+OpWVkOzcaOzK/7pT+4jTE935hxvt/XWW53Z\nxrfFlujk50fWQgnHmjXORPHXv7rjJUtcy7M2O/eePa4SueYa9//xx7uejreiDsf//uc+WlWnEB58\n0CmtYL2c2ti716V1662u0vT1eak6RdK7t6tcQ9n6o6Gw0JkigzU2IsWrBF57LXy4srLgvZeKClfe\nGRnub13e2fJyp3T+9jd/h26we11zTXQ+hkjZuNH1HlescKa+bt3cc3z++epfYI+6PlRWut72kCHO\np+VVHrUNVkkUmoSyALoDaZ7/u3pGPd0cEOZZ4Pue/yfg9vKO2MHd2Hz2mbO9X3CBa+GcfXbw1khL\nYNMmpyhpxYBtAAAgAElEQVQffNC1nGtzzHs5dMiZBdu3dy3HSHw0qi5cRoYzV/Xp43p79VHSu3c7\nM1Wgj0LVtbanT0+sZ+u1c9dllJ6XTZtcxVefNBKBRx5xivOiiyIz/caC0lI3onDixNgqo4amvsqi\nUbZVFZGjPcpgBM4ktVBVTxaRH3syMFNEVgLtcf6KI7ilSH6iqksC0tLGkDkSSkrg4YfdHAjvXIiW\nyo4dbtz6oEHw9ttub+pIKC2FTz+FyZOju9+iRdCmDRxzTOT3qk2ONm1ik5bReFRWwtdfw8iR8ZYk\n8RERNNH34PbM3n4O8M7IbgNcpKprfMJ8AozB+SzA7XlxSSIrC8Ofw4ddZZuaWntYwzAal/oqi8Za\nSNC7vPYyEekAZAM/wC2f62UfbrLeNwBEZB0hpp9Pnz696v/JkyczOdpmqdEgdOgQbwkMw/CyYMEC\nFixYELP0Gqtn0R0oU9WDItIO5+D+g6r+zifML4DbVLWXiBwHPKqqxwVJy3oWhmEYUdJUeha9gWc9\nGxil4EZE/dnXZ4EbWttJREpwcy2uaiTZDMMwjFpolJ5F1c2cCWoBcL+qvhnkWqWqFnn2tHhMVYcF\nScN6FoZhGFESt56FiLQBpgRW+iHCZgD/Ar6BGzY7IDCM1r75URXmszAMwwhPXH0WItIROBe4GLcW\n+lxVfSCCeOnA33EjnX6F20L1QlVd5xPmMmCaqp4rIlcBf1XVlCBpWc/CMAwjShp86KyIdAG+h1tm\nPBXXM7hPPftBRyjk8cB/gZW4ZXIHAQ/itvjzzrP4H67HkYebZ9EbOE5VcwPSMmVhGIYRJY1hhnoC\n+BC4WlVzRSQFmCEiZaq6OpKbqOonOKc2IjIQ57d4XFUP+wTLB+5U1U894T7AbX4U561gDMMwjEiU\nxfWqmu89UNUSEbkH+KOIfKWq/6otARF5GjgP2Itb8+lmX0UhIifjNjYaJiKFwGvh0jOfhWEYRngS\nap6FiHxXVV+JINwJONPSAuCXqvpYwPWTccuB3KmqL3nOrQNONjOUYRhG/amvGapee3BHoig84RYC\n9wDFgYrCh1zcrG48k/LyAxWFYRiGER8aZVKex8F9EVAqIktxTu57cA5tBb4GMoHBHjPUduDyxpDN\nMAzDqJ3GmsF9NbAf6KCqxwRe9EzI+w9wBm7V2XaBCwj6Yj4LwzCM8CSUzyLimzifRUfgNVVtF+T6\nOcDPPHMsJgIfA32CTcgzn4VhGEb0xNVnESken0VBmCCX4pYwB2eWaoVbxtwwDMNIABqrZzELOA3o\nAewAsoBkotz4yJOW9SwMwzCipEmsOquql4vIAGCOqo4JEiQb+F3AhDzDMAwjQWgsBzfAybhJd+uB\np1X19z7XdgJnisjbwGbgKNwSI0Gd3ObgNgzDCE9TdXC3wimBYuBo4AvgUu9CgiIyBbgPN8P7AUJs\nfOQJa2YowzCMKGkqe3C/i+tZCG7y3SLcWlHveDY+QkReA84CNgJXhRo6a8rCMAwjeprEaChgJvCc\nqqaoan/gTWCHV1F4eAzn3FbgfhE5qpFka/IUlhby/Irn2VO4J96iGIbRTGlMn0WGZ72nVjhfRGDN\n9hXVE/NGAPMIskkSROezKK0o5Z2N73BC/xPo2q5rPcT3p6CkgOz8bLYc2IKIMLjLYAZ3GUz7Nu1r\nhC0pL2HXoV10TOlI9/bd/a4Vlxfz5ro3qdRKMrtmMrjLYLq164aIfwNAVdlRsIOu7bqSmpwKOJvk\nod6HuHHejQzoPIAb593I90d9n9sm3cbQbkOjys++on0UlRXRP62/3/lKrWTRjkVUaAXH9j6Wdm1q\nTJOJmpLyEgpKCujevnuNfPpSVFaEqlblNxwLFiyo1Xe15cAW8o7kcWzvY8Pe15fSilJ2H9rNmr1r\nWL13NWv2rqFd63aMTR/LuPRxjOg+go7JHcOmd+DIAb7c9SXH9jm2xjvo7SVHKk8kBJZFcXkxq/as\nYtWeVRyTfgxj08cGjaeqFJYVkl+cT3JSMp3bdiY5KbnOcqgq7256l/vm30fH5I7cc+I9nDbotDrl\nNb84n+0Ht7Pr0C52H96NqtIjtQc9U3vSvX130lLS6JTSiTZJ/iPuI3kv4kFFZQXbDm6juLyY4d2H\n00qCt9sLSgrImp/FnPVzmNRvEicPOJnjMo5jT+EeVu9Zzeq9qxncZTA/PvbHpLVN84vbVH0Wk4AP\ncI7rXcBWYJaq3u4TJqKJeSKib657kyNlRzhUeoi1e9eyaq/7EDK7ZHLt+Gv57lHfJTkpmVkrZ5G1\nIIuu7bqyJX8LPzn2J9wy6Ra6t+/O9oPb+XzH56zas4pDpYc4VHKII+VH6NOxD0O6DmFo16F0SulE\nUVkRRWVF7Cvax7KcZSzLXcbynOUUlhUysPNABnUehKJsPrCZ7PxsUtuk0q5NO9q0akObpDbkFeVR\nUFJA7469yS/OZ1LGJC4bfRkT+k7gXyv+xVNLnmJs+li6tO3CpgOb2LR/ExVawYC0AQzoPIDu7buz\nIW8Dq/ason2b9hwqPcTonqM5od8JvPX3t9DJyl/O/QunDjqVPYV7eHzx4/zly7/QtV1XBqQNYGDn\ngQzuMpjRPUczqsco+qf1Z33eehbtXMSiHYtYuWcl6/ato7yynOSkZNLapnHOkHM4vt/xfLr9U15d\n+yqd23YmNTmVVXtWMbL7SIZ3H86hkkPsP7KfQ6WHyOySyfje4xnfezxd2nZhX9E+9hXto7SilKN6\nHMWYXmNIa5vG8pzlPL30aWatnEWlVnKk/Ah9O/alT8c+dG3XlS7tutApuRPbCraxas8qdhbsJKlV\nEhP6TuDszLM5acBJtJJWFJcXU1JRQptWbWjfpj2pyan845F/8H+//T+/iqigpIClu5fy7qZ3mbN+\nDnsK99AppROVWslloy/jOyO/Q68OvWjbui0pSSls2L+BT7d/yifbP2Fl7kp2H97NweKD9EztyYju\nIxjdczRH9TiKorIiluUsY3nuctbnraesooyu7brSI7UHpw06jW+P/DbH9zuePYV7eOTzR3h66dMM\n6zaM1XtW0z+tPxP7TiS/JJ8NeRvYdGATbVq1YUyvMYzpNYZ+nfqxJX8LX+d9zab9m0hOSqZb+250\na9eNbu270bWtK6cubbuQmpxKaptUOiR3oG+nvgzqPIgu7brwy1/9khN/cCLvbHyHD7d8yIb9GxjW\nbRgju4/kk+2f0KdjH64dfy2je45m4baF/Hfrf1m8czF5R/LcO5CSRlllWZXSGNp1KFOPnsrUMVPp\n07EPAHsK97B091LKK8vpkdqDHu170DGlI6UVpZSUl7CjYAczPp7B9oLtPHDqAxSVFfG7hb+jY0pH\nfnTMj0hLSSOldQrJScm0btWaVtKKJEmipMI1JA4WH2T34d0sy1nGkt1LyDuSR/+0/vTp2IfeHXrT\nSlqxt2gvewr3sLdwLwUlBRSUFJDSOoU+HfswIG0A/dP6k/d2Hv954j+ktK6xj5ofJeUlFJYVUljq\nlOX6vPV8nfc16/PW0zG5I0O6DmFI1yEoyuc7PuezHZ+xas8q0lLS6Jnak56pPTmqx1FMypjEcRnH\n0bVdV/KO5JGdn012fjbbDm6r+m3cv5GN+zfSrX032rRqQ0FJAScPPJmT+p/EUT2OYmi3ofTr1I+X\nVr/EHe/fwTlDzuH6b17PV7u+4uOtH7No5yLSO6Qzuod7HxftXMS8jfO4dvy1/HDsDympKOHAkQMc\nKD7AmZln0iG5g7fubBI+i28B7+OUxW7cUNlAZfEs8LaqviQiE4CFQL9gq86eN+s82rVuR4fkDgzv\nNrzqI16as5SZX83ky11f0r19d7q178bvTvsdJw04iez8bB5c+CAvr3mZdq3bUVJRwqSMSYztNZbO\nbTvTIbkD7dq0Y2fBTjbu38iG/Rs4XHqY1ORU2rdpT+e2nRnTcwzH9D6Gcenj6N2hd40WUqVWkleU\nR0lFCaUVpZRWlNK1XVe6t+9OK2lFYWkhc9bP4d+r/s2iHYv4/qjvc8OEGxjWzX+r8fzifLbmbyU7\nP5u9RXsZ2nUoo3uOplv7bhwpO8IXu75g4baFzH9mPnNnzq3xIRSXF7P5wGa25m9l68GtbNy/kdV7\nV7Nqzyp2H9rtKqyMiUzsO5Fj0o9hRPcR9EztiaKsyF3BvA3zWLh9IRP7TuR7R32PkT1GAnCk7AhL\nc5ayIW8Dndt2pmu7rnRI7sD6vPUs2b2Er3Z/xaHSQ/Ro34Nu7bvRWlpX3Tc1OZWUpBSuPuZqrhp3\nFQM6D6CwtJAdBTvYdWgXB4oPkF+cz8Hig2R0ymB0z9EM7TaU4vJi5m+Zz7yN8/h8x+cktUoiJSmF\nlNYplFeWU1RWRGFpIVve2ELyackM6zaM3h16s3bfWnYf2s3onqM5ddCpnD/sfCb0nUAracXSnKXM\nWjmLuevncrDkIMXlxRSXFzMgbQDH9zueb/X7FuN7j6dPxz50b9+dpFZJYd/vkvISDhQfYGfBTuZt\nnMdra19jR8EOyivLuWLMFdz2rdvon9af8spylucsZ/HOxXRt15Wh3YYypOsQSitKWZG7guU5y9lR\nsIPBXQYzvPtwMrtkUl5ZTt6RPPKK8sg7kseBIwfYf2Q/+4/sd5VbWSGHSw+zo2AHWw5sIalVEkXv\nFTFhmlOwZ2SewdheY6vekYrKCt7Z+A4zl8wkOz+bE/ufyEkDTmJSxiR6dejl15NQVYrKivhy15f8\na8W/eG3tawzrNowdBTsoKivimN7HkJKUUlVpHy49TEqSUwAdUzpywzdv4EfH/KiqtV+plbyx7g3e\n/No19korSimpKKGisoIKraCisoKU1ilVvYSeqT0Z22ss43uPJ7NrZsjWt6+8hWWF7CzYydaD7vt5\n6LcPUXR8ETdPvJlpY6axJX8LS3cvZVnOMrIPZrOzYCc7D+2kuLyYDskdSG2TSlrbNIZ0HcKIbiMY\n2m0oh0sPs3H/RjYd2ESlVjKx70SOyziOMb3GUFhayJ7CPeQW5rI8Zzmf7fiMxTsXU6EVpCSlMLDz\nQAZ0HlClvPqn9SezSybDug2r6jHvKNjBx9kfs3DbQtblrWND3gb2FO5hVM9RPDnlSSb1mxQ23wDZ\n+dn84dM/MHfDXDqldKJz2850aduFx6c8TkanDKDpKIvvAD/CLRbYCmdy2gOsJsqJeZE4uLPzs9ma\nv5WTBpxUo0LfdWgXRWVFZHbJjGnXPx5Mnz7dzyQXCSXlJbW2smJNRWUFOwp2kNEpo9aKt65Mnz6d\nG39xIxv2b2DXoV2M6D6CYd2G0bpVY1paq9mav5XU5NQaZseGRFXJO5LHww88zIO/eTDm6R8pO8Kn\n2z9lcJfBDOw8sEl8P9OnT+fin1zMw589zOyvZzO823DG9x7PuPRxZHbJpG8n17NNS0mLWX4qKiso\nLCukU0qnOqdRVFZE29Zta1WQ0dCUlMVZqnqd53gaMEFVb/IJM4eaE/N+EUxZNLjAhmEYzZCEn8GN\nm3Tn6znN8JwLDNOvljD1yqxhGIZRNxpr6OwXwBARGSAiybiFA2cHhJmNbX5kGIaRkDTW2lAVIvIz\n4D2cgnpaVdeKyI/x+CxU9W0RmSIiG4FC4KrGkM0wDMOonUbz/qnqO8BwABF5WkRygVzvwoIi8hBu\nZdrDwCbcTG7DMAwjAWgsM1Qgz+CW9vDlPWCUqo4DNgB3B0YSkbNFZJ2IrBeROxtBzoRBRDJE5CMR\nWS0iK0XkJs/5LiLynoh8LSLvikhabWk1B0SklYgsEZHZnuOWWg5pIvKyiKz1vBsTW3BZ3CIiq0Rk\nhYi8ICLJLaksvI1wEVnhcy5k/kXkbhHZ4Hl3zqwt/bgoC89mSAcCzn2gqpWew89xDu4qPIsRPo5T\nMqOAy0RkRCOImyiUA7eq6ihgEnCDJ/93AR+o6nDgI4Io2WbKzcAan+OWWg6P4eYnjQTGAutogWUh\nIn2AG4HxHmtFa+AyWlZZBGuEB82/ZzmlS4CRwDnAk1LL2OF49Sxq42rcch++TAA2qOpWVS0DXgQu\nbHTJ4oSq5qjqMs//h4G1OIV6IfCsJ9izwEXxkbDxEJEMYArwlM/pllgOnYATVfUZAFUtV9WDtMCy\n8JAEpIpIa6AdbjRliymLYI1wQuf/AuBFzzuTjbPmTAiXfnxmLIVBRO4FylR1VsClvsD2wHkWIhL7\n2UdNB6+yzPFtFLSguSing19+W2Q5hMhniywLwHd5oPM9f1tcWQTk0Tf/XmXSF/jMJ8xOz7mQJFTP\nQkSuxLUYLw8XTlWb7S8rKyvuMlj+LH8tLW8tIX/e6rOu9XNclIWIPA18CQz1OfddnE+iJzA7iCMq\ncGKfYRiGER3e1b4jmgTtS7x6Fv1xGi5FRLaJyFW4PS/KgEPAaGB+QJwvgCGNKqURMenpbq2g9PSB\n8RbFMIzQvOn5Oxu41DNibBCubl0cLmK8RkOdAXwTWKWq/dU56PYAI1R1PHAsblFB3zgVwM8aXdhG\nJhHX3o+E3NytgHr+hqap5i9SmnP+mnPeoPnnz8ODAKq6Brd/0BrgbeB69bFVBaNRFhIMemORAcAc\nrZ6Ut19Vu/pc9zv2OV9bnow44BxoCgj2fAwj8ajvqrMJNxrKh5A1TjQ75RnxIT19ILm5W+nVawA5\nOdnxFscwWhxNcqe8oDeu2bNYC0xW1VwRSQfmq5toFBjPehYJSGDPwnoahpFYNOWehXh+XmYDT3mc\nLd2A/SKSrKqlcZHOMAzDqCJeQ2dnAZ8Cw3xGQ/0Tt5BgG2Al8DVuKXPDMAwjzsRrNNTlqtpHVVN8\nRkMdBPYCx+Em5rUBdsVDPiM8NkzWMFoeCTODW1V3AX8AtuEmh+Sr6gfxlcoIRqTDZA3DaD4kjLIQ\nkc64tY4GAH2ADiISdtkPwzAMo3FIpKGzpwObVXU/gIi8BnwLCFxQ0IbOxpHoTU8ppKcPtOGzhtHI\nNJuhs4GIyATgadzM7hLc2uxfqOoTAeFs6GwcqV690n9YbLihs4ANnzWMONOUh876oaqLReQVYDnQ\ny3P6myLypaouiqNohmEYLZ6E8VkAqOoM3Brrt6hqZ9zOX2vjK5VRG5GapmwUlWE0XRLGDAVVO38t\nVdXMMGHMDBVHgpmhQp/zN0PZrG7DiB/1NUMlVM8CGATsE5FnRGSJiMwUkXbxFsowDKOlkzA+Cw+t\ngfHADar6pYg8ittwPMs3kI2GMgzDCE+zHQ0FICK9gM9UdbDn+ATgTlU93yeMmaHiiJmhDKNp0mxG\nQwF4VpzdLiKrgU3AV7jNOQzDMIw4kmg+C3Bb+/UCTsaNhnogvuIYUD2SyTCMlklCKQsRyQDGAd8F\nPlbVb6vqwTiLZVC9HpRhGC2ThFIWwCPAHVitZBiGkVAkjLIQkXOBXFVdRs2NkQzDMIw4kkgO7uOB\nC0RkCtAO6Cgiz6nqDwID2tDZRCEliB8j0nMO26vbMBqGZj101ouInAzcpqoXBLlmQ2fjQOBQ2Or/\nw50LfS3UAoSGYTQMzWoGt4hkiMhHuC1WJ4vITXEWyTAMwyDBlAVQDtyqqoNwGyDdICIj4ixTiyM9\nfSBJSamISNXfhrqPYRhNg4Q0Q3kRkTeAP6vqhz7nzAzVwATOyI7M5BRd+FAzvw3DaBia1QxuX0Rk\nIG7Ohe1l0QgsW7aMefPm0bNnz3iLYhhGApKQykJEOgCvADer6uHA6zYaKvbce++DzJt3hFat3muE\nu4UeMWWjogwjNjT70VAi0hqYC8xT1ceCXDczVAMwZcqlzJt3EW3b3khx8T7P2YYzQ0UyUsowjNjR\n3EZDnQ3kAROBtnEWJy7EsiWQmCyItwANSnN+fs05b9D881dfEkZZiEgr4CmgA7AVyBKR1R4F0mJo\n/i/sglqupzTprVeb8/NrznmD5p+/+pJIPosJwEpVzQAQkbsAVdV34iuW0biUAEpurq32YhiJRML0\nLIC+wHaf4x2ec0ajUYRqZbyFMAwjAUkYB7eIfAc4S1Wv8xxPAyao6k0B4RJDYMMwjCZGc5lnsRPo\n73Oc4TnnR30yaxiGYdSNRDJDfQEMEZEBIpIMXArMjrNMhmEYBgnUs1DVChH5GfAeTok9rapr4yyW\nYRiGQZx6FiLytIjkisgK3/OekU8zgUzPX8MwDCMBiJcZ6hngrMCTnj24z8DNs6iBiJwtIutEZL2I\n3NnAMjYq3uXZPXNLVjbH5dlFpJWILBGRZmdeFJE0EXlZRNZ6nuHEeMsUS0TkFhFZJSIrROQFj6m4\nyRKswSoiXUTkPRH5WkTeFZG0eMpYH0Lk7yHP+7lMRF4VkU7RpBkXZaGqC4EDQS559+CugWfS3uM4\nJTMKuKyZLV/uXZ59FDCJ5rk8+83AmngL0UA8BrytqiOBsUCzMaGKSB/gRmC8qo7Bma8vja9U9SZY\ng/Uu4ANVHQ58BNzd6FLFjmD5ew8YparjgA1Emb+EcXCLyAXAdlVdGSLIBGCDqm5V1TLgReDCRhOw\ngVHVHM/+43gWT1xLM5pn4uk1TsHN0m9WeFpoJ6rqMwCqWq6qBXEWK9YkAametdvaA7viLE+9CNFg\nvRB41vP/s8BFjSpUDAmWP1X9QKsnUn2OG3EaMQnh4BaRdsA9OBNU1emAYH2B7YHzLETkwQYWL55c\n2FAbD8WR06F5zpdpjnkKYL/P/+c1l3cz4Lnl+OarOTxT3zz4TD24Gtfgjpi4ObiBL4GhnlOZwGgg\nT0RKgAHAEhEJurmCqgb9hbvWVH5ZWVlxl8HyZ/lraXlrCfnzqX/vBcpUdVY09XY8Hdw/9B6o6ipc\nl6+tqqYABcBrqrrHJ07gpD3DMAwjCkTkSpw5+PJo48ZLWVyPUxgpIrJNRK5Sf3taMW4Pbl++AIY0\nppCGYRjNjDuAC1S1JNqI8RoNdTnOYb1KVfurxzHow2LgzYA4FcDPGknEuNHcd/2z/DVdmnPeoPnn\nz0MH4H3PEPYno4kYt4UERWQAMEfdUDzf8/fihuh9J0Q8DSWzZyeomMtqGIbR1KnvTnkJMRrKi489\n7dRw4Wrbg3v69Ol+YQzDMFoazWIPbs9oqAuADqraznPuu8A/gd3AZuASVT0YJG6tPQvrYRiGYfjT\nVPfg7g8oPg5u3FpQZcAh3DDa+dEkmJ6RXn2QFHBsGIZh1IuE8VmIyDrgZFXNFZF0YIGq1ljuIlTP\nwjuRxtuz8P5vGIZhNN2eRTB6qmouuKUvgKAT8gzDMIzGJ6Ec3AGE7BaEdHCb+ckwDANoJg5uCGqG\nWgvMAr6P6/H0BXqoamlAvLBmKF/MDGUYhuFoymYowX+xwI+AnwPjcSs+7qLpL4NsGIbRLIjXQoKz\ngE+BYT6joR7HLX28GrcyaTZNfBlkwzCM5kLczFDB8OwO91ugCHhPVa8IEsbMUIZhGFHSlM1QfohI\nZ9zmIwNwiwh2EJGoV0Y0DMMwYk8ijYY6HdisqvsBROQ14Fs4p7cftS33YRiG0dJpNqOhAhGRCcDT\nwGnAE7hd84qAi1V1kU84M0MZhmFESbNZSFBVF4vIK8B63JIfbwI/BpLjKphhGIaROD0LqNr4fqmq\nZoYJYz0LwzCMKGk2Dm4Pg4B9IvKMZ3OOmSLSLt5CGYZhtHQSxgzloTVuUt4NqvqliDwK3AVk+QYK\n6+BOAioaXlDDMIxEptk6uAFEpBfwGZAP7AAeAu5U1fN9wpgZyjAMI0qalRnKs+psJbDdc+o0YE38\nJDIMwzAgwcxQIpIB7AWGA72BcuCqOiXmWYE2Z0dO7AQ0DMNooSSaGepl3HIfacBtqnpBkDARm6HA\nTFGGYRjQjMxQInIukKuqy6i5Iq1hGIYRRxLJDHU8cIGITAHaAR1F5DlV/UFgQN/RUH964k8cLDjY\naEIahmE0BZr7aKgM4DncfItuwC9V9U8BYfzMUKHMT14SKX+GYRjxor5mqERTFulAOs5n8QtgCHCh\nqq7zCWPKwjAMI0qazdpQAKqaA3iHL30sIm/gtlddFzqWYRiG0dAkjIM7EBEZCIwDFoUPaRiGYTQ0\nCdWz8CIiHYBXgJtV9XDgdV8Hd1hsroVhGC2U5u7gPht4FOgHfBjJPIvafBZgfgvDMIzmNM+iFfA4\nsBr4O9BfREaECv+XmX/h6WeebizxGo1YtgQSEctf06U55w2af/7qS8IoC2ACbqmPi4DJQA/gQ09v\nowbX//h6rrn6mtpTTYKklCTSM9JjJ2kD0txfWMtf06U55w2af/7qSyL5LPoCK1V1EoCITAMmqOo7\noSK0SW5DWWlZ+FQroLKiktydubGU1TAMo0WRSD2LqOjZuyf9BvaLtxiGYRgtgoRxcIvIccB0VT3b\nc3wXoKr6+4BwiSGwYRhGE6NZzOAWkSTga9weFruBxcBlqro2roIZhmEYieOzUNUKEfkZ8B7OPPa0\nKQrDMIzEoMF8FiLytIjkisgKn3NdROQ9EflaRN4VkTSfa3cDf8btlHeDqj7YULIZhmEY0dGQDu5n\ngLMCzt0FfKCqw4GPgLsBROQo4BJgJHAO8KQEmW0nImeLyDoRWS8idzag7I2OiGSIyEcislpEVorI\nTfGWKdaISCsRWSIis+MtS6wRkTQReVlE1nqe4cR4yxRLROQWEVklIitE5AURSY63TPUh2sZsUyNE\n/h7yvJ/LRORVEekUTZoNpixUdSFwIOD0hcCznv+fxc2pALgAeFFVy1U1G9iAm3dRhc+kvbOAUcBl\n4SbtNUHKgVtVdRQwCbihmeUP4Gaa757qjwFvq+pIYCzQbEyoItIHuBEYr6pjcObrS+MrVb2JuDHb\nRAmWv/eAUao6DlfHRpW/xh4621NVc6FqhdmenvN9ge0+4XZ6zvkyAdigqltVtQx4Ead8mgWqmuPZ\nJRDPelhrqVkGTRbPXiVTgKfiLUus8bTQTlTVZwA8jZ6COIsVa5KAVBFpDbQHdsVZnnoRZWO2yREs\nf6r6gapWeg4/BzKiSTPeDu5ohmL1BbYHDp0Vkebs27gwkrWvmhinQ/McAt0c8xTAfp//z2su72bA\ncxKowycAAAjESURBVMvxzVdzeKa+efAZOns1rsEdMY3ds8gVkV5QtdHRHs/5nbjFA71keM4FRVVb\n7C8rKyvuMiTSLxHK49133+err76KuxyJUh6J8rOy8P95EZF7gTJVnRVN5d3QykI8Py+zgSs9//8Q\neNPn/KUikiwig3A75C0OSGsn0L/hRDWMunHWWWdw6aVN3YRvtARE5EqcOfjyaOM2mBlKRGbhFgTs\nJiLbgCzgQeBlEbka2IobAYWqrhGR/+Ccn2XA9eqrCh1f4JSIYSQcR44cibcIhhEJdwAnqWpJtBEb\ncjTU5araR1VTVLW/qj6jqgdU9XTgr0Af4L8+w/D+Cmz2yHRb4LA1Va0AftZQ8jYVJk+eHG8REgor\nD3+sPKqxsghKB+B9zxD2J6OJ2OjLfXiG4S0ERqhqqYi8BLwNHAXkqepDnjkUXVT1riDxg3Q6DCN+\niAgZGRls37699sCGESea6uZHvsPw2uH8Ec1m2JphGEZzo9GVharuAv4AbMMpiYOq+gHQS4PPwTAM\nwzDiTKPPsxCRzrhexADgIM7hPZWacy5C2pqmT59e9f/kyZPNNmkYhhHAggULYrr7Xzx8Ft8FzlLV\naz3HVwDHAacCk1U11zMHY766pRMC45vPwkgozGdhNAWaos9iG3CciLT1LBZ4Gm7IbKg5GIZhGEac\naXQzlKouFpE5uNnbyUAxMBM3IupLEbkfKADGN7ZshmEYRnDiNRqqD3CzqrYFugOrgZ8Av1fVZOD/\ngOvjJJthGIYRQKMrixArdB7Ehs4ahmEkLPHoWQwC9onIM55ZhDNFpD02dNYwDCNhiYeyaI3zRzyh\nquOBQtymIxEPnTUMwzAal3jsZ7ED2K6qX3qOX8Upi1wR6eUzdHZPqARsnoVhGEZ4mvw8CwAR+Ri4\nVlXXi0gWbuctgP2q+ntbG8poStg8C6MpUN95FvHaKe8m4AURaYNbafYq3HpR/wlcvtwwDMOIP3EZ\nOquqy4GJQCXQ2jMaSj3HeP637oNhGEaCEK95FgA342Zue7kL+EBVhwMfAXfHRSrDMAyjBnFRFiKS\ngdva7ymf0zbPwjAMI0GJV8/iEdz2fr6mJptnYRiGkaDEY4nyc4FcVV0mIpPDBLUlyg3DMOpIkx86\nKyIPANOActwueR2B14FvYEuUG00QGzprNAWa3BLlqnqPqvZX1cHApcBHqnoFMAdbotwwDCMhiedo\nqEAeBM4Qka9xe1w8GGd5DMMwDA/x8FlkAM8BvXDzKv7uuWTzLAzDMBKUePQsyoFbVXUUMAm4QURG\nYPMsDMMwEpZ4+CxyVHWZ5//DwFogA5tnYRiGkbDE1WchIgOBccDn2DwLwzCMhCVeCwkiIh2AV3Db\nqx4WkYj3s7B5FoZhGOFp8vMsAESkNTAXmKeqj3nOrcXmWRhNEJtnYTQFmtw8Cw//ANZ4FYWH2dg8\nC8MwjISk0ZWFiBwPTAVOFZGlnn24zwZ+D1wqIiXAr4AjjS2bYRiGEZx4jIb6RFWTVHWcqh6jquNV\n9R0gH7f0xzCgM3CxZ0it4UMsbZDNASsPf6w8qrGyiC2JNIN7ArBBVbeqahnwIm44reGDfQD+WHn4\nY+VRjZVFbEkkZdEX8PUQ7vCcM4wmhe9oPcNoLiSSsjCMJsmgQZmMHTu26njGjBlxlMYwGoa4DJ0N\nhogcB0xX1bM9x3cBqqq/DwiXGAIbhmE0MeozdDaRlEUS4F1xdjewGLhMVdfGVTDDMAwjfjO4A1HV\nChH5GfAezjz2tCkKwzCMxCBhehaGYRhG4mIO7gRGRDJE5CMRWS0iK0XkJs/5LiLynoh8LSLvikha\nvGVtLESklWci52zPcUsuizQReVlE1nrekYktvDxuEZFVIrJCRF4QkeSWVB4i8rSI5IrICp9zIfMv\nIneLyAbP+3NmbembskhsbO+PmtwMrPE5bsll8RjwtmcNtbHAOlpoeYhIH+BGYLyqjsGZ2C+jZZXH\nM8BZAeeC5l9EjgIuAUYC5wBPikhY57cpiwTG9v7wx7PL4hTgKZ/TLbUsOgEnquozAKparqoH/7+9\nO3atswrjOP590Gptl2JBi7SNkdKpUOgkdVCom5ClUArFQeruIA5m6V8gHcSli0OgLa2i6CJSOoij\npUNaFyGtLZFGRM3gICU8DucEkui9B+xLbpLz/Uz3vrlcznnywu8973vPOXRaj+opYG9dqPQ5YJGO\n6pGZ3wN/bDg8qv8zwNV63twHfqJMjB7JsNgm3PsDgIvAB6xfvr7XWkwDv0XEp/W23KWI2EOn9cjM\nX4CPgAeUkFjOzBt0Wo81XhjR/42ToBdpTII2LLaBjXt/8O+9Pnb8rxQi4i1gqY60xg2Xd3wtqqeB\nE8AnmXkC+Ityy6G7cwMgIvZRrqKngJcoI4xzdFqPMf53/w2LLa4OqT8D5jJzddn2pYh4sf79APDr\npNq3iV4DZiJiAbhCWbV4DnjUYS2gLIfzMDN/qO8/p4RHj+cGwJvAQmb+npkrwBfASfqtx6pR/V8E\nDq353MF6bCTDYutz7w8gM2cz83BmvgKcBW5m5tvA13RWC4B6a+FhRByth04Bd+nw3KgeAK9GxO76\noPYU5YcQvdUjWD/yHtX/ryhbQjwTEdPAEcpE6NFf7DyLravu/fEdME8ZPiYwS/mnXqNcGfwMnMnM\nPyfVzs0WEa8D72fmTEQ8T6e1iIjjlIf9u4AF4B3KQ95e63GBciHxGLgNvEvZ9qCLekTEZeANYD+w\nBFwAvgSu8x/9j4gPgfOUer2Xmd+O/X7DQpLU4m0oSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbD\nQpLUZFhIkpoMC+kJRcRURPxYV369ExHfRMSzk26XNCTDQhrGEeDjzDwGLAOnJ9weaVCGhTSMe5k5\nX1/fAl6eYFukwRkW0jD+XvN6hbLfhLRjGBbSMMbuXyxtd4aFNAyXb9aO5hLlkqQmRxaSpCbDQpLU\nZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNf0DAdE7KCP28WcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116700790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make one big subplots and put everything in it.\n",
    "f, (ax1,ax2,ax3,ax4,ax5)=plt.subplots(5,1)\n",
    "# Plot the data\n",
    "ax1.stem(range(N),x,linefmt='b-', markerfmt='bo')\n",
    "ax1.plot(range(N),lambdas,'r--')\n",
    "ax1.set_ylabel('Counts')\n",
    "\n",
    "# Gibbs sampler\n",
    "E=110\n",
    "BURN_IN=10\n",
    "\n",
    "# Initialize the chain\n",
    "n=int(round(uniform.rvs()*N))\n",
    "lambda1=gamma.rvs(a,scale=1./b)\n",
    "lambda2=gamma.rvs(a,scale=1./b)\n",
    "\n",
    "# Store the samples\n",
    "chain_n=numpy.array([0.]*(E-BURN_IN))\n",
    "chain_lambda1=numpy.array([0.]*(E-BURN_IN))\n",
    "chain_lambda2=numpy.array([0.]*(E-BURN_IN))\n",
    "\n",
    "for e in range(E):\n",
    "\tprint \"At iteration \"+str(e)\n",
    "\t# sample lambda1 and lambda2 from their posterior conditionals, Equation 8 and Equation 9, respectively.\n",
    "\tlambda1=gamma.rvs(a+sum(x[0:n]), scale=1./(n+b))\n",
    "\tlambda2=gamma.rvs(a+sum(x[n:N]), scale=1./(N-n+b))\n",
    "\t\n",
    "\t# sample n, Equation 10\n",
    "\tf_mult_n=numpy.array([0]*N)\n",
    "\tfor i in range(N):\n",
    "\t\tf_mult_n[i]=sum(x[0:i])*log(lambda1)-i*lambda1+sum(x[i:N])*log(lambda2)-(N-i)*lambda2\n",
    "\tmult_n=exp(f_mult_n-max(f_mult_n))\n",
    "\tn=numpy.where(multinomial(1,mult_n/sum(mult_n),size=1)==1)[1][0]\n",
    "\t\n",
    "\t# store\n",
    "\tif e>=BURN_IN:\n",
    "\t\tchain_n[e-BURN_IN]=n\n",
    "\t\tchain_lambda1[e-BURN_IN]=lambda1\n",
    "\t\tchain_lambda2[e-BURN_IN]=lambda2\n",
    "\n",
    "ax2.plot(chain_lambda1,'b',chain_lambda2,'g')\n",
    "ax2.set_ylabel('$\\lambda$')\n",
    "ax3.hist(chain_lambda1,20)\n",
    "ax3.set_xlabel('$\\lambda_1$')\n",
    "ax3.set_xlim([0,12])\n",
    "ax4.hist(chain_lambda2,20,color='g')\n",
    "ax4.set_xlim([0,12])\n",
    "ax4.set_xlabel('$\\lambda_2$')\n",
    "ax5.hist(chain_n,50)\n",
    "ax5.set_xlabel('n')\n",
    "ax5.set_xlim([1,100])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10f3f0b10>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XnclXP+x/HXJ0kmiWGUKSkSKamQXbetZCtbskwyRaMm\n/fySkWVq8LOELENhLFmjMoOoVHTbm1IouqdlTI2iNINQ0l335/fH90rnTnfd3We5zvJ+Ph7n0Tnf\n65xzfc7luD/nu5u7IyIisl61uAMQEZHsosQgIiLlKDGIiEg5SgwiIlKOEoOIiJSjxCAiIuVUT/cJ\nzGwhsAIoA0rdva2Z7Qw8B+wJLAS6uPuKdMciIiJblokaQxlQ5O6t3b1tVHY1MNnd9wVeBwZmIA4R\nEamETCQG28R5OgGPR/cfBzpnIA4REamETCQGByaZ2XQz6xmV1XX3ZQDuvhTYLQNxiIhIJaS9jwE4\n0t2/MLNfARPNbC4hWSTSuhwiIlki7YnB3b+I/l1uZi8AbYFlZlbX3ZeZWT3gy0291syUMEREqsDd\nraqvTWtTkpn9wsx2iO7XAtoDs4GXgO7R0y4CXqzoPdxdN3cGDRoUewzZctO10LXQtdj8LVnprjHU\nBf4W/fKvDjzt7hPN7H1glJn9FlgEdElzHCIiUklpTQzu/i+g1SbKvwJOSOe5RUSkajTzOUcUFRXF\nHULW0LXYQNdiA12L1LFUtEeli5l5NscnIpKNzAzP1s5nERHJPUoMIiJSjhKDyFZ64AFo3RpGj4ay\nsrijEUk9JQaRrVBaCjffDN26we23Q8uWMGYMqCtM8ok6n0W2wpNPwogR8NprIRmMHw9/+AM0agQP\nPgi//nXcEYqo81kkY9xhyBC46qrw2AxOPhlmzAhNS61awRNPqPYguU+JQaSSJkyAatWgffvy5TVq\nwA03hONDh8KRR8LEiUoQkruUGEQqaX1twSqooLdpE2oPfftCv34hQRQXZzREkZRQH4NIJUybBuec\nAwsWwLbbbvn569bBc8+F/odTTw1JpXbt9McpAupjEEm7sjIYNAiuuKJySQFgm23g/PNh9mz48Uc4\n4ACYPDm9cYqkimoMIltwzTXw5pvhD3vNmlV7j/Hj4bLLoFmzkGQOOyy1MYokUo1BJI0efRRGjYIX\nXqh6UgDo2BHmzoXOneHcc8Pj2bNTF6dIKikxiFTgtddg4EB45RXYddfk32+77aBXL5g3L/Q7HHcc\n/N//hUlzItlETUkim/DWW3DWWWHZi3bt0nOOf/8bLr0Uli8PNZMDD0zPeaTwqClJJMUefRTOPhue\nfjp9SQGgYcPQ99C7N3ToAF27wiefpO98IpWlxCASWbcO+veHW26BN96AE09M/znNoEePMAz2oIPg\n+ONDgliyJP3nFqmIEoNI5MorwwS1v/8d9tsvs+feYQcYMCAkiH33DUtsPP64Zk9LPGLrYzCzk4C7\nCcnpEXe/bRPPUR+DZMRbb4XRQrNnwy67xB0NfPghdO8O9evD/feHRfpEKisn+xjMrBpwH9ABaA6c\nZ2YZ/o0mEqxaBRdfDMOHZ0dSgLAg37RpcPjhoYnp0kth4cK4o5JCEVdTUltgvrsvcvdS4FmgU0yx\nSIG75ho49FDolGXfwBo14LrrwvDWX/1KCUIyJ67EUB/4LOHx4qhMJKPeeitMYLv33rgjqdguu4T5\nDokJolcvWLQo7sgkX6nzWQrWBx+E9YyyqQlpczZOEG3aKEFIelSP6bxLgIYJjxtEZT8zePDgn+4X\nFRVRVFSUzrikQPz1r+GP6vDh2deEtCW77AI33RQW9bvrrpAgzjorrMFUX/XuglRcXExxCtd4j2VU\nkpltA8wFjge+AKYB57l7yUbP06gkSSn3sGfzAw+E9Y8OOijuiJL33//CHXfAI4/AbbeF0UwV7Rkh\nhSHZUUlxD1e9hw3DVW/dxHOUGCRlfvgBevYMTTEvvph/+zPPmhWSQr168NBD0KBB3BFJXHJyuCqA\nu09w933dfZ9NJQWRVFq6FI49NsxufvPN/EsKAC1bhsl5hx8e7g8YAF9+GXdUkovU+Sx5rbQ0rI7a\nti2ccgqMHAnbbx93VOmz7bZw/fWh9rB6dZjBPWAArFgRd2SSS5QYJC+9805oNtp999BRe++94Q9m\nobS9N2gAf/5zSBBffx1qEBMnxh2V5Aotuy1559VX4aKLwn7LZ50VVjEtdJMmhUTZvn3oqK5TJ+6I\nJJ1yto9BJB0+/hh+8xsYMyYM51RSCE48MawDVb06NGkSalHffht3VJKtlBgkbyxdGnZGu+ceOOqo\nuKPJPjvuGOZtvPNOGJm1995w663aQU5+Tk1Jkhd++CFsqnPaaaEvQbZs3rxQq/r8cxgxQjvI5RM1\nJYkQNtdp0CAsOieV07QpvPwy/M//hKamP/0pjGQSUY1Bct6CBXDYYWEPA03qqpolS6Bv3zAP4uqr\n4ZJLoGbNuKOSqlKNQQqaO1x+eRiBpKRQdfXrh/Wjxo6FyZND/8Ojj2oHuUKlGoPktBdfDL9wP/oo\n7F8gqfH++2GRwV/9Cv7yF9hjj7gjkq2hGoMUrFWrQvv4ffcpKaTawQfD1Klw9NFh9dYHH4S1a+OO\nSjJFNQbJSaWlYTezVavguefijia/zZ4d+h+WLAkjvs4/P8yHkOyVs6urVoYSg2zK11/DOeeEWsLI\nkZrFmynFxWHPhyVLoEMHaN483A46CHbYIe7oJJESgxSUuXPDXIVTT4Xbb4dttok7osIzdWq4zZkD\nn3wS5kP06xcGAey4Y9zRCSgxSAGZNAkuvDBstNOjR9zRyHpz54YlNiZMCBPmrrxSfT5xU2KQgnD/\n/XDjjaE/oV27uKORTZk7NySFRYvgscfyY3e8XKXEIHnLHT77LGxXOWVKGGO/995xRyWb4w5PPQX9\n+4fBAT17hoUMq2n8Y0ZpuKrkFXcYOhQOPTR0Kh9+eNhk5r33lBRygVlY3fajj+Bf/wqLGdapEzZK\nuvpq7SiXK1RjkKyxfk/m+fPhzjuhRQvYeee4o5JkffNN6KgeORKefjr0Dw0YALvtFndk+Strawxm\nNsjMFpvZzOh2UsKxgWY238xKzKx9umKQ3LF+T+ayMnjjjTCxSkkhP+y0ExxxxIYd5Vavhv33D0uA\nl5XFHZ1sStpqDGY2CPjO3YduVN4MeAY4BGgATAb22VTVQDWGwrB8eWhq+O1vw+qohbL9ZiErKYHu\n3aFWLXjkEWjcOO6I8kuyNYZ0z1/cVGCdgGfdfS2w0MzmA22Bv6c5FslCZWXQrRuce672USgkzZqF\nDYOGDoVDDglLgK/XvHnYQGiXXeKLr9ClOzH83sx+A7wP9Hf3FUB94L2E5yyJyqQADRkC330XhqJK\nYaleHa66Cs47L4w+W2/UKDjggNDU1KlTfPEVsqQSg5lNAuomFgEOXAsMA25wdzezm4A7gZ5be47B\ngwf/dL+oqIiioqIkIpZs8tZbcPfdMH06bLtt3NFIXPbYo/zqrUccAWedBRdfHDqrTzop1CKaNdPM\n6ooUFxdTXFycsvfLyKgkM9sTGOvuLc3sasDd/bbo2ARgkLv/rClJfQz56z//gdatw6qdJ58cdzSS\njVatCt+PDz4Io5r+8Y8wKGHQoLAxk1Qsaye4mVk9d18a3b8COMTdzzez/YGngUMJTUiTUOdzwenR\nIyy8ds89cUciuWLNmjCj+uabw6imvn2hVSvYfXcNWNhYNieGJ4BWQBmwEOjl7suiYwOBHkAp0M/d\nJ1bwHkoMeWjq1NBUUFKipgHZeusTxMiRoSZRWhoSRfPmG/5N7LjeffdwKyRZmxhSQYkh/6xbF0ah\n9O8PF1wQdzSSD5YvD6u8rl/t9ZNP4NtvNxz/97/DrOsrriic1XiVGCSnDBsWFsIrLlb1XzLjX/8K\nc2RWrw41jf32izui9Mvamc8iG1u+HAYPDiulKilIpjRuDK+9FtZwOuqosB3sF1/EHVV2U2KQtFu6\nNNQUOnQI/3O2aBF3RFJoqlWD3r3h44/D/ebNQ4LQon6bpsQgabF0aagZFBWFqvu774ZhhrfdFndk\nUsjq1QuzrefMCSv5tmwJzz4b7ssG6mOQlProo7DN40cfwSmnhL2ZO3SAmjXjjkzk56ZPD2s27bcf\n3HVXWPAPQid1rVqxhpYU9TFI1li8OCSD884LbbhPPRWWNFBSkGx1yCEwYwbsuy8ceCA0aBBuu+wC\nv/89fP993BHGQ4lBUuK77+DUU0NtoVcvJQPJHTVrhklzX38dhrl++234YfPdd6GpKYUrTeQMNSVJ\n0tauDTWD+vXDEgYacST54uWX4Xe/g+OPD6v/NmkSd0SVo6YkiVVpaaghlJZqGKrkn1NPDRPm9tor\nrM908cXwz3/GHVX6KTFIlX39NXTsGKrdY8ZohVTJT3XqhBF1CxZAo0ZhP/Jbbw015XylxCBVMm9e\n+AXVsiWMHas1jyT/7bRTSBDvvw+TJ4flwefMiTuq9FBikK32+edwzDFhQ/ehQwtn/RkRCLWGSZOg\nZ09o1y4/Z1IrMchWu/LKsGx2z63edkkkP5jBpZfC7Nnhh1GLFiFBLFsWd2SpocQgW2XKlLBX7zXX\nxB2JSPzq1YM77wwd1GahafWpp3J/JrWGq0qllZaGjVFuvBHOPDPuaESyz8yZYSZ148bwwAPx7QOh\n4aqSNh07hv1233svPL733rA37xlnxBuXSLZq0yZ0Th94IBxwAPzxj2H0Xq5RYpBNKikJe+126gRd\nu4b1jm65JSQHzVUQqViNGnDDDTBtWhiosc8+IUHk0vIaakqSTfrf/4XttgvJYM0aGDECqlcPG56I\nSOV9+mkY5vrOO/DII3Dssek/p3Zwk5RbvTo0Gf3972HGp4gkb9y4sErAaaeF5edr107fuWLtYzCz\ns83sYzNbZ2ZtNjo20Mzmm1mJmbVPKG9jZrPMbJ6Z3Z3M+SU9nn8eWrdWUhBJpZNPDsNbV68OzUt3\n3gkrV244XlISZlTffXeYJ/H55/GNbkq2j2E2cAbwRmKhmTUDugDNgI7AMLOfWqaHAz3cvSnQ1Mw6\nJBmDpNhDD4VfNiKSWjvtBI8+ChMnhkEde+8NffuGjuoTTghL18+fH1Z7bdUqdGZ/+GHm46yezIvd\nfS5Awh/99ToBz7r7WmChmc0H2prZIqC2u0+PnvcE0Bl4NZk4pOr++9+wvWGzZuFxSQnMnQunnx5v\nXCL5rGXLsL7YrFnwt7+Foa2HHx62HV3PHZ54Atq3hz59YODA0LGdCUklhs2oD7yX8HhJVLYWWJxQ\nvjgqlxisXBlGG82fDyeeGEZOjBgRVpDUgngi6deyZbhtihlcdFGoSfTqFX681a+/4dgFF8All6Rn\nlOAWE4OZTQLqJhYBDlzr7mNTH1J5gwcP/ul+UVERRUVF6T5lQVi3LnyxWrSAN94I+yh06BDGXOfr\nwmAiuah+/bBQ5fvvww8/hLKVK8P+EGPGwMMPw6efFlOcwh2FUjIqycymAP3dfWb0+GrA3f226PEE\nYBCwCJji7s2i8q5AO3e/rIL31aikNOnfP8zSfPXVDdXTVavCNodHHx1vbCKyZWvXwpAhYa/qO++E\nbt02HMummc+JQbwEdDWzGmbWGGgCTHP3pcAKM2sb9Ut0A15MYQxSCQ8+CK+8EkYfJbZZ/uIXSgoi\nuaJ69bBm2ZQpqR/6mlSNwcw6A38GdgW+AT50947RsYFAD6AU6OfuE6Pyg4ARQE1gnLv328z7q8aQ\nYgsXwsEHhxER++wTdzQikg6a4CZb5Ywz4KCD4Lrr4o5ERNIl2cSQrlFJkoXGj4ePP4aRI+OORESy\nmRbRKxCrV8Pll4dF8GrWjDsaEclmSgwF4o47wtDUjh3jjkREsp2akvJYWRlMnQqjR4cZlDNmxB2R\niOQC1Rjy1F//Cg0bhpmRO+0URiE1ahR3VCKSCzQqKQ+VlMAxx8ALL8CRR8YdjYhkWjZNcJMssGoV\ndOkSNthRUhCRqlCNIc/07BlGID35pLbgFClUmscgQFii9+GH4e23w2JbSgoiUlVKDDlu2jQYNSqs\nslizZhiBtMMOcUclIrlMfQw56ssv4Zxz4NxzYfvtw7K8JSVhJygRkWQoMeSgUaPC5h6NG4e9E268\nMSQENR+JSCqoKSnHDBkS9ox94QU47LC4oxGRfKRRSTlk5cpQS3jrLdh337ijEZFspXkMBeSxx+Co\no5QURCS9VGPIEWvXho11Ro5UE5KIbJ5qDAVi9GjYYw8lBRFJPyWGHOAeOp3/8Ie4IxGRQqDEkAMm\nT4bSUu2lICKZkVRiMLOzzexjM1tnZm0Syvc0s1VmNjO6DUs41sbMZpnZPDO7O5nzF4JVq2DQIBgw\nAKopjYtIBiT7p2Y2cAbwxiaOLXD3NtGtd0L5cKCHuzcFmppZhyRjyFtLloTls5s0gfPPjzsaESkU\nSSUGd5/r7vOBTfV+/6zMzOoBtd19elT0BNA5mRjy1fTpcOihYdmLxx+HbbeNOyIRKRTpnPncyMxm\nAiuA6939baA+sDjhOYujMkkwcyacfHJYLbVTp7ijEZFCs8XEYGaTgLqJRYAD17r72Ape9jnQ0N2/\njvoeXjCz/asS4ODBg3+6X1RURFFRUVXeJmeUlcFll4VRSEoKIlIZxcXFFBcXp+z9UjLBzcymAP3d\nfebmjhMSxhR3bxaVdwXauftlFbyu4Ca4PfxwWAvp7bfV2SwiVZNNE9x+CsLMdjWzatH9vYAmwKfu\nvhRYYWZtzcyAbsCLKYwhp331FVx7Ldx/v5KCiMQn2eGqnc3sM+Aw4GUzGx8dOgaYFfUxjAJ6ufs3\n0bE+wCPAPGC+u09IJoZ8cu21cPbZ0Lp13JGISCHTWklZYsYMOOWUsNnOzjvHHY2I5LJsakqSKior\ng9694eablRREJH5KDFng0UdDn0L37nFHIiKipqTYffUVNGsG48dDmzZbfr6IyJYk25SkxBCzyy6D\nbbaB++6LOxIRyRfJJgbt+RyjGTPC3s1z5sQdiYjIBupjiElpaagtqMNZRLKNEkNMrrsOdt0VLroo\n7khERMpTU1IMXnkl7N08c6ZmOItI9lFiyLDPPoMePWDMmFBjEBHJNvq9mkHffw9du8IVV8BRR8Ud\njYjIpikxpNn338Ozz8JZZ0H9+tC0adimU0QkW2keQxp8/z28/DKMHg2TJ8MRR4Sd2Dp3hl/+Mu7o\nRCTfaYJbllm+HFq2hFatlAxEJB5KDFmmRw+oUweGDo07EhEpVJr5nEXeew8mTAhLZ4uI5Cp1PqfI\nunXQp0/Yq3nHHeOORkSk6pQYUuTBB6F2bTj//LgjERFJjvoYUmD5cmjeHF5/HVq0iDsaESl02sEt\nC9x1V5inoKQgIvkgqcRgZkPMrMTMPjSz581sx4RjA81sfnS8fUJ5GzObZWbzzOzuZM6fDb79Fh56\nSJPWRCR/JFtjmAg0d/dWwHxgIICZ7Q90AZoBHYFhZra+WjMc6OHuTYGmZtYhyRhi9Ze/wAknwF57\nxR2JiEhqJJUY3H2yu5dFD6cCDaL7pwPPuvtad19ISBptzaweUNvdp0fPewLonEwMcVqzJjQjXXVV\n3JGIiKROKvsYfguMi+7XBz5LOLYkKqsPLE4oXxyV5aSRI8N+zdqrWUTyyRYnuJnZJKBuYhHgwLXu\nPjZ6zrVAqbuPTHWAgwcP/ul+UVERRUVFqT5FlZSVwe23hxqDiEiciouLKS4uTtn7JT1c1cy6A5cA\nx7n7j1HZ1YC7+23R4wnAIGARMMXdm0XlXYF27n5ZBe+dtcNVX34Zrr8+bLZjVR4UJiKSerEOVzWz\nk4ABwOnrk0LkJaCrmdUws8ZAE2Cauy8FVphZ26gzuhvwYjIxxOE//wl7Kvzxj0oKIpJ/kl0r6c9A\nDWBSNOhoqrv3dvc5ZjYKmAOUAr0Tfvr3AUYANYFx7j4hyRgy6scf4YwzwryFM86IOxoRkdTTzOet\n4A4XXhhGIz33nPZrFpHspNVVM+hPf4IFC6C4WElBRPKXagyVNHcuHH00zJ4Ndetu+fkiInHRWkkZ\ncscdYVltJQURyXeqMVTCF1+E1VPnzYNdd407GhGRzVONIQPuuSd0OispiEghUI1hC779Fho3hhkz\noFGjWEMREakU1RjS7KGHoEMHJQURKRyqMWzGmjVhOe2xY6F169jCEBHZKqoxpNHDD4dOZyUFESkk\nqjFUYP0+zq+9BgccEEsIIiJVkmyNQYmhAj17Qu3aWlZbRHKPlsRIg6lTYdw4KCmJOxIRkcxTH8NG\n1q0LM5yHDIE6deKORkQk81RjSPDdd2Hpi1q14IIL4o5GRCQeBV9jcIfRo+HMM6FBA5g+PYxG0gY8\nIlKoCrrz+YcfQifznDnQrx906gQ775y204mIZIQ6n6to6VLo3DnMaH73Xdh++7gjEhHJDgXZlDR/\nPrRtC6ecAiNHKimIiCQqyKakrl2hZUu45pqUv7WISOxiXRLDzIaYWYmZfWhmz5vZjlH5nma2ysxm\nRrdhCa9pY2azzGyemd2dzPmr4tNPYfJk6Ns302cWEckNyTYlTQSau3srYD4wMOHYAndvE916J5QP\nB3q4e1OgqZl1SDKGrTJ0KPTqFWY1i4jIzyXV+ezukxMeTgXOSnj8s2qMmdUDarv79KjoCaAz8Goy\ncVTW8uXwzDNhFJKIiGxaKjuffwuMT3jcKGpGmmJmR0Vl9YHFCc9ZHJVlxH33wTnnQL16mTqjiEju\n2WKNwcwmAXUTiwAHrnX3sdFzrgVK3f2Z6DmfAw3d/WszawO8YGb7VyXAwYMH/3S/qKiIoqKiqrwN\nK1fC8OHw9ttVermISNYqLi6muLg4Ze+X9KgkM+sOXAIc5+4/VvCcKUB/QsKY4u7NovKuQDt3v6yC\n16VsVNI998Cbb8Lzz6fk7UREslasE9zM7CRgAHBMYlIws12Br9y9zMz2ApoAn7r7N2a2wszaAtOB\nbsC9ycSwJT/+CDfdFLbofDUjPRkiIrkt2ZnPfwZqAJMsLC40NRqBdAxwg5mtAcqAXu7+TfSaPsAI\noCYwzt0nJBlDhWbMgO7dw/acH34Iu++erjOJiOSPvJ3g9sAD8Mc/ho12zj9fi+KJSOHQWkkbWbsW\nrrgibMn57rvQpEncEYmI5Ja8SgwrVkCXLlCtGrz3njbaERGpirxqSrrwQthuO3jwQaieVylPRKTy\n1JQUeeONMBy1pERJQUQkGXmx7HZpKfz+96GjuVatuKMREclteZEY7rsvDEU988y4IxERyX0538fw\nxRdhb4W334Z9981QYCIiWSzZPoacTwwXXggNG8LNN2coKBGRLFfQnc9vvrmhw1lERFIjZ/sYSkuh\nTx91OIuIpFrOJgZ1OIuIpEdO9jGow1lEpGIF2fl84YWwxx5wyy0xBCUikuUKrvN53jyYPBn++c+4\nIxERyU8518cwbhycdpo6nEVE0iXnEsP48dCxY9xRiIjkr5zqY1i5EurVgyVLYMcdYwxMRCSLJdvH\nkFM1hilT4OCDlRRERNIpqcRgZjeY2Udm9oGZTTCzegnHBprZfDMrMbP2CeVtzGyWmc0zs7u35nxq\nRhIRSb9kawxD3P1Ad28NvAIMAjCz/YEuQDOgIzDM7Kddl4cDPdy9KdDUzDpU5kTuSgwiIpmQVGJw\n9+8THtYCyqL7pwPPuvtad18IzAfaRjWK2u4+PXreE0Dnypxr3jxYswZatEgmYhER2ZKk5zGY2U1A\nN+Ab4NiouD7wXsLTlkRla4HFCeWLo/ItWl9bsCp3p4iISGVsscZgZpOiPoH1t9nRv6cBuPt17t4Q\neBrom65A1YwkIpIZW6wxuPuJlXyvZwj9DIMJNYQ9Eo41iMoqKq/Q4MGDWbMmjEi6/PIioKiS4YiI\nFIbi4mKKi4tT9n5JzWMwsybuviC63xc42t27RJ3PTwOHEpqKJgH7uLub2VTgcmA6IZHc6+4TKnh/\nd3deeiksrz1lSpVDFREpGHGvlXSrmTUldDovAn4H4O5zzGwUMAcoBXonzFTrA4wAagLjKkoK602Z\nApdeCsOGJRmpiIhUStbPfN5tN2fkSDjuuLijERHJDXm/7Pa8ec4++8QdiYhI7sj7xJDN8YmIZKOC\nWitJRETST4lBRETKUWIQEZFylBhERKQcJQYRESlHiUFERMpRYhARkXKUGEREpBwlBhERKUeJQURE\nylFiEBGRcpQYRESkHCUGEREpR4lBRETKUWIQEZFylBhERKQcJQYRESknqcRgZjeY2Udm9oGZTTCz\nelH5nma2ysxmRrdhCa9pY2azzGyemd2d7AcQEZHUSrbGMMTdD3T31sArwKCEYwvcvU10651QPhzo\n4e5NgaZm1iHJGApCcXFx3CFkDV2LDXQtNtC1SJ2kEoO7f5/wsBZQlvD4Z/uNRjWK2u4+PSp6Auic\nTAyFQl/6DXQtNtC12EDXInWqJ/sGZnYT0A34Bjg24VAjM5sJrACud/e3gfrA4oTnLI7KREQkS2yx\nxmBmk6I+gfW32dG/pwG4+3Xu3hB4GugbvewLoKG7twH6A8+Y2Q7p+hAiIpI65u6peSOzPYBx7n7A\nJo5NISSIz4Ep7t4sKu8KtHP3yyp4z9QEJyJSYNz9Z835lZVUU5KZNXH3BdHDzkBJVL4r8JW7l5nZ\nXkAT4FN3/8bMVphZW2A6oQnq3oreP5kPJiIiVZNsH8OtZtaU0Om8CPhdVH4McIOZrYmO9XL3b6Jj\nfYARQE1CDWNCkjGIiEgKpawpSURE8kNWznw2s5PM7B/RJLg/xB1PJplZAzN73cw+iTr6L4/Kdzaz\niWY218xeNbM6cceaKWZWLZoo+VL0uCCvhZnVMbPRZlYSfT8OLeBrcYWZfRwNhHnazGoUyrUws0fM\nbJmZzUooq/Czm9lAM5sffW/aV+YcWZcYzKwacB/QAWgOnGdm+8UbVUatBf7X3ZsDhwN9os9/NTDZ\n3fcFXgcGxhhjpvUD5iQ8LtRrcQ+h+bUZcCDwDwrwWpjZrwkjINu4e0tCk/h5FM61eIzw9zHRJj+7\nme0PdAGaAR2BYWa2xb7brEsMQFtgvrsvcvdS4FmgU8wxZYy7L3X3D6P73xM69BsQrsHj0dMep0Am\nBppZA+DGSnmeAAACYElEQVRk4OGE4oK7Fma2I3C0uz8G4O5r3X0FBXgtItsAtcysOrA9sIQCuRbR\nnLCvNyqu6LOfDjwbfV8WAvMJf2M3KxsTQ33gs4THBTsJzswaAa2AqUBdd18GIXkAu8UXWUbdBQwA\nEjvDCvFaNAb+Y2aPRc1qD5nZLyjAa+HunwN3Av8mJIQV7j6ZArwWCXar4LNv/Pd0CZX4e5qNiUGA\naELgGKBfVHPYeJRA3o8aMLNTgGVRDWpz1d+8vxaE5pI2wP3RxNGVhOaDQvxe7ET4hbwn8GtCzeEC\nCvBabEZSnz0bE8MSoGHC4wZRWcGIqsdjgCfd/cWoeJmZ1Y2O1wO+jCu+DDoSON3MPgVGAseZ2ZPA\n0gK8FouBz9z9/ejx84REUYjfixMI86K+cvd1wN+AIyjMa7FeRZ99CbBHwvMq9fc0GxPDdKBJtHR3\nDaAr8FLMMWXao8Acd78noewloHt0/yLgxY1flG/c/Rp3b+juexG+B6+7+2+AsRTetVgGfBbNGwI4\nHviEAvxeEJqQDjOzmlFH6vGEwQmFdC2M8rXoij77S0DXaNRWY8Jk42lbfPNsnMdgZicRRmBUAx5x\n91tjDiljzOxI4E1gNqE66MA1hP+YowjZfxHQJWHSYN4zs3ZAf3c/3cx+SQFeCzM7kNAJvy3wKXAx\noRO2EK/FIMKPhVLgA6AnUJsCuBZm9gxQBOwCLCNsd/ACMJpNfHYzGwj0IFyrfu4+cYvnyMbEICIi\n8cnGpiQREYmREoOIiJSjxCAiIuUoMYiISDlKDCIiUo4Sg4iIlKPEICIi5SgxiIhIOf8PXDWmvXjF\n2a0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11673e090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(f_mult_n)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}