{ "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, 1st Edition](https://www.crcpress.com/Introduction-to-Probability/Blitzstein-Hwang/p/book/9781466575578), 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 104 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": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAFRCAYAAACVJc8eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHuRJREFUeJzt3XuYZHV95/H3JwwIiAjIaBBGxwtLNCYCGQ1Zk2gAo4CC+sgGYyK6JOyuJN42q6PrRjZrNvg8iRfW9YKiAeIFRaIoGoMiMeRRYECiXFQQRmcYhMkCIhcdwe/+cU5r0VPdXd1d1dXV5/16nnr61Klz+f5OndOfOpc6lapCkiR1wy+MuwBJkrR0DH5JkjrE4JckqUMMfkmSOsTglySpQwx+SZI6xOCXJKlDDH5JkjrE4BcASTYmObzn+dVJnjHGkuZkjZMhyQFJvpbkh0leMe56BjXq9y7J3yZ587DnNX1a07ftYU9/KSxmHUpycpKTR1TaRDL4h6DdsLYl2Xta/yuTVJK1A05jaBvnYlXVL1fVRfMZp18bkrw0ycWLrafftBdS4zAN0t5Ba1xu7/+QvRa4qKoeUlWnLsUMh7E8l3L9GuZ6Msy6l9F2N+c6lOSgJP+S5J4klyZ5VJ9h9mw/PDx9Wv+zkpybJIspMsnfJbk5yZ1Jvp3kjxYzvVEx+IfnRuBFU0+S/Aqwy7AmnmTVsKalblkG686jgavHXMPAFru8xrW8l8H7PEqzrkNJ9gM+C7wFeBhwA/DG6cNV1e3AacCre8b9H8ATgT+oxd/D/q+AtVW1O3A08OYkv7bIaQ6dwT88ZwEv6Xl+PHBm7wBJHpnkE0m2Jrlx6pBVkrOARwGfTnJXkte2/TcmeV2SrwN3J1mV5AlJLkpyR3vI7eie6W9M8vok1yS5PckHk+zc8/qM407X+0m/reGm9pPyt5IcttCFlGR9ku+007omyfOnvb7dvOZYPof3dP9Zkq8n+UGSs6e1/eCeQ4Ufb19/80zzXGj7+rR3+imU+bRvtvd6xvb0zHf6ujPjsm+H/2/t8rs7yelJHpHkc+3wX0iy5wxtnK3OC4HfAd7Ztu3fzbCM+q63c62z81yefbe/WZZX7/o1Vx3bjd+nnQcluaKt9Wxg52njz7q99WvXXHX3eMoMy7eSPL6njt7TD4Nsd3P9P5pxm+yzfPpOa5B1CPgb4H1VdV5V3Qt8FHjKDLN6K/CsJI9LcixwIvDcqrpnptoGVVVXV9WPp562j8ctdrpDV1U+FvkANgKHA98CngDsAGyi+ZRawFqaD1mXA38O7AQ8luZT6bN6p9FnulcCa2iOHuwIXA+8oZ3GocAPgQN6hr+qHX4v4F+AN7evDTLu4X3adEDblke2/dcCj5ttOUzr91Lg4p7nxwKPbJfH7wF3A/u0r804r1mWz+E93Ze2094LuBb4z+1rOwHfBV7ZLocXANuAN8/VPuBdwLsW0d7eGgdu32zv12ztmWndGWDZbwS+CjwC2Be4FbgCOAh4EHAh8KY+y2DW9aod5iLgj+bYfrZbb+ea9jyX5yDb3/TltZFmGxikjduNP62NU+/Zq9vpvRD4CT/fPqfmNdf6OL1dM9Y91/JtXyvg8T3D/i3br0d9t7sB3p+NzLBNznc9YpZ1CNgd+BHwqJ5+LwC+0nafDJw8bZwPAP8AbAUOnmG6nwHumOHxmVnW53cB97TL9gpgt5mGHdfDPf7hmtrrfybwTeCmnteeAqyuqr+oqm1VdQPwPuC4OaZ5alVtquZT7CHAbsAp7TQupFk5X9Qz/Dvb4W8D/rLntUHG7ed+mn/8T0yyY1VtrKrvzDL8J9tP7HckuYNmI/iZqvp4VW2pqp9W1dnAdcBTFziv6U5tp30b8GngwJ62r2pf/0lVnUvzD2nOeVbVy6vq5Qtt7zTzad9s79ds7Zm+PKbWnbmWPcD/qapbquom4J+BS6rqa9Xswfw9zYeA+dQ5H/3W27mmPZ/lOcj294DltYA2zjT+1DR2BN7evmfnAJf1GW4h28Bs850y0/+FxRhkucy0TS5kWjM5jGbZfr1nO/wQzQetmbwVeBZwUlVd0W+AqnpOVe0xw+M5M024/X/xEOC3gHOBH8807LgY/MN1FvD7NHt9Z0577dHAI6eFxBto9rBms6mn+5HApqr6aU+/79LsofUb/rvtOIOOu52quh54Fc2n5luTfLQ9ZPri9rDbXUk+1zPK83o3EOABoZnkJWkuepxaBk8C9p5tXrPVN833e7rvoflHMtX2m6qq9/zdpiHNc9b29prnvGZ7v2ZszzQP6Dfbsm/d0tN9b5/nu7G9Ba1Xc9Q6td7OOu15Ls9Btr9+y5C56hhg/KlpTH/PtgumBa6Ps8233zC9/xcWY5DlMtM2uZBpzWQtcN607fBLNHv0M9mJJpDPHWD681ZV91fVxcB+wH8ZxTwWw+Afoqr6Ls1Ffkey/Qq1Cbhx2qfGh1TVkVOjzzTZnu4twJokve/bo3jgkYU1017bMo9xZ2rXh6vqN/n5qYu3VNWHqmq39nHEXNMASPJomr2sPwEe1m6gVwE/u5K237ymXhpkHjO4Gdg3ecAVuz9bTrPMc+jm0b7Z3q9Z29M7u6mOQZb9Ai14vZqm33o757TnsTzn2v76jTNl0DbOto72e8+2u+oc5lwf+81jkG1jpv8L9wC79rz2i/OY9rDe+8VO60E07QAgyWOAdcB5s4zzZOCqqrpvpgHSXN9y1wyPz8003jSrWIbn+A3+4TsBOLSq7p7W/1LgzvZCnF2S7JDkSUmmLkC5hea842wuoTkv+9okO6b5Lu1zaS5kmXJSkv2S7EWzR3P2PMbdTprvzx6a5EE059HupTkcuRAPpvlHsrWd9sto9joHmdcgy2cmX2mn8yftxU/H0B7iHnL7ZjXP9s32fs3YnlnMuuwXYUHrVR/91ttZpz3P5TnX9jfqNn4FuA94RfuevYA+79kA6+NCt4OZ/i9cCfx+uzyeDTx92nizzW9Y7/1ip3UZ8PT2SOQa4MPAf29PL8zkQJq2z6iqjujZuZn+2G5nJ8nDkxyXZLd2eT6L5lTFhQO0YUkZ/ENWVd+pqg19+t9PsyIfSHNU4N+A9wMPbQf5K+CN7WHIP5th2ttoviJyRDv+u4CXVNU3ewb7MPCPNBcu3UBzkdSg4/bzIOCUdpzvAw+n+ccxb1V1Dc3Vt1+h+YfyKzQXGg0yrzmXzyzz3UZzsc8JNBfm/AHN+cMfzzFPkrwnyXvm1dCZDdy+2d6vOdrT1wDLfkEWsV5Nt916O8C057M859r+RtrGnvfspcDtNBdX9jvMPNf2ttDtoO//BZoLRJ9Lsx69GPjktPFmnN8Q3/vFTutCmusHvg1cDJxVVe+bY5wnM0fwL0DRHNbfTPMe/zXwqqr61JDns2h54CknTbIkG2mufP3CuGtZ7pJcArynqj447lqGYZLb43qrUUp7176qOnm8lSwf7vGrE5I8PckvtodZjwd+ldkv/lnWVlp7JC2dlXynJ6nXAcDHaK4q/g7wwqq6ebwlLcpKa480KheNu4DlxkP9kiR1iIf6JUnqEINfkqQOWZHn+Pfee+9au3btuMuQJGnJXH755f9WVavnGm5FBv/atWvZsGG7r9JLkrRiJZnt9wl+xkP9kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdciKvFe/pO2tXX/+ksxn4ylHLcl8JC2Me/ySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKH+CM90jKwVD+gI0nu8UuS1CEGvyRJHWLwS5LUISML/iQfSHJrkqt6+u2V5IIk17V/92z7J8mpSa5P8vUkB/eMc3w7/HVJjh9VvZIkdcEo9/j/Fnj2tH7rgS9W1f7AF9vnAEcA+7ePE4F3Q/NBAXgT8OvAU4E3TX1YkCRJ8zey4K+qLwO3Tet9DHBG230G8Lye/mdW46vAHkn2AZ4FXFBVt1XV7cAFbP9hQpIkDWipz/E/oqpuBmj/Prztvy+wqWe4zW2/mfpLkqQFWC4X96VPv5ql//YTSE5MsiHJhq1btw61OEmSVoqlDv5b2kP4tH9vbftvBtb0DLcfsGWW/tupqtOqal1VrVu9evXQC5ckaSVY6uA/D5i6Mv944FM9/V/SXt1/CPCD9lTA54HfTbJne1Hf77b9JEnSAozslr1JPgI8A9g7yWaaq/NPAT6W5ATge8Cx7eCfBY4ErgfuAV4GUFW3JflfwGXtcH9RVdMvGJQkSQMaWfBX1YtmeOmwPsMWcNIM0/kA8IEhliZJUmctl4v7JEnSEjD4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqkLEEf5JXJ7k6yVVJPpJk5ySPSXJJkuuSnJ1kp3bYB7XPr29fXzuOmiVJWgmWPPiT7Au8AlhXVU8CdgCOA94CvK2q9gduB05oRzkBuL2qHg+8rR1OkiQtwLgO9a8CdkmyCtgVuBk4FDinff0M4Hlt9zHtc9rXD0uSJaxVkqQVY8mDv6puAv4a+B5N4P8AuBy4o6ruawfbDOzbdu8LbGrHva8d/mFLWbMkSSvFOA7170mzF/8Y4JHAg4Ej+gxaU6PM8lrvdE9MsiHJhq1btw6rXEmSVpRxHOo/HLixqrZW1U+Ac4F/D+zRHvoH2A/Y0nZvBtYAtK8/FLht+kSr6rSqWldV61avXj3qNkiSNJHGEfzfAw5Jsmt7rv4w4BrgS8AL22GOBz7Vdp/XPqd9/cKq2m6PX5IkzW0c5/gvoblI7wrgG20NpwGvA16T5Hqac/int6OcDjys7f8aYP1S1yxJ0kqxau5Bhq+q3gS8aVrvG4Cn9hn2R8CxS1GXJEkrnXfukySpQwx+SZI6xOCXJKlDDH5JkjrE4JckqUPGclW/pJVr7frzRz6PjaccNfJ5SCuVwS/NYilCTJKWkof6JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDhko+JM8adSFSJKk0Rt0j/89SS5N8vIke4y0IkmSNDIDBX9V/SbwYmANsCHJh5M8c6SVSZKkoRv4HH9VXQe8EXgd8HTg1CTfTPKCURUnSZKGa9Bz/L+a5G3AtcChwHOr6glt99tGWJ8kSRqiVQMO907gfcAbqureqZ5VtSXJG0dSmSRJGrpBg/9I4N6quh8gyS8AO1fVPVV11siqkyRJQzXoOf4vALv0PN+17SdJkibIoMG/c1XdNfWk7d51NCVJkqRRGTT4705y8NSTJL8G3DvL8JIkaRka9Bz/q4CPJ9nSPt8H+L3RlCRJkkZloOCvqsuS/BJwABDgm1X1k5FWJkmShm7QPX6ApwBr23EOSkJVnTmSqiRJ0kgMFPxJzgIeB1wJ3N/2LsDglyRpggy6x78OeGJV1TBm2v7Qz/uBJ9F8gPiPwLeAs2mOKmwE/kNV3Z4kwDto7iVwD/DSqrpiGHVIktQ1g17VfxXwi0Oc7zuAf6iqXwKeTHMr4PXAF6tqf+CL7XOAI4D928eJwLuHWIckSZ0y6B7/3sA1SS4FfjzVs6qOnu8Mk+wO/Dbw0nYa24BtSY4BntEOdgZwEc0PAh0DnNkebfhqkj2S7FNVN8933pIkdd2gwX/yEOf5WGAr8MEkTwYuB14JPGIqzKvq5iQPb4ffF9jUM/7mtp/BL0nSPA10qL+q/onmvPuObfdlwELPs68CDgbeXVUHAXfz88P6/aRfSdsNlJyYZEOSDVu3bl1gaZIkrWyD/izvHwPnAO9te+0LfHKB89wMbK6qS9rn59B8ELglyT7t/PYBbu0Zfk3P+PsBW5imqk6rqnVVtW716tULLE2SpJVt0Iv7TgKeBtwJUFXXAQ+fdYwZVNX3gU1JDmh7HQZcA5wHHN/2Ox74VNt9HvCSNA4BfuD5fUmSFmbQc/w/rqptzTfrIMkq+hxun4c/BT6UZCfgBuBlNB9CPpbkBOB7wLHtsJ+l+Srf9TRf53vZIuYrSVKnDRr8/5TkDcAuSZ4JvBz49EJnWlVX0twbYLrD+gxbNEccJEnSIg16qH89zZX43wD+E81e+BtHVZQkSRqNQX+k56fA+9qHJEmaUIPeq/9G+pzTr6rHDr0iSZI0MvO5V/+UnWkuvNtr+OVIkqRRGvQGPv+v53FTVb0dOHTEtUmSpCEb9FD/wT1Pf4HmCMBDRlKRJEkamUEP9f9NT/d9tD+bO/RqJEnSSA16Vf/vjLoQSZI0eoMe6n/NbK9X1VuHU44kSRql+VzV/xSa++YDPBf4Mg/8uVxJkrTMDRr8ewMHV9UPAZKcDHy8qv5oVIVJkqThG/SWvY8CtvU83wasHXo1kiRppAbd4z8LuDTJ39Pcwe/5wJkjq0qSJI3EoFf1/2WSzwG/1fZ6WVV9bXRlSZKkURj0UD/ArsCdVfUOYHOSx4yoJkmSNCIDBX+SNwGvA17f9toR+LtRFSVJkkZj0D3+5wNHA3cDVNUWvGWvJEkTZ9Dg31ZVRfvTvEkePLqSJEnSqAwa/B9L8l5gjyR/DHwBeN/oypIkSaMw6FX9f53kmcCdwAHAn1fVBSOtTJIkDd2cwZ9kB+DzVXU4YNhLkjTB5jzUX1X3A/ckeegS1CNJkkZo0Dv3/Qj4RpILaK/sB6iqV4ykKkmSNBKDBv/57UOSJE2wWYM/yaOq6ntVdcZSFSRJkkZnrnP8n5zqSPKJEdciSZJGbK7gT0/3Y0dZiCRJGr25gr9m6JYkSRNorov7npzkTpo9/13abtrnVVW7j7Q6SZI0VLMGf1XtsFSFSNKg1q4f/ZeMNp5y1MjnIY3DoPfqlyRJK4DBL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1yNiCP8kOSb6W5DPt88ckuSTJdUnOTrJT2/9B7fPr29fXjqtmSZIm3Vz36h+lVwLXAlP3+38L8Laq+miS9wAnAO9u/95eVY9Pclw73O+No2AtL0tx21ZJWmnGssefZD/gKOD97fMAhwLntIOcATyv7T6mfU77+mHt8JIkaZ7Gdaj/7cBrgZ+2zx8G3FFV97XPNwP7tt37ApsA2td/0A4vSZLmacmDP8lzgFur6vLe3n0GrQFe653uiUk2JNmwdevWIVQqSdLKM449/qcBRyfZCHyU5hD/24E9kkxdc7AfsKXt3gysAWhffyhw2/SJVtVpVbWuqtatXr16tC2QJGlCLXnwV9Xrq2q/qloLHAdcWFUvBr4EvLAd7HjgU233ee1z2tcvrKrt9vglSdLcltP3+F8HvCbJ9TTn8E9v+58OPKzt/xpg/ZjqkyRp4o3z63xU1UXARW33DcBT+wzzI+DYJS1MkqQVajnt8UuSpBEz+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6pAlD/4ka5J8Kcm1Sa5O8sq2/15JLkhyXft3z7Z/kpya5PokX09y8FLXLEnSSjGOPf77gP9aVU8ADgFOSvJEYD3wxaraH/hi+xzgCGD/9nEi8O6lL1mSpJVh1VLPsKpuBm5uu3+Y5FpgX+AY4BntYGcAFwGva/ufWVUFfDXJHkn2aacjSSOxdv35SzKfjacctSTzkaaM9Rx/krXAQcAlwCOmwrz9+/B2sH2BTT2jbW77SZKkeRpb8CfZDfgE8KqqunO2Qfv0qz7TOzHJhiQbtm7dOqwyJUlaUcYS/El2pAn9D1XVuW3vW5Ls076+D3Br238zsKZn9P2ALdOnWVWnVdW6qlq3evXq0RUvSdIEG8dV/QFOB66tqrf2vHQecHzbfTzwqZ7+L2mv7j8E+IHn9yVJWpglv7gPeBrwh8A3klzZ9nsDcArwsSQnAN8Djm1f+yxwJHA9cA/wsqUtV5KklWMcV/VfTP/z9gCH9Rm+gJNGWpQkSR3hnfskSeoQg1+SpA4x+CVJ6hCDX5KkDjH4JUnqEINfkqQOMfglSeoQg1+SpA4x+CVJ6hCDX5KkDhnHvfrVAWvXnz/uEiRJfbjHL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKHGPySJHWIwS9JUocY/JIkdYjBL0lShxj8kiR1iMEvSVKHrBp3AZLUZWvXnz/yeWw85aiRz0OTwz1+SZI6xOCXJKlDDH5JkjrE4JckqUMMfkmSOsSr+jtmKa4gliQtX+7xS5LUIQa/JEkdYvBLktQhBr8kSR3ixX2StMJ5W2D1co9fkqQOmZjgT/LsJN9Kcn2S9eOuR5KkSTQRwZ9kB+D/AkcATwRelOSJ461KkqTJMynn+J8KXF9VNwAk+ShwDHDNUhXgOTJJ0kowKcG/L7Cp5/lm4NfHVMvIeFc9SZNqqf5/LcUO0krf0ZuU4E+ffvWAAZITgRPbp3cl+dbIq9re3sC/jWG+o2BblqeV1BZYWe2xLUsgb5n3KMuyLQtoB8zdlkcPMpFJCf7NwJqe5/sBW3oHqKrTgNOWsqjpkmyoqnXjrGFYbMvytJLaAiurPbZlebIt25uIi/uAy4D9kzwmyU7AccB5Y65JkqSJMxF7/FV1X5I/AT4P7AB8oKquHnNZkiRNnIkIfoCq+izw2XHXMYexnmoYMtuyPK2ktsDKao9tWZ5syzSpqrmHkiRJK8KknOOXJElDYPAPWZJjk1yd5KdJJvJK0pVye+QkH0hya5Krxl3LYiVZk+RLSa5t169XjrumhUqyc5JLk/xr25b/Oe6aFivJDkm+luQz465lMZJsTPKNJFcm2TDuehYjyR5JzknyzXa7+Y1x17QQSQ5o34+px51JXrWoaXqof7iSPAH4KfBe4M+qaqI2nvb2yN8GnknzNcrLgBdV1ZLdJXFYkvw2cBdwZlU9adz1LEaSfYB9quqKJA8BLgeeN6HvS4AHV9VdSXYELgZeWVVfHXNpC5bkNcA6YPeqes6461moJBuBdVW17L73Pl9JzgD+uare334bbNequmPcdS1G+//5JuDXq+q7C52Oe/xDVlXXVtU4bh40LD+7PXJVbQOmbo88carqy8Bt465jGKrq5qq6ou3+IXAtzR0tJ0417mqf7tg+JnYPJMl+wFHA+8ddixpJdgd+GzgdoKq2TXrotw4DvrOY0AeDX9vrd3vkiQyYlSrJWuAg4JLxVrJw7aHxK4FbgQuqamLbArwdeC3Nkb5JV8A/Jrm8vRvqpHossBX4YHsK5v1JHjzuoobgOOAji52Iwb8ASb6Q5Ko+j4ncM55mztsja3yS7AZ8AnhVVd057noWqqrur6oDae7C+dQkE3kqJslzgFur6vJx1zIkT6uqg2l+CfWk9nTZJFoFHAy8u6oOAu4GJvZ6JYD2dMXRwMcXO62J+R7/clJVh4+7hhGa8/bIGo/2fPgngA9V1bnjrmcYquqOJBcBzwYm8SLMpwFHJzkS2BnYPcnfVdUfjLmuBamqLe3fW5P8Pc2pvy+Pt6oF2Qxs7jmSdA4THvw0H8auqKpbFjsh9/g1nbdHXobaC+JOB66tqreOu57FSLI6yR5t9y7A4cA3x1vVwlTV66tqv6paS7OtXDipoZ/kwe2Fo7SHxX+XyfwwRlV9H9iU5IC212Es4c+4j8iLGMJhfjD4hy7J85NsBn4DOD/J58dd03xU1X3A1O2RrwU+Nqm3R07yEeArwAFJNic5Ydw1LcLTgD8EDu35Ws+R4y5qgfYBvpTk6zQfNC+oqon+GtwK8Qjg4iT/ClwKnF9V/zDmmhbjT4EPtevZgcD/HnM9C5ZkV5pvWg3lSJ9f55MkqUPc45ckqUMMfkmSOsTglySpQwx+SZI6xOCXJKlDDH5JkjrE4JckqUMMfklD0f7wzjuSXN3+pvtjx12TpO0Z/JKG5fXADVX1y8CpwMvHXI+kPvyRHkmL1t7b/flV9WttrxtpfqNe0jJj8EsahsOBNUmubJ/vBXxhjPVImoGH+iUNw4HAn1fVgVV1IPCPwJVzjCNpDAx+ScOwJ3APQJJVND/p+umxViSpL4Nf0jB8Gzik7X41zU+63jjGeiTNwJ/llbRoSfYEPgfsDXwFOLGq7h1vVZL6MfglSeoQD/VLktQhBr8kSR1i8EuS1CEGvyRJHWLwS5LUIQa/JEkdYvBLktQhBr8kSR3y/wFw7OOB5yBZ7wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "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 104 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/rSwjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X2YJVV57/3vTwbwBeKADASHwdFIDEYikNGQmJMoaCIQHcwDCcYocpGQ54iJRhOZ+OQ6auI5Z0xOQkISXzAYwaMBxBeIqBFBjieJoKMgImgYcYQBZEZl8AUFwfv5o1ZD0dPds3um99493d/Pde1rV61aVfuuqt29+u61qipVhSRJkiSp87BxByBJkiRJ84lJkiRJkiT1mCRJkiRJUo9JkiRJkiT1mCRJkiRJUo9JkiRJkiT1mCRJkiRJUo9JkiRJkiT1mCRJGookX0zyzDa9IcmzZ6g74/Jh6sepByV5UpKrk3wnyR+MO55BDft8JnlnkjcO47Nm8zOzI9te6PrnSJK2l0mSpO2S5MQkVyX5XpJNbfplSQJQVT9dVVeMMb6t/shM8tIk/9YvGyTOcSZxY/Qa4Iqq2rOqzhzFB87FcR7l927Qzxp0v+Yq9qk+b9w/j/NFkr2SVJJPTSp/W5IzxhWXpPnHJEnSrCV5NfC3wF8CPw7sB/y/wDOA3cYY2k4nyZJxxzCNxwFfHHcQg9rR4zjO8zCPvwML0aHA14EnJ9l/Uvk14wlJ0nxkkiRpVpI8Gvgz4GVVdWFVfac6V1fVi6rqnlZv8n+zn5bk+iR3JvmnJA+ftOkplyc5PcmtbdjXl5McNcf780CcU31WkncBBwL/kuS7SV7T6h6c5IokW9pQpuf3tnl4b6jae5Oc3x/+0z7z9CTXAt9LsiTJmiRfaetcn+QFk+r/cZJrW8/d2Un2S/KRVv/jSfbq1R/omE23D0kuB54F/H3b55+c5rj9yXTndBvHZzbH+bFJ3pdkc5Kvpjf0b5rj+JDv3Tbi2Gr9KfbzsCSfa7GeDzx80vr9z5rNfm0zdmb4mUnXG/LE3vw7k7xxG5838T2f9pj06v5R+77d1b6/s/p5TPKiJP/R1v16kluSHD1V3VZ/W9//KePZ1jmawqHAOuBSYOL7vgtwCHD1DOtJWmyqypcvX74GfgHPBe4Dlmyj3gbg2b3p64AVwN7AvwNvnFR3q+XAk4BbgMe2eiuBn+it92bgzdv6/F7ZS4F/m6reTJ81eVvArsB64LV0PWdHAt9p29gN+Brwilbv14F7p9jfa9r+PqKVnQA8lu6fV78JfA/Yv1f/Sroeu+XAJuBzwGHA7sDlwOta3RmP2SD70JZfAfzONs7vlOd0G8dnNsf5YcBngf/WtvME4CbgV2c4jg9sY4B93Gr9Sfs4cS7/sG3reOCHvf3sf9bA+zVg7NMe37a8gCf25t85VVxTfM9nPCa9up+m+z7uDdxA11M80HerLVsLfB/4jfaZfwR8bYbv07a+/1vFM8g5muJzzgX+HPht4COt7KeBe4Bdx/V71ZcvX/PvZU+SpNnaB/hGVd03UdD+Y7wlyfeT/NI06/19Vd1SVd8C/jvwwgGW30+XBDw5ya5VtaGqvjKxQlW9rKpeNkOsH2xxbUmyhS6pms6MnzXJEcAewNqqureqLgc+1GI+AlgCnFlVP6yq99P9gTfZmW1/v9/25b1VdVtV/aiqzgduBJ7eq/93VXVHVd0K/F/gqup67+4BPkCXMM1mP2bah0FNd05n2vZsjvPTgGVV9WdtOzcBbwdO7NV5yHHcjn3c1vq7An/TzuWFwGemiXU2+zXIZ8O2f2a2x6Dn/cz2ffwW8C90PTCz2cdDgDOq6oKq+iFdcnJgtu5BBgb6/k8Vz8T+DHqO4MFhdZcA/yXJnq3suhanJAEOt5M0e98E9ukPTaqqX6iqpW3ZdL9XbulNf43uv8IzLq+q9cArgdcDm5Kcl2TyejM5rqqWTryAaROqWX7WY4FbqupHk2Je3pbdWlU1zb5NWZbkJUmu6SV0T6FLSCfc0Zv+/hTze8xyP2bah0FNd06n3fYsj/PjgMdOSnRfS9ejNlUMkw2yj9taf/K5/NpUFbfzuzrTZ09ePtXPzPYY9Lx/vTd9N7DHLPfxEODC3vy+wHer6gdTVR7g+79VPL39GegcJdkdOBi4pqrupPvnxdF0/2DweiRJD2GSJGm2PkU3NGX1LNdb0Zs+ELhtkOVV9Z6q+kW6P5gLeNMsP3dgM3xWTap6G7AiSf936IHArcDtwPKku8tf09+3Bz5uYiLJ4+h6SF4OPKYldNcBmWK9HdmPQfdhUNOd0xm3PYvjfAvw1X6iW93d9o7p1Zm8Tt8g+zjT+lOdywOnqzyL/Rrks2Hmn5m7gUf25n98wO3u0Hkf5LuVZGmLfXOv+HjgI1Ntcwe//7M5R0+h+4fCTW3+g8BxdEmS1yNJegiTJEmzUlVbgDcAb05yfJI9kjwsyaHAo2ZY9bQkByTZm6434PxtLU/3rJ4j23+Af0D3B879c79XDzwXaLrPuoPuepgJV9FdM/GaJLume/7M84Dz6JLI+4GXt4vxV/PQYUNTeRTdH5ybWywn0/1BN9f70TfTPgxqunM67bZneZw/DXy73SzgEUl2SfKUJE8bML4d3cdP0V1/9wftXP4605zLWe7XoGb6mbkG+K12TJ4L/HJv2Uyft93HZBbfrUNa+W+143YsXS/u66fZ9I58/wc+R3TJ0Od7vU4XAcdgT5KkKZgkSZq1qvoL4FV0z9LZRPdH2duA04H/mGa19wAfo/sv7k10N2bY1vLd6S4A/wbdcJt96f5YBCDJW5O8dU52aubP+p/An7ahQH9UVffS3Rnr6Fb/zcBLqupLbdmvA6cAW+guEP8QXe/blKrqeuCv6P7gu4Puj8x/H8J+9D9z2n2YxWdNeU63se3ZHOf76f6APxT4alvnH4FHDxLcju5j71y+FLiT7oYC75+m+sD7NchnNzP9zLyC7thsAV5E1ysyYdrP28FjMtB3i+77+27g5+mO2xvohr5eP9VGd+T7P8tz9JDbfFfV1+huCrEU+Pwgnydp8chDh/FKkuZakquAt1bVP407lrmSZAPd3e8+Pu5YNL8keQvwn1Xlw1kl7bTsSZKkOZbkl5P8eBv+cxLwM8BHxx2XNCKH0N2mW5J2Wj7lW5Lm3pOAC+juwPUV4Piqun28IUkj8xRgNsM2JWnecbidJEmSJPU43E6SJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekyStCgl+WKSZ447jnFJ8qQkVyf5TpI/GHc8kqSHsp2yndJ4mSRpwUmyIcmzJ5W9NMm/TcxX1U9X1RWz3c4C8hrgiqras6rOHHcwkrSY2E4NZMZ2KsleSSrJpyaVvy3JGSOLUguWSZI0BkmWjDmExwFfHHMMkqR5aidopw4Fvg48Ocn+k8qvGWZgWhxMkrQo9f/7luT0JLe2Lv0vJzkqybuAA4F/SfLdJK9pdQ9OckWSLW0oxPN72zy8NzTgvUnOT/LGSZ95epJrge8lWZJkTZKvtHWuT/KCSfX/OMm1Sb6X5Owk+yX5SKv/8SR7zbCPU8aa5HLgWcDft337ySnWfVGS/2j78PUktyQ5eocPvCRpILZTM7dTdMnQOuBSYGK9XYBDgKu387BLDzBJ0qKW5EnAy4GnVdWewK8CG6rqxcDNwPOqao+q+oskuwL/AnwM2Bf4feDd6cZN7wZ8AHgnsDfwz8ALtvpAeCFwLLC0qu4DvgL8F+DRwBuA/z3pP2L/D/Ac4CeB5wEfAV4L7EP38zvlOO2ZYq2qI4H/C7y87dt/TrGJQ4DDgPcBK4C/Bd467YGUJA2F7dS07dRhdD1GHwSOa2U/BewC3DDVZ0qzYZKkheqD7T9TW5JsAd48Tb37gd3puut3raoNVfWVaeoeAewBrK2qe6vqcuBDdA3KEcAS4Myq+mFVvR/49BTbOLOqbqmq7wNU1Xur6raq+lFVnQ/cCDy9V//vquqOqrqVrsG4qqqurqp76Bq7w7Yj1kEcApxRVRdU1Q+Bc4EDkzx8wPUlSTOzndqxdmpiWN0lwH9Jsmcru661W9IOMUnSQnVcVS2deAEvm6pSVa0HXgm8HtiU5Lwkj51mm48FbqmqH/XKvgYsb8turarqLbtlim08pCzJS5Jc02skn0L337cJd/Smvz/F/B7bEesgDgEu7M3vC3y3qn4w4PqSpJnZTm1nO5Vkd+Bg4JqqupMu2TuaB3uXpB1mkqRFr6reU1W/SHeRaAFvmlg0qeptwIok/Z+bA4FbgduB5UnSW7Ziqo+bmEjyOODtdMMoHtMayeuATLHebM0U64ySLKWLfXOv+Hi6IRSSpBGzndrKU+gSsJva/MSQu8PweiTNEZMkLWptnPaR7b9SP6D7pXt/W3wH8IRe9auA7wGvSbJruudXPA84D/hUW+/l7ULX1Tx0OMJUHkXXGG1usZxM94t/LswU67YcQrcvv9X25Vi6/3C+fo5ikyQNyHZqSocBn+/1il0EHIM9SZpDJkla7HYH1gLfoLuV6L50F5wC/E/gT9sQgz+qqnvp7qBzdKv/ZuAlVfWltuzXgVOALcBv042tvme6D66q64G/omu47qBLTv59LnZqplgHWP0Q4N3AzwN30l2oe1yLV5I0WrZTW3vIbb6r6mvABmAp8Pm5iE/KQ4emSporSa4C3lpV/zTuWGYjyVuA/6wqH8YnSQvYztpOSaNgT5I0R5L8cpIfb8MYTgJ+BvjouOPaDofg7VMlacFZQO2UNHTjfpqytJA8CbiA7k4+XwGOr6rbxxvSdnkKMMhwB0nSzmWhtFPS0DncTpIkSZJ6HG4nSZIkST0mSZIkSZLUs1Nfk7TPPvvUypUrxx2GJC1qn/3sZ79RVcvGHcd8ZDslSeO3Pe3UTp0krVy5knXr1o07DEla1JJ8bdwxzFe2U5I0ftvTTg11uF2SDUm+kOSaJOta2d5JLk1yY3vfq5UnyZlJ1ie5Nsnhw4xNkiRJkqYyimuSnlVVh1bVqja/Brisqg4CLmvz0D1x+aD2OhV4ywhikyRJkqSHGMeNG1YD57Tpc4DjeuXnVudKYGmS/ccQnyRJkqRFbNhJUgEfS/LZJKe2sv0mHlzW3vdt5cuBW3rrbmxlD5Hk1CTrkqzbvHnzEEOXJEmStBgN+8YNz6iq25LsC1ya5Esz1M0UZVs96baqzgLOAli1apVPwpUkSZI0p4bak1RVt7X3TcAHgKcDd0wMo2vvm1r1jcCK3uoHALcNMz5JkiRJmmxoSVKSRyXZc2Ia+BXgOuBi4KRW7STgojZ9MfCSdpe7I4C7JoblSZIkSdKoDHO43X7AB5JMfM57quqjST4DXJDkFOBm4IRW/8PAMcB64G7g5CHGJkmSJElTGlqSVFU3AU+dovybwFFTlBdw2rDikSRJkqRBjOMW4JIkSZI0b5kkSZIkSVKPSZIkSZIk9ZgkSZIkSVLPsB8mK0mztnLNJXOynQ1rj52T7Uha+Gb7e8ffL9LCZk+SJEmSJPWYJEmSJElSj0mSJEmSJPWYJEmSJElSj0mSJEmSJPWYJEmSJElSj0mSJEmSJPWYJEmSJElSj0mSJEmSJPWYJEmSJElSz5JxByBJw7JyzSVzsp0Na4+dk+1IkqSdgz1JkiRJktRjT5KkOTNXPTfSqCT5Q+B3gAK+AJwM7A+cB+wNfA54cVXdm2R34FzgZ4FvAr9ZVRvGEbckabjsSZIkLUpJlgN/AKyqqqcAuwAnAm8Czqiqg4A7gVPaKqcAd1bVE4EzWj1J0gJkkiRJWsyWAI9IsgR4JHA7cCRwYVt+DnBcm17d5mnLj0qSEcYqSRoRkyRJ0qJUVbcC/wu4mS45ugv4LLClqu5r1TYCy9v0cuCWtu59rf5jJm83yalJ1iVZt3nz5uHuhCRpKLwmSZK0KCXZi6536PHAFuC9wNFTVK2JVWZY9mBB1VnAWQCrVq3aarlGw2skJe0Ie5IkSYvVs4GvVtXmqvoh8H7gF4ClbfgdwAHAbW16I7ACoC1/NPCt0YYsSRoFkyRJ0mJ1M3BEkke2a4uOAq4HPgEc3+qcBFzUpi9u87Tll1eVPUWStACZJEmSFqWquoruBgyfo7v998PohsmdDrwqyXq6a47ObqucDTymlb8KWDPyoCVJI+E1SZKkRauqXge8blLxTcDTp6j7A+CEUcQlSRove5IkSZIkqcckSZIkSZJ6TJIkSZIkqcdrkiRpG+bqeSsb1h47J9uRJEnDZU+SJEmSJPWYJEmSJElSj8PtJM3ZcDJJkqSFwJ4kSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSepZMu4AJG2/lWsuGXcIkiRJC449SZIkSZLUY5IkSZIkST0mSZIkSZLU4zVJkiRp3vMaTEmjZE+SJEmSJPWYJEmSJElSj0mSJEmSJPWYJEmSJElSj0mSJEmSJPWYJEmSFqUkT0pyTe/17SSvTLJ3kkuT3Nje92r1k+TMJOuTXJvk8HHvgyRpOIZ+C/AkuwDrgFur6teSPB44D9gb+Bzw4qq6N8nuwLnAzwLfBH6zqjYMOz5J0uJUVV8GDoUH2qpbgQ8Aa4DLqmptkjVt/nTgaOCg9vo54C3tXYvQbG5JvmHtsUOMRNIwjOI5Sa8AbgB+rM2/CTijqs5L8lbgFLqG5hTgzqp6YpITW73fHEF8kjQSc/WcF//gGoqjgK9U1deSrAae2crPAa6gS5JWA+dWVQFXJlmaZP+qun0cAUuShmeow+2SHAAcC/xjmw9wJHBhq3IOcFybXt3macuPavUlSRq2E4F/btP7TSQ+7X3fVr4cuKW3zsZW9hBJTk2yLsm6zZs3DzFkSdKwDPuapL8BXgP8qM0/BthSVfe1+X4D80Dj05bf1eo/hI2PJGkuJdkNeD7w3m1VnaKstiqoOquqVlXVqmXLls1FiJKkERtakpTk14BNVfXZfvEUVWuAZQ8W2PhIkubW0cDnquqONn9Hkv0B2vumVr4RWNFb7wDgtpFFKUkamWH2JD0DeH6SDXQ3ajiSrmdpaZKJa6H6DcwDjU9b/mjgW0OMT5IkgBfy4FA7gIuBk9r0ScBFvfKXtLvcHQHc5fVIkrQwDS1Jqqo/qaoDqmol3Vjvy6vqRcAngONbtcmNz0SjdHyrv1VPkiRJcyXJI4HnAO/vFa8FnpPkxrZsbSv/MHATsB54O/CyEYYqSRqhUdzdbrLTgfOSvBG4Gji7lZ8NvCvJeroepBPHEJskaRGpqruZdP1rVX2T7m53k+sWcNqIQpMkjdFIkqSquoLuFqpU1U3A06eo8wPghFHEI0mSJEnTGfbd7SRJkiRpp2KSJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1GOSJEmSJEk9JkmSJEmS1LNk3AFIkmZn5ZpL5mQ7G9YeOyfbkSRpobEnSZIkSZJ6TJIkSZIkqcckSZIkSZJ6TJIkSZIkqcckSZIkSZJ6TJIkSZIkqcckSZIkSZJ6TJIkSZIkqcckSZK0aCVZmuTCJF9KckOSn0+yd5JLk9zY3vdqdZPkzCTrk1yb5PBxxy9JGg6TJEnSYva3wEer6qeApwI3AGuAy6rqIOCyNg9wNHBQe50KvGX04UqSRmHJuAOQFqOVay4ZdwjSopfkx4BfAl4KUFX3AvcmWQ08s1U7B7gCOB1YDZxbVQVc2Xqh9q+q20ccuiRpyOxJkiQtVk8ANgP/lOTqJP+Y5FHAfhOJT3vft9VfDtzSW39jK5MkLTAmSZKkxWoJcDjwlqo6DPgeDw6tm0qmKKutKiWnJlmXZN3mzZvnJlJJ0kg53E6aBYfJSQvKRmBjVV3V5i+kS5LumBhGl2R/YFOv/ore+gcAt03eaFWdBZwFsGrVqq2SKEnS/GdPkiRpUaqqrwO3JHlSKzoKuB64GDiplZ0EXNSmLwZe0u5ydwRwl9cjSdLCZE+SJGkx+33g3Ul2A24CTqb7B+IFSU4BbgZOaHU/DBwDrAfubnWlbZrtKIQNa48dUiSSBmWSJElatKrqGmDVFIuOmqJuAacNPShJ0tg53E6SJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKlnaElSkocn+XSSzyf5YpI3tPLHJ7kqyY1Jzk+yWyvfvc2vb8tXDis2SZIkSZrOMHuS7gGOrKqnAocCz01yBPAm4IyqOgi4Ezil1T8FuLOqngic0epJkiRJ0kgNLUmqznfb7K7tVcCRwIWt/BzguDa9us3Tlh+VJMOKT5IkSZKmsmSYG0+yC/BZ4InAPwBfAbZU1X2tykZgeZteDtwCUFX3JbkLeAzwjUnbPBU4FeDAAw8cZviStKCtXHPJuEOQJGleGuqNG6rq/qo6FDgAeDpw8FTV2vtUvUa1VUHVWVW1qqpWLVu2bO6ClSRJkiRGdHe7qtoCXAEcASxNMtGDdQBwW5veCKwAaMsfDXxrFPFJkiRJ0oRh3t1uWZKlbfoRwLOBG4BPAMe3aicBF7Xpi9s8bfnlVbVVT5IkSZIkDdMwr0naHzinXZf0MOCCqvpQkuuB85K8EbgaOLvVPxt4V5L1dD1IJw4xNkmSJEma0tCSpKq6FjhsivKb6K5Pmlz+A+CEYcUjSdJkSTYA3wHuB+6rqlVJ9gbOB1YCG4DfqKo72x1X/xY4BrgbeGlVfW4ccUuShmug4XZJnjLsQCRJ2l472E49q6oOrapVbX4NcFl7nt9lbR7gaOCg9joVeMsOfKYkaR4b9Jqktyb5dJKXTVxnJEnSPDKX7VT/uX2Tn+d3bnsO4JV0NyLafwc/S5I0Dw2UJFXVLwIvorv73Lok70nynKFGJknSgHagnSrgY0k+257DB7BfVd3etns7sG8rf+B5fk3/WX+SpAVk4GuSqurGJH8KrAPOBA5r47NfW1XvH1aAkiQNYjvbqWdU1W1J9gUuTfKlGT5ioOf5+dBzSdr5DXpN0s8kOYPuFt5HAs+rqoPb9BlDjE+SpG3a3naqqm5r75uAD9DdWOiOiWF07X1Tq/7A8/ya/rP++tv0oeeStJMb9Jqkvwc+Bzy1qk6buJtPa1z+dFjBSZI0oFm3U0kelWTPiWngV4DreOhz+yY/z+8l6RwB3DUxLE+StLAMOtzuGOD7VXU/QJKHAQ+vqrur6l1Di06SpMFsTzu1H/CBbkQeS4D3VNVHk3wGuCDJKcDNPPh4ig+3z1lPdwvwk4e2N5KksRo0Sfo48Gzgu23+kcDHgF8YRlCSJM3SrNup9ty+p05R/k3gqCnKCzhtLoKVJM1vgw63e3hVTTQ8tOlHDickSZJmzXZKkjRnBk2Svpfk8ImZJD8LfH84IUmSNGu2U5KkOTPocLtXAu9NMnEXn/2B3xxOSJIkzZrtlCRpzgyUJFXVZ5L8FPAkuudEfKmqfjjUyCRJGpDtlCRpLg38MFngacDKts5hSaiqc4cSlSRJs2c7JUmaEwMlSUneBfwEcA1wfysuwMZHkjR2tlOSpLk0aE/SKuDJ7fankiTNN7ZTkqQ5M+jd7a4DfnyYgUiStANspyRJc2bQnqR9gOuTfBq4Z6Kwqp4/lKgkSZod2ylJ0pwZNEl6/TCDkCRpB71+3AFIkhaOQW8B/n+SPA44qKo+nuSRwC7DDU2SpMHYTkmS5tJA1yQl+V3gQuBtrWg58MFhBSVJ0mzYTkmS5tKgN244DXgG8G2AqroR2HdYQUmSNEu2U5KkOTPoNUn3VNW9SQBIsoTu+ROSJM0HtlNaMFauuWRW9TesPXZIkUiL16A9Sf8nyWuBRyR5DvBe4F+GF5YkSbNiOyVJmjODJklrgM3AF4DfAz4M/OmwgpIkaZZspyRJc2bQu9v9CHh7e0mSNK/YTkmS5tJASVKSrzLF2O6qesKcRyRJ0izZTkmS5tKgN25Y1Zt+OHACsPfchyNJ0naxnZIkzZmBrkmqqm/2XrdW1d8ARw45NkmSBmI7JUmaS4MOtzu8N/swuv/Y7TmUiCRJmiXbKUnSXBp0uN1f9abvAzYAvzHn0UiStH1spyRJc2bQu9s9a9iBSJK0vWynJElzadDhdq+aaXlV/fXchCNJ0uzZTkmS5tJs7m73NODiNv884JPALcMISpKkWbKdkiTNmUGTpH2Aw6vqOwBJXg+8t6p+Z1iBSZI0C7ZTkqQ5M9AtwIEDgXt78/cCK+c8GkmSto/tlCRpzgzak/Qu4NNJPkD3RPMXAOcOLSpJkmbHdkqSNGcGfZjsfwdOBu4EtgAnV9X/GGZgkiQNakfaqSS7JLk6yYfa/OOTXJXkxiTnJ9mtle/e5te35SuHszeSpHEbdLgdwCOBb1fV3wIbkzx+SDFJkrQ9tredegVwQ2/+TcAZVXUQXdJ1Sis/Bbizqp4InNHqSZIWoIGSpCSvA04H/qQV7Qr872EFJUnSbGxvO5XkAOBY4B/bfIAjgQtblXOA49r06jZPW35Uqy9JWmAG7Ul6AfB84HsAVXUbsOewgpIkaZa2t536G+A1wI/a/GOALVV1X5vfCCxv08tptxRvy+9q9R8iyalJ1iVZt3nz5u3bG0nSWA2aJN1bVUV3MSxJHjW8kCRJmrVZt1NJfg3YVFWf7RdPUbUGWPZgQdVZVbWqqlYtW7Zs25FLkuadQZOkC5K8DVia5HeBjwNvH15YkiTNyva0U88Anp9kA3Ae3TC7v2nbmLj76wHAbW16I7ACoC1/NPCtudwJSdL8MOjd7f78OkI1AAATpklEQVQX3fjr9wFPAv5bVf3dMAOTJGlQ29NOVdWfVNUBVbUSOBG4vKpeBHwCOL5VOwm4qE1f3OZpyy9vvVeSpAVmm89JSrIL8K9V9Wzg0uGHJEnS4IbQTp0OnJfkjcDVwNmt/GzgXUnW0/UgnTgHnyVJmoe2mSRV1f1J7k7y6Kq6axRBSZI0qLlop6rqCuCKNn0T8PQp6vwAOGEHQpUk7SS2mSQ1PwC+kORS2p2DAKrqD4YSlSRJs2M7JUmaM4MmSZe0lyRJ85HtlCRpzsyYJCU5sKpurqpzZqonSdI42E5JkoZhW3e3++DERJL3DTkWSZJmy3ZKkjTntpUk9R+c94RhBiJJ0nawnZIkzbltJUk1zbQkSfOB7ZQkac5t68YNT03ybbr/1D2iTdPmq6p+bKjRSZI0M9spSdKcmzFJqqpdRhWIJEmzZTslSRqGbQ23kyRJkqRFxSRJkiRJknqGliQlWZHkE0luSPLFJK9o5XsnuTTJje19r1aeJGcmWZ/k2iSHDys2SZIkSZrOMHuS7gNeXVUHA0cApyV5MrAGuKyqDgIua/MARwMHtdepwFuGGJskSZIkTWloSVJV3V5Vn2vT3wFuAJYDq4GJJ6OfAxzXplcD51bnSmBpkv2HFZ8kSZIkTWUk1yQlWQkcBlwF7FdVt0OXSAH7tmrLgVt6q21sZZO3dWqSdUnWbd68eZhhS5IkSVqEhp4kJdkDeB/wyqr69kxVpyjb6sGAVXVWVa2qqlXLli2bqzAlSZIkCRhykpRkV7oE6d1V9f5WfMfEMLr2vqmVbwRW9FY/ALhtmPFJkiRJ0mQzPkx2RyQJcDZwQ1X9dW/RxcBJwNr2flGv/OVJzgN+DrhrYlietKNWrrlk3CFIkiRpJzG0JAl4BvBi4AtJrmllr6VLji5IcgpwM3BCW/Zh4BhgPXA3cPIQY5MkSZKkKQ0tSaqqf2Pq64wAjpqifgGnDSseSZIkSRrESO5uJ0mSJEk7C5MkSZIkSeoZ5jVJkiRJGrLZ3pxow9pjhxSJtHDYkyRJkiRJPSZJkiRJktRjkiRJkiRJPSZJkiRJktRjkiRJkiRJPSZJkiRJktRjkiRJWpSSPDzJp5N8PskXk7yhlT8+yVVJbkxyfpLdWvnubX59W75ynPFLkobHJEmStFjdAxxZVU8FDgWem+QI4E3AGVV1EHAncEqrfwpwZ1U9ETij1ZMkLUAmSZKkRak6322zu7ZXAUcCF7byc4Dj2vTqNk9bflSSjChcSdIImSRJkhatJLskuQbYBFwKfAXYUlX3tSobgeVtejlwC0BbfhfwmCm2eWqSdUnWbd68edi7IEkaApMkSdKiVVX3V9WhwAHA04GDp6rW3qfqNaqtCqrOqqpVVbVq2bJlcxesJGlkTJIkSYteVW0BrgCOAJYmWdIWHQDc1qY3AisA2vJHA98abaSSpFEwSZIkLUpJliVZ2qYfATwbuAH4BHB8q3YScFGbvrjN05ZfXlVb9SRJknZ+S7ZdRZKkBWl/4Jwku9D90/CCqvpQkuuB85K8EbgaOLvVPxt4V5L1dD1IJ44jaEnS8JkkSZIWpaq6FjhsivKb6K5Pmlz+A+CEEYQmSRozh9tJkiRJUo9JkiRJkiT1mCRJkiRJUo/XJGleW7nmknGHIEmSpEXGniRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqQekyRJkiRJ6jFJkiRJkqSeJeMOQJIkSaOzcs0ls6q/Ye2xQ4pEmr/sSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSdKilGRFkk8kuSHJF5O8opXvneTSJDe2971aeZKcmWR9kmuTHD7ePZAkDYtJkiRpsboPeHVVHQwcAZyW5MnAGuCyqjoIuKzNAxwNHNRepwJvGX3IkqRRMEmSJC1KVXV7VX2uTX8HuAFYDqwGzmnVzgGOa9OrgXOrcyWwNMn+Iw5bkjQCJkmSpEUvyUrgMOAqYL+quh26RArYt1VbDtzSW21jK5u8rVOTrEuybvPmzcMMW5I0JCZJkqRFLckewPuAV1bVt2eqOkVZbVVQdVZVraqqVcuWLZurMCVJI2SSJElatJLsSpcgvbuq3t+K75gYRtfeN7XyjcCK3uoHALeNKlZJ0uiYJEmSFqUkAc4Gbqiqv+4tuhg4qU2fBFzUK39Ju8vdEcBdE8PyJEkLy5JxByBJ0pg8A3gx8IUk17Sy1wJrgQuSnALcDJzQln0YOAZYD9wNnDzacCVJo2KSJElalKrq35j6OiOAo6aoX8BpQw1KkjQvDG24XZJ3JNmU5LpemQ/okyRJkjSvDbMn6Z3A3wPn9somHtC3NsmaNn86D31A38/RPaDv54YYm4Zs5ZpLxh2CJEmStF2G1pNUVZ8EvjWp2Af0SZIkSZrXRn13ux16QJ8kSZIkDdt8uQX4QA/oA59kLkmSJGm4Rp0k7fAD+nySuSRJkqRhGnWS5AP6JEmSJM1rQ7u7XZJ/Bp4J7JNkI/A6fECfJEmSpHluaElSVb1wmkU+oE+SJEnSvDVfbtwgSZIkSfOCSZIkSZIk9QxtuJ0kSZJ2fivXXDKr+hvWHjukSKTRsSdJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJkiRJknpMkiRJkiSpxyRJkrQoJXlHkk1JruuV7Z3k0iQ3tve9WnmSnJlkfZJrkxw+vsglScNmkiRJWqzeCTx3Utka4LKqOgi4rM0DHA0c1F6nAm8ZUYySpDEwSZIkLUpV9UngW5OKVwPntOlzgON65edW50pgaZL9RxOpJGnUTJIkSXrQflV1O0B737eVLwdu6dXb2Mq2kuTUJOuSrNu8efNQg5UkDYdJkiRJ25YpymqqilV1VlWtqqpVy5YtG3JYkqRhMEmSJOlBd0wMo2vvm1r5RmBFr94BwG0jjk2SNCImSZIkPehi4KQ2fRJwUa/8Je0ud0cAd00My5MkLTxLxh2AJEnjkOSfgWcC+yTZCLwOWAtckOQU4GbghFb9w8AxwHrgbuDkkQcsSRoZkyRJ0qJUVS+cZtFRU9Qt4LThRiRJmi8cbidJkiRJPSZJkiRJktRjkiRJkiRJPV6TJEmSpDmzcs0lA9fdsPbYIUYibT+TJD3EbH6xSZIkSQuRw+0kSZIkqcckSZIkSZJ6TJIkSZIkqcckSZIkSZJ6TJIkSZIkqcckSZIkSZJ6vAW4JEmSxmK2jx7xuUoaFXuSJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnHJEmSJEmSekySJEmSJKnH5yQtELN9zoAkSdLOxucqaVTsSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHpMkSZIkSeoxSZIkSZKkHp+TJEmSpAXJ5yppe82rnqQkz03y5STrk6wZdzySJE1mWyVJC9+86UlKsgvwD8BzgI3AZ5JcXFXXjzcySZI6tlXSwmbPkybMmyQJeDqwvqpuAkhyHrAasOGRJM0XtlWSHjDbpGo2TMDGaz4lScuBW3rzG4Gfm1wpyanAqW32niTXjSC2ndE+wDfGHcQ85bGZnsdmeh6b6T1p3AGM0DbbqkXcTi22n5HFtL/u6xjkTUP/iHmzryMw63ZqPiVJmaKstiqoOgs4CyDJuqpaNezAdkYem+l5bKbnsZmex2Z6SdaNO4YR2mZbtVjbqcW0r7C49td9XZgW277Odp35dOOGjcCK3vwBwG1jikWSpKnYVknSIjCfkqTPAAcleXyS3YATgYvHHJMkSX22VZK0CMyb4XZVdV+SlwP/CuwCvKOqvriN1c4afmQ7LY/N9Dw20/PYTM9jM71Fc2y2o61aNMeGxbWvsLj2131dmNzXGaRqq8t+JEmSJGnRmk/D7SRJkiRp7EySJEmSJKlnp0iSkjw3yZeTrE+yZorluyc5vy2/KsnK0Uc5HgMcm1cluT7JtUkuS/K4ccQ5Dts6Nr16xyepJIviNpgw2LFJ8hvtu/PFJO8ZdYzjMsDP1IFJPpHk6vZzdcw44hy1JO9Ismm6Z/6kc2Y7btcmOXzUMc43g/4OWgiSbEjyhSTXLLRbwk/13U+yd5JLk9zY3vcaZ4xzZZp9fX2SW9u5vWah/M5LsqL9Lr+htXOvaOUL7tzOsK8L9dw+PMmnk3y+7e8bWvnjW55wY8sbdptxQ1U1r190F8Z+BXgCsBvweeDJk+q8DHhrmz4ROH/ccc+jY/Ms4JFt+r96bLaqtyfwSeBKYNW4454vxwY4CLga2KvN7zvuuOfRsTkL+K9t+snAhnHHPaJj80vA4cB10yw/BvgI3XOEjgCuGnfMYz5eA/0OWigvYAOwz7jjGNK+bfXdB/4CWNOm1wBvGnecQ9zX1wN/NO7YhrCv+wOHt+k9gf9sv9MX3LmdYV8X6rkNsEeb3hW4qrVLFwAntvK3TrTl0712hp6kpwPrq+qmqroXOA9YPanOauCcNn0hcFSSqR74t9Bs89hU1Seq6u42eyXdMz0Wg0G+NwB/TvcL8QejDG7MBjk2vwv8Q1XdCVBVm0Yc47gMcmwK+LE2/WgWyTNyquqTwLdmqLIaOLc6VwJLk+w/mujmpUF/B2mem+a73/+74xzguJEGNSQD/JwvGFV1e1V9rk1/B7gBWM4CPLcz7OuC1Nqh77bZXdurgCPp8gQY4NzuDEnScuCW3vxGtj6xD9SpqvuAu4DHjCS68Rrk2PSdQvef3sVgm8cmyWHAiqr60CgDmwcG+d78JPCTSf49yZVJnjuy6MZrkGPzeuC3k2wEPgz8/mhCm/dm+/tooVtsx6OAjyX5bJJTxx3MCOxXVbdD9wcosO+Y4xm2l7dhtO9YCMPPJmuXaRxG1+OwoM/tpH2FBXpuk+yS5BpgE3ApXc/+lpYnwAC/k3eGJGmqHqHJ9y0fpM5CNPB+J/ltYBXwl0ONaP6Y8dgkeRhwBvDqkUU0fwzyvVlCN+TumcALgX9MsnTIcc0HgxybFwLvrKoD6IaYvat9nxa7xfp7eDqL7Xg8o6oOB44GTkvyS+MOSHPmLcBPAIcCtwN/Nd5w5laSPYD3Aa+sqm+PO55hmmJfF+y5rar7q+pQuhFUTwcOnqraTNvYGRr2jcCK3vwBbD285YE6SZbQDYFZDN3Fgxwbkjwb+P+A51fVPSOKbdy2dWz2BJ4CXJFkA91Y1YsXyc0bBv2ZuqiqflhVXwW+TJc0LXSDHJtT6MY1U1WfAh4O7DOS6Oa3gX4fLSKL6nhU1W3tfRPwAbo/ShayOyaGk7b3BTskuaruaH9w/gh4Owvo3CbZlS5peHdVvb8VL8hzO9W+LuRzO6GqtgBX0P2dt7TlCTDA7+SdIUn6DHBQuyPFbnQ3Zrh4Up2LgZPa9PHA5dWuylrgtnls2pCyt9ElSAviB31AMx6bqrqrqvapqpVVtZLueq3nV9WCuivTNAb5mfog3U0/SLIP3fC7m0Ya5XgMcmxuBo4CSHIwXZK0eaRRzk8XAy9pd7k7ArhrYsjKIjXId2lBSPKoJHtOTAO/Akx5F8QFpP93x0nARWOMZagmXVv4AhbIuW3Xrp8N3FBVf91btODO7XT7uoDP7bKJ0S9JHgE8m+46rE/Q5QkwwLldMtPC+aCq7kvycuBf6e4W9I6q+mKSPwPWVdXFdCf+XUnW0/UgnTi+iEdnwGPzl8AewHvbvSxurqrnjy3oERnw2CxKAx6bfwV+Jcn1wP3AH1fVN8cX9WgMeGxeDbw9yR/SddW/dDH8UybJP9MNv9ynXY/1OrqLYamqt9Jdn3UMsB64Gzh5PJHOD9N9l8Yc1rDsB3ygtTFLgPdU1UfHG9Lcmea7vxa4IMkpdP84OWF8Ec6dafb1mUkOpft9twH4vbEFOLeeAbwY+EK7dgXgtSzMczvdvr5wgZ7b/YFzkuxC1yF0QVV9qP1Nc16SN9LdwffsmTaSRdC2S5IkSdLAdobhdpIkSZI0MiZJkiRJktRjkiRJkiRJPSZJkiRJktRjkiRJkiRJPSZJkiRJktRjkiRJkiRJPSZJ0hgk+USS57TpNyY5c9wxSZI0IcnvJakkB/fKbkiycnxRSaOzZNwBSIvU64A/S7IvcBjw/DHHI0lS388A1wDHAjck2R3YD/jaWKOSRsSeJGkMquqTQIBXASdW1f1jDkmSpL5DgLV0SRLATwM3VFWNLyRpdEySpDFIcgiwP3BPVX1n3PFIkjTJk4GLgX2TPJouafrCeEOSRsckSRqxJPsD7wZWA99L8qtjDkmSpAckWQF8s6q+D1wK/Crd8LtrxxqYNEImSdIIJXkk8H7g1VV1A/DnwOvHGpQkSQ/1MzzYa/RhuiF39iRpUYlDSyVJkjQhyRpg96p6Q7thww3AI4Gfqqot441OGg17kiRJktT3QK9RVd3Tpu81QdJiYk+SJEmSJPXYkyRJkiRJPSZJkiRJktRjkiRJkiRJPSZJkiRJktRjkiRJkiRJPSZJkiRJktRjkiRJkiRJPf8/Muzrvwkebf4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "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", "© Blitzstein, Joseph K.; Hwang, Jessica. Introduction to Probability (Chapman & Hall/CRC Texts in Statistical Science)." ] } ], "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.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }