exercise_04_linear_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Programming Task: Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from sklearn.datasets import load_boston\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": [
    "## Exporting the results to PDF\n",
    "Once you complete the assignments, export the entire notebook as PDF and attach it to your homework solutions. \n",
    "The best way of doing that is\n",
    "1. Run all the cells of the notebook (`Kernel -> Restart & Run All`)\n",
    "2. Export/download the notebook as PDF (`File -> Download as -> PDF via LaTeX (.pdf)`)\n",
    "3. Concatenate your solutions for other tasks with the output of Step 2. On Linux you can simply use `pdfunite`, there are similar tools for other platforms too. You can only upload a single PDF file to Moodle.\n",
    "\n",
    "**Make sure** you are using `nbconvert` **Version 5.5 or later** by running `jupyter nbconvert --version`. Older versions clip lines that exceed page width, which makes your code harder to grade."
   ]
  },
  {
   "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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\anaconda3\\envs\\AI_AIMA\\lib\\site-packages\\sklearn\\utils\\deprecation.py:87: FutureWarning: Function load_boston is deprecated; `load_boston` is deprecated in 1.0 and will be removed in 1.2.\n",
      "\n",
      "    The Boston housing prices dataset has an ethical problem. You can refer to\n",
      "    the documentation of this function for further details.\n",
      "\n",
      "    The scikit-learn maintainers therefore strongly discourage the use of this\n",
      "    dataset unless the purpose of the code is to study and educate about\n",
      "    ethical issues in data science and machine learning.\n",
      "\n",
      "    In this special case, you can fetch the dataset from the original\n",
      "    source::\n",
      "\n",
      "        import pandas as pd\n",
      "        import numpy as np\n",
      "\n",
      "\n",
      "        data_url = \"http://lib.stat.cmu.edu/datasets/boston\"\n",
      "        raw_df = pd.read_csv(data_url, sep=\"\\s+\", skiprows=22, header=None)\n",
      "        data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])\n",
      "        target = raw_df.values[1::2, 2]\n",
      "\n",
      "    Alternative datasets include the California housing dataset (i.e.\n",
      "    :func:`~sklearn.datasets.fetch_california_housing`) and the Ames housing\n",
      "    dataset. You can load the datasets as follows::\n",
      "\n",
      "        from sklearn.datasets import fetch_california_housing\n",
      "        housing = fetch_california_housing()\n",
      "\n",
      "    for the California housing dataset and::\n",
      "\n",
      "        from sklearn.datasets import fetch_openml\n",
      "        housing = fetch_openml(name=\"house_prices\", as_frame=True)\n",
      "\n",
      "    for the Ames housing dataset.\n",
      "    \n",
      "  warnings.warn(msg, category=FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "X , y = load_boston(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.2\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1: Fit standard linear regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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",
    "    w=(np.linalg.inv(np.dot(np.transpose(X),X))).dot(np.transpose(X)).dot(y)\n",
    "    return w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2: Fit ridge regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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",
    "    E=np.eye(X.shape[1])\n",
    "    w=(np.linalg.inv(np.transpose(X).dot(X)+reg_strength*E)).dot(np.transpose(X)).dot(y)\n",
    "    return w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3: Generate predictions for new data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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_pred=np.dot(X,w)\n",
    "    return y_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4: Mean squared error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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=np.size(y_pred)\n",
    "    mse=np.sum((y_pred-y_true)**2)/n\n",
    "    return mse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare the two models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reference implementation produces \n",
    "* MSE for Least squares $\\approx$ **23.96**\n",
    "* MSE for Ridge regression $\\approx$ **21.03**\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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE for Least squares = 23.964571384957637\n",
      "MSE for Ridge regression = 21.03493121591776\n"
     ]
    }
   ],
   "source": [
    "# Load the data\n",
    "np.random.seed(1234)\n",
    "X , y = load_boston(return_X_y=True)\n",
    "X = np.hstack([np.ones([X.shape[0], 1]), X])\n",
    "test_size = 0.2\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 = 1\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",
   "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"
  },
  "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
}