phys512/metropolis_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "438a07af",
   "metadata": {},
   "source": [
    "# Metropolis-Hastings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3eade768",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.integrate\n",
    "\n",
    "seed = 239\n",
    "rng = np.random.default_rng(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "087a3409",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acceptance fraction = 0.3567\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGdCAYAAAA1/PiZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABpQ0lEQVR4nO3deVxU9f7H8dcAAwwgoCAqiCtquSSmaS65pOXSZpapmWmZbViabbbYcm9d76/dirIsU7uV3krLsmvlkpWZa5SpqZgZggqogOzb+f0xMoKAAgKHGd7Px4OHcuacM5/DMvPm+/2e79diGIaBiIiIiBNwM7sAERERkYpScBERERGnoeAiIiIiTkPBRURERJyGgouIiIg4DQUXERERcRoKLiIiIuI0FFxERETEaXiYXUB1KywsJCEhgQYNGmCxWMwuR0RERCrAMAxOnDhBaGgobm7lt6u4XHBJSEggPDzc7DJERESkCuLi4mjevHm5j7tMcImOjiY6Opr8/HzAfuH+/v4mVyUiIiIVkZaWRnh4OA0aNDjjfhZXW6soLS2NgIAAUlNTFVxEREScREXfvzU4V0RERJyGgouIiIg4DQUXERERcRoKLiIiIuI0FFxERETEabhMcImOjqZjx45cdNFFZpciIiIiNUS3Q4uIiIjpdDu0iIiIuBwFFxEREXEaCi4iIiLiNBRcRERExGkouIiIiIjTcJnVoUVEziQ+JYvjGbmOzxv6ehIWaDOxIhGpCpcJLtHR0URHR1NQUGB2KSJSx8SnZDHkxXVk5Z16fbBZ3Vl1/wCFFxEn4zLBJSoqiqioKMd94CJSPVyhpeJ4Ri5ZeQW8MiaSiBA/YhPTmb4khuMZuU53LSL1ncsEFxGpfq7WUhER4kfnMP1hI+LMFFxEpFxqqRCRukbBRUTOSi0VIlJXKLiISI0ze5yMNT2eTpb9eCcHgMUP7+R0QkmutecXkeqj4CIiNcr0cTIpcbT7+FJWeGXBMvumCGCVlxdx6RcBAY79yDxa8lifIAgMr/kaRaTCFFxEpEZVZJzM2VpkTn+8rH3KlXkUt/wspuXezT1jriCisR9xe2MIXzsN9+xj9n1S4iC6J+RlljzW6gNRm+pUeCnra1GcM971JVIZCi4iUivKGydzthaZsh4vsQ/JFWopiTXCyA7uAqEB5CSll9w/86g9tIyaB8Ht7duS98DSKfbHKhhcKtUldnoLTwVad8r6WoSSTEPLCcfnWR6BvH//dQov4rJcJrhoAjoR5xL/1x68k3OIO5ZFm/xYHhjagdDQMHZlBpRokTm9xQZwtNqkH9kPnwyuvpaS4PYQGlm166lMl1hZLTwVqPn0r4U1PZ52H0/GLT/LsU+m4UXckQshsGOVrkOkrnOZ4KIJ6EScQ0NfT9pYj3PJN7fgY8khAhjkBazD/uZ9/eoyjyurxcY9+1i1tJRUh/K6xDbvP8bxk4HL0QJzegtPJWt2fC0S9kN+luM8pbrARFyQywQXEXEOYYE2Fk9oh8+HOcQNmkNOYAT+Nish2Qdg6ZSSb7opcXgnHyhxRxA+QYB/yZOW1VKSvKfkv7WkKFQ09PXEZnVn+pIYR3eOt4cbb07obr/W8uquipPnKdUFJuKCFFxEpNaF+HkBEN4u0vHGbcR7kpVnUFhYCNhvYU54dSCHkzJ4wWYhYpm7/WCrD6v9HiL994OcSL/ccc5V6zezes1aWoU0YIqHDbelUwBYEJNLar6VxNaGY98///qbbzfn4hH0ExFd+zm2L/n8a1LdN3H11VfT9OS2Q0eS2LL1Cxo3bkyr0CDSjx8B4O+DCfjbrHTu2Bnvpu1KXWNYoI1V9w8g/cj+kt05H+K4DnsIq1nnNLBZpA5ScBGRGpOdnU3cgb8wjFOhYenSpXz24XwG5OYy+eS2goICbK17kpeXz/pBxwAr7tnHeH9rGjNX5xDWpRffvf8CEZYEWDqFJ2Y9QdqJdJKmdAcgMT2Hm555nyMro/Fp34eF1/6fY8DqmjWzyD6RQqvJFhr6egLw6/adPPRVNv1TlzLlnocctT3xwpvs+fMAHTt2pGkbe/fOT2tXcv39L9EzsiNrr06iqSUHgDFvp7PtUCGf3xTI1a/9BviTfXAnVw24j+6RF/DJJ5+c7BbK5YsdacR1mMSnyWE8fl1PwhvZKPBuhB/BhFXi63n6fDRna00668BmhRdxQgouIk7E7IncwH4Xi3fydvsb50l/JWeyeXc8oaGh9O3bF4Dc3Fx8fHwwDIPm937o2Hf37t28/+kK3COtjuDi7u6Oj82b1Lx0srJzACsATfwstGrRnGO2k3cEnXzOi3t25/t9aXh72Vtu0rLyMPyb0veyq+jT52JunjLS8Xz/jP+F1JTjPDH1UsfXKrRpE0ad70F4184lrm1wv56c3yWSRo0agU8DsPoQEDOXnmFuRLrtBQL5qffbhDQNI3/Z3Xgk7SLQI/fk3UH+FGam8Ne+vYQ18oWEGPtJk/fwzA85bFr8JmGjn+CWr3OBXHKP/ELq2neZNHII0Q/ddPYvfBnz0QBnbLk508BmLdsgzkrBRcRJmD6RG+CW+jezEqbxxUPZ3HWRJ94eFgD+82MBs1ZncPX1Y3m2VWdHoAoODiYt7QSFGamOc1x++eV4ZCYS+dfbJc79x7pl+C4ZxTHfDDpZkvBK8WNSpCf9nv6IIR+lltg3+sVnGPJRKgVe9haUuGNZ2Fp3Y+5L95QawLtk0bulrqPXRd0Yf4MPsddOKrH9jdmPlhxzErWJIZlHGQLEJqUz5KN9vN15GBFhAXzw1Q88/Pr7XOL5KCTvwdtIZ1ALC3dNboivZS+8PcBxmouae+PdsguP3zOCZm3OA+C9RX/x0oJf2fizL3AquAwbNoykpCRef/11evfufaqWMuajASp0G7WWbBBXouAi4iRqZcHD0+YWSd35F+FJmxxdE75HtnHjx6mkZEP/++bT/YKOHPv7dyL33E5gaCvWHfHgytd+dASq2NhYDqQVctXr6x3n7N69O92bucPb80s8ddNWHcDPlwbr7mOFF7AWsPpQ4N0IKBlc/G1WbFZ3XvhmN4O84IVvdmOzRji6gioq7lgW2fGpeCenE1HWDoHhjlCQbaSScFodx40GGFYfLEunEAEsCYTCYBtuY/8LPsGO/V6fXjpcDBs8gIVX3s/N117o2GYYBhs3biQlJQUfHx/H9v/973/84/GHuT4oh9jup+ajEamPFFxE6pCKdAWd/tdzbGL6GfevcPfSaXOLfLg9j/FLs+gb7k5EeDQALYGR53sTH9wPtybnQ2gkx5LSubK9lW/HPkxgy87EHcvisa8T7IEqLADLidTSz1WWwHCI2kTsgQNMWxzDnLGRRLRsSV6GP7C/xK4h2QdYN6ElOYf8YC3MGRuJrWX3Cgc4f5u9K+qFb3az4+tcOln2s8LLPlYmpGLVApBAMHtHr6FDg1xik9KZtjiGl8ZdSoeIs8+h0iwsHL9Og7h0aL8S17dp0ya2b9/Oeeed59j2/fff8/O27XTqZi1xjhdffJFu3brRr18/PD0rF9qKqwtdkCIVpeAiUkdUtiuo+O225e1f0XMuXbqUBXPncFfjVIY/vACC29P773jcPruaHQXh7LlmIe2b+BOblM7u7H28fc9IR3gq8G5EpuFFj20Pw7Zy1gGqqMBwsjP82WGk2lsVAgOgWDcTPkH2MR1Lp5wKGFYfIlq2hEq80Rbd1TR3mB85gQGkHPCAbfaxMsWDS/E39OIBsbg8vzAIDSDbSGWHkWr/vIosR/fSLrg97S5uA0d3ObqB7r77bjqF+tH6l3/xr6La4uN54IEHsFgsHDp0iCZNmgCUGAhdEXWhC1KkMhRcROqIynYFFd1uW/yNtcT+KXFkHThQ5qy0P+6Kp1OLxoA9AH377bd88e33NO9hpbt3SxKN1hDemkVrt/PYV/vJDenqeHM+vbskzy+MITnPs2hc2xLrACUkxJPnF1buG36VnGyVqexU+aWcDEDha6c5NmUaXie7pezKe0M/vTuq6PrO6TqLBbISrD4kX/EuORZ/BlwQRvhRD7Df1EROTg4TJkwgIyPDEVoA7vgym33
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADfpklEQVR4nOx9d7wdRfn+s+fcmuSm914g9BJ67x0pFkARBBUUO4KF/FTUrwUVC4IoNopKFwTpvQkJoQVIICEJhISE9HLTbjtnf3/smd13ZmdmZ8sp9955Pp/knrNndmZ2d3bmnedtjuu6LiwsLCwsLCwsqoBctTtgYWFhYWFh0XthBRELCwsLCwuLqsEKIhYWFhYWFhZVgxVELCwsLCwsLKoGK4hYWFhYWFhYVA1WELGwsLCwsLCoGqwgYmFhYWFhYVE1WEHEwsLCwsLComqoq3YHdCgWi1i+fDlaWlrgOE61u2NhYWFhYWFhANd1sWnTJowePRq5nJ7zqGlBZPny5Rg3bly1u2FhYWFhYWGRAEuXLsXYsWO1ZWpaEGlpaQHgXUj//v2r3BsLCwsLCwsLE7S2tmLcuHH+Oq5DTQsiTB3Tv39/K4hYWFhYWFh0M5iYVVhjVQsLCwsLC4uqwQoiFhYWFhYWFlWDFUQsLCwsLCwsqgYriFhYWFhYWFhUDVYQsbCwsLCwsKgarCBiYWFhYWFhUTVYQcTCwsLCwsKiarCCiIWFhYWFhUXVYAURCwsLCwsLi6rBCiIWFhYWFhYWVYMVRCwsLCwsLCyqBiuIWFhYWFhYWFQNVhCxKBtmvrsWt85aUu1uWFhYWFjUMGo6+65F98Yn/zITADB1RAv2njCoyr2xsIiHdVs6MOu9dZg6oh8mD+tX7e5YWPRYWEHEouz4cOM2AFYQseheuPAfL+OV99cDAOb8+Hj0a7TTpYVFOWBVMxZlR13OqXYXLCxi48MN2/zPG7d1VrEnFhY9G1YQsSg76vN2mFlUBqta29BZKGZSl0s+F4uuspyFhUU62BXCoiwokIm7zgoiFhXAGx9swH4/fwJf+termdTnuvLPFhYW2cKuEBZlAd2VWtWMRSVw3TOLAACPv70yk/pcwokUrSRiYVE2WEHEoizoIIKIVc1YVAJL122LLhQDVBtjBRELi/LBmoFblAWdXYEgkreMiEUZ8cw7q/GVm1/F5vauTOt1rSBiYVER2K2qRVnQWbATt0Vl8PzCNZkLIR6oaqYM1VtYWACwgohFmcB7LthZvLfj5hffx9+ee7csdTOPlp1G9QcATB7aN5N6LSNSm3h39WZ85Jrn8NCbH1a7KxYZwQoiFlqs2dyOe15bhrte+QBvf9hqfF47Uc3Y3WTvRmehiO/9Zw5++sDbWNnalnn9bHi1NHma5kJGQgPvvptJlRYZ4Ou3vYY5y1rxpZuz8Y6yqD6sjYiFFt+683U8PX81AKC5Po9XfnAM+jREDxvKiNjNZO8GdeXeUgYVChtfeccJtZcGlAXpaYxIe1cB/3jhfaze3A4AcAAct8vIbpGK4d3VW6rdBYuMYQURCy1Wtrb7n7d1FrC5vStSEHFdFw/NWcF9t+i9KPciztxs6/KeIJJVc6ZxRApFF13FIhryOThO9zDMfmb+avzswbe5Y0/MW4XHLzk8dl2Footn31mNtVs60NJUhyN3GI6GuvKR7Vs7CmWr26I6sIKIhRZiREmTSf6FRWtx9RMLgjok58x8dy0G9WnADiNb0nbRosbRRQZAOUQSnxHJZcuIuAaMyNsftuKsP89Aa1sX9p80GLd94YBuIYy0tnnM1LjBzdh34mDc/eoybG5LxlY98fZKfOGfr/jff/bRXfHp/Sdk0k8RHV1WR9YTYW1EegFeX7oB37jtNTw1b1Xsc8UJ2EQQeXPZRv4cYflZum4rPvmXmTj+qmdj98ei+6Hc4dGZwMAC55XFRkRR56tL1vuL+ovvres2u/VCyehlhxEt+NzBkwAkZ65WbWrnv7e2K0qmxwuL1vifBzTXl60di8rCCiK9AJffOwf3zl6Oz9/0UuxzxUndZLIa1EeYIIRTFq+1Ot7ehK5yCyKlv4wRyUwVyHnNyIuEGMNsWi472DPJ5xzkSgxO0r6L9zsrRkqGbUTQmzwsG+8oi+rDCiK9AIvXbgWQzHtFnNNNqhjQ3MB9F9t1UPvUtUV2oIt1OdgRNkbrct50Vg5jVZVwIzbVXeyhCkQQYZqkpF0X73dnsYhrn1qIY377TCIWVgfaVJeNVdRjYAWRXoA0kU3FScZkIRFV5KJqphuo0C0yBGVEysGOMIEhcxsRrg15GbGt7rI0BoJIjggiyXov3ptCwcWVj8zHwlWbcfWTC+QnJQSdS7LKsmxRfZRVEHn22WdxyimnYPTo0XAcB/fcc085m7NQIJdi5U+iNxaFlTAjEobrurjiwbdx96sfxG7PorZBF+ty0Pasxrpc+bxmVO9BEhuqtGjrLGDhqs2p2Bf2HOoyUM2I94AKm1vbs7WZoU31VkHkmXdW4//9501OTdXdUVZBZMuWLdhjjz1w7bXXlrMZiwikSfUSFiqipytx12syYb6waC3+/Oy7uOSO1+N10KLmURAYkS3tXfjto/Nxy4tLMqk/5DWTmbFqtNdMWBBRt71k7VZMv/sN3Dt7Wap+nX7t8zjmt8/ghucXJ66DvaM5J1CUJhVsxNPo886lmHy6JIIGbarctke1ivOun4VbXlyCv5YpUnE1UFb33RNPPBEnnnhiOZvoUSgUXTw8ZwVWb/KiT+4/eYgftjoNUqlmEuz4QpOzWEDSnXVbOuJ1zKLboItjRIp4eM4KXP3kQgDA0TsNx4j+TSlbKJNqhjIiis132EZEXd+fnlmIW2ctxa2zluK0Pcck7te8FZsAeB4knztkUqI6KCPiq2YS9ifMiAQ3qy7h3HPTC4vxf/e/hX99fn8cOGWIf5wKS5293JX3w43ZZpuuJmoqjkh7ezva2wPXr9ZW85DiPQHPLliNr9wShC0e2q8RL3//mNT1plPNiN8NGJGCfpdojVXLj5nvrsXvH1+Aj+41BmfuM66qfaFjpqvgYlNbp/89i0irbN3LV0E1E8dGZPWm2hG2fRuRvOPHPUlqSCyeRt//pIzID/87FwBwyR2zMWP60f5xTjVTg4zIui0deOytFfjI7qPRt7G8y2uaeb3WUFPGqldccQUGDBjg/xs3rroTaKWxbrM3UfUv5cxYs7k9Eyt8FSNy4/Pv4ey/ztRmLk3inhjFotD3p7t4GcTF0nVbcftLS6oWgOmXD8/DjHfX4jv/fiOT+tZt6cCF/3gZj7+1UluurbOA22Yt4XLK0IWpUHS5zMxZRF31I6tWQTUjjl/deK6lMPFdlBEpHcuKEaHCWVJGRAX6TGSqm2rj8ze9hO/e9Sa+9583y95W1ve2mqgpQWT69OnYuHGj/2/p0qXV7lJFwV7oycP6kWPp61UJIj+67y28sGgt/vT0wsg+MZgIDqLw8uqS9Tj/hln49ythQ9Qsrs91Xdz/xnK8t6Z24pMc+eun8d273qyaHnfD1s7oQjHw8wffxmNvrcQF/3hZW+6XD8/DZXe/iTP/PMM/RsdQwXXRUcg2IWJgI5Kt+65JiHdxLdS1nLWhbhq5hgU0yzkBI5JUEhHnhHZyU/IZ79ppU9V23y0WXbzy/no8v3ANNmz1NpGvLdkAALhn9vKyt5/G/qbWUFOqmcbGRjQ2Nla7G1UDe8lonobOQhH5XD5VvVHjdX5J5yxDiHo2ePdFI7Jrn1oEAHh6/mp8Yu+xnGLGm8TSvVBPzV+Fr97yGgBg8S9OTlVXVmD3YMaitfjKkdtVvP2GfLZ7DFN99P1veKnZ3y/FrgHC7rvU2yGLxdn3mskH48h13dSh1nn3XTNjVR3rUUuMCHsEnteM9zlp/8RH2N4ZeHPkMt7q0rY6qsyIXP/8e/jpA16+npbGOryUgRo9DrIW8qqJmmJEejvYREAXkSwm6ihd4tJ16kVGnJtMuhM1oVFJPotN4jyNIFUr+N1j7+CfMxZXrL3GevmrvbK1Dadd+zw+e8OsWO6PKmNNEW0Sl0LOfbeQvSAixhHJql4
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAJLCAYAAADKNQzQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC9dklEQVR4nO2deXwUVbr3f9WRYEAICYuEECAEIexLBAybbA7biCjOiMwdlXEbAXEdCXNHRWdGQJ07CqLjNi7vXJC5w+KKMwoqQhAwrMomIZAQgoIJDYZIMF3vH6dPd3XVObV1VXd1cr6fjwLdtZyqPvXUc55VkmVZhkAgEAg8hy/eAxAIBAIBGyGgBQKBwKMIAS0QCAQeRQhogUAg8ChCQAsEAoFHEQJaIBAIPIoQ0AKBQOBRhIAWCAQCjyIEtEAgEHgUIaAFAoHAo7gqoBcsWICBAweiWbNmaNOmDaZMmYIDBw64eUqBQCCoN7gqoD/77DPMmjULX3zxBT766CNcuHABP/vZz1BdXe3maQUCgaBeIMWyWNLJkyfRpk0bfPbZZxgxYoTh9oFAAMePH0ezZs0gSVIMRigQCATuIssyzp49i3bt2sHn09eRL4rRmAAAfr8fAJCens78/vz58zh//nzo3+Xl5ejRo0dMxiYQCASxpKysDO3bt9fdJmYadCAQwOTJk3H69Gls3LiRuc38+fPx2GOPaT4vKytD8+bN3R6iQCAQuM6ZM2eQlZWF06dPIzU1VXfbmAnou+66C2vXrsXGjRu5bw21Bk0vxO/3CwEtEAjqBWfOnEFqaqopuRYTE8fs2bPx3nvvYcOGDboqfePGjdG4ceNYDEkgEAg8j6sCWpZl3H333Vi9ejU+/fRTZGdnu3k6gUAgqFe4KqBnzZqFZcuW4e2330azZs1w4sQJAEBqaipSUlLcPLVAIBAkPK7aoHmhca+99hpuueUWw/2t2GoEAoEgEfCMDVr0oxUIBAL7iFocAoFA4FGEgBYIBAKPIgS0QCAQeBQhoAUCgcCjCAEtEAgEHkUIaIFAIPAoQkALBAKBRxECWiAQCDyKENACgUDgUYSAFggEAo8iBLRAIBB4FCGgBQKBwKMIAS0QCAQeRQhogUAg8ChCQAscocJfg8LiU6jw18R7KA0W8Ru4R7zubUx6EgrqNyu2lWLeqj0IyIBPAhZc1xs3DOwQ72E1KMRv4B7xvLdCgxZERYW/JjR5ASAgA79f9ZXQ4mKI+A3cI973VghoQVSUnKoOTV5KnSzjyKlz8RlQA8QLv0G0JgCvmWfoeIqOVsX13goThyAqsls1hU9CxCROkiR0atUkfoNqYMT7N4jWBBBLE0KFvwYlp6qR3aopMlLZjatf/KwYC9fuhwxACv6nlNGxvLdCg3YZr2kGTpORmoIF1/VGUrBBcJIk4YnrenEnv8B54vkbRGsCiKUJYcW2UgxduB7TX96CoQvXY8W2Us02L24oxoKgcAbCgtkX7H8d6/ktNGgXaSiOmxsGdsCIrq1x5NQ5dGrVRAjnOBCv30DPvGJmDHb3N6MJq7dnvQhGdG0d2r/CX4OFa/dr9pUBLJnWHy0vaRzz+S0EtEuYmRDxGJOVSW2FjNQUIZjjTDx+g2jNK3b2t6P4GL0IdpVV4Z9flkGWtfv6AOR1SovL/BYmDpfwguNGiZnlnUBglWjNK1b3N2MSYZkV6YtACX0RPPDPnbhmaSH+d0sZ85xzJ+TGTfkQGrRLxNtxo8SL2rwaN7X7eBCL6+GdI9b3MlrzipX9jTRhnnZNXwS/X/UV6mQ59CJ4cMVObDpcyTyXTyLC+c4ROZaux0mEgHYJ3oSIh/CJ1k7oNrGy1cdKcFm9Hjvj4p0jXn6PaM0rZvfXU3yMFBH1i2DIgvVgWDQAAL++ogNmjuoS9+dDCGgX8YrzzEvavJpYafexfAlYuR474+KdI7dtM8+vlKJFT/EpLD5lqIjQF8HsfxRxhTMAXJ/X3hP3TAhol/GC8yxe2rwZzTAW2n0sTTxWrsfuuHjn2HaEn1QR7znoJDzFx4wiQufkp9+c5B5fAtCm+cVuDd8SQkA3EFiT2s0lv1nN0GjJ6oSNNZYmHiurFbvj4p1jYKc0z66UnIal+BgpIiu2lWLuyj2Gx5YBz7zUhIBuQCgntZtLfiuaIe+h2nDwpGM21liaeKysVuyOi3eOvllpnvF7xAuedl3hrzElnAFvvdQkWWZF/nmDM2fOIDU1FX6/H82bN4/3cOoNFf4aDF24XiMYNhaMcuRhLiw+hekvb9F8vvz2K5Cf05I7JvpQAWCOb9XMfFz7fKGtca/YVqoRXG46z5TXoze2aMbFO4fZczck7ntrO1bvrNB8ntE8GVXnLuDHn8ikisXcsCLXhAbdAHFiya9nZrCjGSq1e56zJxoba6wdtmZ9D9GMi3cOL/g99Kjw1+DLI5WQJAl5Hd1NANlVVoWtRyqxmRNKd3Gji7D/T1d59qUmBHQDJNolv5GZIVqnpFs2Vq8JLuVLjreyqG+s2FaKgpV7QhEUEoCFU92JqJn5v0X4YM8J3W1+3icDgPfmBkVkEnqAWBdUiib7y2xxmxsGdsDGglFYfvsV2FgwyvIDeOuwbE2BGmpj9WJhJqu/YUPM7Kzw10QIZ4A45Oat2uP43J9lQjg3TU7CA+NyHT2v0wgNOs7EK7HA7tLainnEjlaivB8SgDtGZGPG0OzQcbwSW67ETmJKLML+3IzSsXPsklPVzNjjgOxs1MS4v36GA9/+wPxuWE5LfHf2R4zr2RbTr+iIwuJTns5eFQI6jsQ7BduOADUyj0QjFNT3Qwbw6udHMGNodtTjdgs7v2Eswv7cfPHbPXZ2q6aa2soAOYZTURPz3/mKK5wB4Hfju6FvVlpoBeP1SpPCxBFHvFZQyQx65pFol+1O3A8rpgarZgnW9nbGrFe4xwncrLEczbEzUlOwcGpvKC9dCgpHdRSKHZPf/31ZitcLj3K/n9i7LfpmpcW9jZUVhAYdR1jaqA9Ak+To3ptu15zgJb1Euxpw23lpd1ve9iO6tkZlda3ljhtuZ3a6qaFHe2w6d4qOVEGSgAGqKA4zvwtrfg9Z8DGO+89zzzupd1ss/VWeI9cQS4SAjiPqBxUAAgCufb7Q9pIrVjZttZnBiUlvRXCpH1LWC2Leyj1o2vgiTSiX1ZcJa/uClXsgBV8myrZIPsCUsHXTlu5mYo4Tx85ITcHP+4YzWelnZn4X9fz+9eCOePOLo7p1NS5t3jgknJ26hlghBHScuWFgB+S2bYYpzxeGioXbtUVHo8VGq3XvKfdrPrMz6c0ILtZLKCu9ieYFEQAwe9kOzYvK6suEtb0MhH6viK9Upgs93LKlu6mhO3Vss7+h8ndhze83vuCbNCinztaiwl8TUTApUTIuhYD2ANW1dZpODnaWXHa12Gi17gp/DRYxWgXdMLC96WMo0RNcvJfQqpn5Gq2Ion5RWdWgWNvzcNrRa/fF6aaGHu2xeb/hSzcNYG5PTX6s+W0G1jPgxWggFsJJ6AH0nEZWHCZ2nE9OOEx4D86yrWWmnIVWrpH3EjpXG4hwXqpROu6sxoGrt/dBX1E249g0c83ROl0zUlOQn9PSNS3d7rF5v2EJ556dqw0AYK/SzMB7Bty8P04hNGgPYLVgkNXjqO2vSo3MCduxnoYZbT1k9Xj1tN/8nJYhB9Sct3boashWNShqitp2pAoDO6Vh/4mzEb4DJXZ66o3o2trQpu5mCGYsu7DYyRTlrdJYJEkAIHnefGEGIaA9glpgAJEFg8w+oHqChycYnHD6qJ2dSpQapRUhxBPeei8h6oCqrv3J0MZoxQbMGsvGglE4cuocdh87jSc/PGBKIDAdjqv2AHLQyWjSHusksegAo4SnSOhV42PVZ2GRmnIRdj06zrO1Nawiqtl5FDsV4fTQq2C34eBJw4pqZh7KCn8Nth+twuxlOzRhZw9N6IZFa/drnEK8a+zUqolmvD4JWD1zSCiW1egBdOoh3VVWhSlLCzXXpKyiZ/ZcvN9VSbSV+5QY/W5WKxs6GSVE50tAlnF5p3QA5AXeNDkJ52oDoXu5q6wKt7y6FVU//qR
      "text/plain": [
       "<Figure size 400x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def f(x):\n",
    "    return np.exp(-np.abs(x**3))\n",
    "\n",
    "N = 10**4 # generate N samples\n",
    "\n",
    "x = np.zeros(N)\n",
    "x[0] = 0  # starting point\n",
    "\n",
    "count = 0\n",
    "\n",
    "for i in range(N-1):\n",
    "    \n",
    "    # Proposal\n",
    "    x_try = x[i] + rng.normal(scale = 2)   \n",
    "\n",
    "    # Accept the move or stay where we are\n",
    "    u = rng.uniform()\n",
    "    if u <= f(x_try) / f(x[i]):\n",
    "        x[i+1] = x_try\n",
    "        count = count + 1\n",
    "    else:\n",
    "        x[i+1] = x[i]\n",
    "\n",
    "print(\"Acceptance fraction = %g\" % (count/N,))\n",
    "\n",
    "# Also do uniform sampling for comparison\n",
    "xmax = 4\n",
    "xu = -xmax + 2*xmax*rng.uniform(size = 10*N)    # Generate x values between +-xmax\n",
    "y = rng.uniform(size = 10*N)\n",
    "xu = xu[np.where(y <= f(xu))]  \n",
    "xu = xu[:N]\n",
    "        \n",
    "# Plot the distribution        \n",
    "counts, bins = np.histogram(x, bins=100, density = True)\n",
    "plt.stairs(counts, bins, label = 'Metropolis')    \n",
    "\n",
    "counts, bins = np.histogram(xu, bins=100, density = True)\n",
    "plt.stairs(counts, bins, label = 'Rejection method')    \n",
    "\n",
    "norm,_ = scipy.integrate.quad(f,-np.inf,np.inf)\n",
    "xx = np.linspace(min(x),max(x),1000)\n",
    "plt.plot(xx, f(xx)/norm, 'k:')\n",
    "\n",
    "plt.yscale('log')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# Plot the first 1000 samples to see the time series\n",
    "plt.clf()\n",
    "t = np.arange(N)\n",
    "plt.plot(t[:1000], x[:1000])\n",
    "plt.plot(t[:1000], xu[:1000] - 3.5)\n",
    "plt.show()\n",
    "\n",
    "# Plot the point-to-point correlation\n",
    "plt.clf()\n",
    "plt.figure(figsize=(4,7))\n",
    "plt.plot(x[:-1], x[1:], '.')\n",
    "plt.plot(xu[:-1], xu[1:] - 3.5, '.')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49118673",
   "metadata": {},
   "source": [
    "Notes:\n",
    "\n",
    "- We started at $x=0$ which is right in the middle of the high probability region. If instead you start at a point that has a low probability (e.g. try $x=-5$ as initial condition in line 7), you will see that the chain has a \"burn in\" period in which it moves towards the higher probability regions. The part of the chain should be discarded. \n",
    "- The time series for the MCMC shows clear correlations, whereas the rejection method samples are uncorrelated. This can be seen in the plot of $x_i$ vs $x_{i+1}$ -- in particular you can see the line of points where the proposed jump was not accepted.\n",
    "- Try changing the `scale` in the normal distribution we're using for the jump and see how the acceptance ratio and the correlations behave. You want to choose a jump size that is large enough to explore the probability distribution, but small enough to get a reasonable acceptance rate for the jumps. An acceptance rate of around $0.2$--$0.3$ is good to aim for."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5503d1cb",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}