Metropolis_Monte_Carlo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 04:  Metropolis Monte Carlo for Classical Systems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "try: plt.style.use('./mc_notebook.mplstyle')\n",
    "except: pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluating Phase Space Integrals\n",
    "\n",
    "We have seen that Monte Carlo methods are an ideal way to sample high-dimensional integrals.  We now ask if this can be applied to the phase space integrals that appear when computing equilbrium observables in classical statistical mechanics for $N$ interacting particles.  \n",
    "\n",
    "Starting with the general $N$-body Hamiltonian:\n",
    "\n",
    "\\begin{equation}\n",
    "H = -\\frac{\\hbar^2}{2m} \\sum_{i=1}^{N} \\nabla_i^2 + \\frac{1}{2}\\sum_{i,j} V_{ij}\n",
    "\\end{equation}\n",
    "\n",
    "we are interested in the expectation value of some observable $\\mathcal{O}$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle \\mathcal{O} \\rangle = \\frac{1}{\\mathcal{Z}} \\int d^{d} x_1 \\cdots \\int d^{d} x_N\\, \\mathcal{O}(\\mathbf{x}) \\mathrm{e}^{-\\beta V(\\mathbf{x})}\n",
    "\\end{equation}\n",
    "\n",
    "where the partition function is given by $(\\beta = 1/k_{\\rm B} T)$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathcal{Z} = \\int d^{d} x_1 \\cdots \\int d^{d} x_N\\,  \\mathrm{e}^{-\\beta V(\\mathbf{x})}.\n",
    "\\end{equation}\n",
    "\n",
    "We can immediatly identify this with our previous importance sampling scheme:\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle \\mathcal{O} \\rangle = \\int d^{d} x_1 \\cdots \\int d^{d} x_N\\, \\mathcal{O}(\\mathbf{x}) \\pi( \\mathbf{x}) \\simeq \\frac{1}{M} \\sum_m \\mathcal{O}(\\mathbf{x}_m)\n",
    "\\end{equation}\n",
    "\n",
    "where \n",
    "\n",
    "\\begin{equation}\n",
    "\\pi(\\mathbf{x}) = \\frac{\\mathrm{e}^{-\\beta V(\\mathbf{x})}}{\\mathcal{Z}}\n",
    "\\end{equation}\n",
    "\n",
    "is the Boltzmann weight of configuraiton $\\mathbf{x}$. \n",
    "\n",
    "### Key Problem: \n",
    "\n",
    "How do we effeciently sample the correct Boltzmann probability distribution?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solved by:\n",
    "\n",
    "[N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, *Equation of State Calculations by Fast Computing Machines*, J. Chem. Phys. **21**, 1087 (1953).](https://aip.scitation.org/doi/abs/10.1063/1.1699114)\n",
    "\n",
    "<img src='./metropolis.png', width=600px>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Markov Chain Monte Carlo\n",
    "\n",
    "Since $\\pi(\\mathbf{x})$ is an extremely complicated $dN$-dimensional probability distribution there is no hope of direct sampling.  However, we can recall our strategy to the *pebble throwing game* when we imagined that the entire area couldn't be reached with a single toss: **Markov chain sampling**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"548pt\" height=\"44pt\"\n",
       " viewBox=\"0.00 0.00 548.00 44.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 40)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-40 544,-40 544,4 -4,4\"/>\n",
       "<!-- x0 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>x0</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"18\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"12.5\" y=\"-15.3\" font-family=\"Times,serif\" font-weight=\"bold\" font-style=\"italic\" font-size=\"14.00\" fill=\"#000000\">x</text>\n",
       "<text text-anchor=\"start\" x=\"18.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">0</text>\n",
       "</g>\n",
       "<!-- x1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>x1</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"90\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"84.5\" y=\"-15.3\" font-family=\"Times,serif\" font-weight=\"bold\" font-style=\"italic\" font-size=\"14.00\" fill=\"#000000\">x</text>\n",
       "<text text-anchor=\"start\" x=\"90.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">1</text>\n",
       "</g>\n",
       "<!-- x0&#45;&gt;x1 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>x0&#45;&gt;x1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M36.1686,-18C43.869,-18 53.0257,-18 61.5834,-18\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"61.5868,-21.5001 71.5867,-18 61.5867,-14.5001 61.5868,-21.5001\"/>\n",
       "</g>\n",
       "<!-- x2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>x2</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"162\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"156.5\" y=\"-15.3\" font-family=\"Times,serif\" font-weight=\"bold\" font-style=\"italic\" font-size=\"14.00\" fill=\"#000000\">x</text>\n",
       "<text text-anchor=\"start\" x=\"162.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">2</text>\n",
       "</g>\n",
       "<!-- x1&#45;&gt;x2 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>x1&#45;&gt;x2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M108.1686,-18C115.869,-18 125.0257,-18 133.5834,-18\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"133.5868,-21.5001 143.5867,-18 133.5867,-14.5001 133.5868,-21.5001\"/>\n",
       "</g>\n",
       "<!-- x3 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>x3</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"234\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"228.5\" y=\"-15.3\" font-family=\"Times,serif\" font-weight=\"bold\" font-style=\"italic\" font-size=\"14.00\" fill=\"#000000\">x</text>\n",
       "<text text-anchor=\"start\" x=\"234.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">3</text>\n",
       "</g>\n",
       "<!-- x2&#45;&gt;x3 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>x2&#45;&gt;x3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M180.1686,-18C187.869,-18 197.0257,-18 205.5834,-18\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"205.5868,-21.5001 215.5867,-18 205.5867,-14.5001 205.5868,-21.5001\"/>\n",
       "</g>\n",
       "<!-- x4 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>x4</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"306\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"300.5\" y=\"-15.3\" font-family=\"Times,serif\" font-weight=\"bold\" font-style=\"italic\" font-size=\"14.00\" fill=\"#000000\">x</text>\n",
       "<text text-anchor=\"start\" x=\"306.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">4</text>\n",
       "</g>\n",
       "<!-- x3&#45;&gt;x4 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>x3&#45;&gt;x4</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M252.1686,-18C259.869,-18 269.0257,-18 277.5834,-18\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"277.5868,-21.5001 287.5867,-18 277.5867,-14.5001 277.5868,-21.5001\"/>\n",
       "</g>\n",
       "<!-- x5 -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>x5</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"378\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"372.5\" y=\"-15.3\" font-family=\"Times,serif\" font-weight=\"bold\" font-style=\"italic\" font-size=\"14.00\" fill=\"#000000\">x</text>\n",
       "<text text-anchor=\"start\" x=\"378.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">5</text>\n",
       "</g>\n",
       "<!-- x4&#45;&gt;x5 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>x4&#45;&gt;x5</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M324.1686,-18C331.869,-18 341.0257,-18 349.5834,-18\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"349.5868,-21.5001 359.5867,-18 349.5867,-14.5001 349.5868,-21.5001\"/>\n",
       "</g>\n",
       "<!-- x6 -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>x6</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"450\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"444.5\" y=\"-15.3\" font-family=\"Times,serif\" font-weight=\"bold\" font-style=\"italic\" font-size=\"14.00\" fill=\"#000000\">x</text>\n",
       "<text text-anchor=\"start\" x=\"450.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">6</text>\n",
       "</g>\n",
       "<!-- x5&#45;&gt;x6 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>x5&#45;&gt;x6</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M396.1686,-18C403.869,-18 413.0257,-18 421.5834,-18\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"421.5868,-21.5001 431.5867,-18 421.5867,-14.5001 421.5868,-21.5001\"/>\n",
       "</g>\n",
       "<!-- x7 -->\n",
       "<g id=\"node8\" class=\"node\">\n",
       "<title>x7</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"522\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"516.5\" y=\"-15.3\" font-family=\"Times,serif\" font-weight=\"bold\" font-style=\"italic\" font-size=\"14.00\" fill=\"#000000\">x</text>\n",
       "<text text-anchor=\"start\" x=\"522.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">7</text>\n",
       "</g>\n",
       "<!-- x6&#45;&gt;x7 -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>x6&#45;&gt;x7</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M468.1686,-18C475.869,-18 485.0257,-18 493.5834,-18\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"493.5868,-21.5001 503.5867,-18 493.5867,-14.5001 493.5868,-21.5001\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x10aff4898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display\n",
    "\n",
    "try:\n",
    "    import graphviz\n",
    "    markov_chain = graphviz.Digraph()\n",
    "    markov_chain.attr(rankdir='LR')\n",
    "    for i in range(2**3):\n",
    "        markov_chain.node(name=\"x%d\"%i, label=\"<<i><b>x</b></i><sub>%d</sub>>\"%i, shape='circle')\n",
    "\n",
    "    for i in range(2**3-1):\n",
    "        markov_chain.edge(\"x%d\"%i, \"x%d\"%(i+1)) \n",
    "                          \n",
    "    display(markov_chain)\n",
    "    \n",
    "except:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal is to build a Markov Chain consisting of a sequence of valid configuraitons $\\mathbf{x}_n$ such that $\\pi(\\mathbf{x}_n)$ reflects the stationary probability of the system being in configuraiotn $\\mathbf{x}_n$.\n",
    "\n",
    "We achieve this via a set of *transition probabilities* $p_{n,m} = p(\\mathbf{x}_n \\to \\mathbf{x}_{m})$ which encode the liklihood of moving from configuration $\\mathbf{x}_n \\to \\mathbf{x}_{m}$.  We require:\n",
    "\n",
    "### 1. Normalization\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{m} p_{n,m} = 1\n",
    "\\end{equation}\n",
    "\n",
    "### 2. Detailed Balance\n",
    "\n",
    "The probability of configuration $\\mathbf{x}_n$, $\\pi_n \\equiv \\pi(\\mathbf{x}_n)$, is generated via all possible ways of obtaining it (i.e. the target distribution is stationary):\n",
    "\n",
    "\\begin{align*}\n",
    "\\pi_n &= \\sum_{m} \\pi_{m} p_{m,n} \\newline\n",
    "\\pi_n(1-p_{n,n}) &= \\sum_{m\\ne n} \\pi_{m} p_{m,n} .\n",
    "\\end{align*}\n",
    "\n",
    "Now utilizing the normalization condition:\n",
    "\\begin{equation}\n",
    "\\sum_{m} p_{n,m} = 1 \\Rightarrow 1-p_{nn} = \\sum_{m\\ne n} p_{n,m}\n",
    "\\end{equation}\n",
    "\n",
    "above we can write:\n",
    "\n",
    "\\begin{equation}\n",
    "\\pi_n \\sum_{m \\ne n} p_{n,m} = \\sum_{m\\ne n} \\pi_{m} p_{m,n}.\n",
    "\\end{equation}\n",
    "\n",
    "One obvious solution to this equation which is sufficient but not necessary is to equate each term individually (this is *detailed balance*):\n",
    "\n",
    "\\begin{equation}\n",
    "\\boxed{\\pi_n  p_{n,m} = \\pi_{m} p_{m,n}}.\n",
    "\\end{equation}\n",
    "\n",
    "### 3. Ergodicity\n",
    "\n",
    "Any configuraiton is reachable from any other, i.e. $\\forall m,n,\\ \\exists a \\, : (p^a)_{m,n} \\ne 0$. The graph corresponding to the Markov chain is connected.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metropolis Algorithm\n",
    "\n",
    "Teller's solution to generating a Markov Chain satisfying these conditions was to use the rejection sampling scheme we have already employed in combination with choosing the probability to move from configuration $\\mathbf{x}_n \\to \\mathbf{x}_m$ to be composite:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathcal{P}(\\mathbf{x}_n \\to \\mathbf{x}_m) = \\mathcal{A}(\\mathbf{x}_n \\to \\mathbf{x}_m) p(\\mathbf{x}_n \\to \\mathbf{x}_m)\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathcal{A}(\\mathbf{x}_n \\to \\mathbf{x}_m)$ is an *a priori* proposal probability and $p(\\mathbf{x}_n \\to \\mathbf{x}_m)$ is the acceptance probability.\n",
    "\n",
    "Our *detailed balance* condition requires:\n",
    "\n",
    "\\begin{equation}\n",
    "\\pi(\\mathbf{x}_n) \\mathcal{P}(\\mathbf{x}_n \\to \\mathbf{x}_m) =\\pi(\\mathbf{x}_m) \\mathcal{P}(\\mathbf{x}_m \\to \\mathbf{x}_n) \n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{p(\\mathbf{x}_n \\to \\mathbf{x}_m)}{p(\\mathbf{x}_m \\to \\mathbf{x}_n)} = \\frac{\\pi(\\mathbf{x}_m)\\mathcal{A}(\\mathbf{x}_m \\to \\mathbf{x}_n)}{\\pi(\\mathbf{x}_n)\\mathcal{A}(\\mathbf{x}_n \\to \\mathbf{x}_m)}\n",
    "\\end{equation}\n",
    "\n",
    "which is satified by the choice:\n",
    "\n",
    "\\begin{equation}\n",
    "p(\\mathbf{x}_n \\to \\mathbf{x}_m) = \\text{min}\\left[1,\\frac{\\pi(\\mathbf{x}_m)\\mathcal{A}(\\mathbf{x}_m \\to \\mathbf{x}_n)}{\\pi(\\mathbf{x}_n)\\mathcal{A}(\\mathbf{x}_n \\to \\mathbf{x}_m)}\\right].\n",
    "\\end{equation}\n",
    "\n",
    "The idea is to :\n",
    "\n",
    "1. Starting from configuration $\\mathbf{x}_n$.\n",
    "2. Propose some change to the configuration with *a priori* proposal probability $\\mathcal{A}(\\mathbf{x}_n \\to \\mathbf{x}_m)$.\n",
    "3. Accept the change with probability $p(\\mathbf{x}_n \\to \\mathbf{x}_m)$ given by the expression above where $\\pi(\\mathbf{x})$ is known (i.e. the Boltzmann weight) and $\\mathcal{A}(\\mathbf{x}_n \\to \\mathbf{x}_m)$ can be chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"298pt\" height=\"46pt\"\n",
       " viewBox=\"0.00 0.00 297.70 45.70\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 41.6967)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-41.6967 293.6967,-41.6967 293.6967,4 -4,4\"/>\n",
       "<!-- x0 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>x0</title>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-15.1483\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">⋯</text>\n",
       "</g>\n",
       "<!-- x1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>x1</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"108\" cy=\"-18.8483\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"102.5\" y=\"-16.1483\" font-family=\"Times,serif\" font-weight=\"bold\" font-style=\"italic\" font-size=\"14.00\" fill=\"#000000\">x</text>\n",
       "<text text-anchor=\"start\" x=\"108.5\" y=\"-16.1483\" font-family=\"Times,serif\" font-style=\"italic\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">n</text>\n",
       "</g>\n",
       "<!-- x0&#45;&gt;x1 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>x0&#45;&gt;x1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M54.4229,-18.8483C62.4647,-18.8483 71.2772,-18.8483 79.3909,-18.8483\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"79.5312,-22.3484 89.5312,-18.8483 79.5311,-15.3484 79.5312,-22.3484\"/>\n",
       "</g>\n",
       "<!-- x2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>x2</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"180.8483\" cy=\"-18.8483\" rx=\"18.6979\" ry=\"18.6979\"/>\n",
       "<text text-anchor=\"start\" x=\"174.3483\" y=\"-16.1483\" font-family=\"Times,serif\" font-weight=\"bold\" font-style=\"italic\" font-size=\"14.00\" fill=\"#000000\">x</text>\n",
       "<text text-anchor=\"start\" x=\"180.3483\" y=\"-16.1483\" font-family=\"Times,serif\" font-style=\"italic\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">m</text>\n",
       "</g>\n",
       "<!-- x1&#45;&gt;x2 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>x1&#45;&gt;x2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M126.0075,-18.8483C133.7357,-18.8483 142.965,-18.8483 151.628,-18.8483\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"151.7752,-22.3484 161.7752,-18.8483 151.7752,-15.3484 151.7752,-22.3484\"/>\n",
       "</g>\n",
       "<!-- x3 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>x3</title>\n",
       "<text text-anchor=\"middle\" x=\"262.6967\" y=\"-15.1483\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">⋯</text>\n",
       "</g>\n",
       "<!-- x2&#45;&gt;x3 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>x2&#45;&gt;x3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M199.8317,-18.8483C207.4348,-18.8483 216.4724,-18.8483 225.2905,-18.8483\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"225.4424,-22.3484 235.4423,-18.8483 225.4423,-15.3484 225.4424,-22.3484\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x10e39b390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "try:\n",
    "    markov_chain_snippet = graphviz.Digraph()\n",
    "    markov_chain_snippet.attr(rankdir='LR')\n",
    "\n",
    "    markov_chain_snippet.node(name=\"x%d\"%0, label=\"⋯\", shape='none')\n",
    "    markov_chain_snippet.node(name=\"x%d\"%1, label=\"<<i><b>x</b></i><sub><i>n</i></sub>>\", shape='circle')\n",
    "    markov_chain_snippet.node(name=\"x%d\"%2, label=\"<<i><b>x</b></i><sub><i>m</i></sub>>\", shape='circle')\n",
    "    markov_chain_snippet.node(name=\"x%d\"%3, label=\"⋯\", shape='none')\n",
    "\n",
    "    for i in range(3):\n",
    "        markov_chain_snippet.edge(\"x%d\"%i, \"x%d\"%(i+1)) \n",
    "\n",
    "    display(markov_chain_snippet)\n",
    "except:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How do we sample the $p$ distribution?\n",
    "\n",
    "For a given configuration $\\mathbf{x}_n$, some update scheme that changes to configuration $\\mathbf{x}_m$, and known functions $\\pi$ and $\\mathcal{A}$, we can sample $p$ via rejection:\n",
    "\n",
    "1. Compute $p$ as a real number using the formula above.\n",
    "2. Generate a uniform random number $r \\in \\mathcal{U}_{[0,1)}$\n",
    "3. If $r < p$ accept the configuration change to $\\mathbf{x}_m$; otherwise reject and return to configuration $\\mathrm{x}_n$.\n",
    "\n",
    "#### Let's test this approach for the 2D classical Ising Model!"
   ]
  }
 ],
 "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}