conformal_regression_ex.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conformal inference: regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.datasets import fetch_california_housing\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler, LabelEncoder\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score\n",
    "from sklearn.linear_model import LinearRegression, QuantileRegressor\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "from mapie.metrics import regression_coverage_score\n",
    "from mapie.regression import MapieQuantileRegressor, MapieRegressor\n",
    "\n",
    "from utils import plot_1d_data, x_sinx, get_1d_data_with_constant_noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating the aleatoric uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate data with DGP: homoscedastic 1D non-linear process f(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 10**3\n",
    "min_x, max_x, noise = -5, 5, 0.5\n",
    "\n",
    "X_train, y_train, X_test, y_test, y_mesh = get_1d_data_with_constant_noise(x_sinx, min_x, max_x, n_samples, noise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X_train, y_train, color=\"C0\")             # measurements\n",
    "_ = plt.plot(X_test, y_mesh, color=\"C1\")              # latent regression function\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"f(x)\")\n",
    "plt.title(r'Regression function $f(x) = x\\cdot \\sin(x)$ with measurements')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "polyn_model = Pipeline(\n",
    "    [\n",
    "        (\"poly\", PolynomialFeatures(degree=10)),\n",
    "        (\"model\", RandomForestRegressor(n_estimators=50, random_state=42))\n",
    "        # (\"linear\", QuantileRegressor(\n",
    "        #         solver=\"highs\",\n",
    "        #         alpha=0,\n",
    "        # ))\n",
    "    ]\n",
    ")\n",
    "\n",
    "mapie = MapieRegressor(polyn_model, verbose=1)\n",
    "\n",
    "mapie.fit(X_train, y_train)\n",
    "\n",
    "y_pred, y_set = mapie.predict(X_test, alpha=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_1d_data(\n",
    "    X_train.ravel(),\n",
    "    y_train.ravel(),\n",
    "    X_test.ravel(),\n",
    "    y_mesh.ravel(),\n",
    "    np.full((X_test.shape[0]), 1.96*noise).ravel(),     # true predictive interval\n",
    "    y_pred.ravel(),\n",
    "    y_set[:, 0, 0].ravel(),\n",
    "    y_set[:, 1, 0].ravel(),\n",
    "    title=r'Regression function with prediction intervals'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: California Housing dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "View response variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fetching the California Housing dataset\n",
    "housing = fetch_california_housing()\n",
    "\n",
    "X, y = housing.data, housing.target\n",
    "X.shape, y.shape\n",
    "\n",
    "pd.Series(y).hist(bins=50, figsize=(5, 3));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Partition data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_rest1, y_train, y_rest1 = train_test_split(X, y, train_size=0.8, random_state=42)\n",
    "X_test, X_rest2, y_test, y_rest2 = train_test_split(X_rest1, y_rest1, train_size=0.8, random_state=42)\n",
    "X_calib, X_new, y_calib, y_new = train_test_split(X_rest2, y_rest2, train_size=0.8, random_state=42)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "\n",
    "X_train = scaler.fit_transform(X_train)\n",
    "X_test = scaler.transform(X_test)\n",
    "X_calib = scaler.transform(X_calib)\n",
    "X_new = scaler.transform(X_new)\n",
    "\n",
    "print(X_train.shape)\n",
    "print(X_test.shape)\n",
    "print(X_calib.shape)\n",
    "print(X_new.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = RandomForestRegressor(n_estimators=100, random_state=42)\n",
    "\n",
    "model.fit(X_train, y_train)\n",
    "\n",
    "y_pred = model.predict(X_test)\n",
    "\n",
    "mae = mean_absolute_error(y_test, y_pred)\n",
    "print(f\"MAE: {mae}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mapie.regression import MapieRegressor\n",
    "\n",
    "mapie_reg = MapieRegressor(estimator=model, cv='prefit', n_jobs=-1, verbose=1)\n",
    "\n",
    "mapie_reg.fit(X_calib, y_calib)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "y_pred, y_set = mapie_reg.predict(X_new, alpha=alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0             # first sample\n",
    "interval = y_set[i].flatten()    \n",
    "\n",
    "print(f\"{(1-alpha)*100:.0f}% predictive interval: [{interval[0]:.2f}, {interval[1]:.2f}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mapie.metrics import regression_coverage_score_v2, regression_mean_width_score\n",
    "\n",
    "coverage = regression_coverage_score_v2(y_new, y_set)\n",
    "\n",
    "print(f\"Coverage: {coverage[0]:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = regression_mean_width_score(y_set[0], y_set[1])\n",
    "\n",
    "print(f\"Avg. interval width: {size:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env_py",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}