linearRegression.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:b92a4438c3cfdabf7bba9979e7066763a40dbfd82e6c5132da39def7231a77ea"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "<center> Linear Regression Lab </center>\n",
      "\n",
      "<center> February 17th 2015 </center>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Regression refers to the technique employed when learning a mapping from vector of input data to a quantitative output. For example, coordinates in a room to temperature, someone's Body Mass Index to their life expectancy or a stock's performance over the last 5 days to it's value tomorrow.\n",
      "\n",
      "There are numerous methods for addressing this problem, each with their own set of assumptions and behaviour. Some are ideal for tackling large volumes of data while others provide more informative probabilistic outputs.\n",
      "\n",
      "In this lab, we'll focus on Linear Regression; a good point of reference for many other regression techniques and still used widely around the world today due to their simplicity and favourable scaling characteristics."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Part I: Maximum Likelihood Linear Regression"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this section, we're going to explore what linear regression models are and investigate how they can be trained using a \"Maximum Likelihood\" approach."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Import related packages\n",
      "import numpy as np\n",
      "import scipy.linalg as linalg\n",
      "import scipy.stats as stats\n",
      "import matplotlib.pyplot as pl\n",
      "\n",
      "#Set generated images to appear underneath the related cell\n",
      "%matplotlib inline "
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First, lets generate some synthetic data. Here's we create a random 1-D vector $\\mathbf{x}$ of observation locations and evaluate an \"unknown\" ground truth function at these locations before adding some noise to create our observations, $\\mathbf{y}$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nObs = 20          # number of observations\n",
      "betaInv = 0.2      # variance of Gaussian noise applied to samples\n",
      "x = np.sort(np.random.random(nObs)*10)[:,np.newaxis]  # observation locations\n",
      "y = 0.5*x + 1 + 1*np.sin(x)+ (np.random.randn(nObs)*betaInv)[:,np.newaxis] # ground truth (we don't know this!) plus some noise"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Lets plot the data first to see what it looks like:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pl.figure(figsize=(15,5))\n",
      "pl.plot(x,y,'b.')\n",
      "pl.axis([0,10,0,7])\n",
      "pl.xlabel('x')\n",
      "pl.ylabel('y')\n",
      "pl.title('Raw Data')\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Linear Model with Polynomial Features"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Okay, now lets try to fit a model to it.\n",
      "\n",
      "One of the most straightforward types of regression is maximum likelihood regression using a linear model, (the output is a weighted sum of the features).\n",
      "For now, lets use a very simple model with just two parameters (or weights), $\\mathbf{w}$.\n",
      "\\begin{equation}\n",
      "f(x) = w_0 +  w_1 x\n",
      "\\end{equation}\n",
      "In matrix form this is $\\mathbf{f}= \\mathbf{\\theta w}$ where $\\mathbf{\\theta}$ is $[x^0,x^1]$ (the feature matrix) and $\\mathbf{w}$ is $[w_0, w_1]^T$. \n",
      "\n",
      "The feature matrix has a dimensionality of $N\\times p$ where $N$ is the number of inputs and $p$ is the number of features. The complexity of the model can be increased by increasing $p$. Polynomial features (increasing powers of $x$) are commonly used in practice. Evaluating a polynomial model at locations $[x_1,...,x_N]$ and weights $[w_0,...,w_p]$ is given as:\n",
      "\n",
      "\\begin{equation}\n",
      "\\begin{bmatrix}\n",
      "f_1 \\\\\n",
      "f_2 \\\\\n",
      "\\vdots \\\\\n",
      "f_N\n",
      "\\end{bmatrix} =\n",
      " \\begin{bmatrix}\n",
      "  x_1^0 & x_1^1 & \\cdots & x_1^{p-1} \\\\\n",
      "  x_2^0 & x_2^1 & \\cdots & x_2^{p-1} \\\\\n",
      "  \\vdots  & \\vdots  & \\ddots & \\vdots  \\\\\n",
      "  x_N^0 & x_N^1 & \\cdots & x_N^{p-1}\n",
      " \\end{bmatrix} \n",
      " \\begin{bmatrix}\n",
      "  w_0 \\\\\n",
      "  w_1 \\\\\n",
      "  \\vdots \\\\ \n",
      "  w_{p-1} \n",
      " \\end{bmatrix}\n",
      "\\end{equation}\n",
      "\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Lets query our model and compare it against the data. To generate a polynomial feature matrix $\\theta$, we'll use a\n",
      "custom function called polyFeatureGen() that takes the observation or query locations ($x$ or $x_{\\text{query}}$) and the number of features ($p$) as inputs. It returns an $n \\times p$ matrix such that:\n",
      "\\begin{equation}\n",
      "\\theta_{i,j} = x_i^{j-1}\n",
      "\\end{equation}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "__Complete the function in the cell below:__"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Define Polynomial basis functions -------------------\n",
      "# x: input vector (n x d) where n is the number of inputs and d is their dimensionality (you can assume d=1 for this lab)\n",
      "# p: number of features (scalar)\n",
      "# theta: the feature matrix (n x p)\n",
      "\n",
      "def polyFeatureGen(x,p):\n",
      "    # Code goes here \n",
      "    #\n",
      "    #\n",
      "    return theta    # n by p sized array"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Test polyFeatureGen\n",
      "# Run this cell. The output should be: \n",
      "# [[1  1  1]\n",
      "#  [1  4 16]\n",
      "#  [1  3  9]\n",
      "#  [1  6 36]]\n",
      "\n",
      "xTest = np.array([1,4,3,6])[:,np.newaxis] # 4 1-D inputs\n",
      "pTest = 3 #This test model will have 3 features\n",
      "thetaTest = polyFeatureGen(xTest,pTest) #the test feature matrix (should be 4 x 3)\n",
      "print(thetaTest)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For now, lets stick with two features. Lets pick two arbitriary values for the weights. Say,"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p = 2 #number of features\n",
      "w_guess = np.array([-0.2, 5])[:,np.newaxis]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Generate the query points and the associated feature matrix"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nQuery = 1000 # number of query inputs\n",
      "xQuery = np.linspace(0,10,nQuery)[:,np.newaxis] #linearly spaced from 0 to 10"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "__Write a line of code to generate the associated feature matrix for the query points:__"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Uncomment and complete the line of code below to generate thetaQuery, (the feature matrix for the query points):\n",
      "#\n",
      "# thetaQuery =                  \n",
      "#\n",
      "\n",
      "\n",
      "print(thetaQuery) #Should be a nQuery x p matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now evaluate the model at the query locations by multiplying the query feature matrix by our guessed weights:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fQuery = thetaQuery.dot(w_guess)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And plot the results:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = pl.figure(figsize=(15,5))\n",
      "pl.plot(x,y,'b.')\n",
      "pl.plot(xQuery,fQuery,'r')\n",
      "pl.xlabel('x')\n",
      "pl.ylabel('f(x)')\n",
      "pl.title('Model prediction using guessed weights')\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It's not a great fit! The Mean Squared Error between the model predictions and the observations is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fModel = polyFeatureGen(x,p).dot(w_guess)\n",
      "MSE = np.mean((fModel - y)**2)\n",
      "print('Mean Squared Error: ' + str(MSE))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Maximum Likeliood Estimation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now lets learn the parameters of our model to see if we can improve the fit.\n",
      "The goal here is to estimate the parameter values (i.e. $\\mathbf{w}$) from the data. A popular way of doing this is by maximising the conditional likelihood of the data given a model, $p(y|x,w,\\beta^{-1})$, with respect to the parameters.\n",
      "\n",
      "A common assumption to make is that the measurements we gathered are subject to Gaussian noise\n",
      "i.e. we assume that our observations, y are drawn from a Gaussian distribution with mean $\\mathbf{w\\theta^T}$ and\n",
      "variance $\\beta^{-1}$ (indicating the level of noise). \n",
      "\n",
      "The probability of acquiring a particular set of observation, given a specific model is\n",
      "\\begin{equation}\n",
      "p(y|x,w,\\beta^{-1}) = \\prod_{n=1}^{N} \\mathcal{N}(\\mathbf{w\\theta^T}, \\beta^{-1})\n",
      "\\end{equation}\n",
      "\n",
      "It is often referred to as the likelihood of the data. In log form:\n",
      "\\begin{equation}\n",
      "\\text{ln}p(\\mathbf{y}|\\mathbf{x},\\mathbf{w},\\beta^{-1}) = \\sum_{n=1}^{N} \\text{ln}\\big( \\mathcal{N}(\\mathbf{w\\theta^T}, \\beta^{-1})\\big) \\\\\n",
      "                 = \\frac{N}{2} \\text{ln}(\\beta) - \\frac{N}{2}\\text{ln}(2\\pi) - \\frac{\\beta}{2} \\sum_{n=1}^{N} (\\mathbf{y}-\\mathbf{w\\theta^T})^2   \n",
      "\\end{equation}\n",
      "\n",
      "\n",
      "The final term is the sum of squares error between the observations and the model.\n",
      "If we vary the parameters of the model to maximise the likelihood of the data, we train a model\n",
      "that best fits the data\n",
      "The gradient of the log likelihood with respect to the weights is:\n",
      "\\begin{equation}\n",
      "\\triangledown \\text{ln}p(\\mathbf{y}|\\mathbf{x},\\mathbf{w},\\beta^{-1})= \\sum_{n=1}^{N} (\\mathbf{y}_n - \\mathbf{w\\theta^T})\\theta\n",
      "\\end{equation}\n",
      "Find the max by setting the gradient to zero and then solve for $\\mathbf{w}$ to find the \"maximum-likelihood parameters\" (w_ml):\n",
      "\\begin{equation}\n",
      "\\mathbf{w}_{ml}^T = (\\mathbf{\\theta}^T\\mathbf{\\theta})^{-1}\\mathbf{\\theta^T y}\n",
      "\\end{equation}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "__Generate the feature matrix using the observations locations.__"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#theta =                   # the shape of this matrix will be N x p"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "__Then find the parameters that maximise the likelihood as described above.__"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Implement a line of code to determine w_ml. Hint: use np.linalg.solve(A,B) \n",
      "#rather than np.linalg.inv(A)B for better numerical stability\n",
      "\n",
      "#w_ml = \n",
      "\n",
      "\n",
      "print(w_ml)  # the shape of w_ml should be p x 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "__Now evaluate the model at the query locations using the new parameters and plot the raw data and the new model's prediction as before:__"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# fQuery =                     # the model's predictions at the query locations. This will be an nQuery x 1 matrix\n",
      "\n",
      "fig = pl.figure(figsize=(15,5))\n",
      "pl.plot(x,y,'b.')\n",
      "pl.plot(xQuery,fQuery,'r')\n",
      "pl.axis([0, 10, 0, 8])\n",
      "pl.title('Max Likelihood weights with ' + str(p) + ' features')\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Much better. __Calculate and print the mean squared error between the model predictions and the observations using these new parameters to quantively compare with our initial set of parameters:__"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# fModel =            # the model's predictions at the observation locations. This will be an N x 1 matrix\n",
      "# MSE =               # scalar\n",
      "\n",
      "print('Mean Squared Error: ' + str(MSE))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This MSE is the best we can do with this model."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Increasing the Model Complexity"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To reduce the MSE further, we will need to increase the model's flexibility by adding higher order polynomial features. \n",
      "__In the cell below, change the number of features in the model. Run the cell to learn the max likelihood weights for that set of features, plot the resulting query outputs as before and evaluate the MSE. \n",
      "What is the number of features that gives the lowest MSE?__"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Set the number of features\n",
      "p = 3 #number of features\n",
      "\n",
      "# Build the design feature matrix\n",
      "theta = polyFeatureGen(x,p)\n",
      "\n",
      "# Determine the maximum likelihood weights\n",
      "w_ml = np.linalg.solve((np.transpose(theta).dot(theta)),(np.transpose(theta))).dot(y)\n",
      "\n",
      "# Build the feature matrix for the query points\n",
      "thetaQuery = polyFeatureGen(xQuery,p)\n",
      "\n",
      "# Evaluate the model at the query points\n",
      "fQuery = thetaQuery.dot(w_ml)\n",
      "\n",
      "fig = pl.figure(figsize=(15,5))\n",
      "pl.plot(x,y,'b.')\n",
      "pl.plot(xQuery,fQuery,'r')\n",
      "pl.axis([0, 10, 0, 8])\n",
      "pl.xlabel('x')\n",
      "pl.ylabel('f(x)')\n",
      "pl.title('Number of Features: '+str(p))\n",
      "pl.show()\n",
      "\n",
      "fModel = theta.dot(w_ml)\n",
      "MSE = np.mean((fModel - y)**2)\n",
      "print('Mean Squared Error: ' + str(MSE))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As you may have noticed, increasing the flexibility and model complexity decreases the mean squared error because the model's ability to fit all of the observation exactly increases as well. However, it also results in \"overfitting\". This is when the model fits the observations very well but generalises poorly to out locations. If a number of observations are withheld from training, the model would not do well at predicting their value as it is more concerned with overfitting the observations it can see.\n",
      "\n",
      "\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Basis Functions"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Up until now, the feature matrix has been generated using polynomial equation:\n",
      "\n",
      "\\begin{equation}\n",
      "\\theta_{i,j} = x_i^{j-1}\n",
      "\\end{equation}\n",
      "This is sometimes referred to as the  polynomial basis function. \n",
      "\n",
      "Run the cell below to visualise the basis functions that were used in the last model you trained:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pl.figure(figsize=(15,5))\n",
      "for i in range(p):\n",
      "    pl.plot(xQuery,thetaQuery[:,(i-1)])\n",
      "    pl.axis([0, 10, 0, 80])\n",
      "    pl.title('Individual basis functions')\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The above graph is simply plotting each column or feature of the feature matrix. **Uncomment and complete the line of code in the cell below to plot the basis functions** _**multipled by their associated learnt weight**_."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pl.figure(figsize=(15,5))\n",
      "for i in range(p):\n",
      "    # pl.plot(Code goes here)\n",
      "    pl.axis([0, 10, -80, 80])\n",
      "    pl.title('Individual weighted basis functions')\n",
      "pl.show()\n",
      "print('Learnt weights: ' + str(w_ml.transpose()))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Interestingly, the model's complicated looking prediction is just a linear combination of these basis functions. Hence the name \"linear\" regression."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There are many possible variants of basis functions each giving the resulting model different characteristics (e.g. smooth, non-stationary, discontinuous etc.). \n",
      "\n",
      "Radial basis functions are a commonly used alternative to the polynomial function:\n",
      "\\begin{equation}\n",
      "\\theta_{i,j} = \\text{exp}\\Bigg(-\\frac{(x_i-\\mu_j)^2}{2s^2}\\Bigg)\n",
      "\\end{equation}\n",
      "Here, $\\mu_j$, is the mean of the $j^{\\text{th}}$ radial basis function and $s$ corresponds to the basis functions width or its region of influence. \n",
      "\n",
      "\n",
      "**Implement the basis function radialFeatureGen()** that takes a vector of inputs ($\\textbf{x}$), a vector of $p$ mean locations ($\\mathbf{\\mu}$) with shape $p \\times 1$ and a basis function spatial width ($s$). It then calculates the feature matrix, $\\theta$, using the radial basis function given above. \n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Define Radial basis functions -------------------\n",
      "# Dimensions:\n",
      "# x: vector of inputs (shape: n x 1)\n",
      "# mu: vector of feature means (shape: p x 1)\n",
      "# s: the spatial width of the radial basis functions (scalar)\n",
      "# theta: n x p\n",
      "\n",
      "def radialFeatureGen(x,mu, s):\n",
      "    # Code goes here\n",
      "    #\n",
      "    #\n",
      "    return theta    #n x p sized array"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Using the values for $\\mu$ and $s$ provided below, **write some code that finds the weights with the maximum likelihood** ($w_{ml}$) for a design matrix generated using the radial basis function. Then **evaluate the model at the query locations and plot the results as before.**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mu = np.array([1,2,3,4,5])[:,np.newaxis] #The mean of the radial basis functions used in this model\n",
      "s = 1 # the spatial width of the basis function\n",
      "\n",
      "# Uncomment and complete the following lines\n",
      "# theta =                 # the training feature matrix (N x p) using the radial basis functions\n",
      "# w_ml =                  # the max likelihood weights \n",
      "# thetaQuery =            # the query feature matrix (nQuery x p)\n",
      "# \n",
      "# Evaluate the model at the query points\n",
      "# fQuery =                # the predicted values at xQuery (nQuery x 1) \n",
      "\n",
      "# Plot the results:\n",
      "fig = pl.figure(figsize=(15,5))\n",
      "pl.plot(x,y,'b.')\n",
      "pl.plot(xQuery,fQuery,'r')\n",
      "pl.axis([0, 10, 0, 8])\n",
      "pl.figure(figsize=(15,5))\n",
      "pl.title('Radial Basis Function')\n",
      "for i in range(mu.shape[0]):\n",
      "    pl.plot(xQuery,thetaQuery[:,(i-1)][:,np.newaxis]*w_ml[i])\n",
      "    pl.title('Individual weighted basis functions')\n",
      "pl.show()\n",
      "print('Learnt weights: ' + str(w_ml.transpose()))\n",
      "\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Rerun the above cell but change the basis functions by altering $\\mu$. e.g. **mu = np.array([1,2,3,4,5,6,7,8,9,10])[:,np.newaxis]**\n",
      "\n",
      "* What effect does changing $s$ have on the model? \n",
      "* How might we find the optimal values for $\\mu$ and $s$?\n",
      "* How does the behaviour of the radial basis function model differ from the polynomial case in regions far from data?\n",
      " "
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Downsides to Maximum Likelihood Regression"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "While this approach to regressing over data is appealling, it does have a few shortcomings.\n",
      "For example, we have no estimate of how confident the model is in its predictions (which can be crucial in decision-making applications). Also the model learns a single hypothesis to explain the data. In reality there may be numerous explanations for the same set of observations. In Part II of today's lab, we will investigate a probabilitic approach to regression."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "% Bayesian Linear Regression Tutorial, Alistair Reid, NICTA 2015\n",
      "\\newcommand{\\W} [0] {\\mathbf{w}}\n",
      "\\newcommand{\\X} [0] {\\mathbf{x}}\n",
      "\\newcommand{\\y} [0] {\\mathbf{y}}\n",
      "\\newcommand{\\Xf} [0] {\\mathbf{\\theta}}\n",
      "\\newcommand{\\norm} [0] {\\mathcal{N}}\n",
      "\\newcommand{\\weightprec} [0] {Q}\n",
      "\\newcommand{\\weightprecsimp} [0] {\\alpha I}\n",
      "\\newcommand{\\alphs} [0] {\\mathbf{\\alpha}}\n",
      "\\newcommand{\\invnoise} [0] {\\beta}\n",
      "\\newcommand{\\noise} [0] {\\beta^{-1}}\n",
      "\\newcommand{\\Ainv} [0] {A^{-1}}\n",
      "\\newcommand{\\fq} [0] {\\mathbf{f_*}}\n",
      "\\newcommand{\\xq} [0] {\\mathbf{x_*}}\n",
      "$$ "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Part 2: Bayesian Linear Regression\n",
      "\n",
      "In the previous sections, we have selected the linear regression weights by maximising the likelihood of the data. The learnt model certainly makes useful predictions, but it provides no notion of uncertainty: we cannot distinguish between a prediction that is tightly constrained by the data, and a prediction that is essentially an unconstrained \\`guess'. Also, if there are multiple likely explanations of the data, we have neglected all but one. Recall that for a given feature mapping $\\X \\to \\Xf$, the predictive behaviour of a linear regression is encoded in the weights:\n",
      "\\begin{equation}\n",
      "f(x) = \\Xf \\W \n",
      "\\end{equation}\n",
      "\n",
      "Where $\\W$ is $p \\times 1$ and $\\Xf$ is $n \\times p$. Instead of using a single *maximum likelihood* weight vector, we can take a Bayesian approach and consider a probability distribution over all possible weight vectors.  The distribution can be inferred from the observations using Bayes rule:\n",
      "\n",
      "\\begin{equation}\n",
      "p(\\W|\\y) \\propto p(\\y|\\W) p(\\W)\n",
      "\\end{equation}\n",
      "\n",
      "The likelihood $p(\\y|\\W)$ is exactly the criterion we were using to select the maximum likelihood weights:\n",
      "\n",
      "\\begin{equation}\n",
      "p(\\y | \\W) = \\norm (y |  \\Xf \\W,  \\noise I )\n",
      "\\end{equation}\n",
      "\n",
      "In addition, we must now specify a prior distribution over $\\W$ encoding our belief about reasonable weight values. Any probability distribution we consider to be sensible for the problem could be used here, but some are easier to work with than others. In this tutorial, we will use a simple Gaussian prior with precision (inverse covariance) $\\weightprecsimp$ that is $p \\times p$ with a single free parameter:\n",
      "\n",
      "\\begin{equation}\n",
      "\\W \\sim \\norm( 0,  \\alphs^{-1}I )\n",
      "\\end{equation}\n",
      "\n",
      "Hence we are assuming the weights are independently distributed & penalised for large magnitudes. While the $\\W$'s are model parameters that directly affect the output, we consider the $\\alphs$'s to be *hyperparameters* because they tune the prior distribution over the parameters $\\W$. \n",
      "\n",
      "Because we have assumed a Gaussian form for both the likelihood and the prior, Bayes rule can be computed analytically in this case to obtain a posterior distribution over the weights. This involves expanding the exponents of the two Gaussians and completing the square (an exercise left to the reader).\n",
      "\n",
      "The posterior over the weights is given by: \n",
      "\\begin{equation}\n",
      "p(\\W|\\Xf,y) \\sim \\norm \\left( \\mu_w, \\Ainv \\right)\n",
      "\\end{equation}\n",
      "Where $A$ is the $p \\times p$ posterior precision given by: \n",
      "\\begin{equation}\n",
      "A = \\invnoise \\Xf^\\top \\Xf + \\weightprecsimp\n",
      "\\end{equation}\n",
      "And the weight means are given by:\n",
      "\\begin{equation}\n",
      "\\mu_w = \\invnoise \\Ainv \\Xf^\\top \\y\n",
      "\\end{equation}\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Exercise: Compute Posterior Weight Distribution\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy.linalg as linalg\n",
      "import matplotlib.pyplot as pl\n",
      "import scipy.optimize # In addition, nlopt would be a great optimisation library to consider\n",
      "from math import log, pi\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      " # Generate Noisy Data ---------------------------------\n",
      "nSamp = 5\n",
      "noiseSTD = 0.5\n",
      "xdomain = 10\n",
      "x = np.sort(np.random.random(nSamp)*xdomain)[:,np.newaxis] # x is n by 1\n",
      "y = x + np.random.randn(nSamp)[:,np.newaxis]*noiseSTD - 5  # y is also n by 1\n",
      "\n",
      "pl.figure(figsize=(5,5) )\n",
      "pl.plot(x,y,'k.')\n",
      "pl.grid()\n",
      "pl.title('Noisy observations of the linear target function')\n",
      "pl.xlabel('x')\n",
      "pl.ylabel('y')\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Add some code to compute $A$ and $\\mu_w$ in the cell below:**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use 2 features - (1, x) to learn intercept and slope -------------------\n",
      "p = 2    # n features...\n",
      "n = len(x)\n",
      "\n",
      "# build a feature vector n*[1, x]\n",
      "theta = polyFeatureGen(x,p)   \n",
      "\n",
      "# Define a prior\n",
      "alpha = 2**-2           # Weight precision -> Larger alpha 'relaxes' the complexity penalty - try it\n",
      "beta  = noiseSTD**-2     # noise inverse standard deviation - cheat and use the real value\n",
      "\n",
      "# Compute the posterior\n",
      "\n",
      "A = ... \n",
      "# Hint -> when computing A^-1 b it is much more efficient to use linalg.solve(A,b) then to compute A.inverse() explicitly\n",
      "mu_w = ... \n",
      "\n",
      "\n",
      "# Plot a countour ellipse of the weight distribution\n",
      "th = np.linspace(0,2*pi,100)\n",
      "xy = mu_w - np.log(0.1/2)*linalg.cholesky(linalg.inv(A),lower=False).dot(np.vstack((np.sin(th), np.cos(th))))\n",
      "\n",
      "postpl = pl.plot(xy[0], xy[1],'b')\n",
      "xy = np.zeros(mu_w.shape) - np.log(0.1/2)*np.vstack((np.sin(th), np.cos(th)))/np.sqrt(alpha)\n",
      "pripl = pl.plot(xy[0], xy[1],'g')\n",
      "\n",
      "\n",
      "pl.xlabel('Intercept')\n",
      "pl.ylabel('Gradient')\n",
      "pl.title('Weights (prior=green, posterior=blue) - contour at 0.1 max probability')\n",
      "pl.axis('equal')\n",
      "pl.show()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Marginalising over weights\n",
      "\n",
      "While it is interesting to see the distribution of weights used by the Bayesian linear regression, the ultimate goal is to predict the outputs $\\fq$ corresponding to inputs $\\xq$. To do this, we *marginalise* out the weights $\\W$:\n",
      "\n",
      "\\begin{equation}\n",
      "p(\\fq | \\xq, \\Xf, \\y) = \\int p(\\fq|\\xq,\\W) p(\\W|\\Xf,\\y) d\\W\n",
      "\\end{equation}\n",
      "\n",
      "We know from above that:\n",
      "\\begin{equation}\n",
      "p(\\W|\\Xf,y) \\sim \\norm ( \\mu_w, \\Ainv )\n",
      "\\end{equation}\n",
      "\n",
      "We have also defined a deterministic mapping to the predictive function:\n",
      "\\begin{equation}\n",
      "\\fq = \\xq \\W \n",
      "\\end{equation}\n",
      "\n",
      "This marginalisation has a closed form solution - we can produce predictive distributions, considering all possible weights, without explicity computing them:\n",
      "\\begin{equation}\n",
      "\\fq \\sim \\norm \\left( \\invnoise \\xq \\Ainv \\Xf^\\top \\y, \\hspace{0.5cm} \\xq \\Ainv \\xq^\\top \\right)\n",
      "\\end{equation}\n",
      "\n",
      "Note that this equation gives a multivariate Gaussian prediction: for $m$ queries, the covariance will be $m \\times m$. This is useful for joint draws. However, if you are only interested in the envelope of uncertainty, then the diagonal can be extracted to give the marginal variance of each prediction. Don't forget to add the noise variance back on if you want to predict the range of possible *observations* rather than *latent functions*!\n",
      "\n",
      "### Exercise:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Complete the unfinished lines of code in the cell below:**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make predictions directly from the data-------------------\n",
      "nQuery = 1000\n",
      "xQuery = np.linspace(-1, 11, nQuery)[:,np.newaxis]\n",
      "thetaQuery = polyFeatureGen(xQuery,p)\n",
      "\n",
      "# lets make this a reusable function\n",
      "def linreg(theta, y, thetaQuery, alpha, beta):\n",
      "    I = np.eye(theta.shape[1])\n",
      "    \n",
      "    A = ...\n",
      "    mu_w = ... \n",
      "    f_mean = ... \n",
      "\n",
      "    # Simple Approach\n",
      "    f_covar = ....\n",
      "    f_std = np.sqrt(np.diag(f_covar))[:,np.newaxis]\n",
      "    \n",
      "    # Optional Advanced challenge: use cholesky factorisation to compute the diagonal of f_covar directly\n",
      "    #L = linalg.cholesky(A, lower=True)  # Factorise A into L L^T\n",
      "    #diagonal = ...\n",
      "    #f_std = np.sqrt( diagonal  + 1./beta)[:,np.newaxis] # try with and without noise variance\n",
      "\n",
      "    #f_std = np.sqrt(np.diag(f_covar))\n",
      "    return f_mean, f_std\n",
      "\n",
      "f_mean, f_std = linreg(theta, y, thetaQuery, alpha, beta)\n",
      "\n",
      "# plot the mean, and +- 2 standard deviations\n",
      "pl.figure(figsize=(5,5) )\n",
      "pl.plot(x,y,'k.')\n",
      "pl.plot(xQuery, f_mean, 'b')\n",
      "pl.plot(xQuery, f_mean+2.5*f_std, 'r')\n",
      "pl.plot(xQuery, f_mean-2.5*f_std, 'r')\n",
      "pl.xlabel('Input x')\n",
      "pl.ylabel('Output y')\n",
      "pl.title('Predictive Envelope with Data and Mean')\n",
      "pl.show()\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Exercise: Non-linear Basis Functions\n",
      "\n",
      "We have already been using simple linear basis functions, but lets try a demo with a non-linear problem. Don't worry about visualising the distribution over weights - with many polynomial features they will be too high dimensional to show, but we will get very interesting predictive marginals:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      " # Generate Noisy Sinusoid Data ---------------------------------\n",
      "nSamp = 30\n",
      "noiseSTD = 0.1\n",
      "xdomain = 10\n",
      "x = np.sort(np.random.random(nSamp)) \n",
      "x = x-np.min(x)\n",
      "x = x/np.max(x) * xdomain\n",
      "y = np.sin(x) + np.random.randn(nSamp)*noiseSTD\n",
      "x = x[:,np.newaxis]\n",
      "y = y[:,np.newaxis]\n",
      "\n",
      "\n",
      "pl.figure(figsize=(15,5) )\n",
      "pl.plot(x,y,'k.')\n",
      "pl.xlabel('x')\n",
      "pl.ylabel('y')\n",
      "pl.title('Noisy observations of the target function')\n",
      "pl.grid()\n",
      "pl.show()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This data clearly has a non-linear pattern. As we saw in part 1 of this tutorial, linear regression techniques can perform well on non-linear data by introducing basis functions. Lets apply the polynomial basis functions you implemented in polyFeatureGen to increase the dimensionality of the inputs to $p=7$ and run the regression again."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Complete the unfinished lines of code in the cell below:**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Bayesian linear regression - prediction with polynomial basis functions\n",
      "# Define prediction locations -------------------\n",
      "nQuery = 1000\n",
      "xQuery = np.linspace(0, 10, nQuery)[:,np.newaxis]\n",
      "\n",
      "# transform into feature space (now n by p)\n",
      "p = 7   # n features...\n",
      "...\n",
      "...\n",
      "\n",
      "# Define a prior\n",
      "alpha = 0.7**-2 # per-weight inverse standard deviation - chosen automatically using strategy in final section\n",
      "beta  = 0.3**-2 # noise inverse standard deviation - using the real value for now\n",
      "\n",
      "f_mean, f_std = linreg(theta, y, thetaQuery, alpha, beta)\n",
      "\n",
      "# plot the mean, and +- 2 standard deviations\n",
      "pl.figure(figsize=(15,5) )\n",
      "pl.plot(x,y,'k.')\n",
      "pl.plot(xQuery, f_mean, 'b')\n",
      "pl.plot(xQuery, f_mean+2.5*f_std, 'r')\n",
      "pl.plot(xQuery, f_mean-2.5*f_std, 'r')\n",
      "pl.xlabel('Input x')\n",
      "pl.ylabel('Output y')\n",
      "pl.title('Predictive Envelope with Data and Mean')\n",
      "pl.show()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now the predictions are able to capture the non-linear structure of this problem! There are two ways of interpreting how this is working.\n",
      "\n",
      "The first interpretation is that we are inferring a probability distribution over the weights, and each weight corresponds to a basis function. With the polynomial basis functions, we are therefore considering distributions over the polynomial coefficients in a polynomial equation. You can see in the predictive distribution:\n",
      "\\begin{equation}\n",
      "\\fq \\sim \\norm \\left( \\invnoise \\xq \\Ainv \\Xf^\\top \\y, \\hspace{0.5cm} \\xq \\Ainv \\xq^\\top \\right)\n",
      "\\end{equation}\n",
      "that our mean function is in fact a linear combination of the basis functions composing $\\Xf$, and the space of models is restricted to the span of the basis functions. \n",
      "\n",
      "A second interpretation is that we have used polyFeatureGen to project the inputs and the queries into 7 dimensions, effectively warping the space on the x-axis. The line in the high dimensional space can then project to a nonlinear function in the original space. To understand this principle, think about how a direct path an airline takes around a [3D sphere](https://www.metabunk.org/data/MetaMirrorCache/contrailscience.com_skitch_Google_Earth_20111211_141946.jpg) can appear as a curved arc when projected onto a [2d map](http://www.gcmap.com/map?P=syd-jnb&MS=wls&MR=900&MX=720x360&PM=*) ). You will see this principle applied in the support vector machine tutorial. \n",
      "\n",
      "You can investigate the choice of basis functions by modifying the above code:\n",
      "* What happens to the prediction as you move away from the input data? (try changing the domain above)\n",
      "* We have been using polynomial basis functions. Can you think of a different class of basis functions that might work for periodic data?\n",
      "\n",
      "## Plotting draws from the posterior\n",
      "\n",
      "We have been looking at the mean and envelope of the predictive distribution, but it is also useful to look at draws from the posterior distribution. To do this we must look at the *joint* multivariate distribution of the posterior. Note that given our particular choice of Gaussian prior and likelihood, the posterior predictions are also joint-Gaussian. For example, we can compute a posterior multivariate Gaussian over the weights and draw random weights out of this distribution. This will sample a set of deterministic models each drawn with probability consistent with our Bayesian analysis.\n",
      "\n",
      "How do we draw values out of a multivariate Gaussian distribution? We can do this by linear operations on a vector of independent univariate Gaussian draws. This can be achieved using the cholesky factorisation: if $L$ is the upper Cholesky factor of the  covariance matrix $A^{-1}$, and $\\mathbf{\\Sigma}$ is a vector of iid Gaussians, we may draw from the distribution according to:\n",
      "\\begin{equation}\n",
      "\\hat{\\W} = \\mu_w + L \\Sigma\n",
      "\\end{equation}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Complete the unfinished lines of code in the cell below:**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plotting draws from the function\n",
      "I = np.eye(theta.shape[1])\n",
      "A = ...\n",
      "mu_w = ...\n",
      "L = linalg.cholesky(linalg.inv(A), lower=True)\n",
      "\n",
      "nPlots = 10\n",
      "fig = pl.figure(figsize=(15,5))\n",
      "pl.plot(x,y,'k.')\n",
      "for i in range(nPlots):\n",
      "    Sigma = np.random.randn(p)[:,np.newaxis]\n",
      "    w_hat = ... # use mu_w,  L, and Sigma\n",
      "    f_hat =  ... # use thetaQuery and w_hat\n",
      "    pl.plot( ... \n",
      "        \n",
      "pl.xlabel('x')\n",
      "pl.ylabel('y')\n",
      "pl.title('Draws from the posterior')\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* If we were to average over a very large number of these draws, the density at each input value would be Gaussian distributed with the marginals from the previous section. \n",
      "\n",
      "\n",
      "### Advanced: Hyperparameter selection\n",
      "\n",
      "Until now we have been manually specifying values for $\\alphs^2$ and $\\invnoise$. Such manual selection is not ideal, as it requires a deep understanding of every dataset and large amounts of manual tuning to select appropriate noise and prior hyperparameters. Instead, there are two main approaches for hyperparameter selection based on the data. Both rely on the concept of *marginal likelihood* - the probability of drawing the observations out of the prior. We can obtain this criterion by marginalising over the actual weight vectors and the noise draws. The log marginal likelihood given all our Gausianity assumptions so far is given by:\n",
      "\n",
      "\\begin{equation}\n",
      "\\log p(y | \\Xf) = \\frac{1}{2} \\left( p \\log \\alphs^2 + n \\log \\invnoise - \\log |A| - n \\log 2\\pi - \\invnoise (y-y^*)^\\top (y-y^*) - \\alphs^2 \\mu^{*\\top} \\mu^* \\right)\n",
      "\\end{equation}\n",
      "\n",
      "\\begin{equation}\n",
      "\\mu^* = \\invnoise \\Ainv \\Xf^T y\n",
      "\\end{equation}\n",
      "\n",
      "\\begin{equation}\n",
      "\\y^* = \\Xf \\mu^*\n",
      "\\end{equation}\n",
      "\n",
      "See \"Evaluation of the evidence function\" in Bishop - Pattern recognition and machine learning for how this is derived. Here $|A|$ is the determinant of $A$, and $\\mu^*$ and $\\y^*$ are the mean weights and mean latent function predictions at $\\X$ respectively using the current hyperparameters. This criterion naturally weights model fit against complexity. If the prior is too tight to fit the data, we cannot draw the observations, but if it is too general then the probability of drawing any particular dataset is also low. It can become very small or large and is best computed in log form.\n",
      "\n",
      "The first hyperparameter strategy is to define a *hyper-prior* over $\\alphs$ and $\\mathbf{\\beta}$. If we carefully choose a conjugate prior (for our scenario this will be an inverse gamma distribution) over these parameters, then it is possible to marginalise over the hyperparameters (in the same way we marginalised the weights out of the predictive distribution) and obtain a closed form solution for the weight distribution. An interested student is directed to the above textbook for a detailed explanation.  Unfortunately, this distribution will be of Student-T form rather than Gaussian, so marginalising out the weights to obtain predictions will then be an exercise in numerical integration.\n",
      "\n",
      "A simpler approach that we will take here is to search for the maximum likelihood hyperparameters - this is also called the evidence approximation or type 2 maximum likelihood. The main drawback of this approach is that we loose a level of Bayesian inference compared to the above alternative. Also we must select these parameters via optimisation strategies and the problem is non-convex, meaning we may end up in a local optimum that is better than its neighbours but not the best set of parameters. However, it is a more principled method of model selection than manual tweaking.\n",
      "\n",
      "**Below we will use a simple gradient-free numerical optimisation procedure to find some hyperparameters:**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compute the log marginal likelihood in a function\n",
      "def nLML(params, theta, y, verbose=True):\n",
      "\n",
      "    alpha_ = params[0]**2\n",
      "    beta_  = ...\n",
      "    p_ = theta.shape[1]\n",
      "    I = np.eye(p_)\n",
      "    \n",
      "    if verbose:\n",
      "        print('Querying:', alpha_, beta_)\n",
      "    \n",
      "    A = ... \n",
      "    L = linalg.cholesky(A, lower=True)\n",
      "    \n",
      "    # The log determinant of a positive semi-definite matrix is given by twice the log-trace of its cholesky factorisation.\n",
      "    mn = ... # mean of weights using these parameters\n",
      "    predY = ... # mean of prediction using these parameters \n",
      "    \n",
      "    LML = ... # see equation in text above\n",
      "    return -LML\n",
      "\n",
      "# Now we will use scipy's optimize tool to find some [locally] optimal parameters\n",
      "# Note we havent defined numerical gradients (which might be a good idea later), so we will use\n",
      "# the COBYLA algorithm that can estimate them automatically\n",
      "params0 = [1., 1.] # Define some initial parameters to search from\n",
      "nLMLcriterion = lambda params: nLML(params, theta, y)\n",
      "newparams = scipy.optimize.minimize(nLMLcriterion, params0, method='COBYLA') # gradient free\n",
      "alpha = newparams.x[0]**2\n",
      "beta  = newparams.x[1]**2\n",
      "print('Local Optimum:\\n alpha = ',alpha, ', beta=', beta)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Now try using these values in the previous sections!"
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Part III: Stopping Distance Assignment"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that you know the basics of linear regression, lets apply it to a real data problem. The cell below loads a dataset of stopping distances (metres) of a car travelling at different speeds (km/h). If you were in the same car travelling at 21 km/h and you hit the brakes at a distance of 135 metres from a wall, **what is the probability of you coming to a stop before hitting the wall?**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Load data\n",
      "data = np.load('cars_stopping_dist.npz')\n",
      "speed= data['speed']\n",
      "stoppingDist= data['dist']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Solution goes here.\n",
      "# Hint: You may need to use stats.norm.cdf() at the end to determine the final probability."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}