Ch12.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 12: Markov Chain Monte Carlo\n",
    " \n",
    "This Jupyter notebook is the Python equivalent of the R code in section 12.4 R, pp. 512 - 517, [Introduction to Probability, 2nd Edition](https://www.crcpress.com/Introduction-to-Probability-Second-Edition/Blitzstein-Hwang/p/book/9781138369917), Blitzstein & Hwang.\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metropolis-Hastings\n",
    "\n",
    "Here's how to implement the Metropolis-Hastings algorithm for Example 12.1.8, the Normal-Normal model. First, we choose our observed value of $Y$ and decide on values for the constants $\\sigma$, $\\mu$, and $\\tau$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = 3\n",
    "sigma = 1\n",
    "mu = 0\n",
    "tau = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to choose the standard deviation of the proposals for step 1 of the algorithm, as explained in Example 12.1.8; for this problem, we let $d = 1$. We set the number of iterations to run, and we allocate a NumPy array `theta` of length 10<sup>4</sup> which we will fill with our simulated draws:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 1\n",
    "niter = 10**4\n",
    "theta = np.zeros(niter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for the main loop. We initialize $\\theta$ to the observed value $y$, then run the algorithm described in Example 12.1.8:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta[0] = y\n",
    "\n",
    "np.random.seed(1134903170)\n",
    "\n",
    "from scipy.stats import binom, norm\n",
    "\n",
    "for i in range(1, niter):\n",
    "    theta_p = theta[i-1] + norm.rvs(loc=mu, scale=2, size=1)[0]\n",
    "    numer = norm.pdf(y, loc=theta_p, scale=sigma) * norm.pdf(theta_p, loc=mu, scale=tau)\n",
    "    denom = norm.pdf(y, loc=theta[i-1], scale=sigma) * norm.pdf(theta[i-1], loc=mu, scale=tau)\n",
    "    r = numer / denom\n",
    "    flip = binom.rvs(1, np.min([r, 1]), size=1)[0]\n",
    "    theta[i] = theta_p if flip==1 else theta[i-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's step through each line inside the loop. The proposed value of $\\theta$ is `theta_p`, which equals the previous value of $\\theta$ plus a Normal random variable with mean 0 and standard deviation $d$ (recall that [`scipy.stats.norm.rvs`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html) function takes as parameter `scale` the standard deviation and _not_ the variance). The ratio `r` is\n",
    "\n",
    "\\begin{align}\n",
    "  \\frac{f_{\\theta|Y}(x^{\\prime}|y)}{f_{\\theta|Y}(x|y)} &= \\frac{e^{-\\frac{1}{2 \\, \\sigma^2}(y-x^{\\prime})^2} \\,\\, e^{-\\frac{1}{2 \\, \\tau^2}(x^{\\prime}-\\mu)^2}}{e^{-\\frac{1}{2 \\, \\sigma^2}(y-x)^2} \\,\\, e^{-\\frac{1}{2 \\, \\tau^2}(x-\\mu)^2}}\n",
    "\\end{align}\n",
    "\n",
    "where `theta_p` is playing the role of $x^{\\prime}$ and `theta[i-1]` is playing the role of $x$. The coin flip to determine whether to accept or reject the proposal is `flip`, which is a coin flip with probability `numpy.min([r, 1])` of Heads (encoding Heads as 1 and Tails as 0). Finally, we set `theta[i]` equal to the proposed value if the coin flip lands Heads, or keep it at the previous value otherwise.\n",
    "\n",
    "The array `theta` now contains all of our simulation draws. We typically discard some of the initial draws to give the chain some time to approach the stationary distribution. The following line of code discards the first half of the draws:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = theta[-int(niter/2):]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see what the remaining draws look like, we can create a histogram using [`matplotlib.axes.Axes.hist`](https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.hist.html) function. We can also compute summary statistics such as `numpy.mean(theta)` and `numpy.var(theta)`, which give us the sample mean and sample variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAFRCAYAAACVJc8eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAHuRJREFUeJzt3XuYZHV95/H3JwwIiAjIaBBGxwtL\nNCYCGQ1Zk2gAo4CC+sgGYyK6JOyuJN42q6PrRjZrNvg8iRfW9YKiAeIFRaIoGoMiMeRRYECiXFQQ\nRmcYhMkCIhcdwe/+cU5r0VPdXd1d1dXV5/16nnr61Klz+f5OndOfOpc6lapCkiR1wy+MuwBJkrR0\nDH5JkjrE4JckqUMMfkmSOsTglySpQwx+SZI6xOCXJKlDDH5JkjrE4BcASTYmObzn+dVJnjHGkuZk\njZMhyQFJvpbkh0leMe56BjXq9y7J3yZ587DnNX1a07ftYU9/KSxmHUpycpKTR1TaRDL4h6DdsLYl\n2Xta/yuTVJK1A05jaBvnYlXVL1fVRfMZp18bkrw0ycWLrafftBdS4zAN0t5Ba1xu7/+QvRa4qKoe\nUlWnLsUMh7E8l3L9GuZ6Msy6l9F2N+c6lOSgJP+S5J4klyZ5VJ9h9mw/PDx9Wv+zkpybJIspMsnf\nJbk5yZ1Jvp3kjxYzvVEx+IfnRuBFU0+S/Aqwy7AmnmTVsKalblkG686jgavHXMPAFru8xrW8l8H7\nPEqzrkNJ9gM+C7wFeBhwA/DG6cNV1e3AacCre8b9H8ATgT+oxd/D/q+AtVW1O3A08OYkv7bIaQ6d\nwT88ZwEv6Xl+PHBm7wBJHpnkE0m2Jrlx6pBVkrOARwGfTnJXkte2/TcmeV2SrwN3J1mV5AlJLkpy\nR3vI7eie6W9M8vok1yS5PckHk+zc8/qM407X+0m/reGm9pPyt5IcttCFlGR9ku+007omyfOnvb7d\nvOZYPof3dP9Zkq8n+UGSs6e1/eCeQ4Ufb19/80zzXGj7+rR3+imU+bRvtvd6xvb0zHf6ujPjsm+H\n/2/t8rs7yelJHpHkc+3wX0iy5wxtnK3OC4HfAd7Ztu3fzbCM+q63c62z81yefbe/WZZX7/o1Vx3b\njd+nnQcluaKt9Wxg52njz7q99WvXXHX3eMoMy7eSPL6njt7TD4Nsd3P9P5pxm+yzfPpOa5B1CPgb\n4H1VdV5V3Qt8FHjKDLN6K/CsJI9LcixwIvDcqrpnptoGVVVXV9WPp562j8ctdrpDV1U+FvkANgKH\nA98CngDsAGyi+ZRawFqaD1mXA38O7AQ8luZT6bN6p9FnulcCa2iOHuwIXA+8oZ3GocAPgQN6hr+q\nHX4v4F+AN7evDTLu4X3adEDblke2/dcCj5ttOUzr91Lg4p7nxwKPbJfH7wF3A/u0r804r1mWz+E9\n3Ze2094LuBb4z+1rOwHfBV7ZLocXANuAN8/VPuBdwLsW0d7eGgdu32zv12ztmWndGWDZbwS+CjwC\n2Be4FbgCOAh4EHAh8KY+y2DW9aod5iLgj+bYfrZbb+ea9jyX5yDb3/TltZFmGxikjduNP62NU+/Z\nq9vpvRD4CT/fPqfmNdf6OL1dM9Y91/JtXyvg8T3D/i3br0d9t7sB3p+NzLBNznc9YpZ1CNgd+BHw\nqJ5+LwC+0nafDJw8bZwPAP8AbAUOnmG6nwHumOHxmVnW53cB97TL9gpgt5mGHdfDPf7hmtrrfybw\nTeCmnteeAqyuqr+oqm1VdQPwPuC4OaZ5alVtquZT7CHAbsAp7TQupFk5X9Qz/Dvb4W8D/rLntUHG\n7ed+mn/8T0yyY1VtrKrvzDL8J9tP7HckuYNmI/iZqvp4VW2pqp9W1dnAdcBTFziv6U5tp30b8Gng\nwJ62r2pf/0lVnUvzD2nOeVbVy6vq5Qtt7zTzad9s79ds7Zm+PKbWnbmWPcD/qapbquom4J+BS6rq\na9Xswfw9zYeA+dQ5H/3W27mmPZ/lOcj294DltYA2zjT+1DR2BN7evmfnAJf1GW4h28Bs850y0/+F\nxRhkucy0TS5kWjM5jGbZfr1nO/wQzQetmbwVeBZwUlVd0W+AqnpOVe0xw+M5M024/X/xEOC3gHOB\nH8807LgY/MN1FvD7NHt9Z0577dHAI6eFxBto9rBms6mn+5HApqr6aU+/79LsofUb/rvtOIOOu52q\nuh54Fc2n5luTfLQ9ZPri9rDbXUk+1zPK83o3EOABoZnkJWkuepxaBk8C9p5tXrPVN833e7rvoflH\nMtX2m6qq9/zdpiHNc9b29prnvGZ7v2ZszzQP6Dfbsm/d0tN9b5/nu7G9Ba1Xc9Q6td7OOu15Ls9B\ntr9+y5C56hhg/KlpTH/PtgumBa6Ps8233zC9/xcWY5DlMtM2uZBpzWQtcN607fBLNHv0M9mJJpDP\nHWD681ZV91fVxcB+wH8ZxTwWw+Afoqr6Ls1Ffkey/Qq1Cbhx2qfGh1TVkVOjzzTZnu4twJokve/b\no3jgkYU1017bMo9xZ2rXh6vqN/n5qYu3VNWHqmq39nHEXNMASPJomr2sPwEe1m6gVwE/u5K237ym\nXhpkHjO4Gdg3ecAVuz9bTrPMc+jm0b7Z3q9Z29M7u6mOQZb9Ai14vZqm33o757TnsTzn2v76jTNl\n0DbOto72e8+2u+oc5lwf+81jkG1jpv8L9wC79rz2i/OY9rDe+8VO60E07QAgyWOAdcB5s4zzZOCq\nqrpvpgHSXN9y1wyPz8003jSrWIbn+A3+4TsBOLSq7p7W/1LgzvZCnF2S7JDkSUmmLkC5hea842wu\noTkv+9okO6b5Lu1zaS5kmXJSkv2S7EWzR3P2PMbdTprvzx6a5EE059HupTkcuRAPpvlHsrWd9sto\n9joHmdcgy2cmX2mn8yftxU/H0B7iHnL7ZjXP9s32fs3YnlnMuuwXYUHrVR/91ttZpz3P5TnX9jfq\nNn4FuA94RfuevYA+79kA6+NCt4OZ/i9cCfx+uzyeDTx92nizzW9Y7/1ip3UZ8PT2SOQa4MPAf29P\nL8zkQJq2z6iqjujZuZn+2G5nJ8nDkxyXZLd2eT6L5lTFhQO0YUkZ/ENWVd+pqg19+t9PsyIfSHNU\n4N+A9wMPbQf5K+CN7WHIP5th2ttoviJyRDv+u4CXVNU3ewb7MPCPNBcu3UBzkdSg4/bzIOCUdpzv\nAw+n+ccxb1V1Dc3Vt1+h+YfyKzQXGg0yrzmXzyzz3UZzsc8JNBfm/AHN+cMfzzFPkrwnyXvm1dCZ\nDdy+2d6vOdrT1wDLfkEWsV5Nt916O8C057M859r+RtrGnvfspcDtNBdX9jvMPNf2ttDtoO//BZoL\nRJ9Lsx69GPjktPFmnN8Q3/vFTutCmusHvg1cDJxVVe+bY5wnM0fwL0DRHNbfTPMe/zXwqqr61JDn\ns2h54CknTbIkG2mufP3CuGtZ7pJcArynqj447lqGYZLb43qrUUp7176qOnm8lSwf7vGrE5I8Pckv\ntodZjwd+ldkv/lnWVlp7JC2dlXynJ6nXAcDHaK4q/g7wwqq6ebwlLcpKa480KheNu4DlxkP9kiR1\niIf6JUnqEINfkqQOWZHn+Pfee+9au3btuMuQJGnJXH755f9WVavnGm5FBv/atWvZsGG7r9JLkrRi\nJZnt9wl+xkP9kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKHGPySJHWI\nwS9JUocY/JIkdciKvFe/pO2tXX/+ksxn4ylHLcl8JC2Me/ySJHWIwS9JUocY/JIkdYjBL0lShxj8\nkiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKH+CM90jKwVD+gI0nu8UuS\n1CEGvyRJHWLwS5LUISML/iQfSHJrkqt6+u2V5IIk17V/92z7J8mpSa5P8vUkB/eMc3w7/HVJjh9V\nvZIkdcEo9/j/Fnj2tH7rgS9W1f7AF9vnAEcA+7ePE4F3Q/NBAXgT8OvAU4E3TX1YkCRJ8zey4K+q\nLwO3Tet9DHBG230G8Lye/mdW46vAHkn2AZ4FXFBVt1XV7cAFbP9hQpIkDWipz/E/oqpuBmj/Przt\nvy+wqWe4zW2/mfpLkqQFWC4X96VPv5ql//YTSE5MsiHJhq1btw61OEmSVoqlDv5b2kP4tH9vbftv\nBtb0DLcfsGWW/tupqtOqal1VrVu9evXQC5ckaSVY6uA/D5i6Mv944FM9/V/SXt1/CPCD9lTA54Hf\nTbJne1Hf77b9JEnSAozslr1JPgI8A9g7yWaaq/NPAT6W5ATge8Cx7eCfBY4ErgfuAV4GUFW3Jflf\nwGXtcH9RVdMvGJQkSQMaWfBX1YtmeOmwPsMWcNIM0/kA8IEhliZJUmctl4v7JEnSEjD4JUnqEINf\nkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5Kk\nDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x\n+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfgl\nSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqkLEEf5JXJ7k6yVVJPpJk5ySPSXJJkuuSnJ1kp3bYB7XP\nr29fXzuOmiVJWgmWPPiT7Au8AlhXVU8CdgCOA94CvK2q9gduB05oRzkBuL2qHg+8rR1OkiQtwLgO\n9a8CdkmyCtgVuBk4FDinff0M4Hlt9zHtc9rXD0uSJaxVkqQVY8mDv6puAv4a+B5N4P8AuBy4o6ru\nawfbDOzbdu8LbGrHva8d/mFLWbMkSSvFOA7170mzF/8Y4JHAg4Ej+gxaU6PM8lrvdE9MsiHJhq1b\ntw6rXEmSVpRxHOo/HLixqrZW1U+Ac4F/D+zRHvoH2A/Y0nZvBtYAtK8/FLht+kSr6rSqWldV61av\nXj3qNkiSNJHGEfzfAw5Jsmt7rv4w4BrgS8AL22GOBz7Vdp/XPqd9/cKq2m6PX5IkzW0c5/gvoblI\n7wrgG20NpwGvA16T5Hqac/int6OcDjys7f8aYP1S1yxJ0kqxau5Bhq+q3gS8aVrvG4Cn9hn2R8Cx\nS1GXJEkrnXfukySpQwx+SZI6xOCXJKlDDH5JkjrE4JckqUPGclW/pJVr7frzRz6PjaccNfJ5SCuV\nwS/NYilCTJKWkof6JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfgl\nSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDhko+JM8\nadSFSJKk0Rt0j/89SS5N8vIke4y0IkmSNDIDBX9V/SbwYmANsCHJh5M8c6SVSZKkoRv4HH9VXQe8\nEXgd8HTg1CTfTPKCURUnSZKGa9Bz/L+a5G3AtcChwHOr6glt99tGWJ8kSRqiVQMO907gfcAbqure\nqZ5VtSXJG0dSmSRJGrpBg/9I4N6quh8gyS8AO1fVPVV11siqkyRJQzXoOf4vALv0PN+17SdJkibI\noMG/c1XdNfWk7d51NCVJkqRRGTT4705y8NSTJL8G3DvL8JIkaRka9Bz/q4CPJ9nSPt8H+L3RlCRJ\nkkZloOCvqsuS/BJwABDgm1X1k5FWJkmShm7QPX6ApwBr23EOSkJVnTmSqiRJ0kgMFPxJzgIeB1wJ\n3N/2LsDglyRpggy6x78OeGJV1TBm2v7Qz/uBJ9F8gPiPwLeAs2mOKmwE/kNV3Z4kwDto7iVwD/DS\nqrpiGHVIktQ1g17VfxXwi0Oc7zuAf6iqXwKeTHMr4PXAF6tqf+CL7XOAI4D928eJwLuHWIckSZ0y\n6B7/3sA1SS4FfjzVs6qOnu8Mk+wO/Dbw0nYa24BtSY4BntEOdgZwEc0PAh0DnNkebfhqkj2S7FNV\nN8933pIkdd2gwX/yEOf5WGAr8MEkTwYuB14JPGIqzKvq5iQPb4ffF9jUM/7mtp/BL0nSPA10qL+q\n/onmvPuObfdlwELPs68CDgbeXVUHAXfz88P6/aRfSdsNlJyYZEOSDVu3bl1gaZIkrWyD/izvHwPn\nAO9te+0LfHKB89wMbK6qS9rn59B8ELglyT7t/PYBbu0Zfk3P+PsBW5imqk6rqnVVtW716tULLE2S\npJVt0Iv7TgKeBtwJUFXXAQ+fdYwZVNX3gU1JDmh7HQZcA5wHHN/2Ox74VNt9HvCSNA4BfuD5fUmS\nFmbQc/w/rqptzTfrIMkq+hxun4c/BT6UZCfgBuBlNB9CPpbkBOB7wLHtsJ+l+Srf9TRf53vZIuYr\nSVKnDRr8/5TkDcAuSZ4JvBz49EJnWlVX0twbYLrD+gxbNEccJEnSIg16qH89zZX43wD+E81e+BtH\nVZQkSRqNQX+k56fA+9qHJEmaUIPeq/9G+pzTr6rHDr0iSZI0MvO5V/+UnWkuvNtr+OVIkqRRGvQG\nPv+v53FTVb0dOHTEtUmSpCEb9FD/wT1Pf4HmCMBDRlKRJEkamUEP9f9NT/d9tD+bO/RqJEnSSA16\nVf/vjLoQSZI0eoMe6n/NbK9X1VuHU44kSRql+VzV/xSa++YDPBf4Mg/8uVxJkrTMDRr8ewMHV9UP\nAZKcDHy8qv5oVIVJkqThG/SWvY8CtvU83wasHXo1kiRppAbd4z8LuDTJ39Pcwe/5wJkjq0qSJI3E\noFf1/2WSzwG/1fZ6WVV9bXRlSZKkURj0UD/ArsCdVfUOYHOSx4yoJkmSNCIDBX+SNwGvA17f9toR\n+LtRFSVJkkZj0D3+5wNHA3cDVNUWvGWvJEkTZ9Dg31ZVRfvTvEkePLqSJEnSqAwa/B9L8l5gjyR/\nDHwBeN/oypIkSaMw6FX9f53kmcCdwAHAn1fVBSOtTJIkDd2cwZ9kB+DzVXU4YNhLkjTB5jzUX1X3\nA/ckeegS1CNJkkZo0Dv3/Qj4RpILaK/sB6iqV4ykKkmSNBKDBv/57UOSJE2wWYM/yaOq6ntVdcZS\nFSRJkkZnrnP8n5zqSPKJEdciSZJGbK7gT0/3Y0dZiCRJGr25gr9m6JYkSRNorov7npzkTpo9/13a\nbtrnVVW7j7Q6SZI0VLMGf1XtsFSFSNKg1q4f/ZeMNp5y1MjnIY3DoPfqlyRJK4DBL0lShxj8kiR1\niMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1yNiCP8kOSb6W5DPt88ckuSTJdUnOTrJT\n2/9B7fPr29fXjqtmSZIm3Vz36h+lVwLXAlP3+38L8Laq+miS9wAnAO9u/95eVY9Pclw73O+No2At\nL0tx21ZJWmnGssefZD/gKOD97fMAhwLntIOcATyv7T6mfU77+mHt8JIkaZ7Gdaj/7cBrgZ+2zx8G\n3FFV97XPNwP7tt37ApsA2td/0A4vSZLmacmDP8lzgFur6vLe3n0GrQFe653uiUk2JNmwdevWIVQq\nSdLKM449/qcBRyfZCHyU5hD/24E9kkxdc7AfsKXt3gysAWhffyhw2/SJVtVpVbWuqtatXr16tC2Q\nJGlCLXnwV9Xrq2q/qloLHAdcWFUvBr4EvLAd7HjgU233ee1z2tcvrKrt9vglSdLcltP3+F8HvCbJ\n9TTn8E9v+58OPKzt/xpg/ZjqkyRp4o3z63xU1UXARW33DcBT+wzzI+DYJS1MkqQVajnt8UuSpBEz\n+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfgl\nSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnq\nEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCD\nX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6pAlD/4ka5J8Kcm1Sa5O8sq2/15JLkhyXft3\nz7Z/kpya5PokX09y8FLXLEnSSjGOPf77gP9aVU8ADgFOSvJEYD3wxaraH/hi+xzgCGD/9nEi8O6l\nL1mSpJVh1VLPsKpuBm5uu3+Y5FpgX+AY4BntYGcAFwGva/ufWVUFfDXJHkn2aacjSSOxdv35SzKf\njacctSTzkaaM9Rx/krXAQcAlwCOmwrz9+/B2sH2BTT2jbW77SZKkeRpb8CfZDfgE8KqqunO2Qfv0\nqz7TOzHJhiQbtm7dOqwyJUlaUcYS/El2pAn9D1XVuW3vW5Ls076+D3Br238zsKZn9P2ALdOnWVWn\nVdW6qlq3evXq0RUvSdIEG8dV/QFOB66tqrf2vHQecHzbfTzwqZ7+L2mv7j8E+IHn9yVJWpglv7gP\neBrwh8A3klzZ9nsDcArwsSQnAN8Djm1f+yxwJHA9cA/wsqUtV5KklWMcV/VfTP/z9gCH9Rm+gJNG\nWpQkSR3hnfskSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCD\nX5KkDhnHvfrVAWvXnz/uEiRJfbjHL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lS\nhxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKHrBp3AZLUZWvXnz/y\neWw85aiRz0OTwz1+SZI6xOCXJKlDDH5JkjrE4JckqUMMfkmSOsSr+jtmKa4gliQtX+7xS5LUIQa/\nJEkdYvBLktQhBr8kSR3ixX2StMJ5W2D1co9fkqQOmZjgT/LsJN9Kcn2S9eOuR5KkSTQRwZ9kB+D/\nAkcATwRelOSJ461KkqTJMynn+J8KXF9VNwAk+ShwDHDNUhXgOTJJ0kowKcG/L7Cp5/lm4NfHVMvI\neFc9SZNqqf5/LcUO0krf0ZuU4E+ffvWAAZITgRPbp3cl+dbIq9re3sC/jWG+o2BblqeV1BZYWe2x\nLUsgb5n3KMuyLQtoB8zdlkcPMpFJCf7NwJqe5/sBW3oHqKrTgNOWsqjpkmyoqnXjrGFYbMvytJLa\nAiurPbZlebIt25uIi/uAy4D9kzwmyU7AccB5Y65JkqSJMxF7/FV1X5I/AT4P7AB8oKquHnNZkiRN\nnIkIfoCq+izw2XHXMYexnmoYMtuyPK2ktsDKao9tWZ5syzSpqrmHkiRJK8KknOOXJElDYPAPWZJj\nk1yd5KdJJvJK0pVye+QkH0hya5Krxl3LYiVZk+RLSa5t169XjrumhUqyc5JLk/xr25b/Oe6aFivJ\nDkm+luQz465lMZJsTPKNJFcm2TDuehYjyR5JzknyzXa7+Y1x17QQSQ5o34+px51JXrWoaXqof7iS\nPAH4KfBe4M+qaqI2nvb2yN8GnknzNcrLgBdV1ZLdJXFYkvw2cBdwZlU9adz1LEaSfYB9quqKJA8B\nLgeeN6HvS4AHV9VdSXYELgZeWVVfHXNpC5bkNcA6YPeqes6461moJBuBdVW17L73Pl9JzgD+uare\n334bbNequmPcdS1G+//5JuDXq+q7C52Oe/xDVlXXVtU4bh40LD+7PXJVbQOmbo88carqy8Bt465j\nGKrq5qq6ou3+IXAtzR0tJ0417mqf7tg+JnYPJMl+wFHA+8ddixpJdgd+GzgdoKq2TXrotw4DvrOY\n0AeDX9vrd3vkiQyYlSrJWuAg4JLxVrJw7aHxK4FbgQuqamLbArwdeC3Nkb5JV8A/Jrm8vRvqpHos\nsBX4YHsK5v1JHjzuoobgOOAji52Iwb8ASb6Q5Ko+j4ncM55mztsja3yS7AZ8AnhVVd057noWqqru\nr6oDae7C+dQkE3kqJslzgFur6vJx1zIkT6uqg2l+CfWk9nTZJFoFHAy8u6oOAu4GJvZ6JYD2dMXR\nwMcXO62J+R7/clJVh4+7hhGa8/bIGo/2fPgngA9V1bnjrmcYquqOJBcBzwYm8SLMpwFHJzkS2BnY\nPcnfVdUfjLmuBamqLe3fW5P8Pc2pvy+Pt6oF2Qxs7jmSdA4THvw0H8auqKpbFjsh9/g1nbdHXoba\nC+JOB66tqreOu57FSLI6yR5t9y7A4cA3x1vVwlTV66tqv6paS7OtXDipoZ/kwe2Fo7SHxX+Xyfww\nRlV9H9iU5IC212Es4c+4j8iLGMJhfjD4hy7J85NsBn4DOD/J58dd03xU1X3A1O2RrwU+Nqm3R07y\nEeArwAFJNic5Ydw1LcLTgD8EDu35Ws+R4y5qgfYBvpTk6zQfNC+oqon+GtwK8Qjg4iT/ClwKnF9V\n/zDmmhbjT4EPtevZgcD/HnM9C5ZkV5pvWg3lSJ9f55MkqUPc45ckqUMMfkmSOsTglySpQwx+SZI6\nxOCXJKlDDH5JkjrE4JckqUMMfklD0f7wzjuSXN3+pvtjx12TpO0Z/JKG5fXADVX1y8CpwMvHXI+k\nPvyRHkmL1t7b/flV9WttrxtpfqNe0jJj8EsahsOBNUmubJ/vBXxhjPVImoGH+iUNw4HAn1fVgVV1\nIPCPwJVzjCNpDAx+ScOwJ3APQJJVND/p+umxViSpL4Nf0jB8Gzik7X41zU+63jjGeiTNwJ/llbRo\nSfYEPgfsDXwFOLGq7h1vVZL6MfglSeoQD/VLktQhBr8kSR1i8EuS1CEGvyRJHWLwS5LUIQa/JEkd\nYvBLktQhBr8kSR3y/wFw7OOB5yBZ7wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_, ax = plt.subplots(figsize=(8, 5))\n",
    "\n",
    "ax.hist(theta, bins=16)\n",
    "\n",
    "ax.set_xlabel(r'$\\theta$')\n",
    "ax.set_ylabel(r'Frequency')\n",
    "ax.set_title(r'Metropolis-Hastings: Histogram of posterior distribution of $\\theta | Y=3$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sample mean = 2.3538455429726173\n",
      "sample var = 0.8442196578144929\n"
     ]
    }
   ],
   "source": [
    "sample_mean = np.mean(theta)\n",
    "print('sample mean = {}'.format(sample_mean))\n",
    "\n",
    "sample_var = np.var(theta, ddof=1)\n",
    "print('sample var = {}'.format(sample_var))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gibbs\n",
    "\n",
    "Now let's implement Gibbs sampling for Example 12.2.6, the chicken-egg problem with unknown hatching probability and invisible unhatched eggs. The first step is to decide on our observed value of $X$, as well as the constants $\\lambda$, $a$, $b$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 7\n",
    "# 'lambda' is a reserved keyword in Python!\n",
    "lambd = 10 \n",
    "a = 1 \n",
    "b = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we decide how many iterations to run, and we allocate space for our results, creating two NumPy arrays `p` and `N` of length 10<sup>4</sup> which we will fill with our simulated draws:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "niter = 10**4 \n",
    "p = np.zeros(niter) \n",
    "N = np.zeros(niter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we're ready to run the Gibbs sampler. We initialize `p` and `N` to the values 0.5 and $2x$, respectively, and then we run the algorithm as explained in Example 12.2.6:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "p[0] = 0.5 \n",
    "N[0] = 2*x\n",
    "\n",
    "np.random.seed(1836311903)\n",
    "\n",
    "from scipy.stats import beta, poisson\n",
    " \n",
    "for i in range(1, niter):\n",
    "    p[i] = beta.rvs(x+b, N[i-1]-x+b, size=1)[0]\n",
    "    N[i] = x + poisson.rvs(lambd*(1-p[i-1]), size=1)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we discard the initial draws:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = p[-int(niter/2):]\n",
    "N = N[-int(niter/2):]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see what the remaining draws look like, we can make histograms using `Axes.hist(p)` and `Axes.hist(N)`, which will result in graphs similar to those R-generated ones in Figure 12.5. We can also compute summary statistics such as `numpy.mean(p)` or `numpy.median(p)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0kAAAGgCAYAAAB/rSwjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3X2YJVV57/3vTwbwBeKADASHwdFI\nDEYikNGQmJMoaCIQHcwDCcYocpGQ54iJRhOZ+OQ6auI5Z0xOQkISXzAYwaMBxBeIqBFBjieJoKMg\nImgYcYQBZEZl8AUFwfv5o1ZD0dPds3um99493d/Pde1rV61aVfuuqt29+u61qipVhSRJkiSp87Bx\nByBJkiRJ84lJkiRJkiT1mCRJkiRJUo9JkiRJkiT1mCRJkiRJUo9JkiRJkiT1mCRJkiRJUo9JkiRJ\nkiT1mCRJGookX0zyzDa9IcmzZ6g74/Jh6sepByV5UpKrk3wnyR+MO55BDft8JnlnkjcO47Nm8zOz\nI9te6PrnSJK2l0mSpO2S5MQkVyX5XpJNbfplSQJQVT9dVVeMMb6t/shM8tIk/9YvGyTOcSZxY/Qa\n4Iqq2rOqzhzFB87FcR7l927Qzxp0v+Yq9qk+b9w/j/NFkr2SVJJPTSp/W5IzxhWXpPnHJEnSrCV5\nNfC3wF8CPw7sB/y/wDOA3cYY2k4nyZJxxzCNxwFfHHcQg9rR4zjO8zCPvwML0aHA14EnJ9l/Uvk1\n4wlJ0nxkkiRpVpI8Gvgz4GVVdWFVfac6V1fVi6rqnlZv8n+zn5bk+iR3JvmnJA+ftOkplyc5Pcmt\nbdjXl5McNcf780CcU31WkncBBwL/kuS7SV7T6h6c5IokW9pQpuf3tnl4b6jae5Oc3x/+0z7z9CTX\nAt9LsiTJmiRfaetcn+QFk+r/cZJrW8/d2Un2S/KRVv/jSfbq1R/omE23D0kuB54F/H3b55+c5rj9\nyXTndBvHZzbH+bFJ3pdkc5Kvpjf0b5rj+JDv3Tbi2Gr9KfbzsCSfa7GeDzx80vr9z5rNfm0zdmb4\nmUnXG/LE3vw7k7xxG5838T2f9pj06v5R+77d1b6/s/p5TPKiJP/R1v16kluSHD1V3VZ/W9//KePZ\n1jmawqHAOuBSYOL7vgtwCHD1DOtJWmyqypcvX74GfgHPBe4Dlmyj3gbg2b3p64AVwN7AvwNvnFR3\nq+XAk4BbgMe2eiuBn+it92bgzdv6/F7ZS4F/m6reTJ81eVvArsB64LV0PWdHAt9p29gN+Brwilbv\n14F7p9jfa9r+PqKVnQA8lu6fV78JfA/Yv1f/Sroeu+XAJuBzwGHA7sDlwOta3RmP2SD70JZfAfzO\nNs7vlOd0G8dnNsf5YcBngf/WtvME4CbgV2c4jg9sY4B93Gr9Sfs4cS7/sG3reOCHvf3sf9bA+zVg\n7NMe37a8gCf25t85VVxTfM9nPCa9up+m+z7uDdxA11M80HerLVsLfB/4jfaZfwR8bYbv07a+/1vF\nM8g5muJzzgX+HPht4COt7KeBe4Bdx/V71ZcvX/PvZU+SpNnaB/hGVd03UdD+Y7wlyfeT/NI06/19\nVd1SVd8C/jvwwgGW30+XBDw5ya5VtaGqvjKxQlW9rKpeNkOsH2xxbUmyhS6pms6MnzXJEcAewNqq\nureqLgc+1GI+AlgCnFlVP6yq99P9gTfZmW1/v9/25b1VdVtV/aiqzgduBJ7eq/93VXVHVd0K/F/g\nqup67+4BPkCXMM1mP2bah0FNd05n2vZsjvPTgGVV9WdtOzcBbwdO7NV5yHHcjn3c1vq7An/TzuWF\nwGemiXU2+zXIZ8O2f2a2x6Dn/cz2ffwW8C90PTCz2cdDgDOq6oKq+iFdcnJgtu5BBgb6/k8Vz8T+\nDHqO4MFhdZcA/yXJnq3suhanJAEOt5M0e98E9ukPTaqqX6iqpW3ZdL9XbulNf43uv8IzLq+q9cAr\ngdcDm5Kcl2TyejM5rqqWTryAaROqWX7WY4FbqupHk2Je3pbdWlU1zb5NWZbkJUmu6SV0T6FLSCfc\n0Zv+/hTze8xyP2bah0FNd06n3fYsj/PjgMdOSnRfS9ejNlUMkw2yj9taf/K5/NpUFbfzuzrTZ09e\nPtXPzPYY9Lx/vTd9N7DHLPfxEODC3vy+wHer6gdTVR7g+79VPL39GegcJdkdOBi4pqrupPvnxdF0\n/2DweiRJD2GSJGm2PkU3NGX1LNdb0Zs+ELhtkOVV9Z6q+kW6P5gLeNMsP3dgM3xWTap6G7AiSf93\n6IHArcDtwPKku8tf09+3Bz5uYiLJ4+h6SF4OPKYldNcBmWK9HdmPQfdhUNOd0xm3PYvjfAvw1X6i\nW93d9o7p1Zm8Tt8g+zjT+lOdywOnqzyL/Rrks2Hmn5m7gUf25n98wO3u0Hkf5LuVZGmLfXOv+Hjg\nI1Ntcwe//7M5R0+h+4fCTW3+g8BxdEmS1yNJegiTJEmzUlVbgDcAb05yfJI9kjwsyaHAo2ZY9bQk\nByTZm6434PxtLU/3rJ4j23+Af0D3B879c79XDzwXaLrPuoPuepgJV9FdM/GaJLume/7M84Dz6JLI\n+4GXt4vxV/PQYUNTeRTdH5ybWywn0/1BN9f70TfTPgxqunM67bZneZw/DXy73SzgEUl2SfKUJE8b\nML4d3cdP0V1/9wftXP4605zLWe7XoGb6mbkG+K12TJ4L/HJv2Uyft93HZBbfrUNa+W+143YsXS/u\n66fZ9I58/wc+R3TJ0Od7vU4XAcdgT5KkKZgkSZq1qvoL4FV0z9LZRPdH2duA04H/mGa19wAfo/sv\n7k10N2bY1vLd6S4A/wbdcJt96f5YBCDJW5O8dU52aubP+p/An7ahQH9UVffS3Rnr6Fb/zcBLqupL\nbdmvA6cAW+guEP8QXe/blKrqeuCv6P7gu4Puj8x/H8J+9D9z2n2YxWdNeU63se3ZHOf76f6APxT4\nalvnH4FHDxLcju5j71y+FLiT7oYC75+m+sD7NchnNzP9zLyC7thsAV5E1ysyYdrP28FjMtB3i+77\n+27g5+mO2xvohr5eP9VGd+T7P8tz9JDbfFfV1+huCrEU+Pwgnydp8chDh/FKkuZakquAt1bVP407\nlrmSZAPd3e8+Pu5YNL8keQvwn1Xlw1kl7bTsSZKkOZbkl5P8eBv+cxLwM8BHxx2XNCKH0N2mW5J2\nWj7lW5Lm3pOAC+juwPUV4Piqun28IUkj8xRgNsM2JWnecbidJEmSJPU43E6SJEmSekySJEmSJKnH\nJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmS\nJKnHJEmSJEmSekyStCgl+WKSZ447jnFJ8qQkVyf5TpI/GHc8kqSHsp2yndJ4mSRpwUmyIcmzJ5W9\nNMm/TcxX1U9X1RWz3c4C8hrgiqras6rOHHcwkrSY2E4NZMZ2KsleSSrJpyaVvy3JGSOLUguWSZI0\nBkmWjDmExwFfHHMMkqR5aidopw4Fvg48Ocn+k8qvGWZgWhxMkrQo9f/7luT0JLe2Lv0vJzkqybuA\nA4F/SfLdJK9pdQ9OckWSLW0oxPN72zy8NzTgvUnOT/LGSZ95epJrge8lWZJkTZKvtHWuT/KCSfX/\nOMm1Sb6X5Owk+yX5SKv/8SR7zbCPU8aa5HLgWcDft337ySnWfVGS/2j78PUktyQ5eocPvCRpILZT\nM7dTdMnQOuBSYGK9XYBDgKu387BLDzBJ0qKW5EnAy4GnVdWewK8CG6rqxcDNwPOqao+q+oskuwL/\nAnwM2Bf4feDd6cZN7wZ8AHgnsDfwz8ALtvpAeCFwLLC0qu4DvgL8F+DRwBuA/z3pP2L/D/Ac4CeB\n5wEfAV4L7EP38zvlOO2ZYq2qI4H/C7y87dt/TrGJQ4DDgPcBK4C/Bd467YGUJA2F7dS07dRhdD1G\nHwSOa2U/BewC3DDVZ0qzYZKkheqD7T9TW5JsAd48Tb37gd3puut3raoNVfWVaeoeAewBrK2qe6vq\ncuBDdA3KEcAS4Myq+mFVvR/49BTbOLOqbqmq7wNU1Xur6raq+lFVnQ/cCDy9V//vquqOqrqVrsG4\nqqqurqp76Bq7w7Yj1kEcApxRVRdU1Q+Bc4EDkzx8wPUlSTOzndqxdmpiWN0lwH9Jsmcru661W9IO\nMUnSQnVcVS2deAEvm6pSVa0HXgm8HtiU5Lwkj51mm48FbqmqH/XKvgYsb8turarqLbtlim08pCzJ\nS5Jc02skn0L337cJd/Smvz/F/B7bEesgDgEu7M3vC3y3qn4w4PqSpJnZTm1nO5Vkd+Bg4JqqupMu\n2TuaB3uXpB1mkqRFr6reU1W/SHeRaAFvmlg0qeptwIok/Z+bA4FbgduB5UnSW7Ziqo+bmEjyOODt\ndMMoHtMayeuATLHebM0U64ySLKWLfXOv+Hi6IRSSpBGzndrKU+gSsJva/MSQu8PweiTNEZMkLWpt\nnPaR7b9SP6D7pXt/W3wH8IRe9auA7wGvSbJruudXPA84D/hUW+/l7ULX1Tx0OMJUHkXXGG1usZxM\n94t/LswU67YcQrcvv9X25Vi6/3C+fo5ikyQNyHZqSocBn+/1il0EHIM9SZpDJkla7HYH1gLfoLuV\n6L50F5wC/E/gT9sQgz+qqnvp7qBzdKv/ZuAlVfWltuzXgVOALcBv042tvme6D66q64G/omu47qBL\nTv59LnZqplgHWP0Q4N3AzwN30l2oe1yLV5I0WrZTW3vIbb6r6mvABmAp8Pm5iE/KQ4emSporSa4C\n3lpV/zTuWGYjyVuA/6wqH8YnSQvYztpOSaNgT5I0R5L8cpIfb8MYTgJ+BvjouOPaDofg7VMlacFZ\nQO2UNHTjfpqytJA8CbiA7k4+XwGOr6rbxxvSdnkKMMhwB0nSzmWhtFPS0DncTpIkSZJ6HG4nSZIk\nST0mSZIkSZLUs1Nfk7TPPvvUypUrxx2GJC1qn/3sZ79RVcvGHcd8ZDslSeO3Pe3UTp0krVy5knXr\n1o07DEla1JJ8bdwxzFe2U5I0ftvTTg11uF2SDUm+kOSaJOta2d5JLk1yY3vfq5UnyZlJ1ie5Nsnh\nw4xNkiRJkqYyimuSnlVVh1bVqja/Brisqg4CLmvz0D1x+aD2OhV4ywhikyRJkqSHGMeNG1YD57Tp\nc4DjeuXnVudKYGmS/ccQnyRJkqRFbNhJUgEfS/LZJKe2sv0mHlzW3vdt5cuBW3rrbmxlD5Hk1CTr\nkqzbvHnzEEOXJEmStBgN+8YNz6iq25LsC1ya5Esz1M0UZVs96baqzgLOAli1apVPwpUkSZI0p4ba\nk1RVt7X3TcAHgKcDd0wMo2vvm1r1jcCK3uoHALcNMz5JkiRJmmxoSVKSRyXZc2Ia+BXgOuBi4KRW\n7STgojZ9MfCSdpe7I4C7JoblSZIkSdKoDHO43X7AB5JMfM57quqjST4DXJDkFOBm4IRW/8PAMcB6\n4G7g5CHGJkmSJElTGlqSVFU3AU+dovybwFFTlBdw2rDikSRJkqRBjOMW4JIkSZI0b5kkSZIkSVKP\nSZIkSZIk9ZgkSZIkSVLPsB8mK0mztnLNJXOynQ1rj52T7Uha+Gb7e8ffL9LCZk+SJEmSJPWYJEmS\nJElSj0mSJEmSJPWYJEmSJElSj0mSJEmSJPWYJEmSJElSj0mSJEmSJPWYJEmSJElSj0mSJEmSJPWY\nJEmSJElSz5JxByBJw7JyzSVzsp0Na4+dk+1IkqSdgz1JkiRJktRjT5KkOTNXPTfSqCT5Q+B3gAK+\nAJwM7A+cB+wNfA54cVXdm2R34FzgZ4FvAr9ZVRvGEbckabjsSZIkLUpJlgN/AKyqqqcAuwAnAm8C\nzqiqg4A7gVPaKqcAd1bVE4EzWj1J0gJkkiRJWsyWAI9IsgR4JHA7cCRwYVt+DnBcm17d5mnLj0qS\nEcYqSRoRkyRJ0qJUVbcC/wu4mS45ugv4LLClqu5r1TYCy9v0cuCWtu59rf5jJm83yalJ1iVZt3nz\n5uHuhCRpKLwmSZK0KCXZi6536PHAFuC9wNFTVK2JVWZY9mBB1VnAWQCrVq3aarlGw2skJe0Ie5Ik\nSYvVs4GvVtXmqvoh8H7gF4ClbfgdwAHAbW16I7ACoC1/NPCt0YYsSRoFkyRJ0mJ1M3BEkke2a4uO\nAq4HPgEc3+qcBFzUpi9u87Tll1eVPUWStACZJEmSFqWquoruBgyfo7v998PohsmdDrwqyXq6a47O\nbqucDTymlb8KWDPyoCVJI+E1SZKkRauqXge8blLxTcDTp6j7A+CEUcQlSRove5IkSZIkqcckSZIk\nSZJ6TJIkSZIkqcdrkiRpG+bqeSsb1h47J9uRJEnDZU+SJEmSJPWYJEmSJElSj8PtJM3ZcDJJkqSF\nwJ4kSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIk\nSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSepZMu4AJG2/lWsuGXcIkiRJC449SZIkSZLUY5IkSZIk\nST0mSZIkSZLU4zVJkiRp3vMaTEmjZE+SJEmSJPWYJEmSJElSj0mSJEmSJPWYJEmSJElSj0mSJEmS\nJPWYJEmSFqUkT0pyTe/17SSvTLJ3kkuT3Nje92r1k+TMJOuTXJvk8HHvgyRpOIZ+C/AkuwDrgFur\n6teSPB44D9gb+Bzw4qq6N8nuwLnAzwLfBH6zqjYMOz5J0uJUVV8GDoUH2qpbgQ8Aa4DLqmptkjVt\n/nTgaOCg9vo54C3tXYvQbG5JvmHtsUOMRNIwjOI5Sa8AbgB+rM2/CTijqs5L8lbgFLqG5hTgzqp6\nYpITW73fHEF8kjQSc/WcF//gGoqjgK9U1deSrAae2crPAa6gS5JWA+dWVQFXJlmaZP+qun0cAUuS\nhmeow+2SHAAcC/xjmw9wJHBhq3IOcFybXt3macuPavUlSRq2E4F/btP7TSQ+7X3fVr4cuKW3zsZW\n9hBJTk2yLsm6zZs3DzFkSdKwDPuapL8BXgP8qM0/BthSVfe1+X4D80Dj05bf1eo/hI2PJGkuJdkN\neD7w3m1VnaKstiqoOquqVlXVqmXLls1FiJKkERtakpTk14BNVfXZfvEUVWuAZQ8W2PhIkubW0cDn\nquqONn9Hkv0B2vumVr4RWNFb7wDgtpFFKUkamWH2JD0DeH6SDXQ3ajiSrmdpaZKJa6H6DcwDjU9b\n/mjgW0OMT5IkgBfy4FA7gIuBk9r0ScBFvfKXtLvcHQHc5fVIkrQwDS1Jqqo/qaoDqmol3Vjvy6vq\nRcAngONbtcmNz0SjdHyrv1VPkiRJcyXJI4HnAO/vFa8FnpPkxrZsbSv/MHATsB54O/CyEYYqSRqh\nUdzdbrLTgfOSvBG4Gji7lZ8NvCvJeroepBPHEJskaRGpqruZdP1rVX2T7m53k+sWcNqIQpMkjdFI\nkqSquoLuFqpU1U3A06eo8wPghFHEI0mSJEnTGfbd7SRJkiRpp2KSJEmSJEk9JkmSJEmS1GOSJEmS\nJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1GOS\nJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1LNk3AFIkmZn5ZpL5mQ7\nG9YeOyfbkSRpobEnSZIkSZJ6TJIkSZIkqcckSZIkSZJ6TJIkSZIkqcckSZIkSZJ6TJIkSZIkqcck\nSZIkSZJ6TJIkSZIkqcckSZK0aCVZmuTCJF9KckOSn0+yd5JLk9zY3vdqdZPkzCTrk1yb5PBxxy9J\nGg6TJEnSYva3wEer6qeApwI3AGuAy6rqIOCyNg9wNHBQe50KvGX04UqSRmHJuAOQFqOVay4ZdwjS\nopfkx4BfAl4KUFX3AvcmWQ08s1U7B7gCOB1YDZxbVQVc2Xqh9q+q20ccuiRpyOxJkiQtVk8ANgP/\nlOTqJP+Y5FHAfhOJT3vft9VfDtzSW39jK5MkLTAmSZKkxWoJcDjwlqo6DPgeDw6tm0qmKKutKiWn\nJlmXZN3mzZvnJlJJ0kg53E6aBYfJSQvKRmBjVV3V5i+kS5LumBhGl2R/YFOv/ore+gcAt03eaFWd\nBZwFsGrVqq2SKEnS/GdPkiRpUaqqrwO3JHlSKzoKuB64GDiplZ0EXNSmLwZe0u5ydwRwl9cjSdLC\nZE+SJGkx+33g3Ul2A24CTqb7B+IFSU4BbgZOaHU/DBwDrAfubnWlbZrtKIQNa48dUiSSBmWSJEla\ntKrqGmDVFIuOmqJuAacNPShJ0tg53E6SJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmS\nekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmS\nJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKln\naElSkocn+XSSzyf5YpI3tPLHJ7kqyY1Jzk+yWyvfvc2vb8tXDis2SZIkSZrOMHuS7gGOrKqnAocC\nz01yBPAm4IyqOgi4Ezil1T8FuLOqngic0epJkiRJ0kgNLUmqznfb7K7tVcCRwIWt/BzguDa9us3T\nlh+VJMOKT5IkSZKmsmSYG0+yC/BZ4InAPwBfAbZU1X2tykZgeZteDtwCUFX3JbkLeAzwjUnbPBU4\nFeDAAw8cZviStKCtXHPJuEOQJGleGuqNG6rq/qo6FDgAeDpw8FTV2vtUvUa1VUHVWVW1qqpWLVu2\nbO6ClSRJkiRGdHe7qtoCXAEcASxNMtGDdQBwW5veCKwAaMsfDXxrFPFJkiRJ0oRh3t1uWZKlbfoR\nwLOBG4BPAMe3aicBF7Xpi9s8bfnlVbVVT5IkSZIkDdMwr0naHzinXZf0MOCCqvpQkuuB85K8Ebga\nOLvVPxt4V5L1dD1IJw4xNkmSJEma0tCSpKq6FjhsivKb6K5Pmlz+A+CEYcUjSdJkSTYA3wHuB+6r\nqlVJ9gbOB1YCG4DfqKo72x1X/xY4BrgbeGlVfW4ccUuShmug4XZJnjLsQCRJ2l472E49q6oOrapV\nbX4NcFl7nt9lbR7gaOCg9joVeMsOfKYkaR4b9Jqktyb5dJKXTVxnJEnSPDKX7VT/uX2Tn+d3bnsO\n4JV0NyLafwc/S5I0Dw2UJFXVLwIvorv73Lok70nynKFGJknSgHagnSrgY0k+257DB7BfVd3etns7\nsG8rf+B5fk3/WX+SpAVk4GuSqurGJH8KrAPOBA5r47NfW1XvH1aAkiQNYjvbqWdU1W1J9gUuTfKl\nGT5ioOf5+dBzSdr5DXpN0s8kOYPuFt5HAs+rqoPb9BlDjE+SpG3a3naqqm5r75uAD9DdWOiOiWF0\n7X1Tq/7A8/ya/rP++tv0oeeStJMb9Jqkvwc+Bzy1qk6buJtPa1z+dFjBSZI0oFm3U0kelWTPiWng\nV4DreOhz+yY/z+8l6RwB3DUxLE+StLAMOtzuGOD7VXU/QJKHAQ+vqrur6l1Di06SpMFsTzu1H/CB\nbkQeS4D3VNVHk3wGuCDJKcDNPPh4ig+3z1lPdwvwk4e2N5KksRo0Sfo48Gzgu23+kcDHgF8YRlCS\nJM3SrNup9ty+p05R/k3gqCnKCzhtLoKVJM1vgw63e3hVTTQ8tOlHDickSZJmzXZKkjRnBk2Svpfk\n8ImZJD8LfH84IUmSNGu2U5KkOTPocLtXAu9NMnEXn/2B3xxOSJIkzZrtlCRpzgyUJFXVZ5L8FPAk\nuudEfKmqfjjUyCRJGpDtlCRpLg38MFngacDKts5hSaiqc4cSlSRJs2c7JUmaEwMlSUneBfwEcA1w\nfysuwMZHkjR2tlOSpLk0aE/SKuDJ7fankiTNN7ZTkqQ5M+jd7a4DfnyYgUiStANspyRJc2bQnqR9\ngOuTfBq4Z6Kwqp4/lKgkSZod2ylJ0pwZNEl6/TCDkCRpB71+3AFIkhaOQW8B/n+SPA44qKo+nuSR\nwC7DDU2SpMHYTkmS5tJA1yQl+V3gQuBtrWg58MFhBSVJ0mzYTkmS5tKgN244DXgG8G2AqroR2HdY\nQUmSNEu2U5KkOTPoNUn3VNW9SQBIsoTu+ROSJM0HtlNaMFauuWRW9TesPXZIkUiL16A9Sf8nyWuB\nRyR5DvBe4F+GF5YkSbNiOyVJmjODJklrgM3AF4DfAz4M/OmwgpIkaZZspyRJc2bQu9v9CHh7e0mS\nNK/YTkmS5tJASVKSrzLF2O6qesKcRyRJ0izZTkmS5tKgN25Y1Zt+OHACsPfchyNJ0naxnZIkzZmB\nrkmqqm/2XrdW1d8ARw45NkmSBmI7JUmaS4MOtzu8N/swuv/Y7TmUiCRJmiXbKUnSXBp0uN1f9abv\nAzYAvzHn0UiStH1spyRJc2bQu9s9a9iBSJK0vWynJElzadDhdq+aaXlV/fXchCNJ0uzZTkmS5tJs\n7m73NODiNv884JPALcMISpKkWbKdkiTNmUGTpH2Aw6vqOwBJXg+8t6p+Z1iBSZI0C7ZTkqQ5M9At\nwIEDgXt78/cCK+c8GkmSto/tlCRpzgzak/Qu4NNJPkD3RPMXAOcOLSpJkmbHdkqSNGcGfZjsfwdO\nBu4EtgAnV9X/GGZgkiQNakfaqSS7JLk6yYfa/OOTXJXkxiTnJ9mtle/e5te35SuHszeSpHEbdLgd\nwCOBb1fV3wIbkzx+SDFJkrQ9tredegVwQ2/+TcAZVXUQXdJ1Sis/Bbizqp4InNHqSZIWoIGSpCSv\nA04H/qQV7Qr872EFJUnSbGxvO5XkAOBY4B/bfIAjgQtblXOA49r06jZPW35Uqy9JWmAG7Ul6AfB8\n4HsAVXUbsOewgpIkaZa2t536G+A1wI/a/GOALVV1X5vfCCxv08tptxRvy+9q9R8iyalJ1iVZt3nz\n5u3bG0nSWA2aJN1bVUV3MSxJHjW8kCRJmrVZt1NJfg3YVFWf7RdPUbUGWPZgQdVZVbWqqlYtW7Zs\n25FLkuadQZOkC5K8DVia5HeBjwNvH15YkiTNyva0U88Anp9kA3Ae3TC7v2nbmLj76wHAbW16I7AC\noC1/NPCtudwJSdL8MOjd7f78OkI1AAATpklEQVQX3fjr9wFPAv5bVf3dMAOTJGlQ29NOVdWfVNUB\nVbUSOBG4vKpeBHwCOL5VOwm4qE1f3OZpyy9vvVeSpAVmm89JSrIL8K9V9Wzg0uGHJEnS4IbQTp0O\nnJfkjcDVwNmt/GzgXUnW0/UgnTgHnyVJmoe2mSRV1f1J7k7y6Kq6axRBSZI0qLlop6rqCuCKNn0T\n8PQp6vwAOGEHQpUk7SS2mSQ1PwC+kORS2p2DAKrqD4YSlSRJs2M7JUmaM4MmSZe0lyRJ85HtlCRp\nzsyYJCU5sKpurqpzZqonSdI42E5JkoZhW3e3++DERJL3DTkWSZJmy3ZKkjTntpUk9R+c94RhBiJJ\n0nawnZIkzbltJUk1zbQkSfOB7ZQkac5t68YNT03ybbr/1D2iTdPmq6p+bKjRSZI0M9spSdKcmzFJ\nqqpdRhWIJEmzZTslSRqGbQ23kyRJkqRFxSRJkiRJknqGliQlWZHkE0luSPLFJK9o5XsnuTTJje19\nr1aeJGcmWZ/k2iSHDys2SZIkSZrOMHuS7gNeXVUHA0cApyV5MrAGuKyqDgIua/MARwMHtdepwFuG\nGJskSZIkTWloSVJV3V5Vn2vT3wFuAJYDq4GJJ6OfAxzXplcD51bnSmBpkv2HFZ8kSZIkTWUk1yQl\nWQkcBlwF7FdVt0OXSAH7tmrLgVt6q21sZZO3dWqSdUnWbd68eZhhS5IkSVqEhp4kJdkDeB/wyqr6\n9kxVpyjb6sGAVXVWVa2qqlXLli2bqzAlSZIkCRhykpRkV7oE6d1V9f5WfMfEMLr2vqmVbwRW9FY/\nALhtmPFJkiRJ0mQzPkx2RyQJcDZwQ1X9dW/RxcBJwNr2flGv/OVJzgN+DrhrYlietKNWrrlk3CFI\nkiRpJzG0JAl4BvBi4AtJrmllr6VLji5IcgpwM3BCW/Zh4BhgPXA3cPIQY5MkSZKkKQ0tSaqqf2Pq\n64wAjpqifgGnDSseSZIkSRrESO5uJ0mSJEk7C5MkSZIkSeoZ5jVJkiRJGrLZ3pxow9pjhxSJtHDY\nkyRJkiRJPSZJkiRJktRjkiRJkiRJPSZJkiRJktRjkiRJkiRJPSZJkiRJktRjkiRJWpSSPDzJp5N8\nPskXk7yhlT8+yVVJbkxyfpLdWvnubX59W75ynPFLkobHJEmStFjdAxxZVU8FDgWem+QI4E3AGVV1\nEHAncEqrfwpwZ1U9ETij1ZMkLUAmSZKkRak6322zu7ZXAUcCF7byc4Dj2vTqNk9bflSSjChcSdII\nmSRJkhatJLskuQbYBFwKfAXYUlX3tSobgeVtejlwC0BbfhfwmCm2eWqSdUnWbd68edi7IEkaApMk\nSdKiVVX3V9WhwAHA04GDp6rW3qfqNaqtCqrOqqpVVbVq2bJlcxesJGlkTJIkSYteVW0BrgCOAJYm\nWdIWHQDc1qY3AisA2vJHA98abaSSpFEwSZIkLUpJliVZ2qYfATwbuAH4BHB8q3YScFGbvrjN05Zf\nXlVb9SRJknZ+S7ZdRZKkBWl/4Jwku9D90/CCqvpQkuuB85K8EbgaOLvVPxt4V5L1dD1IJ44jaEnS\n8JkkSZIWpaq6FjhsivKb6K5Pmlz+A+CEEYQmSRozh9tJkiRJUo9JkiRJkiT1mCRJkiRJUo/XJGle\nW7nmknGHIEmSpEXGniRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQe\nkyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJ\nkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqSeJeMOQJIkSaOzcs0ls6q/Ye2xQ4pEmr/s\nSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSdKilGRFkk8kuSHJF5O8opXv\nneTSJDe2971aeZKcmWR9kmuTHD7ePZAkDYtJkiRpsboPeHVVHQwcAZyW5MnAGuCyqjoIuKzNAxwN\nHNRepwJvGX3IkqRRMEmSJC1KVXV7VX2uTX8HuAFYDqwGzmnVzgGOa9OrgXOrcyWwNMn+Iw5bkjQC\nJkmSpEUvyUrgMOAqYL+quh26RArYt1VbDtzSW21jK5u8rVOTrEuybvPmzcMMW5I0JCZJkqRFLcke\nwPuAV1bVt2eqOkVZbVVQdVZVraqqVcuWLZurMCVJI2SSJElatJLsSpcgvbuq3t+K75gYRtfeN7Xy\njcCK3uoHALeNKlZJ0uiYJEmSFqUkAc4Gbqiqv+4tuhg4qU2fBFzUK39Ju8vdEcBdE8PyJEkLy5Jx\nByBJ0pg8A3gx8IUk17Sy1wJrgQuSnALcDJzQln0YOAZYD9wNnDzacCVJo2KSJElalKrq35j6OiOA\no6aoX8BpQw1KkjQvDG24XZJ3JNmU5LpemQ/okyRJkjSvDbMn6Z3A3wPn9somHtC3NsmaNn86D31A\n38/RPaDv54YYm4Zs5ZpLxh2CJEmStF2G1pNUVZ8EvjWp2Af0SZIkSZrXRn13ux16QJ8kSZIkDdt8\nuQX4QA/oA59kLkmSJGm4Rp0k7fAD+nySuSRJkqRhGnWS5AP6JEmSJM1rQ7u7XZJ/Bp4J7JNkI/A6\nfECfJEmSpHluaElSVb1wmkU+oE+SJEnSvDVfbtwgSZIkSfOCSZIkSZIk9QxtuJ0kSZJ2fivXXDKr\n+hvWHjukSKTRsSdJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJ\nkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSp\nxyRJkiRJknpMkiRJkiSpxyRJkrQoJXlHkk1JruuV7Z3k0iQ3tve9WnmSnJlkfZJrkxw+vsglScNm\nkiRJWqzeCTx3Utka4LKqOgi4rM0DHA0c1F6nAm8ZUYySpDEwSZIkLUpV9UngW5OKVwPntOlzgON6\n5edW50pgaZL9RxOpJGnUTJIkSXrQflV1O0B737eVLwdu6dXb2Mq2kuTUJOuSrNu8efNQg5UkDYdJ\nkiRJ25YpymqqilV1VlWtqqpVy5YtG3JYkqRhMEmSJOlBd0wMo2vvm1r5RmBFr94BwG0jjk2SNCIm\nSZIkPehi4KQ2fRJwUa/8Je0ud0cAd00My5MkLTxLxh2AJEnjkOSfgWcC+yTZCLwOWAtckOQU4Gbg\nhFb9w8AxwHrgbuDkkQcsSRoZkyRJ0qJUVS+cZtFRU9Qt4LThRiRJmi8cbidJkiRJPSZJkiRJktRj\nkiRJkiRJPV6TJEmSpDmzcs0lA9fdsPbYIUYibT+TJD3EbH6xSZIkSQuRw+0kSZIkqcckSZIkSZJ6\nTJIkSZIkqcckSZIkSZJ6TJIkSZIkqcckSZIkSZJ6vAW4JEmSxmK2jx7xuUoaFXuSJEmSJKnHJEmS\nJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnH5yQtELN9zoAkSdLOxucq\naVTsSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHp+T\nJEmSpAXJ5yppe82rnqQkz03y5STrk6wZdzySJE1mWyVJC9+86UlKsgvwD8BzgI3AZ5JcXFXXjzcy\nSZI6tlXSwmbPkybMmyQJeDqwvqpuAkhyHrAasOGRJM0XtlWSHjDbpGo2TMDGaz4lScuBW3rzG4Gf\nm1wpyanAqW32niTXjSC2ndE+wDfGHcQ85bGZnsdmeh6b6T1p3AGM0DbbqkXcTi22n5HFtL/u6xjk\nTUP/iHmzryMw63ZqPiVJmaKstiqoOgs4CyDJuqpaNezAdkYem+l5bKbnsZmex2Z6SdaNO4YR2mZb\ntVjbqcW0r7C49td9XZgW277Odp35dOOGjcCK3vwBwG1jikWSpKnYVknSIjCfkqTPAAcleXyS3YAT\ngYvHHJMkSX22VZK0CMyb4XZVdV+SlwP/CuwCvKOqvriN1c4afmQ7LY/N9Dw20/PYTM9jM71Fc2y2\no61aNMeGxbWvsLj2131dmNzXGaRqq8t+JEmSJGnRmk/D7SRJkiRp7EySJEmSJKlnp0iSkjw3yZeT\nrE+yZorluyc5vy2/KsnK0Uc5HgMcm1cluT7JtUkuS/K4ccQ5Dts6Nr16xyepJIviNpgw2LFJ8hvt\nu/PFJO8ZdYzjMsDP1IFJPpHk6vZzdcw44hy1JO9Ismm6Z/6kc2Y7btcmOXzUMc43g/4OWgiSbEjy\nhSTXLLRbwk/13U+yd5JLk9zY3vcaZ4xzZZp9fX2SW9u5vWah/M5LsqL9Lr+htXOvaOUL7tzOsK8L\n9dw+PMmnk3y+7e8bWvnjW55wY8sbdptxQ1U1r190F8Z+BXgCsBvweeDJk+q8DHhrmz4ROH/ccc+j\nY/Ms4JFt+r96bLaqtyfwSeBKYNW4454vxwY4CLga2KvN7zvuuOfRsTkL+K9t+snAhnHHPaJj80vA\n4cB10yw/BvgI3XOEjgCuGnfMYz5eA/0OWigvYAOwz7jjGNK+bfXdB/4CWNOm1wBvGnecQ9zX1wN/\nNO7YhrCv+wOHt+k9gf9sv9MX3LmdYV8X6rkNsEeb3hW4qrVLFwAntvK3TrTl0712hp6kpwPrq+qm\nqroXOA9YPanOauCcNn0hcFSSqR74t9Bs89hU1Seq6u42eyXdMz0Wg0G+NwB/TvcL8QejDG7MBjk2\nvwv8Q1XdCVBVm0Yc47gMcmwK+LE2/WgWyTNyquqTwLdmqLIaOLc6VwJLk+w/mujmpUF/B2mem+a7\n3/+74xzguJEGNSQD/JwvGFV1e1V9rk1/B7gBWM4CPLcz7OuC1Nqh77bZXdurgCPp8gQY4NzuDEnS\ncuCW3vxGtj6xD9SpqvuAu4DHjCS68Rrk2PSdQvef3sVgm8cmyWHAiqr60CgDmwcG+d78JPCTSf49\nyZVJnjuy6MZrkGPzeuC3k2wEPgz8/mhCm/dm+/tooVtsx6OAjyX5bJJTxx3MCOxXVbdD9wcosO+Y\n4xm2l7dhtO9YCMPPJmuXaRxG1+OwoM/tpH2FBXpuk+yS5BpgE3ApXc/+lpYnwAC/k3eGJGmqHqHJ\n9y0fpM5CNPB+J/ltYBXwl0ONaP6Y8dgkeRhwBvDqkUU0fwzyvVlCN+TumcALgX9MsnTIcc0Hgxyb\nFwLvrKoD6IaYvat9nxa7xfp7eDqL7Xg8o6oOB44GTkvyS+MOSHPmLcBPAIcCtwN/Nd5w5laSPYD3\nAa+sqm+PO55hmmJfF+y5rar7q+pQuhFUTwcOnqraTNvYGRr2jcCK3vwBbD285YE6SZbQDYFZDN3F\ngxwbkjwb+P+A51fVPSOKbdy2dWz2BJ4CXJFkA91Y1YsXyc0bBv2ZuqiqflhVXwW+TJc0LXSDHJtT\n6MY1U1WfAh4O7DOS6Oa3gX4fLSKL6nhU1W3tfRPwAbo/ShayOyaGk7b3BTskuaruaH9w/gh4Owvo\n3CbZlS5peHdVvb8VL8hzO9W+LuRzO6GqtgBX0P2dt7TlCTDA7+SdIUn6DHBQuyPFbnQ3Zrh4Up2L\ngZPa9PHA5dWuylrgtnls2pCyt9ElSAviB31AMx6bqrqrqvapqpVVtZLueq3nV9WCuivTNAb5mfog\n3U0/SLIP3fC7m0Ya5XgMcmxuBo4CSHIwXZK0eaRRzk8XAy9pd7k7ArhrYsjKIjXId2lBSPKoJHtO\nTAO/Akx5F8QFpP93x0nARWOMZagmXVv4AhbIuW3Xrp8N3FBVf91btODO7XT7uoDP7bKJ0S9JHgE8\nm+46rE/Q5QkwwLldMtPC+aCq7kvycuBf6e4W9I6q+mKSPwPWVdXFdCf+XUnW0/UgnTi+iEdnwGPz\nl8AewHvbvSxurqrnjy3oERnw2CxKAx6bfwV+Jcn1wP3AH1fVN8cX9WgMeGxeDbw9yR/SddW/dDH8\nUybJP9MNv9ynXY/1OrqLYamqt9Jdn3UMsB64Gzh5PJHOD9N9l8Yc1rDsB3ygtTFLgPdU1UfHG9Lc\nmea7vxa4IMkpdP84OWF8Ec6dafb1mUkOpft9twH4vbEFOLeeAbwY+EK7dgXgtSzMczvdvr5wgZ7b\n/YFzkuxC1yF0QVV9qP1Nc16SN9LdwffsmTaSRdC2S5IkSdLAdobhdpIkSZI0MiZJkiRJktRjkiRJ\nkiRJPSZJkiRJktRjkiRJkiRJPSZJkiRJktRjkiRJkiRJPSZJ0hgk+USS57TpNyY5c9wxSZI0Icnv\nJakkB/fKbkiycnxRSaOzZNwBSIvU64A/S7IvcBjw/DHHI0lS388A1wDHAjck2R3YD/jaWKOSRsSe\nJGkMquqTQIBXASdW1f1jDkmSpL5DgLV0SRLATwM3VFWNLyRpdEySpDFIcgiwP3BPVX1n3PFIkjTJ\nk4GLgX2TPJouafrCeEOSRsckSRqxJPsD7wZWA99L8qtjDkmSpAckWQF8s6q+D1wK/Crd8LtrxxqY\nNEImSdIIJXkk8H7g1VV1A/DnwOvHGpQkSQ/1MzzYa/RhuiF39iRpUYlDSyVJkjQhyRpg96p6Q7th\nww3AI4Gfqqot441OGg17kiRJktT3QK9RVd3Tpu81QdJiYk+SJEmSJPXYkyRJkiRJPSZJkiRJktRj\nkiRJkiRJPSZJkiRJktRjkiRJkiRJPSZJkiRJktRjkiRJkiRJPf8/Muzrvwkebf4AAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<Figure size 1008x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))\n",
    "\n",
    "# graph for hist(p)\n",
    "ax1.hist(p, bins=16)\n",
    "ax1.set_xlim((0,1))\n",
    "ax1.set_xlabel(r'$x$')\n",
    "ax1.set_ylabel('Frequency')\n",
    "ax1.set_title(r'Histogram of $p$')\n",
    "\n",
    "# graph for hist(N)\n",
    "ax2.hist(N, bins=16)\n",
    "ax2.set_xticks(range(0,35,5))\n",
    "ax2.set_xlabel(r'$N$')\n",
    "ax2.set_ylabel('Frequency')\n",
    "ax2.set_title(r'Histogram of $N$')\n",
    "\n",
    "plt.suptitle(r'Gibbs: Histograms of posterior distributions $p$ and $N$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean of p = 0.6822187872541091\n",
      "median of p = 0.6921336843178136\n"
     ]
    }
   ],
   "source": [
    "mean_p = np.mean(p)\n",
    "print('mean of p = {}'.format(mean_p))\n",
    "\n",
    "med_p = np.median(p)\n",
    "print('median of p = {}'.format(med_p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "Joseph K. Blitzstein and Jessica Hwang, Harvard University and Stanford University, &copy; 2019 by Taylor and Francis Group, LLC"
   ]
  }
 ],
 "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.3"
  },
  "notebook_info": {
   "author": "Joseph K. Blitzstein, Jessica Hwang",
   "chapter": "12. Markov Chain Monte Carlo",
   "creator": "buruzaemon",
   "section": "12.4",
   "title": "Introduction to Probability, 1st Edition"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}