{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "972a028c",
   "metadata": {},
   "source": [
    "<div style=\"width: 100%; overflow: hidden;\">\n",
    "    <div style=\"width: 150px; float: left;\"> <img src=\"data/D4Sci_logo_ball.png\" alt=\"Data For Science, Inc\" align=\"left\" border=\"0\"> </div>\n",
    "    <div style=\"float: left; margin-left: 10px;\"> <h1>Machine Learning From the Ground Up</h1>\n",
    "<h1>Post 13: Linear models from least squares to quantiles</h1>\n",
    "        <p>Bruno Gon\u00e7alves<br/>\n",
    "        <a href=\"http://www.data4sci.com/\">www.data4sci.com</a><br/>\n",
    "            @bgoncalves, @data4sci</p></div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3020bae0",
   "metadata": {},
   "source": [
    "**Post 13** in the 48-part Data For Science series.\n",
    "\n",
    "**Learning goals**\n",
    "- Fit and compare least-squares, ridge, lasso, elastic net, and quantile regression on the same housing problem.\n",
    "- Build linear classifiers and connect them to generalized linear models and link functions.\n",
    "- Replace squared loss with robust and structured losses for outliers, quantiles, and multiple targets.\n",
    "\n",
    "**Background posts**: Posts 2 to 4, 11, and 12.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2903fad7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:06:52.776701Z",
     "iopub.status.busy": "2026-08-31T11:06:52.776606Z",
     "iopub.status.idle": "2026-08-31T11:06:53.442184Z",
     "shell.execute_reply": "2026-08-31T11:06:53.441517Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Setup complete. Python 3.12.3\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import random\n",
    "import shutil\n",
    "import sys\n",
    "from pathlib import Path\n",
    "import requests\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import watermark\n",
    "from tqdm import tqdm\n",
    "from sklearn.datasets import fetch_california_housing\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler, PolynomialFeatures\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet, LogisticRegression, Perceptron, HuberRegressor, QuantileRegressor, RANSACRegressor\n",
    "from sklearn.isotonic import IsotonicRegression\n",
    "from sklearn.metrics import mean_squared_error, mean_absolute_error, accuracy_score, hamming_loss, label_ranking_loss\n",
    "from sklearn.multioutput import MultiOutputRegressor, MultiOutputClassifier\n",
    "\n",
    "# Seed every random stream so the notebook reproduces exactly.\n",
    "random.seed(42)\n",
    "np.random.seed(42)\n",
    "os.environ['PYTHONHASHSEED'] = '42'\n",
    "print('Setup complete. Python', sys.version.split()[0])\n",
    "\n",
    "# Needed by the template watermark and style cells\n",
    "%load_ext watermark\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed894855",
   "metadata": {},
   "source": [
    "We start by print out the versions of the libraries we're using for future reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b9909df3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:06:53.444059Z",
     "iopub.status.busy": "2026-08-31T11:06:53.443924Z",
     "iopub.status.idle": "2026-08-31T11:06:53.477492Z",
     "shell.execute_reply": "2026-08-31T11:06:53.476774Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python implementation: CPython\n",
      "Python version       : 3.12.3\n",
      "IPython version      : 9.15.0\n",
      "\n",
      "Compiler    : GCC 13.3.0\n",
      "OS          : Linux\n",
      "Release     : 6.17.0-1026-nvidia\n",
      "Machine     : aarch64\n",
      "Processor   : aarch64\n",
      "CPU cores   : 20\n",
      "Architecture: 64bit\n",
      "\n",
      "matplotlib: 3.11.1\n",
      "numpy     : 2.5.1\n",
      "pandas    : 2.3.3\n",
      "requests  : 2.34.2\n",
      "sklearn   : 1.9.0\n",
      "tqdm      : 4.69.0\n",
      "watermark : 2.6.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%watermark -n -v -m -iv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b0642dd",
   "metadata": {},
   "source": [
    "Load default figure style"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c6e34009",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:06:53.479500Z",
     "iopub.status.busy": "2026-08-31T11:06:53.479430Z",
     "iopub.status.idle": "2026-08-31T11:06:53.483258Z",
     "shell.execute_reply": "2026-08-31T11:06:53.482071Z"
    }
   },
   "outputs": [],
   "source": [
    "plt.style.use('d4sci.mplstyle')\n",
    "colors = plt.rcParams['axes.prop_cycle'].by_key()['color']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48ab6e22",
   "metadata": {},
   "source": [
    "Detect the available hardware acceleration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "faa4596c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:06:53.484914Z",
     "iopub.status.busy": "2026-08-31T11:06:53.484842Z",
     "iopub.status.idle": "2026-08-31T11:06:54.240323Z",
     "shell.execute_reply": "2026-08-31T11:06:54.239277Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using device: cuda\n"
     ]
    }
   ],
   "source": [
    "# Prefer CUDA, then Apple MPS, then the CPU. Cells that train models read the\n",
    "# device variable instead of naming hardware. torch is optional on purpose.\n",
    "try:\n",
    "    import torch\n",
    "    if torch.cuda.is_available():\n",
    "        device = torch.device(\"cuda\")\n",
    "    elif torch.backends.mps.is_available():\n",
    "        device = torch.device(\"mps\")\n",
    "    else:\n",
    "        device = torch.device(\"cpu\")\n",
    "except ImportError:\n",
    "    device = \"cpu\"\n",
    "print(f\"Using device: {device}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "556e009f",
   "metadata": {},
   "source": [
    "### Dataset card\n",
    "\n",
    "- **Source**: California Housing, fetched through scikit-learn from the StatLib repository.\n",
    "- **Size**: 20,640 rows, 8 predictive features, 1 target (`MedHouseVal` in $100k).\n",
    "- **License**: Public research data from StatLib, redistributed by scikit-learn.\n",
    "- **Why this dataset**: Continuous target with a heavy right tail and spatial structure. It exposes the differences between squared, Huber, and quantile losses.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e674ea17",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:06:54.242782Z",
     "iopub.status.busy": "2026-08-31T11:06:54.242608Z",
     "iopub.status.idle": "2026-08-31T11:06:54.254963Z",
     "shell.execute_reply": "2026-08-31T11:06:54.253056Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape: (20640, 9)\n",
      "Columns: ['MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup', 'Latitude', 'Longitude', 'MedHouseVal']\n"
     ]
    }
   ],
   "source": [
    "# Download to a local cache, load, and validate.\n",
    "data_dir = Path('data')\n",
    "data_dir.mkdir(exist_ok=True)\n",
    "\n",
    "def load_california_housing():\n",
    "    try:\n",
    "        data = fetch_california_housing(data_home=str(data_dir), download_if_missing=True)\n",
    "        frame = pd.DataFrame(data.data, columns=data.feature_names)\n",
    "        frame['MedHouseVal'] = data.target\n",
    "        assert frame.shape[1] == 9\n",
    "        return frame\n",
    "    except Exception:\n",
    "        # A broken cache must not survive; delete and download once more.\n",
    "        shutil.rmtree(data_dir, ignore_errors=True)\n",
    "        data_dir.mkdir(exist_ok=True)\n",
    "        data = fetch_california_housing(data_home=str(data_dir), download_if_missing=True)\n",
    "        frame = pd.DataFrame(data.data, columns=data.feature_names)\n",
    "        frame['MedHouseVal'] = data.target\n",
    "        return frame\n",
    "\n",
    "df = load_california_housing()\n",
    "required_cols = ['MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup', 'Latitude', 'Longitude', 'MedHouseVal']\n",
    "assert set(required_cols).issubset(df.columns)\n",
    "print('Shape:', df.shape)\n",
    "print('Columns:', list(df.columns))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "828c0a66",
   "metadata": {},
   "source": [
    "### Exploratory data analysis\n",
    "\n",
    "We inspect the raw data before any modeling. Transformations can hide missing values, duplicates, and long tails.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fd30ff71",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:06:54.257918Z",
     "iopub.status.busy": "2026-08-31T11:06:54.257810Z",
     "iopub.status.idle": "2026-08-31T11:06:55.186050Z",
     "shell.execute_reply": "2026-08-31T11:06:55.185307Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rows: 20640 Columns: 9\n",
      "Memory (MB): 1.49\n",
      "dtypes:\n",
      "MedInc         float64\n",
      "HouseAge       float64\n",
      "AveRooms       float64\n",
      "AveBedrms      float64\n",
      "Population     float64\n",
      "AveOccup       float64\n",
      "Latitude       float64\n",
      "Longitude      float64\n",
      "MedHouseVal    float64\n",
      "dtype: object\n",
      "\n",
      "First five raw rows before any transformation:\n",
      "   MedInc  HouseAge  AveRooms  AveBedrms  Population  AveOccup  Latitude  Longitude  MedHouseVal\n",
      "0  8.3252      41.0  6.984127   1.023810       322.0  2.555556     37.88    -122.23        4.526\n",
      "1  8.3014      21.0  6.238137   0.971880      2401.0  2.109842     37.86    -122.22        3.585\n",
      "2  7.2574      52.0  8.288136   1.073446       496.0  2.802260     37.85    -122.24        3.521\n",
      "3  5.6431      52.0  5.817352   1.073059       558.0  2.547945     37.85    -122.25        3.413\n",
      "4  3.8462      52.0  6.281853   1.081081       565.0  2.181467     37.85    -122.25        3.422\n",
      "\n",
      "Missing values per column:\n",
      "MedInc         0\n",
      "HouseAge       0\n",
      "AveRooms       0\n",
      "AveBedrms      0\n",
      "Population     0\n",
      "AveOccup       0\n",
      "Latitude       0\n",
      "Longitude      0\n",
      "MedHouseVal    0\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 118628 characters]",
      "text/plain": [
       "<Figure size 1800x1200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean target (regression base rate): 2.069\n",
      "Median target: 1.797\n",
      "\n",
      "Duplicate rows: 0\n",
      "Rows after dropping duplicates: 20640\n",
      "Leakage candidates: Latitude and Longitude encode neighborhoods and are kept as legitimate features.\n",
      "\n",
      "Summary statistics:\n",
      "               count         mean          std         min         25%          50%          75%           max\n",
      "MedInc       20640.0     3.870671     1.899822    0.499900    2.563400     3.534800     4.743250     15.000100\n",
      "HouseAge     20640.0    28.639486    12.585558    1.000000   18.000000    29.000000    37.000000     52.000000\n",
      "AveRooms     20640.0     5.429000     2.474173    0.846154    4.440716     5.229129     6.052381    141.909091\n",
      "AveBedrms    20640.0     1.096675     0.473911    0.333333    1.006079     1.048780     1.099526     34.066667\n",
      "Population   20640.0  1425.476744  1132.462122    3.000000  787.000000  1166.000000  1725.000000  35682.000000\n",
      "AveOccup     20640.0     3.070655    10.386050    0.692308    2.429741     2.818116     3.282261   1243.333333\n",
      "Latitude     20640.0    35.631861     2.135952   32.540000   33.930000    34.260000    37.710000     41.950000\n",
      "Longitude    20640.0  -119.569704     2.003532 -124.350000 -121.800000  -118.490000  -118.010000   -114.310000\n",
      "MedHouseVal  20640.0     2.068558     1.153956    0.149990    1.196000     1.797000     2.647250      5.000010\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 402424 characters]",
      "text/plain": [
       "<Figure size 3600x3000 with 9 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 481332 characters]",
      "text/plain": [
       "<Figure size 2700x2400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Outlier counts by field (1.5 IQR fences):\n",
      "MedInc          681\n",
      "HouseAge          0\n",
      "AveRooms        511\n",
      "AveBedrms      1424\n",
      "Population     1196\n",
      "AveOccup        711\n",
      "Latitude          0\n",
      "Longitude         0\n",
      "MedHouseVal    1071\n"
     ]
    }
   ],
   "source": [
    "# Shape, memory, and schema first.\n",
    "print('Rows:', df.shape[0], 'Columns:', df.shape[1])\n",
    "print('Memory (MB):', round(df.memory_usage(deep=True).sum() / 1e6, 2))\n",
    "print('dtypes:')\n",
    "print(df.dtypes)\n",
    "\n",
    "print('\\nFirst five raw rows before any transformation:')\n",
    "print(df.head().to_string())\n",
    "\n",
    "missing = df.isna().sum()\n",
    "print('\\nMissing values per column:')\n",
    "print(missing.to_string())\n",
    "assert missing.sum() == 0\n",
    "\n",
    "# Target balance and base rate.\n",
    "plt.figure(figsize=(6, 4))\n",
    "plt.hist(df['MedHouseVal'], bins=50, edgecolor='white')\n",
    "plt.xlabel('Median house value ($100k)')\n",
    "plt.ylabel('Number of districts')\n",
    "plt.title('Target distribution before modeling')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "print('Mean target (regression base rate):', round(df['MedHouseVal'].mean(), 3))\n",
    "print('Median target:', round(df['MedHouseVal'].median(), 3))\n",
    "\n",
    "# Duplicates and leakage candidates.\n",
    "duplicate_count = df.duplicated().sum()\n",
    "print('\\nDuplicate rows:', duplicate_count)\n",
    "df = df.drop_duplicates().reset_index(drop=True)\n",
    "print('Rows after dropping duplicates:', len(df))\n",
    "print('Leakage candidates: Latitude and Longitude encode neighborhoods and are kept as legitimate features.')\n",
    "\n",
    "# Summary statistics for the numeric fields.\n",
    "print('\\nSummary statistics:')\n",
    "print(df.describe().T.to_string())\n",
    "\n",
    "# Histograms for all numeric columns.\n",
    "fig, axes = plt.subplots(3, 3, figsize=(12, 10))\n",
    "for ax, col in zip(axes.ravel(), df.columns):\n",
    "    ax.hist(df[col], bins=50, color='steelblue', edgecolor='white')\n",
    "    ax.set_title(col)\n",
    "    ax.set_xlabel('Value')\n",
    "    ax.set_ylabel('Count')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# No categorical columns exist, so the category bar chart step is skipped.\n",
    "\n",
    "# Correlation heatmap.\n",
    "corr = df.corr()\n",
    "fig, ax = plt.subplots(figsize=(9, 8))\n",
    "mat = ax.imshow(corr, cmap='coolwarm', vmin=-1, vmax=1)\n",
    "ax.set_xticks(range(len(corr)))\n",
    "ax.set_xticklabels(corr.columns, rotation=90)\n",
    "ax.set_yticks(range(len(corr)))\n",
    "ax.set_yticklabels(corr.columns)\n",
    "for i in range(len(corr)):\n",
    "    for j in range(len(corr)):\n",
    "        ax.text(j, i, f'{corr.iloc[i, j]:.2f}', ha='center', va='center', fontsize=8)\n",
    "plt.colorbar(mat, ax=ax, fraction=0.046)\n",
    "plt.title('Pearson correlation of California Housing features')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# Outlier fences using the 1.5 IQR rule.\n",
    "def count_outliers(series):\n",
    "    q1 = series.quantile(0.25)\n",
    "    q3 = series.quantile(0.75)\n",
    "    iqr = q3 - q1\n",
    "    return int(((series < q1 - 1.5 * iqr) | (series > q3 + 1.5 * iqr)).sum())\n",
    "\n",
    "outlier_series = pd.Series({col: count_outliers(df[col]) for col in df.columns}, name='outliers')\n",
    "print('\\nOutlier counts by field (1.5 IQR fences):')\n",
    "print(outlier_series.to_string())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac9e6b87",
   "metadata": {},
   "source": [
    "**EDA findings**\n",
    "\n",
    "- `MedHouseVal` has a long right tail and a hard cap near 5.0, so squared loss will be pulled by top-coded districts.\n",
    "- `MedInc` has the strongest linear correlation with the target.\n",
    "- No missing values exist; no imputation is required.\n",
    "- `Latitude` and `Longitude` are spatial leakage candidates and are kept as legitimate features.\n",
    "- `AveOccup`, `Population`, and `AveRooms` contain extreme outliers; scaling before L1 and L2 models is mandatory.\n",
    "- No categorical columns exist, so one-hot encoding is skipped.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7e32a1a",
   "metadata": {},
   "source": [
    "# Linear regression family\n",
    "\n",
    "We start with a naive mean baseline, then solve the same least-squares problem two ways. **Ordinary Least Squares (OLS)** minimizes squared error and has a closed-form solution through the **normal equations**. **Ridge regression** adds an **L2 regularization** penalty, **Lasso regression** adds an **L1 regularization and sparsity** penalty, and **Elastic Net** mixes both. A **basis expansion** adds polynomial features so a linear model can fit curvature.\n",
    "\n",
    "Ridge objective: `min ||y - Xw||^2 + alpha * ||w||^2`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c2c225d4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:06:55.188478Z",
     "iopub.status.busy": "2026-08-31T11:06:55.188421Z",
     "iopub.status.idle": "2026-08-31T11:06:55.319466Z",
     "shell.execute_reply": "2026-08-31T11:06:55.319247Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train rows: 16512 Test rows: 4128\n",
      "Naive mean baseline RMSE: 1.1449\n",
      "Normal equations RMSE: 0.7456\n",
      "Sklearn OLS RMSE: 0.7456\n",
      "Normal equations and sklearn OLS agree: True\n",
      "Ridge RMSE: 0.7456\n",
      "Lasso RMSE: 0.7404\n",
      "Lasso zero coefficients: 1\n",
      "Elastic Net RMSE: 0.7416\n",
      "Polynomial Ridge RMSE: 0.6877\n",
      "baseline   RMSE=1.1449\n",
      "OLS        RMSE=0.7456\n",
      "Ridge      RMSE=0.7456\n",
      "Lasso      RMSE=0.7404\n",
      "ElasticNet RMSE=0.7416\n",
      "PolyRidge  RMSE=0.6877\n",
      "\n",
      "Scaled coefficients:\n",
      "               OLS   Ridge   Lasso  ElasticNet\n",
      "MedInc      0.8544  0.8543  0.8010      0.8240\n",
      "HouseAge    0.1225  0.1226  0.1271      0.1300\n",
      "AveRooms   -0.2944 -0.2942 -0.1628     -0.2161\n",
      "AveBedrms   0.3393  0.3390  0.2062      0.2563\n",
      "Population -0.0023 -0.0023 -0.0000     -0.0000\n",
      "AveOccup   -0.0408 -0.0408 -0.0306     -0.0359\n",
      "Latitude   -0.8969 -0.8962 -0.7901     -0.7884\n",
      "Longitude  -0.8698 -0.8691 -0.7557     -0.7570\n"
     ]
    }
   ],
   "source": [
    "def rmse(y_true, y_pred):\n",
    "    return float(np.sqrt(mean_squared_error(y_true, y_pred)))\n",
    "\n",
    "X = df.drop(columns=['MedHouseVal'])\n",
    "y = df['MedHouseVal']\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n",
    "assert not set(X_train.index).intersection(set(X_test.index))\n",
    "print('Train rows:', X_train.shape[0], 'Test rows:', X_test.shape[0])\n",
    "\n",
    "scaler = StandardScaler()\n",
    "scaler.fit(X_train)\n",
    "X_train_s = scaler.transform(X_train)\n",
    "X_test_s = scaler.transform(X_test)\n",
    "\n",
    "baseline_pred = np.full_like(y_test, y_train.mean())\n",
    "print('Naive mean baseline RMSE:', round(rmse(y_test, baseline_pred), 4))\n",
    "\n",
    "# Normal equations closed-form solution.\n",
    "X_b = np.hstack([np.ones((X_train_s.shape[0], 1)), X_train_s])\n",
    "theta, _, _, _ = np.linalg.lstsq(X_b, y_train, rcond=None)\n",
    "y_pred_ne = np.hstack([np.ones((X_test_s.shape[0], 1)), X_test_s]) @ theta\n",
    "print('Normal equations RMSE:', round(rmse(y_test, y_pred_ne), 4))\n",
    "\n",
    "# Ordinary Least Squares through sklearn.\n",
    "ols = LinearRegression()\n",
    "ols.fit(X_train_s, y_train)\n",
    "y_pred_ols = ols.predict(X_test_s)\n",
    "print('Sklearn OLS RMSE:', round(rmse(y_test, y_pred_ols), 4))\n",
    "print('Normal equations and sklearn OLS agree:', np.allclose(ols.coef_, theta[1:], atol=1e-6))\n",
    "assert rmse(y_test, y_pred_ols) < rmse(y_test, baseline_pred)\n",
    "\n",
    "# Ridge.\n",
    "ridge = Ridge(alpha=1.0)\n",
    "ridge.fit(X_train_s, y_train)\n",
    "y_pred_ridge = ridge.predict(X_test_s)\n",
    "print('Ridge RMSE:', round(rmse(y_test, y_pred_ridge), 4))\n",
    "\n",
    "# Lasso.\n",
    "lasso = Lasso(alpha=0.01, max_iter=100000)\n",
    "lasso.fit(X_train_s, y_train)\n",
    "y_pred_lasso = lasso.predict(X_test_s)\n",
    "print('Lasso RMSE:', round(rmse(y_test, y_pred_lasso), 4))\n",
    "print('Lasso zero coefficients:', int(np.sum(np.abs(lasso.coef_) < 1e-10)))\n",
    "\n",
    "# Elastic Net.\n",
    "elastic_net = ElasticNet(alpha=0.01, l1_ratio=0.5, max_iter=100000)\n",
    "elastic_net.fit(X_train_s, y_train)\n",
    "y_pred_en = elastic_net.predict(X_test_s)\n",
    "print('Elastic Net RMSE:', round(rmse(y_test, y_pred_en), 4))\n",
    "\n",
    "# Basis expansion.\n",
    "poly_pipe = make_pipeline(PolynomialFeatures(degree=2, include_bias=False), StandardScaler(), Ridge(alpha=1.0))\n",
    "poly_pipe.fit(X_train, y_train)\n",
    "y_pred_poly = poly_pipe.predict(X_test)\n",
    "print('Polynomial Ridge RMSE:', round(rmse(y_test, y_pred_poly), 4))\n",
    "\n",
    "model_rmse = {\n",
    "    'baseline': rmse(y_test, baseline_pred),\n",
    "    'OLS': rmse(y_test, y_pred_ols),\n",
    "    'Ridge': rmse(y_test, y_pred_ridge),\n",
    "    'Lasso': rmse(y_test, y_pred_lasso),\n",
    "    'ElasticNet': rmse(y_test, y_pred_en),\n",
    "    'PolyRidge': rmse(y_test, y_pred_poly)\n",
    "}\n",
    "for name, value in model_rmse.items():\n",
    "    print(f'{name:10s} RMSE={value:.4f}')\n",
    "\n",
    "print('\\nScaled coefficients:')\n",
    "coef_table = pd.DataFrame({\n",
    "    'OLS': ols.coef_,\n",
    "    'Ridge': ridge.coef_,\n",
    "    'Lasso': lasso.coef_,\n",
    "    'ElasticNet': elastic_net.coef_\n",
    "}, index=X.columns)\n",
    "print(coef_table.round(4).to_string())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67309b48",
   "metadata": {},
   "source": [
    "All linear models beat the mean baseline. Ridge and Lasso shrink the OLS coefficients, and Lasso drives some coefficients exactly to zero. The polynomial basis expansion lowers RMSE further, confirming the curvature seen in the EDA histograms.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d6853e4",
   "metadata": {},
   "source": [
    "# Linear classifiers\n",
    "\n",
    "**Logistic regression** is a **Generalized Linear Model (GLM)**. A **link function** maps the linear predictor to the target distribution. For binary outcomes the logit link gives `P(y=1) = 1 / (1 + exp(-x w))`. **Softmax regression** extends the same idea to K classes. The **perceptron convergence proof** bounds the number of mistakes on linearly separable data by a function of the margin.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "27f9f1f2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:06:55.324994Z",
     "iopub.status.busy": "2026-08-31T11:06:55.324908Z",
     "iopub.status.idle": "2026-08-31T11:06:55.462119Z",
     "shell.execute_reply": "2026-08-31T11:06:55.461333Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary base rate: 0.503\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Logistic regression accuracy: 0.8263\n",
      "Softmax regression accuracy: 0.6059\n",
      "Perceptron updates until convergence: 2\n",
      "Perceptron separates all points: 1\n",
      "Sklearn Perceptron accuracy on separable data: 1.0\n"
     ]
    }
   ],
   "source": [
    "# Binary target: above the training median.\n",
    "bin_median = np.median(y_train)\n",
    "y_bin_train = (y_train > bin_median).astype(int)\n",
    "y_bin_test = (y_test > bin_median).astype(int)\n",
    "bin_base = max(y_bin_test.mean(), 1 - y_bin_test.mean())\n",
    "print('Binary base rate:', round(bin_base, 3))\n",
    "\n",
    "log_reg = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000, random_state=42))\n",
    "log_reg.fit(X_train, y_bin_train)\n",
    "y_bin_pred = log_reg.predict(X_test)\n",
    "log_acc = accuracy_score(y_bin_test, y_bin_pred)\n",
    "print('Logistic regression accuracy:', round(log_acc, 4))\n",
    "assert log_acc > bin_base\n",
    "\n",
    "# Multiclass target: quartiles of the training target.\n",
    "train_bins = np.quantile(y_train, [0.25, 0.5, 0.75])\n",
    "y_multi_train = np.clip(np.digitize(y_train, train_bins), 0, 3)\n",
    "y_multi_test = np.clip(np.digitize(y_test, train_bins), 0, 3)\n",
    "\n",
    "soft_reg = make_pipeline(StandardScaler(), LogisticRegression(solver='lbfgs', max_iter=1000, random_state=42))\n",
    "soft_reg.fit(X_train, y_multi_train)\n",
    "y_multi_pred = soft_reg.predict(X_test)\n",
    "soft_acc = accuracy_score(y_multi_test, y_multi_pred)\n",
    "print('Softmax regression accuracy:', round(soft_acc, 4))\n",
    "\n",
    "# Perceptron convergence on a separable synthetic set.\n",
    "rng = np.random.default_rng(42)\n",
    "X_per = np.vstack([rng.uniform(0, 1, (30, 2)), rng.uniform(-1, 0, (30, 2))])\n",
    "y_per = np.array([1] * 30 + [-1] * 30)\n",
    "X_per_b = np.hstack([np.ones((len(X_per), 1)), X_per])\n",
    "w = np.zeros(3)\n",
    "updates = 0\n",
    "for _ in range(100):\n",
    "    mistakes = 0\n",
    "    for xi, yi in zip(X_per_b, y_per):\n",
    "        if yi * (w @ xi) <= 0:\n",
    "            w = w + yi * xi\n",
    "            updates += 1\n",
    "            mistakes += 1\n",
    "    if mistakes == 0:\n",
    "        break\n",
    "print('Perceptron updates until convergence:', updates)\n",
    "print('Perceptron separates all points:', int(np.all(y_per * (X_per_b @ w) > 0)))\n",
    "\n",
    "# Compare with the sklearn Perceptron implementation.\n",
    "perceptron_sk = Perceptron(random_state=42)\n",
    "perceptron_sk.fit(X_per, y_per)\n",
    "print('Sklearn Perceptron accuracy on separable data:', round(accuracy_score(y_per, perceptron_sk.predict(X_per)), 4))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bbe83ef",
   "metadata": {},
   "source": [
    "Logistic regression beats the binary base rate. Softmax regression on quartiles shows the multinomial extension. The perceptron converges on the synthetic separable set in a finite number of updates, matching the perceptron convergence proof.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e477338",
   "metadata": {},
   "source": [
    "# Losses for outliers and special targets\n",
    "\n",
    "Squared loss overweights outliers. **Huber loss** is a robust **M-estimator**. **Quantile regression** minimizes **pinball loss** and estimates conditional quantiles. **RANSAC** fits a model to inliers. **Isotonic regression** learns a monotone curve. **Ordinal regression** respects ordered categories, and **multi-label learning** and **multi-output learning** predict several targets at once.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5200596a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:06:55.463970Z",
     "iopub.status.busy": "2026-08-31T11:06:55.463836Z",
     "iopub.status.idle": "2026-08-31T11:07:15.016120Z",
     "shell.execute_reply": "2026-08-31T11:07:15.015330Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Huber RMSE: 0.7584\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "  0%|          | 0/3 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      " 33%|\u2588\u2588\u2588\u258e      | 1/3 [00:06<00:12,  6.37s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quantile 0.1: RMSE=1.0471, pinball=0.0903, MAE=0.7592\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      " 67%|\u2588\u2588\u2588\u2588\u2588\u2588\u258b   | 2/3 [00:12<00:06,  6.51s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quantile 0.5: RMSE=0.7636, pinball=0.2560, MAE=0.5120\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 3/3 [00:19<00:00,  6.46s/it]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 3/3 [00:19<00:00,  6.46s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quantile 0.9: RMSE=1.2778, pinball=0.1540, MAE=1.0315\n",
      "Mean baseline pinball at q=0.1: 0.4598\n",
      "Mean baseline pinball at q=0.5: 0.4530\n",
      "Mean baseline pinball at q=0.9: 0.4463\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Robust and quantile losses on the same regression split.\n",
    "huber = make_pipeline(StandardScaler(), HuberRegressor(epsilon=1.35, max_iter=1000))\n",
    "huber.fit(X_train, y_train)\n",
    "y_pred_huber = huber.predict(X_test)\n",
    "print('Huber RMSE:', round(rmse(y_test, y_pred_huber), 4))\n",
    "\n",
    "def pinball_loss(y_true, y_pred, quantile):\n",
    "    residual = y_true - y_pred\n",
    "    return float(np.mean(np.where(residual >= 0, quantile * residual, (quantile - 1) * residual)))\n",
    "\n",
    "qr_models = {}\n",
    "for q in tqdm([0.1, 0.5, 0.9]):\n",
    "    qr = make_pipeline(StandardScaler(), QuantileRegressor(quantile=q, alpha=0.0001, solver='highs'))\n",
    "    qr.fit(X_train, y_train)\n",
    "    y_pred_q = qr.predict(X_test)\n",
    "    qr_models[q] = y_pred_q\n",
    "    print(f'Quantile {q}: RMSE={rmse(y_test, y_pred_q):.4f}, pinball={pinball_loss(y_test, y_pred_q, q):.4f}, MAE={mean_absolute_error(y_test, y_pred_q):.4f}')\n",
    "\n",
    "mean_pred = np.full_like(y_test, y_train.mean())\n",
    "for q in [0.1, 0.5, 0.9]:\n",
    "    print(f'Mean baseline pinball at q={q}: {pinball_loss(y_test, mean_pred, q):.4f}')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7a1f04a9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:07:15.017800Z",
     "iopub.status.busy": "2026-08-31T11:07:15.017744Z",
     "iopub.status.idle": "2026-08-31T11:07:15.398212Z",
     "shell.execute_reply": "2026-08-31T11:07:15.397190Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RANSAC RMSE: 0.8169\n",
      "RANSAC inliers: 12664\n",
      "Isotonic (MedInc only) RMSE: 0.8336\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 756744 characters]",
      "text/plain": [
       "<Figure size 1800x1200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# RANSAC fits only inliers and resists label outliers.\n",
    "ransac = RANSACRegressor(estimator=LinearRegression(), random_state=42, min_samples=0.5)\n",
    "ransac.fit(X_train, y_train)\n",
    "y_pred_ransac = ransac.predict(X_test)\n",
    "print('RANSAC RMSE:', round(rmse(y_test, y_pred_ransac), 4))\n",
    "print('RANSAC inliers:', int(np.sum(ransac.inlier_mask_)))\n",
    "\n",
    "# Isotonic regression learns a monotone curve for MedInc.\n",
    "iso_model = IsotonicRegression(out_of_bounds='clip')\n",
    "iso_model.fit(X_train['MedInc'].values, y_train)\n",
    "y_pred_iso = iso_model.predict(X_test['MedInc'].values)\n",
    "print('Isotonic (MedInc only) RMSE:', round(rmse(y_test, y_pred_iso), 4))\n",
    "\n",
    "plt.figure(figsize=(6, 4))\n",
    "sort_idx = np.argsort(X_train['MedInc'].values)\n",
    "plt.plot(X_train['MedInc'].values[sort_idx], iso_model.predict(X_train['MedInc'].values[sort_idx]), linewidth=2, label='Isotonic fit')\n",
    "plt.scatter(X_train['MedInc'].values, y_train, s=4, alpha=0.2, label='Train points')\n",
    "plt.xlabel('MedInc ($10k)')\n",
    "plt.ylabel('MedHouseVal ($100k)')\n",
    "plt.title('Isotonic regression on MedInc')\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b05b1c9d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:07:15.401526Z",
     "iopub.status.busy": "2026-08-31T11:07:15.401431Z",
     "iopub.status.idle": "2026-08-31T11:07:15.514590Z",
     "shell.execute_reply": "2026-08-31T11:07:15.513771Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ordinal regression accuracy: 0.5102\n",
      "Ordinal mean absolute error: 0.5264\n",
      "Multi-label hamming loss: 0.0587\n",
      "Multi-label ranking loss: 0.0033\n",
      "Multi-output MedHouseVal RMSE: 0.7462\n",
      "Multi-output AveOccup RMSE: 2.0053\n"
     ]
    }
   ],
   "source": [
    "# Ordinal regression treats quartile labels as ordered numbers.\n",
    "ordinal = make_pipeline(StandardScaler(), LinearRegression())\n",
    "ordinal.fit(X_train, y_multi_train)\n",
    "y_ordinal_pred = np.clip(np.round(ordinal.predict(X_test)), 0, 3).astype(int)\n",
    "print('Ordinal regression accuracy:', round(accuracy_score(y_multi_test, y_ordinal_pred), 4))\n",
    "print('Ordinal mean absolute error:', round(mean_absolute_error(y_multi_test, y_ordinal_pred), 4))\n",
    "\n",
    "# Multi-label learning predicts three binary properties at once.\n",
    "ml_labels = pd.DataFrame({\n",
    "    'high_value': (df['MedHouseVal'] > df['MedHouseVal'].median()).astype(int),\n",
    "    'high_income': (df['MedInc'] > df['MedInc'].median()).astype(int),\n",
    "    'old_housing': (df['HouseAge'] > df['HouseAge'].median()).astype(int),\n",
    "})\n",
    "X_ml = df.drop(columns=['MedHouseVal'])\n",
    "X_ml_train, X_ml_test, Y_ml_train, Y_ml_test = train_test_split(X_ml, ml_labels, test_size=0.2, random_state=42)\n",
    "ml_clf = make_pipeline(StandardScaler(), MultiOutputClassifier(LogisticRegression(max_iter=1000, random_state=42)))\n",
    "ml_clf.fit(X_ml_train, Y_ml_train)\n",
    "Y_ml_pred = ml_clf.predict(X_ml_test)\n",
    "print('Multi-label hamming loss:', round(hamming_loss(Y_ml_test, Y_ml_pred), 4))\n",
    "ml_score = np.column_stack([p[:, 1] for p in ml_clf.predict_proba(X_ml_test)])\n",
    "print('Multi-label ranking loss:', round(label_ranking_loss(Y_ml_test.to_numpy(), ml_score), 4))\n",
    "\n",
    "# Multi-output learning predicts two continuous targets.\n",
    "X_mo = df.drop(columns=['MedHouseVal', 'AveOccup'])\n",
    "Y_mo = df[['MedHouseVal', 'AveOccup']]\n",
    "X_mo_train, X_mo_test, Y_mo_train, Y_mo_test = train_test_split(X_mo, Y_mo, test_size=0.2, random_state=42)\n",
    "multi_out = make_pipeline(StandardScaler(), MultiOutputRegressor(LinearRegression()))\n",
    "multi_out.fit(X_mo_train, Y_mo_train)\n",
    "Y_mo_pred = multi_out.predict(X_mo_test)\n",
    "for i, col in enumerate(Y_mo.columns):\n",
    "    print(f'Multi-output {col} RMSE:', round(rmse(Y_mo_test.iloc[:, i], Y_mo_pred[:, i]), 4))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cb6d49c",
   "metadata": {},
   "source": [
    "Huber and quantile regression reduce the influence of top-coded prices. RANSAC identifies an inlier set and fits clean. Isotonic regression improves the single-feature `MedInc` fit with a monotone curve. The ordinal, multi-label, and multi-output variants show how the same linear core handles structured targets.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2d75570",
   "metadata": {},
   "source": [
    "## Recap\n",
    "\n",
    "EDA found a right-tailed target, no missing values, a strong `MedInc` signal, and spatial leakage candidates. The linear family confirmed that `MedInc` dominates the coefficients and that L1 sparsity zeroes redundant features. The polynomial basis expansion captured the curvature suggested by the EDA histograms. Robust, quantile, and structured-target losses then gave more stable models for the outliers that the first plots exposed.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3287d09d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-31T11:07:15.516804Z",
     "iopub.status.busy": "2026-08-31T11:07:15.516709Z",
     "iopub.status.idle": "2026-08-31T11:07:15.519470Z",
     "shell.execute_reply": "2026-08-31T11:07:15.519004Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Headline median quantile regression RMSE: 0.7636\n",
      "Quantile regression is a robust alternative to squared-loss regression.\n"
     ]
    }
   ],
   "source": [
    "# Headline rerun: median quantile regression RMSE.\n",
    "headline_rmse = rmse(y_test, qr_models[0.5])\n",
    "print(f'Headline median quantile regression RMSE: {headline_rmse:.4f}')\n",
    "print('Quantile regression is a robust alternative to squared-loss regression.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fca856f8",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "1. Tune `alpha` in Ridge with a grid from 0.01 to 10 and plot the test RMSE curve.\n",
    "2. Run `cross_val_score` for Huber and quantile 0.5 on the same folds and compare the mean RMSE.\n",
    "3. Add `AveBedrms` to the multi-output target and report the new RMSE.\n",
    "4. Verify that pinball loss at q=0.5 equals half the mean absolute error.\n",
    "5. Explain why the perceptron is guaranteed to converge on a separable set and why it can loop forever on a nonseparable set.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eddd395",
   "metadata": {},
   "source": [
    "## Next post\n",
    "\n",
    "Next post: nonlinear models. We introduce kernels and the kernel trick. Further reading: Hastie, Tibshirani, and Friedman, The Elements of Statistical Learning, Chapter 3; Huber, 1964; Koenker and Bassett, 1978.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3852800b",
   "metadata": {},
   "source": [
    "<center>\n",
    "     <img src=\"data/D4Sci_logo_full.png\" alt=\"Data For Science, Inc\" align=\"center\" border=\"0\" width=300px> \n",
    "</center>"
   ]
  }
 ],
 "metadata": {
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}