HW4.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:90a2d2af48b4e3c63543de21e3b8ff702201ee0df603133d9559f20b05bcdc93"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "HW 4"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "_ _ _ _ _\n",
      "\n",
      "Pavlos Protopapas <br>\n",
      "Handed out: Thursday, February 27th, 2014<br>\n",
      "Due: 11.59 P.M. Wednesday March 5th, 2014\n",
      "\n",
      "**Instructions**:\n",
      "\n",
      "+ Upload your answers in an ipython notebook to the dropbox.\n",
      "\n",
      "+ Your individual submissions use the following filenames: AM207_YOURNAME_HM4.ipynb\n",
      "\n",
      "+ Your code should be in code cells as part of your ipython notebook. Do not use a different language (or format) unless you get permission from the TFs. If you use any special libraries you must include them with your code (program should run as is). \n",
      "\n",
      "+ If you have multiple files (e.g. you've added code files and images) create a tarball for all files in a single file and name it: AM207_YOURNAME_HM4.tar or AM207_YOURNAME_HM4.zip\n",
      "\n",
      "_ _ _ _ _"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "\\newcommand{\\fw}[1]{\\tilde{f}(#1)}\n",
      "\\newcommand{\\hw}[1]{\\tilde{h}(#1)}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Problem 1. Ising, you sing, we all sing (but only when close)\n",
      "\n",
      "The true power of the Monte Carlo method becomes apparent when we deal with integrals over many\n",
      "dimensions, such as those often found in statistical physics. In this exercise we consider the 2-D Ising\n",
      "model, which is a very idealized representation of a ferromagnetic metal. The model consists of an LxL lattice of spins \n",
      "$s_{ij}$ (with $i, j = 1, . . . , L$) each of which can take only two values, up (s$_{ij}$ = 1) and down (s$_{ij}$ = -1). The total energy\n",
      "for the system is taken to be\n",
      "\n",
      "$$\n",
      "E = -\\sum_{ij} \\frac{1}{2} s_{ij}  \\left[ s_{\\rm{up}(ij)} + s_{\\rm{down}(ij)} + s_{\\rm{left}(ij)} + s_{\\rm{right}(ij)} \\right]  \n",
      "$$\n",
      "\n",
      "where s$_{\\rm{up}(ij)}$ denotes the spin immediately above s$_{ij}$ in the grid, s$_{\\rm{right}(ij)}$ the spin immediately to its\n",
      "right, and so on. We see that the configurations where nearby spins are aligned are energetically preferred.\n",
      "When the system is in thermal equilibrium at the temperature T, the probability of finding it in a given\n",
      "configuration $\\{s\\}=(s_{11}, s_{12}, \\ldots)$ of the spins is given by the Boltzmann factor,\n",
      "\n",
      "$$W\\{{s} \\} = Z^{-1} e^{ -\\beta E[ \\{ s \\} ] } $$\n",
      "\n",
      "where Z is given by the sum of the exponential over all possible spin configurations, and where $\\beta$ is the\n",
      "inverse temperature (measured in units of the Boltzmann constant $k_B$). The values at thermal equilibrium\n",
      "of all the macroscopic quantities for the system (such as the magnetization $M =\\sum_{ij} s_{ij}$ ) are found by\n",
      "averaging over all spin configurations, weighted by W.  For instance,\n",
      "\n",
      "\\begin{equation} \n",
      "\\left< M \\right> = \\sum_{s_{11}=\\pm 1}   \\sum_{s_{12}=\\pm 1} \\ldots \\sum_{s_{LL}=\\pm 1} \\left[ W\\{s\\}    \\times  \\sum_{ij} s_{ij} \\right]  \n",
      "\\end{equation} \n",
      "\n",
      "  where $\\left<M\\right>$ denotes such an ensemble average. The above  equation is effectively \n",
      "  is  completely out of the reach of standard techniques, but not of Monte Carlo.\n",
      "\n",
      "The straightforward application of the basic Monte Carlo method is not however the best way to compute\n",
      "these averages, because it samples all spin configurations uniformly;  The Metropolis algorithm\n",
      "generates a sequence of configurations that have approximately the\n",
      "required probability distribution. This is how it works:\n",
      "\n",
      "Write a program that generates a Metropolis sequence of spin configurations\n",
      "\n",
      "* We start with a completely random spin configuration $\\{s \\}$ (each spin chosen to be up or down with\n",
      "50% probability).\n",
      "\n",
      "* To generate the next configuration in the sequence, we select a point $(i, j$) in the lattice, and we\n",
      "compute the energy change $\\Delta E$ by flipping the corresponding spin, s$_{ij}$ ,\n",
      "\n",
      "\\begin{equation}\n",
      "\\Delta E = 2s_{ij}  \\left[ s_{\\rm{up}(ij)} + s_{\\rm{down}(ij)} + s_{\\rm{left}(ij)} + s_{\\rm{right}(ij)} \\right]  \n",
      "\\end{equation}\n",
      "\n",
      "if E is zero or negative, we flip the spin (i.e., we always accept changes that lower the total energy).\n",
      "If E is positive, we flip the spin with probability $P (\\Delta E) = \\exp(-\\beta \\Delta E)$.\n",
      "We repeat this Metropolis step L $\\times$ L times, once for each location $(i, j)$ in the lattice. (The repeated Metropolis steps constitute a Metropolis sweep.) The updated points are used as they become available.\n",
      "\n",
      "* The Metropolis  estimate of the ensemble average of a macroscopic quantity (say M) is then\n",
      "\n",
      "\\begin{equation} \n",
      "\\left< M \\right>  = \\frac{1}{N}  \\sum_{k=1}^{N} \\sum_{ij} s_{ij}^{(k)} \n",
      "\\end{equation}\n",
      "\n",
      "where $ s_{ij}^{(k)}$  is the sequence of the spin configurations obtained by repeatedly applying Metropolis\n",
      "iteration to the initial configuration. It is usually best to omit a certain number of configurations at\n",
      "the beginning of the sequence; these are not distributed with the proper probabilities, because the\n",
      "system has not yet reached thermal equilibrium. (How long the system takes to thermalize depends\n",
      "on several factors, and the length of the required warm-up period is usually estimated empirically.\n",
      "\n",
      "\n",
      "(a) **Make the standard assumption of periodic boundary conditions, write a program that generates a Metropolis sequence of spin configurations for the  $100 \\times 100$ L 2-D Ising model (note, this will take a long time, perhaps an hour or more), and that visualizes them on the screen as they are computed.**\n",
      "\n",
      "(b) **Use Metropolis-Monte Carlo estimates to evaluate the energy per particle and magnetization perparticle of the 2-D Ising model as a function of temperature, then plot them.**\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Problem 2. Gibbs was a Gaussian\n",
      "\n",
      "The multivariate normal distribution can be drawn using a Gibbs sampler. Consider a 3-dimensional multivariate normal with mean vector (1, 2, 3), variances each equal to 4, and an intraclass correlation\n",
      "coefficient of 1.2; that is,\n",
      "\n",
      "$$ \n",
      "Y \\sim  N \\left(  \n",
      "\\left[\n",
      "\\begin{array}{c}\n",
      "    1  \\\\ \n",
      "    2  \\\\ \n",
      "    3  \\\\ \n",
      "  \\end{array} \n",
      "\\right],\n",
      "\\left[\n",
      "  \\begin{array}{ccc}\n",
      "    4 & 1.2 & 1.2 \\\\ \n",
      "    1.2 & 4 & 1.2 \\\\ \n",
      "    1.2 & 1.2 & 4 \\\\  \n",
      "  \\end{array} \n",
      "  \\right] \\right)\n",
      "$$ \n",
      "\n",
      "\n",
      "Solve for the conditional distributions of each component of Y with respect to the other two.\n",
      "The  covariance matrices can be expressed as \n",
      "\n",
      "$$ \\Sigma = \\left[ \\begin{array}{cc} \\Sigma_{11} & \\Sigma_{12} \\\\ \\Sigma_{21} & \\Sigma_{22}  \\end{array} \\right] $$\n",
      "\n",
      "where for the above example \n",
      "\n",
      "$$ \\Sigma_{12} = \\left[ 1.2 \\quad 1.2 \\right] $$\n",
      "\n",
      "and \n",
      "\n",
      "$$ \\Sigma_{22}= \\left[ \\begin{array}{cc} 4 & 1.2 \\\\ 1.2 & 4 \\end{array} \\right] $$\n",
      "\n",
      "You may use the standard decomposition of the multivariate normal, given below. The conditional distributions are univariate Normals with mean and standard deviations given\n",
      "\n",
      "$$ \\hat{x} = \\mu_x + \\Sigma_{12} \\Sigma^{-1}_{22}  \\left[   \\begin{array}{cc}  y-\\mu_y \\\\ z-\\mu_z \\end{array} \\right] $$\n",
      "$$ \\hat{y} = \\mu_y + \\Sigma_{12} \\Sigma^{-1}_{22}  \\left[   \\begin{array}{cc}  x-\\mu_x \\\\ z-\\mu_z \\end{array} \\right] $$\n",
      "$$ \\hat{z} = \\mu_z + \\Sigma_{12} \\Sigma^{-1}_{22}  \\left[   \\begin{array}{cc}  x-\\mu_x \\\\ y-\\mu_y \\end{array} \\right] $$\n",
      "\n",
      "where $\\mu_x, \\mu_y, \\mu_z$ are $1,2,3$ and $\\Sigma_{ij}$ represents the elements of the covariance matrix above.\n",
      "The similarity of the expressions comes from the fact that the matrix is very symmetric. The variance\n",
      "of the conditional distributions is given by\n",
      "\n",
      "$$ \\hat{\\Sigma} = \\Sigma_{11}-\\Sigma_{12} \\Sigma_{22}^{-1} \\Sigma_{21} $$\n",
      "\n",
      "(a) **Write the Gibbs sampler for this distribution by sampling sequentially from each of these conditional distributions. Choose a thinning parameter, burnin factor and total number of iterations that allow you to take 1000 non-autocorrelated draws. **\n",
      "\n",
      "(b) **Run this sampler 10 times from various starting positions and save the outcome draws. Check the mixing potential of this\n",
      "distribution using the Gelman-Rubin diagnostic described in the lecture.** \n",
      "\n",
      "(c) **Verify one of the ten chains against a true draw from the multivariate normal pdf by making a Q-Q plot (see `scipy.stats.probplot`), which displays a quantile-quantile plot of two samples. If the samples do come from the same distribution, the plot will be linear.**\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}