exercise-linear-regression-MAP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Exercise - Maximum a posterior for the linear regression model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Table of Contents\n",
    "* [Introduction](#Introduction)\n",
    "* [Requirements](#Requirements) \n",
    "  * [Python Modules](#Python-Modules)\n",
    "  * [Knowledge](#Knowledge)\n",
    "  * [Modules](#Python-Modules)\n",
    "* [Recap](#Recap)\n",
    "  * [Linear Regression](#Linear-Regression)\n",
    "  * [Normal Distribution](#Normal-Distribution)\n",
    "* [Exercises](#Exercises)\n",
    "  * [Exercise - Proof](#Exercise---Covariance-Matrix)\n",
    "* [Literature](#Literature)\n",
    "* [Licenses](#Licenses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Introduction\n",
    "\n",
    "One of the first (or even the first) machine learning model students learn is the linear regression. Typically it is described from the point of view of linear algebra. In this notebook you will see, that it is also possible to describe it with a probabilistic viewpoint, regularization included."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Requirements\n",
    "\n",
    "### Knowledge\n",
    "\n",
    "To complete this exercise notebook, you should possess knowledge about the following topics.\n",
    "- Linear regression\n",
    "- Multivariate gaussian\n",
    "- Variance\n",
    "- Covariance\n",
    "- Covariance matrix\n",
    "\n",
    "The following material can help you to acquire this knowledge:\n",
    "* Gaussian, variance, covariance, covariance matrix:\n",
    " * Chapter 3 of the [Deep Learning Book](http://www.deeplearningbook.org/) [GOO16]\n",
    " * Chapter 1 of the book Pattern Recognition and Machine Learning by Christopher M. Bishop [BIS07]\n",
    "* Univariate gaussian:\n",
    " * [Video](https://www.khanacademy.org/math/statistics-probability/modeling-distributions-of-data/more-on-normal-distributions/v/introduction-to-the-normal-distribution) and the follwoing of Khan Academy [KHA18]\n",
    "* Multivariate gaussian:\n",
    " * Video PP 6.1 and following in the playlist [Probability Primer](https://www.youtube.com/watch?v=TC0ZAX3DA88&t=9s&index=35&list=PL17567A1A3F5DB5E4) of the youtube user *mathematicalmonk* [MAT18]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Python Modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# External Modules\n",
    "import daft\n",
    "from matplotlib import rc\n",
    "rc(\"font\", family=\"serif\", size=18)\n",
    "rc(\"text\", usetex=False)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_linear_regression_model():\n",
    "    \n",
    "    pgm = daft.PGM([5.9, 4.5], origin=[-1.0, 0.0], label_params={'fontsize':18})\n",
    "\n",
    "    pgm.add_node(daft.Node(\"sigma\", r\"$\\sigma$\", 4, 1, fixed=True))\n",
    "    pgm.add_node(daft.Node(\"gamma\", r\"$\\gamma$\", 0, 3, fixed=True))\n",
    "    \n",
    "    pgm.add_node(daft.Node(\"y\", r\"$y^{(i)}$\", 3, 1, scale=2))\n",
    "    pgm.add_node(daft.Node(\"w\", r\"$w$\", 1, 3, scale=2))\n",
    "\n",
    "    # Data.\n",
    "    pgm.add_node(daft.Node(\"x\", r\"$x^{(i)}$\", 3, 3, observed=True, scale=2))\n",
    "\n",
    "    # Add in the edges.\n",
    "    pgm.add_edge(\"sigma\", \"y\")\n",
    "    pgm.add_edge(\"gamma\", \"w\")\n",
    "    pgm.add_edge(\"w\", \"y\")  \n",
    "    pgm.add_edge(\"x\", \"y\")\n",
    "\n",
    "    # And a plate.\n",
    "    pgm.add_plate(daft.Plate([2.2, 0.2, 1.5, 3.5], label=r\"$i = 1, \\cdots, n$\", shift=-0.1))\n",
    "\n",
    "    # Render\n",
    "    pgm.render()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recap\n",
    "\n",
    "This chapter gives a short recap over the linear regression with regularization and the multivariate gaussian distribution. It is only meant for refreshment and is not a substitute for a deep involvement with the topic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Regression\n",
    "\n",
    "In linear regression we want to train a model to fit our data $D_{train} = \\{(x^{(1)},y^{(1)}) \\ldots (x^{(m)},y^{(m)})\\}$. The target $y$ is a single number and the features $x$ can be anything, e.g. also a single number or a $p$-dimensional vector.\n",
    "\n",
    "Considering polynomial regression for example the model could be defined with:\n",
    "$$\n",
    "h_w(x^{(i)}) = w_0 + w_1 x{(i)} + w_2 x^{(i)2} + w_3 x^{(i)3}\n",
    "$$\n",
    "\n",
    "And the loss for the whole training data as mean squared error:\n",
    "$$\n",
    "J_D(W)=\\frac{1}{2m}\\sum_{i=1}^{m}{(h_w(x^{(i)})-y^{(i)})^2}\n",
    "$$\n",
    "\n",
    "To prevent overfitting we can use $L^2$-regularization, adding penalty for too large weights:\n",
    "$$\n",
    "J_D(W)=\\frac{1}{2m}\\sum_{i=1}^{m}{(h_w(x^{(i)})-y^{(i)})^2} + \\lambda \\sum_{i=1}^3 w_i^2\n",
    "$$\n",
    "\n",
    "For simplification, it is possible to omit the normalizing term:\n",
    "$$\n",
    "J_D(W)=\\sum_{i=1}^{m}{(h_w(x^{(i)})-y^{(i)})^2} + \\lambda \\sum_{i=1}^3 w_i^2\n",
    "$$\n",
    "\n",
    "Finally, the task is to find the weights, which minimize the costs:\n",
    "$$\n",
    "\\text{arg}\\min_w \\left(\\sum_{i=1}^{m}{(h_w(x^{(i)})-y^{(i)})^2} + \\lambda \\sum_{i=1}^3 w_i^2 \\right)\n",
    "$$\n",
    "\n",
    "This is typically done with gradient descent, but that is not necessary for the context of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normal Distribution\n",
    "\n",
    "In the general case, a random variable $\\bf z$ is said to be normally distributed, if its density is:\n",
    "\n",
    "$$\n",
    "\\mathcal N({\\bf z} \\mid \\mu, \\Sigma) =  \n",
    "\\frac{1}{ \\left( (2\\pi)^p \\text{det}(\\Sigma) \\right)^{1/2}  }\n",
    "\\exp\\left( - \\frac{1}{2} ({\\bf z}-\\mu)^T \\Sigma^{-1} ({\\bf z}-\\mu)  \\right)\n",
    "$$\n",
    "\n",
    "with:\n",
    "- $\\bf z$ the column vector\n",
    "- $p$ the dimesion of the vectors $\\bf z$ \n",
    "- $\\mu$ the expected vector (or empirical mean of $n$ examples if $n \\rightarrow \\infty$)\n",
    "- $\\Sigma$ the covariance matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Regression Model\n",
    "\n",
    "For the probabilistic view of linear regression, first take a look at the following plot by executing the code cell below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_linear_regression_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two things here might be unfamiliar: The *sigma* ($\\sigma$) and the *gamma* ($\\gamma$).\n",
    "\n",
    "The *sigma* ($\\sigma$) depicts, we assume that the data is noisy, which is given by a Gaussian (motivated by the central limit theorem). So the model can written as:\n",
    "\n",
    "$$\n",
    " p(y \\mid X, w) = \\mathcal N (y \\mid X w, \\sigma^2 \\mathbb{1})\n",
    "$$\n",
    "\n",
    "with\n",
    "- the (training) data matrix $X$\n",
    "- identity matrix $\\mathbb{1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ h_w(X) = \n",
    "X w = \\begin{pmatrix}\n",
    "1 & x_1^{(1)} & x_2^{(1)} & \\dots & x_n^{(1)} \\\\\n",
    "1 &x_1^{(2)} & x_2^{(2)} &\\dots  & x_n^{(2)} \\\\\n",
    "\\vdots & \\vdots & \\vdots  & \\ddots & \\vdots \\\\\n",
    "1 & x_1^{(m)} & x_2^{(m)} & \\dots  & x_n^{(m)}\n",
    "\\end{pmatrix} \n",
    "\\begin{pmatrix}\n",
    "w_0  \\\\\n",
    "w_1  \\\\\n",
    "w_2  \\\\\n",
    "\\vdots \\\\\n",
    "w_n\n",
    "\\end{pmatrix} \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "h_w(x^{(i)}) = w^T x'^{(i)}\n",
    "$$\n",
    "\n",
    "with $$\n",
    "x'^{(i)}= \\begin{pmatrix}1 \\\\ x_1^{(i)} \\\\ x_2^{(i)} \\\\ \\vdots \\\\ x_n^{(i)}\\end{pmatrix} \n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *sigma* ($\\sigma$):\n",
    "\n",
    "For the weights we also assume Gaussian, so the prior for $w$ is:\n",
    "\n",
    "$$\n",
    "p(w) = \\mathcal N (w \\mid 0, \\gamma^2 \\mathbb 1)\n",
    "$$\n",
    "\n",
    "$w$ is the parameter vector: \n",
    "$$\n",
    "w^T = (w_0, w_1, \\dots, w_n)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MAP (Maximum a posterior)\n",
    "\n",
    "The MAP for $w$ given the data $(X, y)$ is\n",
    "$$\n",
    "\\text{arg}\\max_w p(w \\mid X, y)\n",
    "$$\n",
    "\n",
    "by Bayes rule we get\n",
    "\n",
    "$$\n",
    "\\text{arg}\\max_w p(w \\mid X, y) = \\text{arg}\\max_w \\frac{p(y \\mid X, w) p(w \\mid X)}{p(y\\mid X)}\n",
    "$$\n",
    "\n",
    "We can read from the graph $w \\perp X $ ($w$ independant of $X$). So we have:\n",
    "\n",
    "$$\n",
    "\\text{arg}\\max_w p(w \\mid X, y) = \\text{arg}\\max_w \\frac{p(y \\mid X, w) p(w)}{p(y\\mid X)}\n",
    "$$\n",
    "\n",
    "The denominator is w.r.t. $w$ a constant:\n",
    "$$\n",
    "\\text{arg}\\max_w p(w \\mid X, y) = \\text{arg}\\max_w p(y \\mid X, w) p(w)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the model assumptions above:\n",
    "$$\n",
    "\\text{arg}\\max_w p(w \\mid X, y) = \\text{arg}\\max_w  \\left( \\mathcal N (y \\mid Xw, \\sigma^2 \\mathbb{1}) \\mathcal N (w \\mid 0, \\gamma^2 \\mathbb 1)\\right)\n",
    "$$    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise - Proof\n",
    "\n",
    "**Task:**\n",
    "\n",
    "Show that the MAP of the model above is equivalent to the minimization of the squared error cost with $L^2$-regularization:\n",
    "$$\n",
    "\\text{arg}\\max_w p(w \\mid X, y) = \\text{arg}\\min_w \\left(\\sum_i^m (w^T x'^{(i)}- y^{(i)})^2 + \\lambda \\sum_i w_i^2 \\right)\n",
    "$$\n",
    "\n",
    "**Hint:**\n",
    "\n",
    "Maximization of $p(w \\mid X, y)$ is equivalent to minimization of $-\\log p(w \\mid X, y)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Literature\n",
    "\n",
    "<table>\n",
    "        <tr>\n",
    "        <td>\n",
    "            <a name=\"BIS07\"></a>[BIS07]\n",
    "        </td>\n",
    "        <td>\n",
    "            Christopher M. Bishop, Pattern recognition and machine learning, 5th Edition. Springer 2007, ISBN 9780387310732.\n",
    "        </td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td>\n",
    "            <a name=\"GOO16\"></a>[GOO16]\n",
    "        </td>\n",
    "        <td>\n",
    "            Goodfellow, Ian, et al. Deep learning. Vol. 1. Cambridge: MIT press, 2016.\n",
    "        </td>\n",
    "    </tr>\n",
    "        <tr>\n",
    "        <td>\n",
    "            <a name=\"MAT18\"></a>[MAT18]\n",
    "        </td>\n",
    "        <td>\n",
    "            mathematicalmonk. (2018, September 30). (ML 15.3) Logistic regression (binary) - intuition. And following in the playlist Machine Learning. Retrieved from https://www.youtube.com/watch?v=-Z2a_mzl9LM&list=PLD0F06AA0D2E8FFBA&t=740s&index=110\n",
    "        </td>\n",
    "    </tr>\n",
    "        </tr>\n",
    "        <tr>\n",
    "        <td>\n",
    "            <a name=\"KHA18\"></a>[KHA18]\n",
    "        </td>\n",
    "        <td>\n",
    "            Khan Academy. (2018, October 05). Deep definition of the normal distribution. And following in the playlist. Retrieved from https://www.khanacademy.org/math/statistics-probability/modeling-distributions-of-data/more-on-normal-distributions/v/introduction-to-the-normal-distribution\n",
    "        </td>\n",
    "    </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Licenses\n",
    "\n",
    "### Notebook License (CC-BY-SA 4.0)\n",
    "\n",
    "*The following license applies to the complete notebook, including code cells. It does however not apply to any referenced external media (e.g., images).*\n",
    "\n",
    "Exercise - Maximum a posterior for the linear regression model <br/>\n",
    "by Christian Herta, Klaus Strohmenger<br/>\n",
    "is licensed under a [Creative Commons Attribution-ShareAlike 4.0 International License](http://creativecommons.org/licenses/by-sa/4.0/).<br/>\n",
    "Based on a work at https://gitlab.com/deep.TEACHING.\n",
    "\n",
    "\n",
    "### Code License (MIT)\n",
    "\n",
    "*The following license only applies to code cells of the notebook.*\n",
    "\n",
    "Copyright 2018 Christian Herta, Klaus Strohmenger\n",
    "\n",
    "Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the \"Software\"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:\n",
    "\n",
    "The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.\n",
    "\n",
    "THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "deep_teaching_kernel",
   "language": "python",
   "name": "deep_teaching_kernel"
  },
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}