exercise_linear_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Programming Task: Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from sklearn.datasets import fetch_california_housing\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your task"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook provides a code skeleton for performing linear regression. \n",
    "Your task is to complete the functions where required. \n",
    "You are only allowed to use built-in Python functions, as well as any `numpy` functions. No other libraries / imports are allowed.\n",
    "\n",
    "In the beginning of every function there is docstring which specifies the input and and expected output.\n",
    "Write your code in a way that adheres to it.\n",
    "You may only use plain python and anything that we imported for you above such as numpy functions (i.e. no scikit-learn classifiers)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load and preprocess the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this assignment we will work with the Boston Housing Dataset.\n",
    "The data consists of 506 samples. Each sample represents a district in the city of Boston and has 13 features, such as crime rate or taxation level. The regression target is the median house price in the given district (in $1000's).\n",
    "\n",
    "More details can be found here: http://lib.stat.cmu.edu/datasets/boston"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X , y = fetch_california_housing(return_X_y=True)\n",
    "\n",
    "# Add a vector of ones to the data matrix to absorb the bias term\n",
    "# (Recall slide #7 from the lecture)\n",
    "X = np.hstack([np.ones([X.shape[0], 1]), X])\n",
    "# From now on, D refers to the number of features in the AUGMENTED dataset (i.e. including the dummy '1' feature for the absorbed bias term)\n",
    "\n",
    "# Split into train and test\n",
    "test_size = 0.9\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20640, 9)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "(20640,)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape # has 9 feature including the absorbed bias term - NxD Data Matrix\n",
    "y.shape # N-dim vector with targets\n",
    "type(X)\n",
    "type(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1: Fit standard linear regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_least_squares(X, y):\n",
    "    \"\"\"Fit ordinary least squares model to the data.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    X : array, shape [N, D]\n",
    "        (Augmented) feature matrix.\n",
    "    y : array, shape [N]\n",
    "        Regression targets.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    w : array, shape [D]\n",
    "        Optimal regression coefficients (w[0] is the bias term). \n",
    "        \n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE ###\n",
    "    A=np.transpose(X)@X\n",
    "    w_temp=np.linalg.inv(A)@np.transpose(X)\n",
    "    w=w_temp@y\n",
    "    \n",
    "    return w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Test\n",
    "# A=np.transpose(X)@X\n",
    "# B=np.linalg.inv(A)\n",
    "# B@np.transpose(X)@y\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2: Fit ridge regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_ridge(X, y, reg_strength):\n",
    "    \"\"\"Fit ridge regression model to the data.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    X : array, shape [N, D]\n",
    "        (Augmented) feature matrix.\n",
    "    y : array, shape [N]\n",
    "        Regression targets.\n",
    "    reg_strength : float\n",
    "        L2 regularization strength (denoted by lambda in the lecture)\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    w : array, shape [D]\n",
    "        Optimal regression coefficients (w[0] is the bias term).\n",
    "    \n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE ###\n",
    "    \n",
    "    n=X.shape[1]\n",
    "    I=np.identity(n,dtype=float)\n",
    "    A=np.transpose(X)@X\n",
    "    B=np.linalg.inv(A+reg_strength*I)\n",
    "    w=B@np.transpose(X)@y\n",
    "    \n",
    "    return w\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 2., 0., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 2., 0., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 2., 0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 2., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 2., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 2., 0., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 0., 2., 0.],\n",
       "       [0., 0., 0., 0., 0., 0., 0., 0., 2.]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n=X.shape[1]\n",
    "2*np.identity(n,dtype=float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3: Generate predictions for new data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict_linear_model(X, w):\n",
    "    \"\"\"Generate predictions for the given samples.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    X : array, shape [N, D]\n",
    "        (Augmented) feature matrix.\n",
    "    w : array, shape [D]\n",
    "        Regression coefficients.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    y_pred : array, shape [N]\n",
    "        Predicted regression targets for the input data.\n",
    "        \n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE ###\n",
    "    y_hat=X@w\n",
    "    \n",
    "    return y_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4: Mean squared error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_squared_error(y_true, y_pred):\n",
    "    \"\"\"Compute mean squared error between true and predicted regression targets.\n",
    "    \n",
    "    Reference: `https://en.wikipedia.org/wiki/Mean_squared_error`\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y_true : array\n",
    "        True regression targets.\n",
    "    y_pred : array\n",
    "        Predicted regression targets.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    mse : float\n",
    "        Mean squared error.\n",
    "        \n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE ###\n",
    "    n=y.shape[0]\n",
    "    err=y_true-y_pred\n",
    "    mse=np.inner(err,err)/n\n",
    "    \n",
    "    return mse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# y.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare the two models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reference implementation produces \n",
    "* MSE for Least squares $\\approx$ **0.5347**\n",
    "* MSE for Ridge regression $\\approx$ **0.5331**\n",
    "\n",
    "You results might be slightly (i.e. $\\pm 1\\%$) different from the reference soultion due to numerical reasons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE for Least squares = 0.48123921834064254\n",
      "MSE for Ridge regression = 0.47981572810065615\n"
     ]
    }
   ],
   "source": [
    "# Load the data\n",
    "np.random.seed(1234)\n",
    "X , y = fetch_california_housing(return_X_y=True)\n",
    "X = np.hstack([np.ones([X.shape[0], 1]), X])\n",
    "test_size = 0.9\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size)\n",
    "\n",
    "# Ordinary least squares regression\n",
    "w_ls = fit_least_squares(X_train, y_train)\n",
    "y_pred_ls = predict_linear_model(X_test, w_ls)\n",
    "mse_ls = mean_squared_error(y_test, y_pred_ls)\n",
    "print('MSE for Least squares = {0}'.format(mse_ls))\n",
    "\n",
    "# Ridge regression\n",
    "reg_strength = 1e-2\n",
    "w_ridge = fit_ridge(X_train, y_train, reg_strength)\n",
    "y_pred_ridge = predict_linear_model(X_test, w_ridge)\n",
    "mse_ridge = mean_squared_error(y_test, y_pred_ridge)\n",
    "print('MSE for Ridge regression = {0}'.format(mse_ridge))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.9.7"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "48px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}