badgradientMCMC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook follows \"MinimizingChiSq\" and sets up some of the basic ingredients of an MCMC algorithm.  (However, this is *not* the implementation you should use in real life, for reasons that will become clear.)\n",
    "\n",
    "We start by loading the same dataset as before, and reminding ourselves what the data and $\\chi^2$ parameter space look like.  Note that the $\\chi^2$ function has been modified to account for 2d errors!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data = np.loadtxt('Lopez-Morales07_table1.txt', comments=';', delimiter=',')\n",
    "Mstar = data[:,0]\n",
    "Merr = data[:,1]\n",
    "Rstar = data[:,2]\n",
    "Rerr = data[:,3]\n",
    "\n",
    "plt.xlabel('Mass (M$_\\odot$)')\n",
    "plt.ylabel('Radius (R$_\\odot$)')\n",
    "plt.errorbar(Mstar, Rstar, xerr=Merr, yerr=Rerr, fmt='o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def chi2(A0, A1, X, D, SIGX, SIGD): \n",
    "    ''' Calculates chi^2 between the data (X, D) and the line described by y = A0 + A1*x\n",
    "    \n",
    "        Input\n",
    "        =====\n",
    "        A0 : y-intercept of the model\n",
    "        A1 : slope of the model\n",
    "        X : array of x-values from the data\n",
    "        D : array of y-values from the data\n",
    "        SIGX : the expected variance of the x-values\n",
    "            (could be a single value, or an array of the same length as X and D)\n",
    "        SIGD : the expected variance of the y-values\n",
    "            (could be a single value, or an array of the same length as X and D)\n",
    "    '''\n",
    "\n",
    "    chi2 = np.sum((A0 + A1*X - D)**2/(SIGD**2+(SIGX*A1)**2))\n",
    "    \n",
    "    return chi2 \n",
    "\n",
    "a0grid = np.linspace(-0.5, 0.5, 100)\n",
    "a1grid = np.linspace(0.5, 1.5, 100)\n",
    "chi2grid = np.zeros((100,100))\n",
    "\n",
    "for i in range(0,100):\n",
    "    for j in range(0,100):\n",
    "        chi2grid[i,j] = chi2(a0grid[i],a1grid[j],Mstar,Rstar,Merr,Rerr)\n",
    "        \n",
    "plt.imshow(np.flipud(np.log10(chi2grid).T), extent=[-0.5, 0.5, 0.5, 1.5])\n",
    "plt.xlabel('a0')\n",
    "plt.ylabel('a1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1**\n",
    "\n",
    "We could use our time much more efficiently if we didn't have to sample that entire grid of $a_0$ and $a_1$.  This is especially true in more complicated cases where the model has more than two parameters, since the volume of \"parameter space\" scales exponentially!  \n",
    "\n",
    "Let's implement a Markov Chain Monte Carlo algorithm that will always move in the direction of increasing likelihood.  Starting at an arbitrary point ($a_0$, $a_1$), the basic steps in the algorithm are as follows:\n",
    "\n",
    "1) Choose small steps $\\Delta a_0$ and $\\Delta a_1$ by which to update your position in parameter space.\n",
    "\n",
    "2) Evaluate the likelihoods at the current position $(a_0, a_1)$ and at the new position $(a_0+\\Delta a_0, a_1 + \\Delta a_1)$.\n",
    "\n",
    "3) Compare the two likelihoods.  If the likelihood is higher at the new position, move there!  If not, stay put.\n",
    "\n",
    "4) Return to step 1.\n",
    "\n",
    "Let's start with step 1.  A good way to choose $\\Delta a_0$ and $\\Delta a_1$ is by drawing from a normal distribution centered around 0 (so that you have an equal chance of drawing a positive or negative change).  It can be tricky to select the best value of $\\sigma_{\\Delta a}$, but for now let's just use 0.1 for both parameters.  You can draw a value from this distribution randomly using np.random.normal() like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(np.random.normal(loc=0.0, scale=0.1))\n",
    "print(np.random.normal(loc=0.0, scale=0.1))\n",
    "print(np.random.normal(loc=0.0, scale=0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**To do:** Try writing some code that random walks in $(a_0, a_1)$ space, saving the locations as it goes along, and plots the path on top of the $\\chi^2$ contours.  A skeleton is below.\n",
    "\n",
    "Try out different initial locations, step sizes, and total step numbers!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Initial location\n",
    "a0 = 0.\n",
    "a1 = 1.\n",
    "# Step sizes (sigmas of normal distributions)\n",
    "da0 = 0.1\n",
    "da1 = 0.1\n",
    "# Total number of steps to take\n",
    "nsteps = 25\n",
    "# Lists for keeping track of location\n",
    "a0list = [a0]\n",
    "a1list = [a1]\n",
    "\n",
    "for i in range(0,nsteps):\n",
    "    # YOUR CODE HERE\n",
    "    # Draw values for delta a0 and delta a1 from a normal distribution\n",
    "    # Update the values of a0 and a1\n",
    "    \n",
    "    # Add the new values to the lists\n",
    "    a0list.append(a0)\n",
    "    a1list.append(a1)\n",
    "\n",
    "# Plotting\n",
    "plt.imshow(np.flipud(np.log10(chi2grid).T), extent=[-0.5, 0.5, 0.5, 1.5])\n",
    "plt.xlabel('a0')\n",
    "plt.ylabel('a1')\n",
    "plt.plot(a0list, a1list, color='black')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's add in the condition that we will only take a step if it increases the likelihood value.  In this case, it is actually easier to stick with $\\chi^2$ rather than evaluating the likelihood itself, because computers don't really like comparing exponentials of large negative numbers.\n",
    "\n",
    "You can keep the same code from above, but now include an \"if\" statement before accepting the new step.\n",
    "\n",
    "As before, mess around with different starting locations, step sizes, and total number of steps.  How do these affect the behavior?  How many steps do you typically need to converge to within 0.01 of the minimum values (-0.0132, 1.1029)?  Within 0.001?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Initial location\n",
    "a0 = -0.4\n",
    "a1 = 0.6\n",
    "# Step sizes (sigmas of normal distributions)\n",
    "da0 = 0.1\n",
    "da1 = 0.1\n",
    "# Total number of steps to take\n",
    "nsteps = 25\n",
    "# Lists for keeping track of location\n",
    "a0list = [a0]\n",
    "a1list = [a1]\n",
    "\n",
    "for i in range(0,nsteps):\n",
    "    # YOUR CODE HERE\n",
    "    # Draw values for delta a0 and delta a1 from a normal distribution\n",
    "    # IF the chi2 value is lower (likelihood is higher) at the new location,\n",
    "    # Update the values of a0 and a1\n",
    "    \n",
    "    # Add the new values to the lists\n",
    "    a0list.append(a0)\n",
    "    a1list.append(a1)\n",
    "\n",
    "plt.imshow(np.flipud(np.log10(chi2grid).T), extent=[-0.5, 0.5, 0.5, 1.5])\n",
    "plt.xlabel('a0')\n",
    "plt.ylabel('a1')\n",
    "plt.plot(a0list, a1list, color='black')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "As you may have deduced, there is a major problem with using this method - what if you have local extrema?  If you wander in, you will get stuck there and will be unable to find the global minimum/maximum.\n",
    "\n",
    "Tomorrow, we will discuss how to get around this issue and introduce ourselves to the Metropolis-Hastings algorithm!"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}