Lab1_Cornelissen_deGroot_Mulder_Westerdiep.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 1: Linear Regression and Overfitting\n",
    "\n",
    "### Machine Learning and Pattern Recognition, September 2015\n",
    "\n",
    "* The lab exercises should be made in groups of two or three people.\n",
    "* The deadline is sunday September 20, 23:59.\n",
    "* Assignment should be sent to Philip Versteeg. (p.j.j.p.versteeg@uva.nl) The subject line of your email should be \"\\#lab\\_lastname1\\_lastname2\\_lastname3\".\n",
    "* Put your and your teammates' names in the body of the email\n",
    "* Attach the .IPYNB (IPython Notebook) file containing your code and answers. Naming of the file follows the same rule as the subject line. For example, if the subject line is \"lab01\\_Kingma\\_Hu\", the attached file should be \"lab01\\_Kingma\\_Hu.ipynb\". Only use underscores (\"\\_\") to connect names, otherwise the files cannot be parsed.\n",
    "* Make sure we can run your notebook / scripts!\n",
    "\n",
    "Notes on implementation:\n",
    "\n",
    "* You should write your code and answers in this IPython Notebook: http://ipython.org/notebook.html. If you have problems, please contact me.\n",
    "* Please write your answers right below the questions.\n",
    "* Among the first lines of your notebook should be \"%pylab inline\". This imports all required modules, and your plots will appear inline.\n",
    "* Refer to last week's lab notes, i.e. http://docs.scipy.org/doc/, if you are unsure about what function to use. There are different correct ways to implement each problem!\n",
    "* For this lab, your regression solutions should be in closed form, i.e., should not perform iterative gradient-based optimization but find the exact optimum directly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\bPhi}{\\mathbf{\\Phi}}$\n",
    "$\\newcommand{\\bx}{\\mathbf{x}}$\n",
    "$\\newcommand{\\bw}{\\mathbf{w}}$\n",
    "$\\newcommand{\\bt}{\\mathbf{t}}$\n",
    "$\\newcommand{\\by}{\\mathbf{y}}$\n",
    "$\\newcommand{\\bm}{\\mathbf{m}}$\n",
    "$\\newcommand{\\bS}{\\mathbf{S}}$\n",
    "$\\newcommand{\\bI}{\\mathbf{I}}$\n",
    "\n",
    "## Part 1: Polynomial Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1. Generate sinusoidal data (5 points)\n",
    "Write a method `gen_sinusoidal(N)` that generates toy data like in fig 1.2 of the MLPR book. The method should have a parameter $N$, and should return $N$-dimensional vectors $\\bx$ and $\\bt$, where $\\bx$ contains evenly spaced values from 0 to (including) 2$\\pi$, and the elements $t_i$ of $\\bt$ are distributed according to:\n",
    "\n",
    "$$t_i \\sim \\mathcal{N}(\\mu_i, \\sigma^2)$$\n",
    "\n",
    "where $x_i$ is the $i$-th elements of $\\bf{x}$, the mean $\\mu_i = sin(x_i)$ and the standard deviation $\\sigma = 0.2$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def gen_sinusoidal(N):\n",
    "    x = linspace(0, 2*pi, N)\n",
    "    t = random.normal(sin(x), .2)\n",
    "    return x, t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x, t = gen_sinusoidal(500)\n",
    "plot(x, t, '+')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Polynomial regression (15 points)\n",
    "\n",
    "Write a method `fit_polynomial(x, t, M)` that finds the maximum-likelihood solution of an _unregularized_ $M$-th order polynomial for some dataset `x`. The error function to minimize w.r.t. $\\bw$ is:\n",
    "\n",
    "$E(\\bw) = \\frac{1}{2} (\\bPhi\\bw - \\bt)^T(\\bPhi\\bw - \\bt)$\n",
    "\n",
    "where $\\bPhi$ is the _feature matrix_ (or _design matrix_) as explained in the MLPR book at section 3.1.1, $\\bt$ is the vector of target values. Your method should return a vector $\\bw$ with the maximum-likelihood parameter estimates. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fit_polynomial(x, t, M):\n",
    "    # Alternative: linalg.inv( Phi.T.dot(Phi) ).dot(Phi.T)\n",
    "    Phi = fliplr(vander(x, M+1)) # M+1 for M'th order polynomial\n",
    "    return linalg.pinv(Phi).dot(t)  # pinv = Moore-Penrose inverse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Plot (5 points)\n",
    "Sample a dataset with $N=9$, and fit four polynomials with $M \\in (0, 1, 3, 9)$.\n",
    "For each value of $M$, plot the prediction function, along with the data and the original sine function. The resulting figure should look similar to fig 1.4 of the MLPR book. Note that you can use matplotlib's `plt.pyplot(.)` functionality for creating grids of figures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def show_regression(x, t, w, xs=False):\n",
    "    if xs == False:\n",
    "        xs = linspace(0, 2*pi, 100)\n",
    "    plot(xs, polyval(w[::-1], xs), '-r', lw=2)\n",
    "    plot(x, t, 'o')\n",
    "    plot(xs, sin(xs), '-g')\n",
    "    show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x, t = gen_sinusoidal(9)\n",
    "    \n",
    "for M in [0, 1, 3, 9]:\n",
    "    w = fit_polynomial(x, t, M)\n",
    "    show_regression(x, t, w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Regularized linear regression (10 points)\n",
    "\n",
    "Write a method `fit_polynomial_reg(x, t, M, lamb)` that fits a _regularized_ $M$-th order polynomial to the sinusoidal data, as discussed in the lectures, where `lamb` is the regularization term _lambda_. (Note that 'lambda' cannot be used as a variable name in Python since it has a special meaning). The error function to minimize w.r.t. $\\bw$:\n",
    "\n",
    "$E(\\bw) = \\frac{1}{2} (\\bPhi\\bw - \\bt)^T(\\bPhi\\bw - \\bt) + \\frac{\\lambda}{2} \\mathbf{w}^T \\mathbf{w}$\n",
    "\n",
    "For background, see section 3.1.4 of the MLPR book."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fit_polynomial_reg(x, t, M, lamb):\n",
    "    Phi = fliplr( vander(x, M+1) )\n",
    "    return linalg.inv( lamb * eye(M+1) + Phi.T.dot(Phi) ).dot(Phi.T).dot(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### 1.5 Model selection by cross-validation (10 points)\n",
    "Use cross-validation to find a good choice of $M$ and $\\lambda$, given a dataset of $N=9$ datapoints generated with `gen_sinusoidal(9)`. You should write a function that tries (loops over) a reasonable range of choices of $M$ and $\\lambda$, and returns the choice with the best cross-validation error. In this case you can use $K=9$ folds, corresponding to _leave-one-out_ crossvalidation.\n",
    "\n",
    "You can let $M \\in (0, 1, ..., 10)$, and let $\\lambda \\in (e^{-10}, e^{-9}, ..., e^{0})$.\n",
    "\n",
    "To get you started, here's a method you can use to generate indices of cross-validation folds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def kfold_indices(N, k):\n",
    "    all_indices = np.arange(N,dtype=int)\n",
    "    np.random.shuffle(all_indices)\n",
    "    idx = np.floor(np.linspace(0,N,k+1))\n",
    "    train_folds = []\n",
    "    valid_folds = []\n",
    "    for fold in range(k):\n",
    "        valid_indices = all_indices[idx[fold]:idx[fold+1]]\n",
    "        valid_folds.append(valid_indices)\n",
    "        train_folds.append(np.setdiff1d(all_indices, valid_indices))\n",
    "    return train_folds, valid_folds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def select_model(x, t):\n",
    "\n",
    "    N = len(x)\n",
    "    train_folds, valid_folds = kfold_indices(N, N)\n",
    "    Ms = range(0,11)\n",
    "    lambdas = [exp(-10 + i) for i in range(0,11)];\n",
    "    errors = zeros((len(Ms), len(lambdas)))\n",
    "\n",
    "    for i, M in enumerate(Ms):\n",
    "        for j, lamb in enumerate(lambdas):\n",
    "            err = []\n",
    "            for k, idx in enumerate(train_folds):\n",
    "\n",
    "                # Regularized regression on training data\n",
    "                train_xs = x[idx]\n",
    "                train_ts = t[idx]\n",
    "                w = fit_polynomial_reg(train_xs, train_ts, M+1, lamb)\n",
    "\n",
    "                # Compare regression with validation data using RMS error\n",
    "                valid_x = x[valid_folds[k]][0]\n",
    "                valid_t = t[valid_folds[k]][0]\n",
    "                estimate_t = polyval(w[::-1], valid_x)\n",
    "\n",
    "                err += [sqrt(mean(square(estimate_t - valid_t)))]\n",
    "            errors[i,j] = mean(err)\n",
    "\n",
    "    # Transpose: Ms vertically and lambdas horizontally\n",
    "    errors = errors.T\n",
    "    min_id = argmin(errors)\n",
    "    M_id = min_id % errors.shape[1] # horizontal position\n",
    "    lamb_id = min_id // errors.shape[1] # vertical position\n",
    "\n",
    "    # Plot this\n",
    "    fig, ax = plt.subplots(figsize=(6,5))\n",
    "    heatmap = ax.pcolor(log(errors), cmap=plt.cm.Blues)\n",
    "\n",
    "    # Highlight optimal choice\n",
    "    ax.add_patch(matplotlib.patches.Rectangle(\n",
    "            (M_id, lamb_id), 1, 1, facecolor=\"none\", edgecolor='red', linewidth=2))\n",
    "\n",
    "    ax.set_xticks(arange(len(Ms))+0.5, minor=False);\n",
    "    ax.set_yticks(arange(len(lambdas))+0.5, minor=False);\n",
    "    ax.set_xticklabels(Ms, minor=False);\n",
    "    ax.set_yticklabels(['$\\exp('+str(int(e))+')$' for e in log(lambdas)], minor=False);\n",
    "    cbar = plt.colorbar(heatmap)\n",
    "    title('$\\log($error$)$ for different values of $\\lambda$ and $M$.\\\n",
    "            \\n$\\lambda=\\exp(%s)$ and $M=%s$ is optimal (red square)\\n'\n",
    "            % (log(lambdas[lamb_id]), Ms[M_id]))\n",
    "    show()\n",
    "    \n",
    "    return Ms[M_id], lambdas[lamb_id]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a comprehensible plot of the cross-validation error for each choice of $M$ and $\\lambda$. Highlight the best choice. \n",
    "\n",
    "_Question_: Explain over-fitting and underfitting, illuminated by your plot. Explain the relationship with model bias and model variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.6 Plot best cross-validated fit (5 points)\n",
    "\n",
    "For some dataset with $N = 9$, plot the model with the optimal $M$ and $\\lambda$ according to the cross-validation error, using the method you just wrote. Let the plot make clear which $M$ and $\\lambda$ were found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x,t = gen_sinusoidal(9)\n",
    "M, lamb =  select_model(x, t)\n",
    "w = fit_polynomial_reg(x, t, M, lamb)\n",
    "\n",
    "title('Regularized least squares regression \\n\\\n",
    "    using $\\lambda=\\exp(%s)$ and $M=%s$' % (log(lamb), M))\n",
    "show_regression(x, t, w)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Bayesian Linear (Polynomial) Regression\n",
    "\n",
    "### 2.1 Sinusoidal data 2 (5 points)\n",
    "\n",
    "Write a function `gen_sinusoidal2(N)` that behaves identically to `gen_sinusoidal(N)` except that the generated values $x_i$ are not linearly spaced, but drawn from a uniform distribution between $0$ and $2 \\pi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def gen_sinusoidal2(N):\n",
    "    x = sort(random.uniform(0, 2*pi, N))\n",
    "    t = random.normal(sin(x),.2)\n",
    "    return x,t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x, t = gen_sinusoidal2(500)\n",
    "plot(x, t, '+')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Compute Posterior (15 points)\n",
    "\n",
    "You're going to implement a Bayesian linear regression model, and fit it to the sinusoidal data. Your regression model has a zero-mean isotropic Gaussian prior over the parameters, governed by a single (scalar) precision parameter $\\alpha$, i.e.:\n",
    "\n",
    "$$p(\\bw \\;|\\; \\alpha) = \\mathcal{N}(\\bw \\;|\\; 0, \\alpha^{-1} \\bI)$$\n",
    "\n",
    "The covariance and mean of the posterior are given by:\n",
    "\n",
    "$$\\bS_N= \\left( \\alpha \\bI + \\beta \\bPhi^T \\bPhi \\right)^{-1} $$\n",
    "$$\\bm_N = \\beta\\; \\bS_N \\bPhi^T \\bt$$\n",
    "\n",
    "where $\\alpha$ is the precision of the predictive distribution, and $\\beta$ is the noise precision. \n",
    "See MLPR chapter 3.3 for background.\n",
    "\n",
    "Write a method `fit_polynomial_bayes(x, t, M, alpha, beta)` that returns the mean $\\bm_N$ and covariance $\\bS_N$ of the posterior for a $M$-th order polynomial, given a dataset, where `x`, `t` and `M` have the same meaning as in question 1.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fit_polynomial_bayes(x, t, M, alpha, beta):\n",
    "    Phi = fliplr(vander(x, M+1))\n",
    "    S_N = linalg.inv(alpha * eye(M+1) + beta * Phi.T.dot(Phi))\n",
    "    m_N = beta * S_N.dot( Phi.T.dot(t) )\n",
    "    return S_N, m_N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Prediction (10 points)\n",
    "\n",
    "The predictive distribution of Bayesian linear regression is:\n",
    "\n",
    "$$ p(t \\;|\\; \\bx, \\bt, \\alpha, \\beta) = \\mathcal{N}(t \\;|\\; \\bm_N^T \\phi(\\bx), \\sigma_N^2(\\bx))$$\n",
    "\n",
    "$$ \\sigma_N^2 = \\frac{1}{\\beta} + \\phi(\\bx)^T \\bS_N \\phi(\\bx) $$\n",
    "\n",
    "where $\\phi(\\bx)$ are the computed features for a new datapoint $\\bx$, and $t$ is the predicted variable for datapoint $\\bx$. \n",
    "\n",
    "Write a function that `predict_polynomial_bayes(x, m, S, beta)` that returns the predictive mean and variance given a new datapoint `x`, posterior mean `m`, posterior variance `S` and a choice of model variance `beta`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def predict_polynomial_bayes(x, m, S, beta):\n",
    "    N = S.shape[0] # Note: N = M+1\n",
    "    phi_x = fliplr(vander([x], N))[0]\n",
    "    sigma = phi_x.T.dot(S).dot(phi_x)\n",
    "    mu = m.T.dot(phi_x)\n",
    "    return mu, sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Plot predictive distribution (10 points)\n",
    "\n",
    "a) (5 points) Generate 7 datapoints with `gen_sinusoidal2(7)`. Compute the posterior mean and covariance for a Bayesian polynomial regression model with $M=5$, $\\alpha=\\frac{1}{2}$ and $\\beta=\\frac{1}{0.2^2}$.\n",
    "Plot the Bayesian predictive distribution, where you plot (for $x$ between 0 and $2 \\pi$) $t$'s predictive mean and a 1-sigma predictive variance using `plt.fill_between(..., alpha=0.1)` (the alpha argument induces transparency).\n",
    "\n",
    "Include the datapoints in your plot.\n",
    "\n",
    "b) (5 points) For a second plot, draw 100 samples from the parameters' posterior distribution. Each of these samples is a certain choice of parameters for 5-th order polynomial regression. \n",
    "Display each of these 100 polynomials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta = 25.0\n",
    "x, t = gen_sinusoidal2(7)\n",
    "S_N, m_N = fit_polynomial_bayes(x, t, 5, .5, beta)\n",
    "\n",
    "xs = linspace(0, 2*pi, 100)\n",
    "# Matrix with mu's and sigma's as columns\n",
    "stats = np.array([predict_polynomial_bayes(d, m_N, S_N, beta) for d in xs])\n",
    "\n",
    "figsize(10,5)\n",
    "plot(x, t, 'o')\n",
    "plot(xs, sin(xs), '--', color='black')\n",
    "plot(xs, stats[:,0], 'r-')\n",
    "ylim(-1.5,1.5)\n",
    "xlim(0,2*pi)\n",
    "fill_between(xs, stats[:,0] + stats[:,1], stats[:,0] - stats[:,1], alpha=.2, color='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def show_regression(x, t, w, xs=False):\n",
    "    if xs == False:\n",
    "        xs = linspace(0, 2*pi, 100)\n",
    "    plot(xs, polyval(w[::-1], xs), '-r', lw=2)\n",
    "    plot(x, t, 'o')\n",
    "    plot(xs, sin(xs), '-g')\n",
    "    show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta = 25.0\n",
    "x, t = gen_sinusoidal2(7)\n",
    "S_N, m_N = fit_polynomial_bayes(x, t, 5, .5, beta)\n",
    "\n",
    "xs = linspace(0, 2*pi, 100)\n",
    "ws = random.multivariate_normal(m_N, S_N, 100)\n",
    "for w in ws :\n",
    "    plot(xs, polyval(w[::-1], xs), 'r-', alpha=.1)\n",
    "\n",
    "figsize(10,5)\n",
    "ylim(-2,2)\n",
    "xlim(0,2*pi)\n",
    "plot(x, t, 'o', color='black')\n",
    "plot(xs, sin(xs), '-', color='black')\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 Additional questions (10 points)\n",
    "\n",
    "a) (5 points) Why is $\\beta=\\frac{1}{0.2^2}$ the best choice of $\\beta$ in section 2.4?\n",
    "\n",
    "b) (5 points) In the case of Bayesian linear regression, both the posterior of the parameters $p(\\bw \\;|\\; \\bt, \\alpha, \\beta)$ and the predictive distribution $p(t \\;|\\; \\bw, \\beta)$ are Gaussian. In consequence (and conveniently), $p(t \\;|\\; \\bt, \\alpha, \\beta)$ is also Gaussian (See MLPR section 3.3.2 and homework 2 question 4). This is actually one of the (rare) cases where we can make Bayesian predictions without resorting to approximative methods.\n",
    "\n",
    "Suppose you have to work with some model $p(t\\;|\\;x,\\bw)$ with parameters $\\bw$, where the posterior distribution $p(\\bw\\;|\\;\\mathcal{D})$ given dataset $\\mathcal{D}$ can not be integrated out when making predictions, but where you can still generate samples from the posterior distribution of the parameters. Explain how you can still make approximate Bayesian predictions using samples from the parameters' posterior distribution.\n",
    "\n",
    "## Answers\n",
    "a) In section 2.4, we considered data points generated by the `gen_sinusoidal2` method. In its implementation, we considered a regression function $h(x) = \\sin(x)$ (or $h(x) = \\sin(2 \\pi x)$, depending on whether you follow Bishop or the lab exercise) with target value deviation $\\sigma = 0.2$. In the book, it is said that the precision $\\beta$ should be set to $\\beta = 1/\\sigma^2$ (see page 24). This is supported by the deduction in page 142, equation (3.21): we see that maximizing the log likelihood of (3.11) wrt $\\beta$ yields a $\\beta_{ML}$ that is of exactly this form. So in conclusion, setting $\\beta = 1/(0.2)^2$ is consistent with the theory.\n",
    "\n",
    "b) We can use a sampling scheme (for example Markov chain Monte Carlo, chapter 11.2) to draw a set $z^{(l)}$ of $L$ independent samples from $p(\\bw\\;|\\;\\mathcal{D})$. Now we can approximate $p(t \\;|\\; \\bt, \\alpha, \\beta)$ with $\\frac{1}{L}\\sum_{l=1}^L p(t \\;|\\; z^{(l)}, \\beta)$. This is analogous to equation (11.2). Therefore, this approximation has the correct expected value."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}