01_linear_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/ascarrambad/ml-21-22/blob/main/01_linear_regression.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pLNkc6gv9x5t"
   },
   "source": [
    "# Machine Learning SP 2021/2022\n",
    "\n",
    "Prof. Cesare Alippi<br>\n",
    "Giorgia Adorni ([`giorgia.adorni@usi.ch`](mailto:giorgia.adorni@usi.ch))<br>\n",
    "Luca Butera ([`luca.butera@usi.ch`](mailto:luca.butera@usi.ch))<br>\n",
    "Matteo Riva ([`matteo.riva@usi.ch`](mailto:matteo.riva@usi.ch))\n",
    "\n",
    "---\n",
    "\n",
    "# Lab 01: Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LSzcDPkmJ6qz"
   },
   "source": [
    "# 01.A) Let's collect some data\n",
    "... or let someone do it for us :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aXrZ4FvtNCDT"
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_boston\n",
    "boston = load_boston()\n",
    "X = boston.data\n",
    "y = boston.target\n",
    "\n",
    "print(boston.DESCR)\n",
    "print(\"X shape:\", X.shape)\n",
    "print(\"y shape:\", y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KBJX9OqvOcZu"
   },
   "source": [
    "Data set of $n=506$ observations $\\{(\\mathbf x_1, y_1), (\\mathbf x_2, y_2) ,\\dots,(\\mathbf x_n, y_n)\\}$, where $\\mathbf x_i\\in\\mathbb R^{d}$, with $d=13$ and $y_i\\in\\mathbb R$. All the observations are stack to form\n",
    "\n",
    "$$\n",
    "X = \\left[ \n",
    "\\begin{array}{c}\n",
    "\\mathbf x_1\\\\\n",
    "\\mathbf x_2\\\\\n",
    "\\vdots \\\\\n",
    "\\mathbf x_n\n",
    "\\end{array}\n",
    "\\right] \n",
    "\\in \\mathbb{R}^{n\\times d},\n",
    "\\qquad \n",
    "Y = \\left[ \n",
    "\\begin{array}{c}\n",
    "y_1 \\\\\n",
    "y_2 \\\\\n",
    "\\vdots \\\\\n",
    "y_n\n",
    "\\end{array}\n",
    "\\right] \n",
    "\\in \\mathbb{R}^{n}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "qRFryiYXNCdf"
   },
   "outputs": [],
   "source": [
    "# Let's consider only the RM index, for now.\n",
    "x = X[:, 5]  \n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.scatter(x, y, c='k', marker='.');\n",
    "plt.xlabel(\"RM\");\n",
    "plt.ylabel(\"thousand dollars\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "FWfVU_hXOlu8"
   },
   "source": [
    "# 01.B) System model\n",
    "\n",
    "We assume there is a function $g(x)$ that links RM index to the house price:\n",
    "$$\n",
    "y = g(x) + \\eta\n",
    "$$\n",
    "where $\\eta \\sim N(0, \\sigma^2_\\eta)$.\n",
    "\n",
    "Every line in the plane (except the vertical ones) can be written in the form\n",
    "$$f(x; \\boldsymbol \\theta) = \\theta_0 + \\theta_1 x$$ \n",
    "with $\\boldsymbol \\theta = (\\theta_0, \\theta_1)$ and $\\theta_0,\\theta_1 \\in\\mathbb R$.\n",
    "\n",
    "We also assume that $g(.)$ is linear, that is, there exists $\\boldsymbol \\theta^o$ so that $ g(x) = f(x; \\boldsymbol \\theta^o)$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9pQjXiWDRkxf"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def my_lin_fun(x, theta):\n",
    "    y = theta[0] + x * theta[1]\n",
    "    return y\n",
    "\n",
    "# some guesses\n",
    "xx = np.array([4., 9.])\n",
    "plt.plot(xx, my_lin_fun(xx, [.1, 3.]))\n",
    "plt.plot(xx, my_lin_fun(xx, [-1., 5.]))\n",
    "plt.scatter(x, y, c='k', marker='.');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2XM0AHt6Rk8T"
   },
   "source": [
    "# 01.C) Model approximation\n",
    "\n",
    "Given a new value of $x$ our best prediction for $y$ is\n",
    "$$\\hat y = E[y] = E[f(x; \\boldsymbol \\theta^o) + \\eta] = f(x; \\boldsymbol \\theta^o) + E[\\eta] = f(x; \\boldsymbol \\theta^o).$$\n",
    "\n",
    "Since ${\\boldsymbol \\theta^o}$ is unknown, we estimate it from the data, by minimising \n",
    "$$ \\hat{\\boldsymbol \\theta} = \\mathop{\\mathrm{arg\\,min}}_{\\boldsymbol \\theta} \\sum_{i=1}^n \\left\\lVert y_i - f(x;\\boldsymbol \\theta) \\right\\rVert^2_2 $$\n",
    "\n",
    "Finally, we predict new house prices with \n",
    "$$\\hat y = f\\left(x; \\hat{\\boldsymbol \\theta}\\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PfdvHYsQiKSb"
   },
   "source": [
    "## 01.D) Parameter estimation\n",
    "\n",
    "__Data in compact form:__ Prepending a '1' to each $\\mathbf x$, then any \n",
    "$f(\\mathbf x,\\boldsymbol \\theta)= \\mathbf x^\\top \\boldsymbol \\theta$. In fact, $\\mathbf x^\\top \\boldsymbol \\theta = \\theta_0 1 +\\theta_1 x$.\n",
    "\n",
    "We showed in class that the solution $\\hat{\\boldsymbol \\theta}$ to $\\mathop{\\mathrm{arg\\,min}}_{\\boldsymbol \\theta} \\sum_{i=1}^n \\left\\lVert y_i - f(x;\\boldsymbol \\theta) \\right\\rVert^2_2$ \n",
    "can be found by solving the linear system\n",
    "$$\n",
    "X^\\top Y = X^\\top X \\boldsymbol \\theta\n",
    "$$\n",
    "with respect to the $\\boldsymbol \\theta$.\n",
    "\n",
    "In our 1-D case,\n",
    "$$\n",
    "X = \\left[ \n",
    "\\begin{array}{c}\n",
    "1, x_1 \\\\\n",
    "1, x_2 \\\\\n",
    "\\vdots \\\\\n",
    "1, x_n\n",
    "\\end{array}\n",
    "\\right] \n",
    "\\in \\mathbb{R}^{n\\times 2},\n",
    "\\qquad \n",
    "Y = \\left[ \n",
    "\\begin{array}{c}\n",
    "y_1 \\\\\n",
    "y_2 \\\\\n",
    "\\vdots \\\\\n",
    "y_n\n",
    "\\end{array}\n",
    "\\right] \n",
    "\\in \\mathbb{R}^{n}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "FfBx7O5sZw6q"
   },
   "outputs": [],
   "source": [
    "# Solving the linear system\n",
    "x_col_vec = x.reshape(-1, 1)\n",
    "ones_col_vec = np.ones(shape=(y.shape[0], 1))\n",
    "X = np.hstack((ones_col_vec, x_col_vec))\n",
    "theta_hat = np.linalg.solve(a=X.T.dot(X), b=X.T.dot(y))  # solves ax=b wrt x\n",
    "print(theta_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "38pagQFna7bv"
   },
   "outputs": [],
   "source": [
    "# plot the result\n",
    "plt.plot(xx, my_lin_fun(xx, theta_hat), 'g')\n",
    "plt.scatter(x, y, c='k', marker='.');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "vX6kECiZbzfI"
   },
   "outputs": [],
   "source": [
    "# There are also libraries for this\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# init the model\n",
    "lr = LinearRegression(fit_intercept=False)  \n",
    "\n",
    "# estimate parameters\n",
    "lr.fit(X, y)\n",
    "theta_hat2 = lr.coef_\n",
    "\n",
    "print('theta1 = {}'.format(theta_hat))\n",
    "print('theta2 = {}'.format(theta_hat2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7FLmt1kEdCn9"
   },
   "outputs": [],
   "source": [
    "# We can also avoid creating a column of ones\n",
    "lr = LinearRegression(fit_intercept=True)  # default is True  \n",
    "lr.fit(x_col_vec, y)\n",
    "theta_hat3 = [lr.intercept_, lr.coef_[0]]\n",
    "print('theta3 = {}'.format(theta_hat3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "M2nZ-Fm0dm1x"
   },
   "source": [
    "# 01.E) More generally \n",
    "\n",
    "## i. Multidimensional data\n",
    "\n",
    "When the input (regressor) $\\mathbf x=[x_1,\\dots,x_d]$ is $d$-dimensional, then\n",
    "$$y = f(\\mathbf x, \\boldsymbol \\theta)  + \\eta = \\mathbf x^\\top \\boldsymbol \\theta = x_1 \\theta_1 + x_2 \\theta_2+ \\dots + x_d\\theta_d + \\eta.$$\n",
    "with $\\mathbf x,\\boldsymbol \\theta \\in\\mathbb R^d$.\n",
    "\n",
    "When also the output (target) $\\mathbf y$ is multidimensional:\n",
    "$$\\mathbf y = f(\\mathbf x, \\Theta) +\\eta= \\mathbf x^\\top \\Theta +\\eta$$\n",
    "with $\\mathbf y\\in\\mathbb R^f$, $\\Theta \\in\\mathbb R^{d\\times f}$ is a matrix and $\\eta\\sim N(0,I\\sigma_\\eta)$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "jrnkOkEhg9CE"
   },
   "outputs": [],
   "source": [
    "# We can also avoid creating a column of ones\n",
    "lr_1d = LinearRegression()  \n",
    "lr_1d.fit(x_col_vec, y)\n",
    "y_pred_1d = lr_1d.predict(x_col_vec)\n",
    "err_1d = ((y_pred_1d - y)**2).sum()\n",
    "\n",
    "# Let's try adding LSTAT...\n",
    "lr_2d = LinearRegression()  \n",
    "X2 = np.vstack((boston.data[:, 5], boston.data[:, 12])).T\n",
    "lr_2d.fit(X2, y)\n",
    "y_pred_2d = lr_2d.predict(X2)\n",
    "err_2d = ((y_pred_2d - y)**2).sum()\n",
    "\n",
    "# ... and finally with all the features \n",
    "lr_md = LinearRegression()  \n",
    "Xall = boston.data\n",
    "lr_md.fit(Xall, y)\n",
    "y_pred_md = lr_md.predict(Xall)\n",
    "err_md = ((y_pred_md - y)**2).sum()\n",
    "\n",
    "print(err_1d)\n",
    "print(err_2d)\n",
    "print(err_md)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Lh0h2IsEimFV"
   },
   "outputs": [],
   "source": [
    "# plot the result\n",
    "%matplotlib inline\n",
    "\n",
    "# 1-d\n",
    "xx = np.array([3., 10.])\n",
    "plt.plot(x_col_vec, y_pred_1d, 'g')\n",
    "plt.scatter(x_col_vec, y, c='k', marker='.');\n",
    "\n",
    "# 2-d\n",
    "from mpl_toolkits.mplot3d import Axes3D  # this is necessary for 3-d plots \n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d');\n",
    "ax.view_init(elev=10., azim=0)\n",
    "ax.set_xlabel('x_5')\n",
    "ax.set_ylabel('x_12')\n",
    "ax.set_zlabel('y')\n",
    "ax.plot_trisurf(X2[:, 0], X2[:, 1], y_pred_2d, alpha=0.3, label='est fun');\n",
    "ax.scatter(X2[:, 0], X2[:, 1], y, c='k', marker='.');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XXOYwKgMdnCZ"
   },
   "source": [
    "## ii. 'Linear' means linear in the parameters\n",
    "\n",
    "The regressor can be of the form \n",
    "$$\\boldsymbol \\phi(\\mathbf x) = [\\phi_1(\\mathbf x), \\phi_2(\\mathbf x), \\dots, \\phi_m(\\mathbf x)]$$ \n",
    "for any collection of functions $\\phi_1,\\dots,\\phi_m$. \n",
    "Function $f$ become\n",
    "$$f(\\mathbf x, \\boldsymbol \\theta) = \\boldsymbol \\theta^\\top\\boldsymbol \\phi(\\mathbf x) =  \\theta_1 \\phi_1(\\mathbf x) + \\theta_2 \\phi_d(\\mathbf x) + \\dots + \\theta_d \\phi_d(\\mathbf x).$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cN5W-BU6oUZT"
   },
   "source": [
    "\n",
    "### Example: Polynomials\n",
    "\n",
    "$$f(x;\\boldsymbol \\theta) = \\theta_0 + x \\theta_1 + x^2 \\theta_2 + \\dots + x^d \\theta_d$$ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "nKYw9sXpmtGn"
   },
   "outputs": [],
   "source": [
    "def pol_fun(x):\n",
    "    return -1 -x - .1 * x**2 + .5*x**3\n",
    "\n",
    "# generate data\n",
    "n = 100\n",
    "X = np.linspace(-1, 2, n).reshape(n,1) \n",
    "interval = np.linspace(-1, 2, 100).reshape(100,1)\n",
    "sigma = 0.3\n",
    "eta = np.random.normal(loc=0, scale=sigma, size=(n,1))\n",
    "Y = pol_fun(X) + eta\n",
    "\n",
    "# create regressor\n",
    "degree = 1\n",
    "\n",
    "# generate polynomial features\n",
    "from sklearn.preprocessing import PolynomialFeatures \n",
    "pol_feat = PolynomialFeatures(degree=degree, include_bias=False) \n",
    "Xpol = pol_feat.fit_transform(X)\n",
    "\n",
    "# estimate parameter\n",
    "lr = LinearRegression()\n",
    "lr.fit(Xpol, Y)\n",
    "Y_est = lr.predict(pol_feat.transform(interval[:,...]))\n",
    "\n",
    "# plot results\n",
    "plt.plot(interval, pol_fun(interval), label='true fun')\n",
    "plt.scatter(X, Y, label='noisy data', c='k', marker='.', alpha=0.5)\n",
    "plt.plot(interval, Y_est, label='est fun')\n",
    "plt.ylim((-2, 0.8))\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QLuQeYuNn5hc"
   },
   "source": [
    "# 01.F) Regularizations\n",
    "\n",
    "## 1. Ridge regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "IQuCKsazC3KC"
   },
   "outputs": [],
   "source": [
    "# generate data [200 100 50 20 10]\n",
    "n = 200\n",
    "X = np.linspace(-1, 2, n).reshape(n,1) \n",
    "sigma = 0.3\n",
    "eps = np.random.normal(loc=0, scale=sigma, size=(n,1))\n",
    "Y = pol_fun(X) + eps\n",
    "\n",
    "# create regressor\n",
    "degree = 5\n",
    "pol_feat = PolynomialFeatures(degree=degree, include_bias=False) \n",
    "Xpol = pol_feat.fit_transform(X)\n",
    "\n",
    "# linear regression\n",
    "lr = LinearRegression()\n",
    "lr.fit(Xpol, Y)\n",
    "Y_est = lr.predict(Xpol)\n",
    "\n",
    "# ridge regression\n",
    "from sklearn.linear_model import Ridge \n",
    "rid = Ridge()\n",
    "rid.fit(Xpol, Y)\n",
    "Y_rid_est = rid.predict(Xpol)\n",
    "\n",
    "# plot results\n",
    "plt.plot(X, pol_fun(X), label='true fun')\n",
    "plt.scatter(X, Y, label='noisy data', c='k', marker='.')\n",
    "plt.plot(X, Y_est, label='lin.reg.')\n",
    "plt.plot(X, Y_rid_est, label='ridge reg.')\n",
    "plt.legend()\n",
    "\n",
    "# estimated theta\n",
    "print('theta_lr  = {}, {}'.format(lr.intercept_, lr.coef_[0]))\n",
    "print('theta_rid = {}, {}'.format(rid.intercept_, rid.coef_[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1lQ_Uy_Anlyo"
   },
   "source": [
    "## ii. Lasso regression\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "P-Y4YxP0y-bL"
   },
   "outputs": [],
   "source": [
    "# Try yourself!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 01.G) Liner regression with gradient descent\n",
    "\n",
    "Consider the following simple regression problem, where:\n",
    "$$g(x) = 2 + {3 \\over 10}  x$$\n",
    "\n",
    "First we define $g(x)$ in Python, we get some data $y_i = g(x_i) + \\eta$ and we use them to fit a model $f(x; \\boldsymbol \\theta) = \\theta_0 + \\theta_1 x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "0Ks9V9Wzf7ZA"
   },
   "outputs": [],
   "source": [
    "def g(x):\n",
    "    y = 2 + 0.3 * x\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "wUZW0MHdiJPe"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "np.random.seed(1233)\n",
    "\n",
    "# create n observations\n",
    "n = 20  # number of data points\n",
    "\n",
    "# regressor\n",
    "x = np.linspace(-1, 1, n)  \n",
    "# noise\n",
    "sigma = 0.1  # std of the noise\n",
    "eta = np.random.normal(loc=0, scale=sigma, size=n)\n",
    "# response\n",
    "y = g(x) + eta\n",
    "\n",
    "# plot\n",
    "plt.plot(x, g(x), label='true fun', color ='red');         # real function\n",
    "plt.scatter(x, y, label='noisy data', color='black');      # data affected by noise\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "GTSR1GbmX-Pl"
   },
   "source": [
    "We can estimate $\\hat{\\boldsymbol \\theta}$ using the tools that we saw in the first lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "K5ngjm0uiJc3"
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "x = x.reshape(-1, 1)\n",
    "y = y.reshape(-1, 1)\n",
    "\n",
    "model = LinearRegression()\n",
    "model.fit(x, y)\n",
    "\n",
    "theta_hat = np.c_[model.intercept_, model.coef_[0]]\n",
    "print('theta_hat = {}'.format(theta_hat.ravel()))\n",
    "# estimated response\n",
    "y_est = model.predict(x)\n",
    "\n",
    "# plot\n",
    "plt.plot(x, g(x), label='true fun', color='red');              # the real function\n",
    "plt.scatter(x, y, label='noisy data', color='black');          # data affected by noise\n",
    "plt.plot(x, y_est, label='est fun (sklearn)', color='green');  # estimate linear function\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "pv12xP-dj7n7"
   },
   "source": [
    "During the lectures we have seen that, given a model $f(x; {\\boldsymbol \\theta})$ it is possible minimize the training error iteratively using gradient descent:\n",
    "\n",
    "\n",
    "$${\\boldsymbol \\theta}^{i+1} \\gets {\\boldsymbol \\theta}^i - \\varepsilon_L \\frac{\\partial V_n({\\boldsymbol \\theta})}{\\partial {\\boldsymbol \\theta}} \\bigg \\rvert_{{\\boldsymbol \\theta} = {\\boldsymbol \\theta}^i}$$\n",
    "\n",
    "Consider the mean squared error, \n",
    "$$V_n({\\boldsymbol \\theta}) = {1 \\over n}\\sum_{i=1}^n\\left(y_i - f(x_i; \\boldsymbol \\theta)\\right)^2$$\n",
    "then, in our case:\n",
    "$$\\frac{\\partial V_n({\\boldsymbol \\theta})}{\\partial {\\boldsymbol \\theta}} =\n",
    "\\left[ \n",
    "\\begin{array}{c}\n",
    "\\frac{\\partial V_n({\\boldsymbol \\theta})}{\\partial {\\theta_0}} \\\\\n",
    "\\frac{\\partial V_n({\\boldsymbol \\theta})}{\\partial {\\theta_1}}\n",
    "\\end{array}\n",
    "\\right] =\n",
    "\\left[ \n",
    "\\begin{array}{c}\n",
    "-{2 \\over n}\\sum_{i=1}^n\\left(y_i - f(x_i; \\boldsymbol \\theta)\\right) \\\\\n",
    "-{2 \\over n}\\sum_{i=1}^n\\left(y_i - f(x_i; \\boldsymbol \\theta)\\right)x_i\n",
    "\\end{array}\n",
    "\\right] = -{2 \\over n}X^T(Y - X\\boldsymbol \\theta)$$\n",
    "\n",
    "Let's do it in numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "nZumYWJcoaX5"
   },
   "outputs": [],
   "source": [
    "X = np.hstack((np.ones((x.shape[0], 1)), x)) # add a column of ones\n",
    "y = y.reshape(-1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "n3qwMDohRx6M"
   },
   "outputs": [],
   "source": [
    "def V(X, theta, Y):\n",
    "  return np.mean((Y - np.dot(X, theta))**2)\n",
    "\n",
    "theta = np.array([[-10.], [-5.]]) # initial value theta, it is random you can change it\n",
    "eps = .3  # step size\n",
    "steps = 100 # number of GD steps\n",
    "history = [theta.ravel().copy()]\n",
    "errs = [V(X,theta,y)]\n",
    "\n",
    "for _ in range(steps):  # (Note: underscore `_` is used as name for useless variables) \n",
    "  grad = - 2 * np.dot(X.T, (y - np.dot(X, theta))) / X.shape[0]\n",
    "  theta = theta - eps * grad\n",
    "  # log theta and loss\n",
    "  history.append(theta.ravel().copy())\n",
    "  errs.append(V(X,theta,y))\n",
    "\n",
    "\n",
    "history = np.array(history)\n",
    "print('theta_hat:', theta.ravel())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JlGLqgJMhHd7"
   },
   "source": [
    "As you can see, the iterative procedure converge to the same values of the closed form solution.\n",
    "\n",
    "Let's visualize the descent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "DsPYR5pRTPyc"
   },
   "outputs": [],
   "source": [
    "# code for plotting\n",
    "\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "r = np.abs(history).max()\n",
    "x_range = np.linspace(-r, r, 100)\n",
    "y_range = np.linspace(-r, r, 100)\n",
    "\n",
    "theta_0, theta_1 = np.meshgrid(x_range, y_range)\n",
    "zs = np.array([V(X, t.reshape(-1,1), y) \n",
    "               for t in np.c_[np.ravel(theta_0), np.ravel(theta_1)]])\n",
    "\n",
    "zs = zs.reshape(theta_0.shape)\n",
    "\n",
    "fig = plt.figure(figsize=(12,5))\n",
    "ax1 = fig.add_subplot(122)\n",
    "ax1.set_xlabel(r'$\\theta_0$')\n",
    "ax1.set_ylabel(r'$\\theta_1$')\n",
    "c = ax1.contour(theta_0, theta_1, zs, levels=25,  cmap='viridis')\n",
    "ax1.plot(history[:,0], history[:,1], '-x',color='red')\n",
    "\n",
    "ax2 = fig.add_subplot(121, projection='3d')\n",
    "ax2.plot_trisurf(theta_0.ravel(), theta_1.ravel(), zs.ravel(), cmap='viridis')\n",
    "ax2.set_xlabel(r'$\\theta_0$')\n",
    "ax2.set_ylabel(r'$\\theta_1$')\n",
    "ax2.set_zlabel(r'$V(\\theta)$')\n",
    "ax2.plot(history[:,0], history[:,1], errs, '-x', color='red', alpha=1.)\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OvAmK69ukSBH"
   },
   "source": [
    "You can try to play with the step size. What happens if you increase/decrease it?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JYBhDMHDzJYy"
   },
   "source": [
    "# 01.H) More on linear regression\n",
    "## i. Confidence intervals for the parameters\n",
    "\n",
    "Assume that $X^\\top X$ is invertible, then\n",
    "\n",
    "$$\n",
    "\\hat \\theta = X^+Y \\sim N\\big(\\theta, \\sigma_\\eta^2 (X^\\top X)^{-1}\\big)\n",
    "$$\n",
    "\n",
    "$$\n",
    "E[\\hat \\theta] = E[X^+Y] = X^+E[Y] = X^+ X\\theta = (X^\\top X)^{-1}X^\\top X \\theta = \\theta \n",
    "$$\n",
    "\n",
    "$$\n",
    "Var[\\hat \\theta] = Var[X^+Y] = X^+Var[Y] (X^+)^\\top = \\sigma_\\eta^2 (X^\\top X)^{-1} X^\\top X (X^\\top X)^{-1} = \\sigma_\\eta^2 (X^\\top X)^{-1} \n",
    "$$\n",
    "\n",
    "A rule of thumb is the following\n",
    "\n",
    "* Extract the diagonal from $\\sigma_\\eta^2 (X^\\top X)^{-1}$, which gives you an idea of the variance of each component of $\\theta$.\n",
    "* For each component $\\theta_i$, check if the interval $(\\theta_i - 2\\sigma_i, \\theta_i + 2\\sigma_i)$ contains the zero; if that is the case, we are not very confident that the $\\theta_i\\ne 0$, thus that $x_i$ is relevant in the model."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "01_linear_regression.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}