Exercise-04.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# California housing dataset with linear and polynomial regression \n",
    "\n",
    "In this notebook, we'll use [linear regression](https://scikit-learn.org/stable/modules/linear_model.html#ordinary-least-squares), [regularized linear regression](https://scikit-learn.org/stable/modules/linear_model.html#ridge-regression), and [polynomial regression](https://scikit-learn.org/stable/modules/linear_model.html#polynomial-regression-extending-linear-models-with-basis-functions) to estimate median house values on Californian housing districts.\n",
    "\n",
    "First, the needed imports. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "from sklearn import datasets, __version__\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.preprocessing import StandardScaler, PolynomialFeatures\n",
    "from sklearn.pipeline import Pipeline\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "\n",
    "Then we load the California housing data. First time we need to download the data, which can take a while."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chd = datasets.fetch_california_housing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first convert the data into a Pandas DataFrame to inspect some basic statistics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame(data=chd.data, columns=chd.feature_names)\n",
    "df['Target'] = pd.Series(chd.target, index=df.index)\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the data consists of 20640 housing districts, each characterized with 8 attributes: *MedInc, HouseAge, AveRooms, AveBedrms, Population, AveOccup, Latitude, Longitude*. We also defined  a target value (median house value) for each housing district.\n",
    " \n",
    "Let's then plot all attributes against the target value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,10))\n",
    "for i in range(8):\n",
    "    plt.subplot(4,2,i+1)\n",
    "    plt.scatter(chd.data[:,i], chd.target, s=2, label=chd.feature_names[i])\n",
    "    plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll now split the data into a training and a test set. Let's use 5000 samples as test data.\n",
    "\n",
    "Let's also select a single attribute to start the analysis with, for example *MedInc*. This way we can plot the regression functions against the target value. Later we will use all attributes in the regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_size = 5000\n",
    "single_attribute = 'MedInc'\n",
    "\n",
    "X_train_all, X_test_all, y_train, y_test = train_test_split(\n",
    "    chd.data, chd.target, test_size=test_size, shuffle=True)\n",
    "\n",
    "attribute_index = chd.feature_names.index(single_attribute)\n",
    "X_train_single = X_train_all[:, attribute_index].reshape(-1, 1)\n",
    "X_test_single = X_test_all[:, attribute_index].reshape(-1, 1)\n",
    "     \n",
    "print()\n",
    "print('California housing data: train:',len(X_train_all),'test:',len(X_test_all))\n",
    "print()\n",
    "print('X_train_all:', X_train_all.shape)\n",
    "print('X_train_single:', X_train_single.shape)\n",
    "print('y_train:', y_train.shape)\n",
    "print()\n",
    "print('X_test_all', X_test_all.shape)\n",
    "print('X_test_single', X_test_single.shape)\n",
    "print('y_test', y_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The training data matrix `X_train_all` is a matrix of size (`n_train`, 8), and `X_train_single` contains only the first attribute (*MedInc* by default) of each housing district. The vector `y_train` contains the target value (median house value) for each housing district in the training set.\n",
    "\n",
    "Let's start our analysis with the single attribute. Later, you can set `only_single_attribute = False` to use all eight attributes in the regression. \n",
    "\n",
    "As the final step, let's scale our input data to zero mean and unit variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "only_single_attribute = True\n",
    "\n",
    "if only_single_attribute:\n",
    "    X_train = X_train_single\n",
    "    X_test = X_test_single\n",
    "else:\n",
    "    X_train = X_train_all\n",
    "    X_test = X_test_all\n",
    "\n",
    "scaler = StandardScaler().fit(X_train)\n",
    "X_train = scaler.transform(X_train)\n",
    "X_test = scaler.transform(X_test)\n",
    "print('X_train: shape:', X_train.shape, 'mean:', X_train.mean(axis=0), 'std:', X_train.std(axis=0))\n",
    "print('X_test: shape:', X_test.shape, 'mean:', X_test.mean(axis=0), 'std:', X_test.std(axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression\n",
    "\n",
    "We begin with linear regression:\n",
    "\n",
    "$$J(w) = \\|y - Xw\\|^2_2$$\n",
    "\n",
    "### Learning\n",
    "\n",
    "The parameters of linear regression can be solved in closed form as:\n",
    "\n",
    "$$\\hat{w} = (X^TX)^{-1}X^Ty$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "lin_reg = LinearRegression()\n",
    "lin_reg.fit(X_train, y_train)\n",
    "print('coefficients:', lin_reg.coef_)\n",
    "print('intercept:', lin_reg.intercept_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize the results if we are using only a single attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if X_train.shape[1] == 1:\n",
    "    plt.figure(figsize=(10, 10))\n",
    "    plt.scatter(X_train, y_train, s=5)\n",
    "    reg_x = np.arange(np.min(X_train), np.max(X_train), 0.01).reshape(-1, 1)\n",
    "    plt.scatter(reg_x, lin_reg.predict(reg_x), s=8, label='linear')\n",
    "    plt.legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inference\n",
    "\n",
    "We use *mean squared error* as the performance measure for our regression algorihm: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "predictions = lin_reg.predict(X_test)\n",
    "print(\"Mean squared error: %.3f\"\n",
    "      % mean_squared_error(y_test, predictions))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regularized linear regression: Ridge\n",
    "\n",
    "Ridge regression adds $L_2$ regularization: \n",
    "\n",
    "$$J(w) = \\|y - Xw\\|^2_2 + \\alpha \\|w\\|^2_2$$\n",
    "\n",
    "where $\\alpha \\ge 0$ is the penalty parameter for the weights. You can experiment with different values of $\\alpha$.\n",
    "\n",
    "You can also try out [`Lasso()`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html#sklearn-linear-model-lasso) or [`ElasticNet()`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html#sklearn-linear-model-elasticnet). Note that Elastic net has also a second parameter `l1_ratio`. \n",
    "\n",
    "### Learning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "alpha = 1e3\n",
    "\n",
    "rdg_reg = Ridge(alpha=alpha)\n",
    "rdg_reg.fit(X_train, y_train)\n",
    "print('coefficients:', rdg_reg.coef_)\n",
    "print('intercept:', rdg_reg.intercept_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if X_train.shape[1] == 1:\n",
    "    plt.figure(figsize=(10, 10))\n",
    "    plt.scatter(X_train, y_train, s=5)\n",
    "    reg_x = np.arange(np.min(X_train), np.max(X_train), 0.01).reshape(-1, 1)\n",
    "    plt.scatter(reg_x, lin_reg.predict(reg_x), s=8, label='linear');\n",
    "    plt.scatter(reg_x, rdg_reg.predict(reg_x), s=8, label=r'ridge, $\\alpha=${}'.format(alpha))\n",
    "    plt.legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "predictions = rdg_reg.predict(X_test)\n",
    "print(\"Mean squared error: %.3f\"\n",
    "      % mean_squared_error(y_test, predictions))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Polynomial regression\n",
    "\n",
    "Polynomial regression can be performed by constructing polynomial features, e.g.:\n",
    "\n",
    "$$z=[1,\\,x_1,\\,x_2,\\,x_1x_2,\\,x_1^2,\\,x_2^2]$$\n",
    "\n",
    "and using a linear model with the new features:\n",
    "\n",
    "$$J(w) = \\|z - X'w\\|^2_2$$\n",
    "\n",
    "### Learning\n",
    "\n",
    "To implement polynomial regression, we use scikit-learn's [Pipeline](https://scikit-learn.org/stable/modules/compose.html#pipeline), a tool for building composite estimators. \n",
    "\n",
    "Note that the polynomial features contain all possible combinations, so the number of features grows quickly especially when using many attributes. Also, you can try using regularized linear regression with polynomial features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "degree = 5\n",
    "\n",
    "poly_model = Pipeline([('poly', PolynomialFeatures(degree=degree)),\n",
    "                      ('linear', LinearRegression(fit_intercept=False))])\n",
    "poly_model.fit(X_train, y_train)\n",
    "print('coefficients:', poly_model.steps[1][1].coef_)\n",
    "print('intercept:', poly_model.steps[1][1].intercept_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if X_train.shape[1] == 1:\n",
    "    plt.figure(figsize=(10, 10))\n",
    "    plt.scatter(X_train, y_train, s=5)\n",
    "    reg_x = np.arange(np.min(X_train), np.max(X_train), 0.01).reshape(-1, 1)\n",
    "    plt.scatter(reg_x, lin_reg.predict(reg_x), s=8, label='linear');\n",
    "    plt.scatter(reg_x, poly_model.predict(reg_x), s=8, label='polynomial')\n",
    "    plt.legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "predictions = poly_model.predict(X_test)\n",
    "print(\"Mean squared error: %.3f\"\n",
    "      % mean_squared_error(y_test, predictions))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model tuning\n",
    "\n",
    "Try to reduce the mean squared error of the regression. Experiment with several single attributes and with using all attributes.\n",
    "\n",
    "To further improve the results, it is possible to replace [StandardScaler](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html), that is scaling the input data to zero mean and unit variance, with more advanced preprocessing.\n",
    "See [Preprocessing data](https://scikit-learn.org/stable/modules/preprocessing.html#preprocessing-data) for more information."
   ]
  },
  {
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}