Metropolis-Hastings_exponential.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two Different Ways to Produce Exponentially Distributed Random Variables\n",
    "\n",
    "We produce exponentially distributed random variables using the Metropolis-Hastings algorithm. Then, we compare the results to the traditional inverse transform method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import RandomState\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code for the Metropolis-Hastings algorithm is inspired by Antti Honkela's lecture material for the course [Computational Statistics I](https://www.cs.helsinki.fi/u/ahonkela/teaching/compstats1/book/markov-chain-monte-carlo-basics.html#implementation-of-the-metropolis-hastings-algorithm)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def metropolis_algorithm(rng, chain_start, n, logtarget, candidate_generating_density):\n",
    "    x = chain_start\n",
    "    chain = np.zeros(n)\n",
    "    accepted = 0\n",
    "    for i in range(n):\n",
    "        candidate = candidate_generating_density(x)\n",
    "        if candidate >= 0:\n",
    "            if np.log(rng.uniform()) < logtarget(candidate) - logtarget(x):\n",
    "                x = candidate\n",
    "                accepted += 1\n",
    "        \n",
    "        chain[i] = x\n",
    "        \n",
    "    print(f\"Acceptance rate {accepted/n}\")\n",
    "    return chain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we set some initial parameters for our experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = RandomState(1729)\n",
    "n = 1_000_000\n",
    "param = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We initialize the logarithmic density of the exponential distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_exponential = lambda x: np.log(param) - param * x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use a normal distribution to explore the state space. The variance of the the standard normal distribution is changed by multiplying it with a suitable constant. Thus, the exploration of the state space can be made more efficient. The constant $ 0.6 $ was determined by trial and error. See the [lecture material of Computational Statistics I](https://www.cs.helsinki.fi/u/ahonkela/teaching/compstats1/book/markov-chain-monte-carlo-basics.html#mcmc-sampling-efficiency) for more detail on this step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "candidate_generating_density = lambda x: x + 0.60 * rng.normal()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can use the Metropolis-Hastings algorithm to generate the observations. We start the Markov chain from the mean of the distribution. Then, we remove the first half of the observations as the burn-in period. Additionally, we take only every tenth value of the remaining observations in order to reduce auto-correlation present in the sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acceptance rate 0.47233\n",
      "(50000,)\n"
     ]
    }
   ],
   "source": [
    "chain = metropolis_algorithm(rng, (1/param), n, log_exponential, candidate_generating_density)\n",
    "chain = chain[chain.size // 2:]\n",
    "chain = chain[::10]\n",
    "print(chain.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the inverse transform method to produce another set of exponentially distributed random variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(50000,)\n"
     ]
    }
   ],
   "source": [
    "rng_2 = RandomState(88)\n",
    "exponential_rvs = -(1 / param) * np.log(1 - rng_2.rand(chain.size))\n",
    "print(exponential_rvs.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the results and compare the two methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3AAAAGrCAYAAACBjHUSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd5ydZZ3//9dnJr0nZICEkgBSpRvALghSRfSHBdBFXITdtay6NvTr2nVt21wLsorYUMCKBkUUVwRBCC2UUEIa6UMmjSSkzfX7477P5CSZyUySM3PfZ87r+XjMY859n/uc8znjw3x43/d1XXeklJAkSZIklV9T0QVIkiRJknrGACdJkiRJdcIAJ0mSJEl1wgAnSZIkSXXCACdJkiRJdcIAJ0mSJEl1wgAn9aKIuCQibq/afjYiDiyypu5ExMsi4vGi65Akqa9FxGcj4pmIWFx0LVJXDHAqtYiYExEbImL8Nvvvj4gUEZO7ef3JETG/N2vcGSmlESmlWTv7uvzvcNo2+7YKh7sq/zs+r6rGv6SUDt3d95Uk9b3O+kU9iYhrIuKzBX32/sD7gSNSSnsXUYPUEwY41YPZwIWVjYg4ChhWqzePiAG1ei9JkhpNX/bRXv6s/YFlKaWlO/tC/1tCfckAp3rwA+Diqu23At+vbETE4Ij4SkTMi4glEXFlRAyNiOHAb4GJ+dDFZyNiYkR8MiJ+GhE/jIhVwCX5e/xXRCzMf/4rIgbn739yRMyPiI/mwyrmRMSbqz5/dER8PyJaI2JuRHwsIjr9/1b11a6IODsiHo2I1RGxICI+sDt/pIi4IiKeyt/v0Yh4XdVzz4uIP0fEyvw7XJfvvy0/5MH87/Omba9a5t/3AxExPX/9dRExpOr5D0XEovzv9vbe/I6SpJ6rjNTIe+TyiJgdEWflz70pIqZtc/z7IuLG/HGnvTV/rtIXP5wPNfxuRIyPiN9ExIqIaIuIv1R6Yd57f5b3ydkR8c9d1Hs58GbgQ3lP+nW+f07+WdOBNRExoJue1+X3rnp+Vv7a2RHx5siuWt7Clv9muCY/9jUR8Uj+vf4vIg6vep/O6poTER/Me+aaiPhOROwVEb/NP+8PETF2d/+3VWMzwKke3AWMiojDI6IZuAD4YdXzXwAOAY4FngfsA3w8pbQGOAtYmA9dHJFSWpi/5jzgp8AY4EfA/wNemL/HMcCJwMeqPmNvYHz+3m8FroqIyjDD/wFGAwcCryALm2/rwff6DvAPKaWRwJHArT37c3TpKeBleS2fAn4YERPy5z4D/B4YC+yb10xK6eX588fkf5/runjvNwJnAgcARwOXAETEmcC/AKeR/e1P3uZ1tf6OkqSdcxLwOFkP+xLwnYgI4NfAoRFxcNWxFwHX5o877a1Vx+4NjAMmAZeTDT2cD7QAewEfBVIe4n4NPJi/x6nAeyPijG0LTSldRdaTv5T3pHOrnr4QOAcYk1LaxI57XpffO7KTu18Fzsp704uBB1JKf2Dr/2a4JCIOAX4MvDf/XjcBv46IQTuoC+B84FX53+9cspPJH83fownoNMBKPWWAU72oXIV7FTADWJDvD7LG8b6UUltKaTXwebKQtyN3ppR+mVJqTymtIzvj9+mU0tKUUitZM/i7bV7zryml9SmlPwNTgTdWBcqPpJRWp5TmAP/eyWs7sxE4IiJGpZSWp5Tu6+b4X+ZnAFdExArgG9VPppRuSCktzL/TdcCTZEG08lmTgIkppedSSjs7d+6r+Xu3kTXiY/P9bwS+m1J6JKW0Fvjkbn5HSVJtzU0p/W9KaTPwPWACsFf+b/avyKco5EHuMODGPOB111vbgU/kfXEd2b/3E4BJKaWN+XzqBJwAtKSUPp1S2pDPA/9fuu/T2/pqSunp/LO663ldfu+q2o+MiKEppUUppUe6+Mw3AVNTSreklDYCXwGGkoW+TuvK/U9KaUlKaQHwF+BvKaX7U0rPAb8AjtvJ7y5txQCnevEDsjODl1A1fJLsbNYw4N6qYPO7fP+OPL3N9kRgbtX23HxfxfL8it62z48HBnby2n26+XzIztCdDcyNbHjjiwDyYRaVIZ9vrjr+tSmlMZUf4B3VbxYRF0fEA1V/hyPz+gA+RBZ2786Hgvx9D+qrVr0a11pgRP54Ilv/Lbf9u3b6HSVJfabj3+88tMGWf8OvZcsc84uAX+bH9KS3tuaBpOLLwEzg9/nwxCvy/ZPIhiVWn4D8KFvCVE9t1V+66Xldfu+8l78J+EdgUURMjYjDuvjMrf7bIKXUntdR3eO37XsAS6oer+tkewTSbnDCpepCSmluRMwmCwOXVj31DNk/hs/Pz3Rt99Ku3nKb7YVkTaZyFm7/fF/F2IgYXhXi9gcezj+/cnXr0arnOqtl2+90D3BeRAwE3gVcD+yXUjprx6/cXkRMIjujeSrZ1cXNEfEAWWgjpbQYuCw/9qXAHyLitpTSzJ39rG0sIhuSWbFf9ZNdfcfd/ExJUm3cArRExLFkQe59+f7ueits00fzq3TvB94fEUcCt0bEPWQBZ3ZK6eBO3qPb9+1sf3c9r9sPSOlm4ObI5vR9Nn+vl3Vy6ELgqKrPDbIeVv036apeqdd4BU715FLgldtcCWsn+4f3PyNiT4CI2KdqbP0SYI+IGN3Ne/8Y+FhEtER2y4KPs/U8O4BPRcSgiHgZ8GrghnxoxvXA5yJiZN5U/qWT124lf583R8TofFjGqvy77KrhZE2kNX//t5Gdjax83hsiohK0lufHVj5vCdn8vV1xPfC2fH7iMOBfqz6z1t9RklRD+b/NN5BdPRtHFugqV5p21Fu3ExGvjmzBrABWApvJ/s2/G1idL/YxNCKaI+LIiDihi7fqSU/aYc/bkXxBkfPyuXDrgWfpujddD5wTEafmJyLfn7/mrz35LKm3GOBUN1JKT6WUpnXy1IfJhm3cFdmqkn8ADs1f8xhZOJuVD7OY2MnrITsDNw2YDjwE3Jfvq1hMFnwWkk2w/sf8vQHeDawBZgG3kw1JuboHX+nvgDl5zf9INg9vl6SUHiWbe3cnWfM7Crij6pATgL9FxLPAjcB7qu5H90nge/nf5407+bm/JZsM/ify/w3yp9bnv2v2HSVJveJasoWobqhahAN20Fu7cHB+zLNkvegbKaU/5Sc6X002d3o22dW9b5MtPtKZ75DNnV4REb/s7IAe9LwdaSI70boQaCNbfOyfuvicx4G3kC389QzZgiTnppQ29PCzpF4R2fxSSV2JiJOBH6aU9u3u2EYX2fLKDwODt/kPAUmSJNWAV+Ak7ZaIeF1k9wsaC3wR+LXhTZIkqXcY4CTtrn8AlpLdk2czXQxFkSRJ0u5zCKUkSZIk1QmvwEmSJElSnSjlfeDGjx+fJk+eXHQZkqRedu+99z6TUmrp/kiB/VGSGklXPbKUAW7y5MlMm9bZavGSpP4kIuYWXUM9sT9KUuPoqkc6hFKSJEmS6oQBTpIkSZLqhAFOkiRJkuqEAU6SJEmS6oQBTpIkSZLqhAFOkiRJkuqEAU6SJEmS6oQBTpIkSZLqhAFOkiRJkuqEAU6SJEmS6oQBTpIkSZLqhAFOkiRJkuqEAU6SJEmS6kS3AS4i9ouIP0XEoxHxSES8p5NjIiK+GhEzI2J6RBxf9dxbI+LJ/Oettf4CkiQVxR4pSeprA3pwzCbg/Sml+yJiJHBvRNySUnq06pizgIPzn5OAbwInRcQ44BPAFCDlr70xpbS8pt9CkqRi2CMlSX2q2ytwKaVFKaX78sergRnAPtscdh7w/ZS5CxgTEROAM4BbUkpteUO6BTizpt9gWwsXwrvfDR/5SK9+jCRJddcjr70WLr0U/vKXXv0YSVLv2ak5cBExGTgO+Ns2T+0DPF21PT/f19X+zt778oiYFhHTWltbd6asrT33HHzta1mTkiSpj/RWj6xZfwS4/Xa4+mp48MHdex9JUmF6HOAiYgTwM+C9KaVVtS4kpXRVSmlKSmlKS0vLrr9R5bW72+QkSeqh3uyRNeuPYI+UpH6gRwEuIgaSNaYfpZR+3skhC4D9qrb3zfd1tb/3jBgBgwfDunWwZk2vfpQkSXXVIw1wklT3erIKZQDfAWaklP6ji8NuBC7OV9p6IbAypbQIuBk4PSLGRsRY4PR8X++JgPHjs8fPPNOrHyVJamx11yPtj5JU93qyCuVLgL8DHoqIB/J9HwX2B0gpXQncBJwNzATWAm/Ln2uLiM8A9+Sv+3RKqa125XehpQUWLMjOME6a1OsfJ0lqWPXVI70CJ0l1r9sAl1K6HYhujknAO7t47mrg6l2qblfZoCRJfaDueqT9UZLq3k6tQlk3bFCSJG3P/ihJda9/BjjH+EuStL099sh+L1sG7e3F1iJJ2iX9M8B5hlGSpO0NGgSjR8PmzbBiRdHVSJJ2gQFOkqRGYo+UpLrWvwOcQyglSdqaAU6S6lr/DHCVOXA2J0mStuY8cUmqa/0zwHl2UZKkztkjJamuGeAkSWok9khJqmv9M8CNHQtNTdkKWxs3Fl2NJEnl4TQDSapr/TPANTfDuHHZ42XLiq1FkqQycaEvSapr/TPAgUNEJEnqjP1RkuqaAU6SpEZif5Skutb/A5xDRCRJ2sI5cJJU1/pvgLNBSZK0PU9wSlJd678BziEikiRtb/hwGDIE1q2DNWuKrkaStJP6f4DzDKMkSVtEeJJTkupY/w9wNidJkrbmNANJqlv9NsC95ZdPAfDXux4ruBJJkspj8hVTuW159viSL91UbDGSpJ3WbwNc27DRAIxbu7LgSiRJKpdllR65zh4pSfWm3wa4ZcNGATBu3aqCK5EkqVw8ySlJ9avfBrjlQ7PmNHbdKkip4GokSSqPtqHZSc491nqSU5LqTb8NcBsGDGT1oKEMbN8MK1YUXY4kSaVRuQI31lEqklR3+m2Agy0NylW2JEnaom1Y5QqcQyglqd707wCXD6P0XnCSJG2xzDlwklS3+nWAqyxk4hU4SZK2WJ7PgXMIpSTVn34d4CoLmRjgJEnaonIFziGUklR/+nWA67gC5xBKSZI6rBwygs3RxOj1a2DjxqLLkSTthH4d4FzERJKk7aVoYvnQkdmGJzklqa4Y4CRJakAu9CVJ9al/B7ihLmIiSVJn2lzoS5LqUv8OcMM8uyhJUmeWOUpFkupSvw5wNidJkjq33FEqklSX+nWAcwilJEmdW+YoFUmqS/06wK0ZNJT1zQNg7drsR5IkAS70JUn1ql8HOCJcZUuSpE44SkWS6lP/DnB4hlGSpM7YHyWpPvX7AOdCJpIkba/jNgKOUJGkujKguwMi4mrg1cDSlNKRnTz/QeDNVe93ONCSUmqLiDnAamAzsCmlNKVWhfdUxypbNihJUo3Vc49cNtQTnJJUj3pyBe4a4MyunkwpfTmldGxK6VjgI8CfU0ptVYeckj/f5+ENvFGpJKlXXUOd9sjl1Vfg2tv7+uMlSbuo2wCXUroNaOvuuNyFwI93q6IacwilJKm31HOP3Ng8kFWDhsHmzbBiRdHlSJJ6qGZz4CJiGNlZyJ9V7U7A7yPi3oi4vJvXXx4R0yJiWmsNw5aTtCVJRdudHtlb/RGqeqTTDCSpbtRyEZNzgTu2GRry0pTS8cBZwDsj4uVdvTildFVKaUpKaUpLS0vNilrmbQQkScXb5R7ZW/0RnGYgSfWolgHuArYZGpJSWpD/Xgr8Ajixhp/XI8ttTpKk4pWyRzrNQJLqT00CXESMBl4B/Kpq3/CIGFl5DJwOPFyLz9sZrrIlSSpSmXvkcm/mLUl1pye3EfgxcDIwPiLmA58ABgKklK7MD3sd8PuU0pqql+4F/CIiKp9zbUrpd7UrvWccHiJJ6i313yOdZiBJ9abbAJdSurAHx1xDtpRy9b5ZwDG7WlitrBwyApqaYPly2LQJBnT7lSVJ6pF675GOUpGk+lPLOXCl1N7UDOPGZRvLlhVbjCRJJeI8cUmqP/0+wAFQWbXLBiVJUgcXMZGk+tNYAc4x/pIkdWirLGJif5SkutFYAc4zjJIkdfAKnCTVn8YIcOPHZ79tUJIkdfA2ApJUfxojwHkFTpKk7awZNJT1zQNh3ToOe//Pii5HktQDjRXgHOMvSdIWER3DKPdYu6rgYiRJPdFYAc4rcJIkbaVyM+9xa1cWXIkkqScaI8A5B06SpE5VVqLcwwAnSXWhIQLc2dc+BsCMh2YVXIkkSeXSlt/Me+w6h1BKUj1oiAC3LG9O42xOkiRtpW2oQyglqZ40RIBbnjensWtXQUoFVyNJUnl0LGKyzgAnSfWgIQLchgEDWT1oKIPaN8FKG5QkSRXLK6NUXIVSkupCQwQ42LLKlrcSkCRpi2WVIZROM5CkutAwAa4yRITFi4stRJKkElk2POuPLc8uL7gSSVJPNEyAWzJij+zBokXFFiJJUolU+mPLmraCK5Ek9UTDBLilI8ZmDwxwkiR1WDo8648ta5a70Jck1YHGCXDDx2UPDHCSJHVYP3AwqwYPZ/DmTdDmVThJKrvGCXAj8gC3cGGxhUiSVDJL7JGSVDcaL8B5BU6SpK04zUCS6kcDBTibkyRJnXGagSTVj8YJcMMdHiJJUmccQilJ9aNhAlzbsFFsbGrOJmivX190OZIklUar0wwkqW40TIBL0URrvlSyN/OWJGmLyq0EDHCSVH4NE+Cgah6cQ0QkSeqwZGR2M2/7oySVX4MFuLxBeYZRkqQOXoGTpPrRWAHOBiVJ0na2utVOSsUWI0naocYKcK6yJUnSdtYOGsrqQUPhuedgxYqiy5Ek7UBDBbglrrIlSVKnltojJakuNFSAszlJktS5VqcZSFJdaKwAN9xVKCVJ6sySEa5EKUn1oLECnKtQSpLUqY5b7dgjJanUGirALRs2CpqaoLUVNm4suhxJkkpj6XCnGUhSPWioANfe1Ax77ZVtLF5cbDGSJJXIkpGu1CxJ9aChAhwAEydmvz3DKElSh1avwElSXWi8ADdhQvbbBiVJUgfnwElSfeg2wEXE1RGxNCIe7uL5kyNiZUQ8kP98vOq5MyPi8YiYGRFX1LLwXVYJcA4RkSTtpv7UI7dahTKlYouRJHWpJ1fgrgHO7OaYv6SUjs1/Pg0QEc3A14GzgCOACyPiiN0ptia8AidJqp1r6Cc98tlBQ2HYMFi7FlavLrIUSdIOdBvgUkq3AW278N4nAjNTSrNSShuAnwDn7cL71JZz4CRJNdKvemSEJzklqQ7Uag7ciyLiwYj4bUQ8P9+3D/B01THz832diojLI2JaRExrbW2tUVmdcAilJKlv7VaP7LP+CFtOctojJam0ahHg7gMmpZSOAf4H+OWuvElK6aqU0pSU0pSWlpYalNUFzy5KkvrObvfIPuuPYI+UpDqw2wEupbQqpfRs/vgmYGBEjAcWAPtVHbpvvq9QJ13zKABLn5hTbCGSpH6v3nqkAU6Sym+3A1xE7B0RkT8+MX/PZcA9wMERcUBEDAIuAG7c3c/bXc8MH0s7wfg1K2DTpqLLkST1Y/XWIx1CKUnlN6C7AyLix8DJwPiImA98AhgIkFK6Eng98E8RsQlYB1yQUkrApoh4F3Az0AxcnVJ6pFe+xU7Y3NTMsmGjaVm7ApYu3dKsJEnaSf2tR3oFTpLKr9sAl1K6sJvnvwZ8rYvnbgJu2rXSes/SEeOyALdokQFOkrTL+l2PNMBJUunVahXKurJ0xNjsgUNEJEnqcNpPngTgqelPFlyJJKkrjRngho/LHniGUZKkDktHZP1xz2d35dZ2kqS+0JABbskIA5wkSdtaNXg465sHMnLDOnj22aLLkSR1oiEDXOUMo0MoJUmqEuFJTkkquYYMcK2VOXA2J0mStrLUACdJpdaQAW7JiD2yBzYnSZK2snS4JzklqcwaMsC5CqUkSZ1bMjI/yWmPlKRSasgA11o5u7hkCWzeXGwxkiSVSKtX4CSp1BoywG1sHsiyoaOy8PbMM0WXI0lSaTgHTpLKrSEDHLgSpSRJnVlif5SkUmvYAOcQEUmStucVOEkqt4YNcDYoSZK25yqUklRuDRvgHCIiSdL2lg8dxYamAbBiBaxbV3Q5kqRtNGyA8wqcJEmdiNhyux17pCSVTuMGOIeISJLUqdbhnuSUpLJq2AC3ZIQ3KpUkqTNLRjrNQJLKqmEDnMNDJEnq3FKvwElSaTVsgGutzIFbvBja24stRpKkEvEkpySVV8MGuPUDBsHYsbBxIyxbVnQ5kiSVhis1S1J5NWyAA2DChOy3ZxglSergIiaSVF4GOLBBSZJUxVvtSFJ5NXaAmzgx++0QEUmSOjiEUpLKq7EDXOUKnA1KkqQObcNGwYAB0NYGzz1XdDmSpCqNHeD23Tf7/fTTxdYhSVKJpGjaMkpl/vxii5EkbaWxA9ykSdnvuXOLrUOSpLKxR0pSKRngwOYkSdK27JGSVEoNHeCO/u5jAKydOYvJH/5NwdVIklQilQA3b16xdUiSttLQAW7VkBGsGjSMYRvXM+a51UWXI0lSeXgFTpJKqaEDHMCC0XsCsM/KpQVXIklSiRjgJKmUDHCjWgDYd5UBTpKkilfeMBuAuffPKLgSSVI1A1zHFbjWgiuRJKk8Kic4J6x6BtrbC65GklTR8AFuYd6g9vEKnCRJHdYPHMwzw0YzqH0TLFpUdDmSpFzDB7gFo/IrcAY4SZK2UumRzoOTpPIwwI1yERNJkjpTGUZpgJOk8mj4ADe/MgdulXPgJEmqVpknboCTpPJo+AD3zPAxrG8ewLh1q2DNmqLLkSSpNBxCKUnl022Ai4irI2JpRDzcxfNvjojpEfFQRPw1Io6pem5Ovv+BiJhWy8JrJUVTx0ImNihJ0s7o7z3SK3CSVD49uQJ3DXDmDp6fDbwipXQU8Bngqm2ePyWldGxKacquldj7PMMoSdpF19CPe6T9UZLKZ0B3B6SUbouIyTt4/q9Vm3cB++5+WX3LBiVJ2hX9vUfOr74ClxJEFFuQJKnmc+AuBX5btZ2A30fEvRFx+Y5eGBGXR8S0iJjW2tq3C4o4RESS1Ad2qUcW2R9XDR7O6kFDszniy5f36WdLkjpXswAXEaeQNacPV+1+aUrpeOAs4J0R8fKuXp9SuiqlNCWlNKWlpaVWZfWIV+AkSb1pd3pkkf2RCHukJJVMTQJcRBwNfBs4L6W0rLI/pbQg/70U+AVwYi0+r9YWjHYRE0lS76j/HmmAk6Qy2e0AFxH7Az8H/i6l9ETV/uERMbLyGDgd6HSVrqLN9+yiJKkX9Ice6RU4SSqXbhcxiYgfAycD4yNiPvAJYCBASulK4OPAHsA3IpvcvClfTWsv4Bf5vgHAtSml3/XCd9hti0eOZ3M00bxwIWzYAIMGFV2SJKkONEKPdJSKJJVLT1ahvLCb598OvL2T/bOAY7Z/Rflsah7AkhHjmLj6GZg/Hw48sOiSJEl1oBF6pFfgJKlcar0KZd2yQUmStD37oySViwEu1zFEZN68YguRJKlE5ruIiSSVigEut3CUY/wlSdrWM8PHZHPDn3kmux+cJKlQBricQ0QkSdpeiibmDNsDgFPf84OCq5EkGeByBjhJkjpXmWaw76qlBVciSTLA5RzjL0lS5yonOfcxwElS4QxwuY4rcPPmQXt7scVIklQiHQFupQFOkopmgMutGzQE9tgju5H3kiVFlyNJUmksyEepTFzVWnAlkiQDXJWHBowB4HUfurbgSiRJKg+HUEpSeRjgqlTOMNqgJEnaojJPfJ+VXoGTpKIZ4Kp4hlGSpO0tHrkH7QR7P7sMNm4suhxJamgGuCpbJml7hlGSpIqNzQNZMmIczakdFiwouhxJamgGuCqV+9x4BU6SpK0t8HY7klQKBrgq810mWZKkTnXcbscAJ0mFMsBV2WoRk5QKrkaSpPKojFIxwElSsQxwVVYMGcmagUMYuWEdrFhRdDmSJJWGV+AkqRwMcNUibFCSJHXC/ihJ5WCA28bCUQ4RkSRpWwvsj5JUCga4bTjGX5Kk7XWsQjlvHrS3F1uMJDUwA9w2OoaIzJtXbCGSJJXI2kFDWT5kJKxfD0tdrVmSimKA28Z873MjSVKntroKJ0kqhAFuGx1X4ObMKbQOSZLKpmMenD1SkgpjgNvGvDETsgczZ3ovOEmSqsyt9Minniq2EElqYAa4bbQOH8Ozg4Zm94FbtqzociRJKo054yZmD554othCJKmBGeC2FcHssXmDevLJYmuRJKlEZo/dJ3tgf5SkwhjgOjF7XN6gPMMoSVKHWV6Bk6TCGeA60XGG0QYlSVKHJSP2gGHDoLUVli8vuhxJakgGuE54hlGSpE5E8MiIvQE47z3fLbgYSWpMBrhOzBlrgJMkqTOVeeIHLF9YcCWS1JgMcJ3omAM3cya0txdbjCRJJVLpkQe0GeAkqQgGuE6sGjICWlpg7VpYaIOSJKli9rjKFbgFBVciSY3JANeFewa1AHDhB79fcCWSJJVHZaGvA9oMcJJUBANcFyrz4A60QUmS1GHLFbiFkFLB1UhS4zHAdaGjQRngJEnqsGLoKJYPGcmIDetg8eKiy5GkhmOA68KsfJL2ZFfZkiRpK5WTnDz5ZLGFSFIDMsB1oWOVLQOcJElb6Vit2dvtSFKfM8B1Yc6YCQDsv2IxbNxYcDWSJJXHrLEGOEkqSo8CXERcHRFLI+LhLp6PiPhqRMyMiOkRcXzVc2+NiCfzn7fWqvDetn7gYBaMbGFg+2aYM6fociRJJdSI/RG2LPRlgJOkvtfTK3DXAGfu4PmzgIPzn8uBbwJExDjgE8BJwInAJyJi7K4W29cc4y9J6sY1NGR/zK/A2R8lqc/1KMCllG4D2nZwyHnA91PmLmBMREwAzgBuSSm1pZSWA7ew40ZXKo7xlyTtSKP2xzljs2kGzJwJmzcXW4wkNZhazYHbB3i6ant+vq+r/duJiMsjYlpETGttba1RWbtntmP8JUm7p1/2x7WDhrJoxB6wYQPMm1d0OZLUUEqziElK6aqU0pSU0pSWlpaiywGqhlAa4CRJBSljfwSYY4+UpELUKsAtAPar2t4339fV/row20nakqTd0y/7IzhKRZKKUqsAdyNwcb7a1guBlSmlRcDNwOkRMTafnH16vq8uzB+9FxubmuHpp2HduqLLkSTVn37ZHwFmuYv+mUQAACAASURBVNCXJBViQE8OiogfAycD4yNiPtnKWQMBUkpXAjcBZwMzgbXA2/Ln2iLiM8A9+Vt9OqW0o8nepbKpeQDzxuzNQW0LsonaRx1VdEmSpBJp1P4ILvQlSUXpUYBLKV3YzfMJeGcXz10NXL3zpZXDnLETswD3xBMGOEnSVhq9PwIGOEnqY6VZxKSsnAcnSdL25o3Zm83RxOY5cznkA78suhxJahgGuG54s1JJkra3sXkgT4/ei+bUzn4rFhddjiQ1DANcN2Y5xl+SpE5Vbrdz4PK6WkBTkuqaAa4bLpMsSVLnKj3ygDYDnCT1FQNcN5aMHAdDh0JrKyxfXnQ5kiSVRuUKnAFOkvqOAa4bKZrg4IOzDefBSZLUoXIF7sDlCwuuRJIahwGuJw45JPttgJMkqUNloa/JBjhJ6jMGuB742oLsz/Tf35hacCWSJJXHwlHjWd88kL2ebYPVq4suR5IaggGuByo3Kz3AVbYkSeqQook5YydkG45SkaQ+YYDrgVmusiVJUqdme7sdSepTBrgeqKyyNXn5Qkip4GokSSqPjtvteAVOkvqEAa4Hlg8dxYohIxi5YR0sXlx0OZIklcas/CQnjz9ebCGS1CAMcD0Rwcw99sseP/JIsbVIklQiT9kfJalPGeB66LGWydmD6dMLrUOSpDJ5fPyk7MGjj8LGjcUWI0kNwADXQ4/teUD2wAAnSVKHNYOHMXfM3rBhgwuZSFIfMMD10AyvwEmS1ClHqUhS3zHA9dDjleb0yCOwaVOhtUiSVCaPtThKRZL6igGuh54dPIx5o/dyiIgkSdt4tDLN4MEHiy1EkhqAAW4nPGaDkiRpO4/tOTl74BU4Sep1Brid4Dw4SZK2N2/M3jBsGCxYAMuWFV2OJPVrBrid4CRtSZK2l6IJjjoq23jooWKLkaR+zgC3E7yVgCRJXTj66Oy3PVKSepUBbifMHbM3DB0K8+dDW1vR5UiSVB4GOEnqEwa4ndDe1AxHHpltOEREkqQOb7xzDQAPTr2t4EokqX8zwO2sY47JfnuGUZKkDpV54oc+Mxc2by62GEnqxwxwO8shIpIkbWfVkBHMH9XCkE0bYObMosuRpH7LALezKgHOe8FJkrQVV2uWpN5ngNtJx/x8EQDr7p/uEBFJkqq4WrMk9T4D3E5aOXQkC0eOZ+im9fDUU0WXI0lSaXgFTpJ6nwFuF9igJEna3oyW/Aqc0wwkqdcY4HbBDIeISJK0nTnjJvLcgEEwdy6sXFl0OZLULxngdoFX4CRJ2t7mpmaeGL9/tuH9UiWpVxjgdoFX4CRJ6pwnOSWpdxngdsHscfuwvnkAzJ4Nq1YVXY4kSaXxWIsnOSWpNxngdsHmpmaeHD8p23CIiCRJHWbsOTl7YICTpF5hgNtFDhGRJGl7W/XH9vZCa5Gk/qhHAS4izoyIxyNiZkRc0cnz/xkRD+Q/T0TEiqrnNlc9d2Mtiy+S8+AkSfbH7S0fNhomToQ1a7KpBpKkmhrQ3QER0Qx8HXgVMB+4JyJuTCk9WjkmpfS+quPfDRxX9RbrUkrH1q7kcpjhFThJamj2xx04+mhYuDDrkQcdVHQ1ktSv9OQK3InAzJTSrJTSBuAnwHk7OP5C4Me1KK7MHqtcgXvoIYeISFJjsj924crlwwH4z3//acGVSFL/05MAtw/wdNX2/HzfdiJiEnAAcGvV7iERMS0i7oqI13b1IRFxeX7ctNbW1h6UVay2YaNh771h9ershqWSpEZjf+xCZZTKYa1zCq1DkvqjWi9icgHw05TS5qp9k1JKU4CLgP+KiE7HUqSUrkopTUkpTWlpaalxWb3k6KOz3w6jlCTtWEP1x8oolcNanQMnSbXWkwC3ANivanvffF9nLmCb4SEppQX571nA/7H1+P+6duWKEQD895euY/IVUwuuRpLUx+yPXZiV3y/1gOWLYOXKosuRpH6lJwHuHuDgiDggIgaRNaHtVsuKiMOAscCdVfvGRsTg/PF44CXAo9u+tl49OOEQAI5d9HjBlUiSCmB/7MLG5oE8umd+QfGee4otRpL6mW4DXEppE/Au4GZgBnB9SumRiPh0RLym6tALgJ+klFLVvsOBaRHxIPAn4AvVq3PVu/snHgbAcQsfJ5ILmUhSI7E/7tj9Ew/NHtx5544PlCTtlG5vIwCQUroJuGmbfR/fZvuTnbzur8BRu1FfqS0eNZ6FI8czcfUzHLisq1EzkqT+yv7YtfsnHgr3AnfdVXQpktSv1HoRk4ZTOcN4/MLHCq5EkqTyuG+fw7MHd90FW118lCTtDgPcbqoEuOMMcJIkdVgwqoWlw8dCWxs8+WTR5UhSv2GA2033TczOMB630IVMJEnqEMF9+2RzxR1GKUm1Y4DbTY/sfRAbmgZwaOvc7KbekiQJqFrIxAAnSTVjgNtN6wcM4tG9DqSJBHffXXQ5kiSVRmW1ZleilKTaMcDVgGcYJUna3vS9nwfNzTB9OqxZU3Q5ktQvGOBqwDOMkiRt77mBQ3ho/GRob+eCy75adDmS1C8Y4Gpgq0naLpUsSVKHyklOF/uSpNowwNXA/FF70jp8DCxbBk89VXQ5kiSVRuUkpwFOkmrDAFcLEQ6jlCSpEx33S13wmKNUJKkGDHA1ct9E73UjSdK25o6ZwLKho2hZuwLmzCm6HEmqewa4GnElSkmSOhFhj5SkGjLA1cj0vQ/Olkp+8EGXSpYkqYrTDCSpdgxwNbJu0BA4+mjYvBnuvbfociRJKg2vwElS7RjgaugHTADgCx//bsGVSJJUHg9OOIR2Au6/H9atK7ocSaprBrgauq/jXjePFVyJJEnlsWbwMB5vmQSbNsF99xVdjiTVNQNcDXUslbzwcZdKliSpisMoJak2DHA1NGfsRNqGjmLPNcth7tyiy5EkqTRcyESSasMAV0sulSxJUqe8X6ok1YYBrsYMcJIkbW/WHvvAmDGwYAHMn190OZJUtwxwNdZxhvGOO4otRJKkEknRBCedlG3YIyVplxngauz+iYexvnlAdi+4traiy5EkqTxOPjn7/Yc/FFqGJNUzA1yNrRs0hPv2OTxbhfLWW4suR5Kk0jh3xhAA5l9/o6s1S9IuMsD1gr9MPi57cMstxRYiSVKJPLLXgSwfMpJ9Vy2FmTOLLkeS6pIBrhfcPvnY7IEBTpKkDu1Nzdwx6Zhswx4pSbvEANcLHt7rIBg7FmbPhlmzii5HkqTS6DjJ6Tw4SdolBrhe0N7UDKeemm14hlGSpA63H5BPM7j1Vti0qdhiJKkOGeB6y2mnZb8NcJIkdZg/ei9mj50AK1fCtGlFlyNJdccA11te9ars9623wubNxdYiSVKJ3O5iX5K0ywxwveXAA7Of5cuze8JJkiQAbp/kYl+StKsMcL2pchXOBiVJUoc7Jx0NTU1w552wenXR5UhSXTHA9ZLJV0zlH5eMyzZcaUuSpA6rhoyAE07IFjG57baiy5GkumKA60V/nXQM7QTccQesWVN0OZIklcZXmyYDcPUnriq2EEmqMwa4XrRqyAimT3gebNzoGUZJkqpUbifw0jkPFFyJJNUXA1wvc6UtSZK2d//EQ1kzcAiHLJsHCxYUXY4k1Q0DXC+7fbIrbUmStK2NzQO5a/+jsg3niktSj/UowEXEmRHxeETMjIgrOnn+kohojYgH8p+3Vz331oh4Mv95ay2Lrwf3TTwchg2Dhx+GRYuKLkeSVGP2yF3nSU5J2nndBriIaAa+DpwFHAFcGBFHdHLodSmlY/Ofb+evHQd8AjgJOBH4RESMrVn1dWDDgIHwildkG55hlKR+xR65ezruB/eHP0BKxRYjSXWiJ1fgTgRmppRmpZQ2AD8Bzuvh+58B3JJSakspLQduAc7ctVLr12eemwjAz750TbGFSJJqzR65G54cvz+LR4yDJUuykSqSpG71JMDtAzxdtT0/37et8yNiekT8NCL228nXEhGXR8S0iJjW2trag7Lqx1YrbXmGUZL6k17vkf25PxLBHQ6jlKSdUqtFTH4NTE4pHU12BvF7O/sGKaWrUkpTUkpTWlpaalRWOTw+fhJLRoxjr2fb4AGXS5akBrNbPbI/90eA2yqrNU+dWmwhklQnehLgFgD7VW3vm+/rkFJallJan29+G3hBT1/bECK45XknZY9/9rNia5Ek1ZI9cjf96aATYMAA+POfob9dYZSkXtCTAHcPcHBEHBARg4ALgBurD4iICVWbrwFm5I9vBk6PiLH5xOzT830N56ZDX5I9uOEGh1FKUv9hj9xNq4aMgNNOg82b4Ve/KrocSSq9bgNcSmkT8C6ypjIDuD6l9EhEfDoiXpMf9s8R8UhEPAj8M3BJ/to24DNkDe4e4NP5vobzt/2Pom3oKHjiCXjkkaLLkSTVgD2yRl7/+uz3T39abB2SVAcG9OSglNJNwE3b7Pt41eOPAB/p4rVXA1fvRo39wuamZm4++IVcOP33WYM68siiS5Ik1YA9sgbOOw/+4R/gj3+EtjYYN67oiiSptGq1iIl64LeVYZSeYZQkaYvx4+GUU2DTJrjxxu6Pl6QGZoDrQ3+ddAyMHZsNoZwxo/sXSJLUACZfMZWPNB8GwK2f+XrB1UhSuRng+tCm5gHw2tdmG16FkySpw+8PfiGboym7Z+qKFUWXI0mlZYDra07UliRpO8uGj+Fv+x3JoPZN8OtfF12OJJWWAa6vnXoqjB4N06dnK1JKkiSg6pY7nuSUpC4Z4Pra4MHwmnxlaW/qLUlSh5sPeRHtBNx8M6xaVXQ5klRKBrg+NvmKqbx9zeRswzOMkiR1aB0xjnv2PQLWr4epU4suR5JKyQBXgL8ccBzPDhoK990Hs2YVXY4kSaXhLXckaccMcAVYP2AQfzzoxGzDYZSSJHX43SEvzh7cdBM8+2yxxUhSCRngCuJEbUmStrd41HjunXgYPPcc73jL54ouR5JKxwBXkD8feDwMGwZ33w1z5xZdjiRJpVE5yXn243cUXIkklY8BriDPDRzCb/Y7HoAvvPWTTL7CydqSJAH8Lg9wr3zqblejlKRtGOAK9PMjXwnAm6bfTKT2gquRJKkcFozek7v2O5JhG9fDtdcWXY4klYoBrkD/d+ALWDhyPAcsX8SL5k4vuhxJkkrj2mPPyh5861uQUrHFSFKJGOAK1N7UzHVHnw7ARQ/8ruBqJEkqj98d8mLaho6CBx6Ae+4puhxJKg0DXMGuO/p0NkcTZzx5JyxZUnQ5kiSVwoYBA/npkadmG1ddVWwxklQiBriCLR41nlsPOoGB7ZvhmmuKLkeSpNL48bFn5g9+DCtXFluMJJWEAa4EflRpUFddBe0uZiJJEsDscfvAKafA2rXwox8VXY4klYIBrgRuO+B45o9qgVmz4I9/LLocSZJK490jpgAw45NfdjETScIAVwrVi5k4zl+SpC1uPuTFLBs6isNb58DddxddjiQVzgBXEtcdfTo0N8MvfwmLFxddjiRJpbBhwEBuOOq0bONb3yq2GEkqAQNcSSwduQecey5s2gTf/W7R5UiSVBo/OeaM/MFPYMWKYouRpIIZ4Mrk8suz3//7vy5mIklSbs64fbhj0tGwbh388IdFlyNJhTLAlciBf9rI/FF7wuzZ8Ic/FF2OJEmlce0xZ2UPrrrKxUwkNTQDXIm0NzXz48owkW98o9hiJEkqkd8f8kLYc0946CG4446iy5GkwhjgSub6o09nffNA+NWvsiYlSZLY2DyQrx50CgD/d/F7C65GkopjgCuZ1hFjubZyY+/Pfa7YYiRJKpGrTziPNQOHcPLse+Gee4ouR5IKYYAroW+deD7rmwfQft31nHrZlUWXI0lSKawYOoofHH9OtvHZzxZbjCQVxABXQotHjeeGo15FE4l33nl90eVIklQa/3vC61g3YDDceCM88EDR5UhSnzPAldQ3X/gGNjY1c96jf4aZM4suR5KkUlg2fAw/qkw18CqcpAZkgCupBaP35OfPfyXNqR0+//miy5EkqTS+ddL52YJfP/sZp1/69aLLkaQ+ZYArsa+/6I1siib4wQ+ye8NJkiRaR4zruO3Ou5xqIKnBGOBKbN7YCfzqiFfApk1c+7p3MPmKqUWXJElSKXzrpPPZ0DSAV8/4Czz2WNHlSFKfMcCV3Ndf9CbaCV7/0B+YsKq16HIkSSqFRaNauOHo02giOdVAUkMxwJXcrD325TeHv4xB7Zv4h7/9rOhyJEkqjcqCX1x7rQt+SWoYBrg68D8vehMAFz54s3PhJEnKzR+9Fz9//ith82b45CeLLkeS+kSPAlxEnBkRj0fEzIi4opPn/yUiHo2I6RHxx4iYVPXc5oh4IP+5sZbFN4onWybxiyNOZvDmjfC+9xVdjiQpZ38s3v+8+E0weDD86Efwl78UXY4k9bpuA1xENANfB84CjgAujIgjtjnsfmBKSulo4KfAl6qeW5dSOjb/eU2N6m44/3by21g9aCj86ldw001FlyNJDc/+WA7zx+wNV+TZ+Z3vhE2bii1IknpZT67AnQjMTCnNSiltAH4CnFd9QErpTymltfnmXcC+tS1TS0fuwX+95KJs45//GZ57rtiCJEn2x7L48IfhgAPgoYfg694XTlL/1pMAtw/wdNX2/HxfVy4Fflu1PSQipkXEXRHx2l2oUbnvveBceP7z4amn4CtfKbocSWp09seSmPypW7n0+LdkGx//OCxeXGxBktSLarqISUS8BZgCfLlq96SU0hTgIuC/IuKgLl57ed7IprW2ulx+ZzY1D4CvfS3b+PznYe7cYguSJPWI/bH3/fF5J/HHg06AVauyK3KS1E/1JMAtAPar2t4337eViDgN+H/Aa1JK6yv7U0oL8t+zgP8DjuvsQ1JKV6WUpqSUprS0tPT4CzSayb9bw42HvxzWrXNBE0kqlv2xZD516uXZgibf/z7cfnvR5UhSr+hJgLsHODgiDoiIQcAFwFarZUXEccC3yJrT0qr9YyNicP54PPAS4NFaFd+oPnfK3/PsoKHwi1/A735XdDmS1KjsjyUzb+wE+NCHsg0XNJHUT3Ub4FJKm4B3ATcDM4DrU0qPRMSnI6KyataXgRHADdssh3w4MC0iHgT+BHwhpWSD2k1LRo7nv198Ybbx7nfD+vU7foEkqebsj+V02JpjmD9qT5g+Hb75zaLLkaSai5RS0TVsZ8qUKWnatGm79R6Tr5hao2rKacDmTcyc+hGYMSMb6/+FLxRdkiTttIi4N58Hph6wP/bM6U/cyVW/+ByMHAn33w8HdTq9UJJKraseWdNFTNR3NjUPgG9/G5qb4YtfhFtuKbokSZJK4fcHv5Cph74EVq+GCy+EDRuKLkmSasYAV8cm37icr+RDKVtf90ZYsqTgiiRJKoEIPnLmu5k/qgXuuYcrX3Zh0RVJUs0Y4OrcN174Bu7c/yha1qyAiy+G9vaiS5IkqXCrhozgn8/9EJuiiX+8++dw881FlyRJNWGAq3PtTc2899Xvp23oKPj97+Hf/73okiRJKoX79j2c/3zpm7ONiy/2Bt+S+gUDXD+wZOR4PnD2e7ONj34U7r672IIkSSqJb77w9dwx6WhYutSRKpL6BQNcP3Hr807k6he8BjZtYu4Z58HKlUWXJElS4dqbmnnfOe+H8eOzBb++/OWiS5Kk3WKA60e+cPLbeHivg5i0YjFccgls3lx0SZIkFW7pyD3ge9/LNj72Mbj11mILkqTdYIDrRzYMGMi7X/MhVg0eDr/8JbznPVDC+/xJktTnzj4bPvhB2LQJXvc6ePDBoiuSpF1igOtnZo/bh8v+v4+xvnkAfP3rfPGUtxVdkiRJhZt8xVQOiJfxm0NfCqtWwVlnwdy5RZclSTvNANcP/W3/o3jfqz9AO8GH//w9+P73iy5JkqTCpWji/a/+F3jFK2DRoizEtbUVXZYk7RQDXD9102Ev5dOnXpZtXHqp97+RJAlYP2AQRx/7Tzw2fhLMmAHnngvr1hVdliT1mAGuH7tmymu48qTzs/H+558P995bdEmSJBVu1ZARXPKGT7Fw5Hj461/hootc+EtS3TDA9XNffMVb4c1vhjVrsqEi991XdEmSJBVu8ajxvPUNn4IxY7KFvy65BDZuLLosSeqWAa6fS9EEV1+dhbfWVjj5ZPjTn4ouS5Kkwj3ZMok3nH0FawYOgR/+MFudcu3aosuSpB0ywDWAyR+/hYOffzm/PuxlsHo1nHkm/OIXRZclSVLh7tnvSN58wedoGzoKpk6F00+H5cuLLkuSumSAaxAbmwfynnM/AO98J2zYAK9/PXz720WXJUlS4R6YeChvuOiLsO++cMcd8PKXw8KFRZclSZ0ywDWQ9qZmJg8/k/98yUXQ3g6XXQZf+II3+5YkNbynxu/Hi879DDPH7QsPP8zTRxwPTz5ZdFmStB0DXKOJ4L9fehEfe9U/QQR85CNw8cXw7LNFVyZJUqEWjWrhDW/+Ig9MOIT9Vi6Bk06C3/ym6LIkaSsGuAb1w+PPgeuug2HDsonbJ5wADz9cdFmSJBVq+bDRXHTB5/jDQSdkc+HOPRc++EFXqJRUGga4Bjb53mGcduFXeGKP/eGxx+DEE7MVKx1SKUlqYGsHDeWy8/+Vz5/8NjZFE3zlK/CKV8C8eUWXJkkGuEY3c/z+nHfxf3DDkafBunVw6aX87KhTHVIpSWpoKZq46qTzedNFX8hu+H3nnXDccQ6plFQ4A5xYN2gIHzznvXzg7PeybsBgzn/kT3DkkXDjjUWXJklSoe7d9wjOueS/ufXAKdDWlg2pfMtbYMmSokuT1KAMcOrw06NO4zUX/weP7HkgzJ0L550Hr32tQ0YkSQ1t+bDRXPr6j/P5k9/GcwMGwY9+BIceCt/8JmzeXHR5khqMAU5bebJlEq9563/yqVMvY/WgofCrX8Hhh8OXvuQEbklSw6oMqTz977/O/x3wAli5Et7xDnjxi+G++4ouT1IDMcBpO5ubmvnulPM47e3f5DeHvQzWroUPf5iZex8IP/lJdg85SZIa0LyxE7jkDZ/kn867gsUjxsHdd2crOf/938OsWUWXJ6kBGODUpSUjx/Ou8z7MxW/4FHPGTOB5bfPhwgvhqKPg+usNcpKkxhTBbw97Kae+/Uq+PeU8NiXgu9/NhlW+/e0we3bRFUrqxwxw6tZtB76AV739G3zkjHfB/vvDo4/Cm94ExxwDN9zg+H9JUkNaM3gYnz31Mk697Ep+euSp2YnN73wHDjkELrvMK3KSeoUBTj2ysXkgPz72TA5+43/z0TPeyYKRLdmNv9/4RjjgAPjc52Dx4qLLlCSpz80dO5EPnPM+Tr70m/zs+aeweXM7fPvb8LznwTnnwNSpnuyUVDMGOO2Ujc0DufbYszjl8qv42OnvYPbYCfD00/Cxj8F++8EFF8Cf/+zNwCVJDWfOuH14/6vfz6su/QY/PfJU1jc1w003watfDQcdBP/2b7B0adFlSqpzBjjtkg0DBvLD487mlZd9i79746e5+eAXZmccr7sOTj4ZDjwQPvzhbGUuw5wkqYHM2mNfPnDO+3jRO77Hv518CfNG75XdnuejH2XT3hP484Ev4INnvxeWLy+6VEl1KFIJ/+N6ypQpadq0abv1HpOvmFqjatRTE1a1csGDN/Om6b9n72fbtjxx0EHZUMvzz4fjjoMmzxtIykTEvSmlKUXXUS/sj/UpUjsvn30/b7l/Kqc8NY0BKV8EbOBAeNWrsh559tnQ0lJsoZJKpaseaYBTzUVq54T5j3LOY3/h7MfvoGXNii1Pjh+fNaszzoDTT4cJE4orVFLhDHA7x/5Y/8auXckZT9zJOY/dzsuenr5lRecIOP74rDeecQa86EUwaFCxxUoqlAFOhWhq38yJ8x/hnMdu55SnprHvqm3G/h95JLz0pdmNUF/84mzoZUQxxUrqcwa4nWN/7F/2WLOCM5/4K2c8cScvX/QorF+/5ckRI7bujyeeCCNHFlespD5ngFPxUmLO3x8MN98MN9/M2lv+yLCN67c+pqUla1QveEF2m4JjjsluXWCok/olA9zOsT/2X0M2PseJTz/Cy2ffx8tn388hy+ZtfUBTExx9NJx0Ehx7bPZz1FEwfHgxBUvqdQY4lc6gTRs5ZtHjHL/wMY5fkP20rF2x/YFjxmRB7sgjs5ukHnpodo+d/fd3Pp1U5wxwO8f+2Dj2Wv0MU+bP4AULZnD8whk8f8ksBrZvfSuCdoI5/3979xtbV13Hcfz9XXvbu623XdfRrmvn1lJExjSABEMQIhLMFoj4QMNMNMaQ8AQU4wMjPvFPfKBP/JNoTMjAgKILATGLiCh/RCXiBohhUKbLtrB2dHeuW9vb9vbeO74+OOey29rSG+7tzj2Hzyv55Z5z7unh+8s2Pv2d+zu/29nLcPcAN33m+nMZefHF0N4eUeUiUi9LZWRzFMWIQLCS5f7N29m/eXtwwJ33nRnjiuOvs+3EYS7JHmFb9jBdZ84EX03w7LPzL5BOB9+xMzAAW7bA1q1B27IF+vuDT/Oams53t0RERGp2IrOBxy65lscuuRYIPqH70NghPvjmf9h28gjbThxm6NQxBk8fZ/D0cfjuc/MvsHFjsIhYZTZu3Rrc/OzrC6ZoikgsaQAnjcOMNzp7eaOzl99een1wzJ3u3DjbskcYOvUGg+PHGRwfYXB8lO7p08GXiR84sPj1mpqCANu0KWg9PcGgrrs7eC23zk5Yvz6YhqKpmiIi0oDyqTT7Nm9nX/mmJ5A6W2To1DE+kD3KDy9bDQcPwsGD5IcPkh4bg7ExeO65xS+YyZzLx97e+dlY3u7qCjKyszNYMVNEGkJVAzgz2wH8GGgCdrv79xa83wo8AHwYOAXc6u5Hw/fuBm4DzgJfdvcn6la9JJ8Z2UwX2UwXf75w/ifImblptpx+k76JLP2TWfonTtA/Ebx258bpmp2E0dGgVaO5+VxQdXQErb39XMtkgjuWbW3BYK/8umYNrF49/zWdy1gXiAAAB3BJREFUDlprqz4FFEk4ZaREpdiUYrh7kOHuQR4tARdeDReC7XyLvsmT83KxfyJL32SWTZMn6cmNk56aenvAV5W2tiAf162bn4/l7cpsrMzHcqvMyXJGplK6cSryLiw7gDOzJuCnwI3ACLDfzPa6+2sVp90GnHb3ITPbBXwfuNXMtgG7gEuBTcCTZvZ+d58/iVvkXZhqXcuBjUMc2Di06PstpSIXTJ+mJ3eK7tw4G6bP0DUzwfrZCbpmJlk/M0Hn7CQd+Rzr8lPBgionTwatnlKpIKhaWoIBXeVrS0vwfio1fzuVCgaUC1tT0/ztd2qrVgWtcrvczJbeL2+b/f/2wgZLv7fY++X9yuNli23X+n6lao8tZSXOjdMvLovVmskEz6O+hykjpRG5rWKko4eRjp4lTnDa56bpzo3TkxunJ3eK9TNBNnbNnHl7uyM/FWZkjqZcDnI5OHasvsWWB3OLZeRi2bhURi7Mx3fKycr8W7i/MB8Xy8jlsnG5Bkvvl7crX5c79m7fr3S+MnKlcjdKS9U5NBTc4FgB1XwCdxVwyN0PA5jZHuAWoDKcbgG+FW4/DPzEzCw8vsfd54AjZnYovN7f61O+yNIKzSlGO7oZ7eiu6vyWUpGOfI72fI5MYYa2uRkyc9O0FWZon5thTWGWtcU8q4t51hbyb++nSwVaS3OsLs6RLhVYXZyjtVSg5WyJ1lKBVcUiFIsr3FuRCNxwAzz5ZNRVRE0ZKfFjxmS6jcl0G4c2vG/5891pK8yyLj9Fez7Ixba5GdoKM2TmghZk4yxryhlZnJ2Xi+lSgXRxjnRpjpazRVpLRVreKkE+HzSRpHn8cdixY0UuXc0Arg+ovN0yAnxkqXPcvWRmE0BXePz5BT/bt9h/xMxuB24Pd3NmVuVn+kvaAPy3xms0GvUpHpLYJ0hmv9SnWjz1VD3ukG6pRykRWvGMVD5WLYn9Up/iQX2Kj/PXr50763GVRTOyYRYxcfd7gHvqdT0zeyFpS1OrT/GQxD5BMvulPkkcKB+rk8R+qU/xoD7FR1L6Vc2XaI0Cmyv2+8Nji55jZs1AB8GD2tX8rIiISFwpI0VE5LyqZgC3H7jIzAbMrIXggeu9C87ZC3wh3P408LQH3xC+F9hlZq1mNgBcBOyrT+kiIiKRU0aKiMh5tewUynC+/p3AEwRLJN/n7q+a2XeAF9x9L3Av8IvwAexxggAjPO8hgoe5S8Ad53F1rbpNN2kg6lM8JLFPkMx+qU9Sk5hmZFL/jiSxX+pTPKhP8ZGIfllwE1BEREREREQaXTVTKEVERERERKQBaAAnIiIiIiISE4kbwJnZDjM7aGaHzOzrUddTD2Z2n5llzexA1LXUi5ltNrNnzOw1M3vVzO6KuqZamVnazPaZ2b/CPn076prqxcyazOyfZva7qGupFzM7amavmNnLZvZC1PXUg5mtM7OHzex1Mxs2s6ujrkkaizKy8SUxH0EZGSfKx8aXqGfgzKwJ+DdwI8EXou4HPuvur0VaWI3M7DogBzzg7tujrqcezKwX6HX3l8wsA7wIfCrOf1ZmZsBad8+ZWQr4G3CXuz+/zI82PDP7KnAl0O7uN0ddTz2Y2VHgSndPzBeVmtn9wF/dfXe4IuIadz8TdV3SGJSR8ZDEfARlZJwoHxtf0j6Buwo45O6H3b0A7AFuibimmrn7XwhWLksMd3/T3V8Kt6eAYaAv2qpq44FcuJsKW+zvkJhZP3ATsDvqWmRpZtYBXEew4iHuXohzOMmKUEbGQBLzEZSREp0k5mPSBnB9wLGK/RES8D+9pDOzrcDlwD+iraR24TSKl4Es8Cd3j32fgB8BXwPeirqQOnPgj2b2opndHnUxdTAAnAR+Hk7l2W1ma6MuShqKMjJmkpSPoIyMEeVjg0vaAE5ixszagEeAr7j7ZNT11Mrdz7r7ZUA/cJWZxXo6j5ndDGTd/cWoa1kBH3X3K4CdwB3hNKw4awauAH7m7pcD00AinnESeS9KWj6CMjJGlI8NLmkDuFFgc8V+f3hMGlA4B/4R4EF3/03U9dRT+NH8M8COqGup0TXAJ8P58HuAj5vZL6MtqT7cfTR8zQKPEkwvi7MRYKTijvbDBIElUqaMjIkk5yMoIxud8rHxJW0Atx+4yMwGwgcUdwF7I65JFhE+zHwvMOzuP4i6nnowswvMbF24vZpgoYDXo62qNu5+t7v3u/tWgn9PT7v75yIuq2ZmtjZcHIBwGsUngFivYOfuY8AxM7s4PHQDEOtFD6TulJExkMR8BGVkXCgf46E56gLqyd1LZnYn8ATQBNzn7q9GXFbNzOzXwMeADWY2AnzT3e+NtqqaXQN8HnglnA8P8A13/32ENdWqF7g/XOltFfCQuydiSeEE6gEeDX5Pohn4lbv/IdqS6uJLwIPhL+eHgS9GXI80EGVkbCQxH0EZGRfKxxhI1NcIiIiIiIiIJFnSplCKiIiIiIgklgZwIiIiIiIiMaEBnIiIiIiISExoACciIiIiIhITGsCJiIiIiIjEhAZwIiIiIiIiMaEBnIiIiIiISEz8DyvgBX1W1/tkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 2)\n",
    "fig.set_figheight(7)\n",
    "fig.set_figwidth(15)\n",
    "width = 2\n",
    "bins = 100\n",
    "x = np.arange(0, max(chain) + 0.5, 0.1)\n",
    "\n",
    "ax[0].hist(chain, bins, density=True)\n",
    "ax[0].plot(x, param * np.exp(-param * x), linewidth=width, color=\"red\")\n",
    "ax[0].set_title(\"Metropolis-Hastings\")\n",
    "\n",
    "ax[1].hist(exponential_rvs, bins, density=True)\n",
    "ax[1].plot(x, param * np.exp(-param * x), linewidth=width, color=\"red\")\n",
    "ax[1].set_title(\"Inverse transform\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results look very similar. Additionally, the results compare rather well against the theoretical density, which is plotted in red."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Statistical Test for the Metropolis-Hastings Algorithm\n",
    "\n",
    "Finally, we test whether our set of observations pass statistical tests for being exponentially distributed. Our null hypothesis is that the set of observations are exponentially distributed.\n",
    "\n",
    "Here **our interest is specifically in the Metropolis-Hastings algorithm**, since the inverse transform method is already mathematically proven to produce correct results. However, to rule out any errors in our implementation of the inverse transform method, we test it as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AndersonResult(statistic=0.5561689854584984, critical_values=array([0.922, 1.078, 1.341, 1.606, 1.957]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))\n",
      "--------------------------------------------------\n",
      "AndersonResult(statistic=0.5178271921322448, critical_values=array([0.922, 1.078, 1.341, 1.606, 1.957]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))\n"
     ]
    }
   ],
   "source": [
    "print(stats.anderson(chain, dist=\"expon\"))\n",
    "print(\"--------------------------------------------------\")\n",
    "print(stats.anderson(exponential_rvs, dist=\"expon\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, for any sensible significance levels we cannot reject the null hypothesis. Thus, our implementation of the Metropolis-Hastings algorithm seems to produce believable results."
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}