lecture-2-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JplqrIfpRq8l"
   },
   "source": [
    "# Lecture 2-2: Random numbers and Markov chain Monte Carlo\n",
    "\n",
    "### Introduction\n",
    "\n",
    "Random numbers can have many uses in programming and in physics. Here, we'll start to explore how random numbers can be used for a method of stochastic simulation called Markov chain Monte Carlo. Let's begin by first looking at how to generate random numbers.\n",
    "\n",
    "\n",
    "### Generating random numbers\n",
    "\n",
    "Python provides a standard library, `random`, for random number generation [here](https://docs.python.org/3/library/random.html). Numpy has a more comprehensive library [here](https://docs.scipy.org/doc/numpy-1.15.0/reference/routines.random.html). For the examples below, we'll use `numpy.random`. \n",
    "\n",
    "Random number generators in Python are not really random. They use values called **seeds** to determine the initial state of the random number generator. Given a specific seed as input, the output of a random number generator is totally deterministic. However, the numbers that appear are very close to random in a statistical sense. (For much more detail, see the related [Wikipedia page](https://en.wikipedia.org/wiki/Random_number_generation).)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 35
    },
    "colab_type": "code",
    "id": "Hvbde_XGRq8n",
    "outputId": "6b54fbe7-ad56-4c9d-851f-beba8658eecf"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as rng\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# First, let's define a new random number generator and set the seed\n",
    "\n",
    "r = rng.RandomState(0)\n",
    "\n",
    "# We can now use it to generate random numbers\n",
    "\n",
    "n = r.rand()\n",
    "print(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The sequence of random numbers generated by the same generator \n",
    "# with the same seed is always the same\n",
    "\n",
    "print(np.array([n] + list(r.rand(5))))\n",
    "\n",
    "r2 = rng.RandomState(0)\n",
    "print(r2.rand(6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "x_QAvAsdRq8q"
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uEdRJuZdRq8u"
   },
   "source": [
    "### Using random numbers to \"measure\" the state of a stochastic system\n",
    "\n",
    "As we just discussed, the energy of a particle with magnetic moment $\\mathbf{\\mu}$ in a magnetic field $\\mathbf{B}$ is given by\n",
    "\n",
    "$$ E = -\\mathbf{\\mu}\\cdot \\mathbf{B}\\,.$$\n",
    "\n",
    "Let's consider a statistical model of a particle in an external magnetic field. To keep things simple, we'll assume that the magnetic field is oriented in a fixed direction, and that the magnetic moment of the particle can be either aligned or anti-aligned with the magnetic field. In statistical physics models, these particles are often called **spins**.\n",
    "\n",
    "Mathematically, we could describe the orientation of the spin with a variable $\\sigma \\in \\{-1, 1\\}$. Let's assume that the magnetic field is oriented in the positive direction. The energy of the spin is then\n",
    "\n",
    "$$ E(\\sigma) = -\\epsilon \\sigma\\,, $$\n",
    "\n",
    "where $\\epsilon = \\mu B$. If we maintain the system with the spin at a constant temperature $T$, then the Gibbs distribution for the spin states $\\sigma$ is\n",
    "\n",
    "$$ P(\\sigma) = \\frac{e^{-E(\\sigma)/T}}{Z} \\,.$$\n",
    "\n",
    "Here we've chosen units such that Boltzmann's constant $k_B=1$. In other words, each time we observe the system, there's a probability $P(+1)$ to find the spin in the $+$ state, and a probability $P(-1)$ to find the spin in the $-$ state.\n",
    "\n",
    "**If we observed the system 10 times, what would the average value of the spin be?**\n",
    "\n",
    "To begin answering this question, let's write down the probabilities $P(+1)$ and $P(-1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 35
    },
    "colab_type": "code",
    "id": "4vc_DkLORq8w",
    "outputId": "d7c358b5-d672-4bd6-be29-259278dceacd"
   },
   "outputs": [],
   "source": [
    "eps = 1 # magnitude of the energy epsilon = mu * B\n",
    "T = 1   # temperature of the system\n",
    "\n",
    "Z =  # fill this in\n",
    "p_plus =  # fill this in\n",
    "p_minus =  # fill this in\n",
    "\n",
    "print('P(+) =', p_plus)\n",
    "print('P(-) =', p_minus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "vnSjhyuLRq8z"
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IkPigVt-Rq81"
   },
   "source": [
    "### Conditional expressions\n",
    "\n",
    "Each time we measure the spin, we know the probability to find it in the plus or the minus state. If we want to use these probabilities to set the value of the spin each time we measure it, we can use a **conditional expression**. \n",
    "\n",
    "The most common conditional expressions are if/else statements, which work in a very intuitive way! Consider the example below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "cuAZQeZIRq82"
   },
   "outputs": [],
   "source": [
    "def check_size(x):\n",
    "    \"\"\" This function checks to see whether or not a number is larger than one.\"\"\"\n",
    "    if x>1:\n",
    "        print('The number is larger than one')\n",
    "    else:\n",
    "        print('The number is not larger than one')\n",
    "        \n",
    "check_size(0)\n",
    "check_size(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "KWybzNzpRq85"
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_WG6B7vFRq87"
   },
   "source": [
    "How does this work? Python takes everything in between the `if` and the `:` and evaluates it as a **boolean** (True/False) statement. If the expression evaluates to `True`, then Python executes the next block of code. If it evaluates to `False`, then Python executes the block of code that follows `else:`.\n",
    "\n",
    "Here are a list of some common comparison operators that you can use to generate True/False expressions:\n",
    "- `>` greater than  \n",
    "- `<` less than  \n",
    "- `==` equal to (**note that there are two equal signs!**)  \n",
    "- `!=` not equal to  \n",
    "- `>=` greater than or equal to  \n",
    "- `<=` less than or equal to  \n",
    "\n",
    "These expressions can also be modified in multiple ways, but we won't go into the details at the moment. Let's practice with a few examples in the space below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "XejD_xxNRq88"
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "X7Yf7gUjRq8_"
   },
   "source": [
    "### Bringing random numbers and conditional expressions together\n",
    "\n",
    "Behold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 283
    },
    "colab_type": "code",
    "id": "5C0lf0s6WrzK",
    "outputId": "e846eb84-5b72-4362-f2ac-7145ee193dc4"
   },
   "outputs": [],
   "source": [
    "x = r.rand()\n",
    "if x<0.5:\n",
    "    print('%lf is <0.5' % x)\n",
    "else:\n",
    "    print('%lf is >=0.5' % x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "i0WebzdYRq9A"
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2IGCgFcHRq9C"
   },
   "source": [
    "We can use conditional expressions and probabilities to **simulate** a number of measurements of the spin, allowing us to estimate the magnetization. For example, let's imagine that we measured the spin 10 times. To find the average magnetization, we could use the following block of code.\n",
    "\n",
    "We can also consider what happens as we increase the number of measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_measure = 10\n",
    "m_avg = 0\n",
    "for i in range(n_measure):\n",
    "    if r.rand()<p_plus:\n",
    "        m_avg = m_avg + 1\n",
    "    else:\n",
    "        m_avg = m_avg - 1\n",
    "\n",
    "m_avg = m_avg / n_measure # normalize by dividing by number of measurements\n",
    "print('The average magnetization is', m_avg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "CBMgFleCRq9D"
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Ising model\n",
    "\n",
    "As a simple model of a ferromagnetic material, we can consider a set of spins arranged on a lattice. The energy of a pair of spins will be lower if they are aligned in the same direction. To make our analysis simple, we'll assume that only **nearest neighbor** spins are coupled together. Spins that aren't nearest neighbors could also interact, but these interactions will be much weaker due to the greater distance between them.\n",
    "\n",
    "Mathematically, we can write the energy of a configuration $\\underline{\\sigma}=\\{\\sigma_1, \\sigma_2, \\ldots, \\sigma_N\\}$ of spins $\\sigma_i\\in\\{-1,1\\}$ in a one-dimensional lattice as\n",
    "\n",
    "$$ \n",
    "E(\\underline{\\sigma}) = -\\sum_{i=1}^{N-1} J \\sigma_i\\,\\sigma_{i+1}\\,.\n",
    "$$\n",
    "\n",
    "The coupling $J$ tells us how strong the interaction is between two neighboring spins.\n",
    "\n",
    "To start, let's use a list to keep track of all the Ising spins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 20      # number of spins in our system\n",
    "spins = []  # the list of spins\n",
    "\n",
    "# Fill in code here to initialize the list of spins\n",
    "# Can we make the spins alternate + - + - ...?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's write a function to compute the energy of a configuration. For the moment, let's set $J=1$ for simplicity. The input to the function should be the configuration of spins, and the output should be the energy of the configuration. **Note: we need to be careful about the boundaries.** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def compute_energy(spins):\n",
    "    \"\"\" This function computes the energy of a spin configuration assuming J=1. \"\"\"\n",
    "    \n",
    "    energy = 0\n",
    "    # Fill in code here to add up the energy\n",
    "    \n",
    "    return energy\n",
    "\n",
    "# Let's check to make sure that the energy function works\n",
    "\n",
    "print('E( -1 -1 ) =', compute_energy([-1, -1]), '\\t(expected -1)')\n",
    "print('E( +1 +1 ) =', compute_energy([ 1,  1]), '\\t(expected -1)')\n",
    "print('E( -1 +1 ) = ', compute_energy([-1,  1]), '\\t(expected  1)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "lecture-3.ipynb",
   "provenance": []
  },
  "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.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}