4 Multiple Linear Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thinking in Tensors, writing in PyTorch\n",
    "\n",
    "A hands-on course by [Piotr Migdał](https://p.migdal.pl) (2019).\n",
    "This notebook prepared by [Weronika Ormaniec](https://github.com/werkaaa).\n",
    "\n",
    "## Notebook 4: Multiple Linear Regression\n",
    "\n",
    "\n",
    "<a href=\"https://colab.research.google.com/github/stared/thinking-in-tensors-writing-in-pytorch/blob/master/4%20Multiple%20Linear%20Regression.ipynb\">\n",
    "    <img src=\"https://colab.research.google.com/assets/colab-badge.svg\"/>\n",
    "</a>\n",
    "\n",
    "Simple linear regression is a useful tool when it comes to predicting an output given single predictor input. However, in practice we often come across problems which are described by more than one predictor. In this case we use Multiple Linear Regression.\n",
    "\n",
    "Instead of fitting several linear equations for each predictor, we will create one equation that will take the form:\n",
    "$$ Y = \\alpha_0 + \\alpha_1 \\cdot X_1 + \\alpha_2\\cdot X_2 + ... + \\alpha_n\\cdot X_n$$\n",
    "where $X_i$ is one of the predictors, $\\alpha_1$ is a coefficient, we want to get to know and $n$ is the number of predictors.\n",
    "\n",
    "The learning process in Multiple Linear Regression is the same as the one in Simple Linear Regression. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.datasets import load_boston\n",
    "import matplotlib.pyplot as plt\n",
    "from livelossplot import PlotLosses\n",
    "\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "from torchvision import transforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we will analyze The Boston Housing Dataset. It contains information about 506 houses in Boston. There are 13 features of the houses, which have grate or little impact on the price of the house. Using PyTorch we will implement a model that will predict the price of the house and then we will try to answer the question, which parameters have the biggest impact on the price of the houses\n",
    "\n",
    "We will take the dataset from scikit learn datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston = load_boston()\n",
    "boston_df = pd.DataFrame(boston.data, columns=boston.feature_names)\n",
    "boston_df[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameters description:\n",
    "    \n",
    "* CRIM: Per capita crime rate by town\n",
    "* ZN: Proportion of residential land zoned for lots over 25,000 sq. ft\n",
    "* INDUS: Proportion of non-retail business acres per town\n",
    "* CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)\n",
    "* NOX: Nitric oxide concentration (parts per 10 million)\n",
    "* RM: Average number of rooms per dwelling\n",
    "* AGE: Proportion of owner-occupied units built prior to 1940\n",
    "* DIS: Weighted distances to five Boston employment centers\n",
    "* RAD: Index of accessibility to radial highways\n",
    "* TAX: Full-value property tax rate per $10,000\n",
    "* PTRATIO: Pupil-teacher ratio by town\n",
    "* B: 1000(Bk — 0.63)², where Bk is the proportion of [people of African American descent] by town\n",
    "* LSTAT: Percentage of lower status of the population\n",
    "\n",
    "Target: Median value of owner-occupied homes in $1000s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all, let's check which parameters have the most linear correlation with the price of the houses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4,figsize=(15,7))\n",
    "ax1.scatter(boston_data_frame['RM'], boston.target)\n",
    "ax1.set_xlabel('RM')\n",
    "ax1.set_ylabel('Target')\n",
    "\n",
    "ax2.scatter(boston_data_frame['CRIM'], boston.target)\n",
    "ax2.set_xlabel('CRIM')\n",
    "ax2.set_ylabel('Target')\n",
    "\n",
    "ax3.scatter(boston_data_frame['PTRATIO'], boston.target)\n",
    "ax3.set_xlabel('PTRATIO')\n",
    "ax3.set_ylabel('Target')\n",
    "\n",
    "ax4.scatter(boston_data_frame['LSTAT'], boston.target)\n",
    "ax4.set_xlabel('LSTAT')\n",
    "ax4.set_ylabel('Target')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can say that the correlation between RM (number of rooms per dwelling) and the price may be linear. The same as correlation between LSTAT (percentage of lower status of the population) and the price. What about CRIM (per capita crime rate by town) and PTRATIO (pupil-teacher ratio by town)? Those relationships are clearly not linear. Let's check how linear model will put up with it!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the data, we can see that some predictors have different orders of magnitude. That can be an obstacle during model training. That is why, we will normalize the data, so they will be in range $[-1,1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = torch.tensor(boston.data, dtype=torch.float32)\n",
    "Y = torch.tensor(boston.target, dtype=torch.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Normalize(data):\n",
    "    data_mean = torch.mean(data, dim=0)\n",
    "    data_max = torch.max(data, dim=0)[0]\n",
    "    data_min = torch.min(data, dim=0)[0]\n",
    "    data = (data-data_mean)/(data_max-data_min)\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_normalized = Normalize(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston_df = pd.DataFrame(np.array(X_normalized), columns=boston.feature_names)\n",
    "boston_df[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time we will divide the data into training and test sets because we will be able to measure how well the model is doing in general, on the examples it has not seen during training process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = X_normalized[:400]\n",
    "Y_train = Y[:400]\n",
    "X_test = X_normalized[401:]\n",
    "Y_test = Y[401:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Linear(nn.Module):\n",
    "    def __init__(self):\n",
    "        super().__init__()\n",
    "        self.linear = nn.Linear(in_features=13, out_features=1)\n",
    "        \n",
    "    def forward(self, x):\n",
    "        return self.linear(x).squeeze(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linear_model = Linear()\n",
    "print(linear_model.linear.weight)\n",
    "print(linear_model.linear.bias)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_predict_train = linear_model(X_train)\n",
    "rmse_train = torch.sqrt(F.mse_loss(Y_train, y_predict_train))\n",
    "\n",
    "y_predict_test = linear_model(X_test)\n",
    "rmse_test = torch.sqrt(F.mse_loss(Y_test, y_predict_test))\n",
    "\n",
    "print(\"The PyTorch model performance:\")\n",
    "print('RMSE_train is {}'.format(rmse_train))\n",
    "print('RMSE_test is {}'.format(rmse_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#optim = torch.optim.SGD(linear_model.parameters(), lr=0.1)\n",
    "optim = torch.optim.Adam(linear_model.parameters(), lr=1.)\n",
    "loss_function = F.mse_loss\n",
    "loss = loss_function(linear_model(X), Y)\n",
    "print(loss)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train(X, Y, model, loss_function, optim, num_epochs):\n",
    "    loss_history = []\n",
    "    liveloss = PlotLosses()\n",
    "\n",
    "\n",
    "    for epoch in range(num_epochs):\n",
    "        \n",
    "        epoch_loss = 0.0\n",
    "        \n",
    "        Y_pred = model(X)\n",
    "        loss = loss_function(Y_pred, Y)\n",
    "        \n",
    "        loss.backward()\n",
    "        optim.step()\n",
    "        optim.zero_grad()\n",
    "        \n",
    "\n",
    "        epoch_loss = loss.data.item()\n",
    "        \n",
    "        avg_loss = epoch_loss / len(X)\n",
    "\n",
    "        liveloss.update({\n",
    "            'loss': avg_loss,\n",
    "        })\n",
    "        liveloss.draw()\n",
    "\n",
    "train(X_train, Y_train, linear_model, loss_function, optim, num_epochs=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_predict_train = linear_model(X_train)\n",
    "rmse_train = torch.sqrt(F.mse_loss(Y_train, y_predict_train))\n",
    "\n",
    "y_predict_test = linear_model(X_test)\n",
    "rmse_test = torch.sqrt(F.mse_loss(Y_test, y_predict_test))\n",
    "\n",
    "print(\"The PyTorch model performance:\")\n",
    "print('RMSE train is {:.3f}'.format(rmse_train))\n",
    "print('RMSE test is {:.3f}'.format(rmse_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A we can see, our model fits the data better after training. \n",
    "\n",
    "We can now compare it with scikit learn linear regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lin_model = LinearRegression()\n",
    "lin_model.fit(X_train.numpy(), Y_train.numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_ptrain = lin_model.predict(X_train)\n",
    "rmse_train = np.sqrt(mean_squared_error(Y_train, y_ptrain))\n",
    "\n",
    "y_ptest = lin_model.predict(X_test)\n",
    "rmse_test = np.sqrt(mean_squared_error(Y_test, y_ptest))\n",
    "\n",
    "print(\"The model performance for training set\")\n",
    "print('RMSE train is {:.3f}'.format(rmse_train))\n",
    "print('RMSE test is {:.3f}'.format(rmse_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our model is not perfect but it has learned some intuition about the data and is able to make predictions even on the data it has not seen during learning process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the coefficients of both models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_groups = 13\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(15, 7))\n",
    "\n",
    "index = np.arange(n_groups)\n",
    "bar_width = 0.35\n",
    "opacity = 0.8\n",
    "\n",
    "rects1 = plt.bar(index, linear_model.linear.weight.detach().squeeze(), bar_width,\n",
    "alpha=opacity,\n",
    "color='b',\n",
    "label='our model')\n",
    "\n",
    "rects2 = plt.bar(index + bar_width, lin_model.coef_, bar_width,\n",
    "alpha=opacity,\n",
    "color='g',\n",
    "label='scikit learn')\n",
    "\n",
    "plt.xlabel('Variable')\n",
    "plt.ylabel('Value')\n",
    "plt.title('Coefficients')\n",
    "plt.xticks(index + bar_width, ('CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',\n",
    "       'TAX', 'PTRATIO', 'B', 'LSTAT'))\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## To do\n",
    "* Proper dataset slicing"
   ]
  }
 ],
 "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.6.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}