TP2-Stein-and-Akaike-criterion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2e1763ba",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "# Shrinkage estimators and Akaike criterion\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa9d8b13",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from itertools import combinations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e0996cd",
   "metadata": {},
   "source": [
    "## Part 1\n",
    "- Generate a dataset according to a Gaussian sequence model\n",
    "- Take the mean to be a cosine function and plot the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4da5dabd",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def generate_dataset(mu, sigma):\n",
    "    \"\"\"Generate a dataset according to the Gaussian sequence model\n",
    "    y = mu + z where z is a Gaussian variable\n",
    "    with variance sigma^2 I\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    mu : array of size n\n",
    "        Mean\n",
    "    sigma : float\n",
    "        Standard deviation\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Y : array of size (n)\n",
    "        Noisy observation\n",
    "    \"\"\"\n",
    "    raise NotImplementedError"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0079477",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def cosine_dataset():\n",
    "    \"\"\"Take a grid of n=100 points between 0 and 10, take\n",
    "    mu[i] = cos(i) where i is the value of each point\n",
    "    and generate a dataset according to the Gaussian sequence\n",
    "    model with sigma^2 = 0.5\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    mu : array of size n\n",
    "        Mean\n",
    "    sigma : float\n",
    "        Standard deviation\n",
    "    y : array of size (n)\n",
    "        Noisy observation\n",
    "    grid : array of size n\n",
    "        The grid used to generate mu\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "def plot_dataset():\n",
    "    \"\"\"\n",
    "    Plot on the same graph mu in function of i and Y in function\n",
    "    of i (taking sigma^2 = 0.5)\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "plot_dataset()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf4a5a5d",
   "metadata": {},
   "source": [
    "## Part 2\n",
    "- Implement a shrinkage estimator\n",
    "- Compute the expected error for a grid of shrinkage coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c3180a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "mu, sigma, _, _ = cosine_dataset()\n",
    "\n",
    "def shrinkage(y, s):\n",
    "    \"\"\" The shrinkage estimator\n",
    "    Observe y and returns s * y\n",
    "    Parameters\n",
    "    -----------\n",
    "    y: array of size n\n",
    "        Observation\n",
    "    s: float between 0 and 1\n",
    "        Shrinkage coefficient\n",
    "    Return\n",
    "    ------\n",
    "    sy: array of size n\n",
    "        Shrinked observations\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "def expected_error(mu_h, mu, sigma, n_trials):\n",
    "    \"\"\"\n",
    "    Compute the expected error between mu_h and mu\n",
    "    for a number (trials) of different datasets\n",
    "    (with the same mean mu)\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    mu_h: function from an array of size n to an array of size n\n",
    "        Estimator of mu\n",
    "    mu: array of size n\n",
    "        True value of mu\n",
    "    n_trials: int\n",
    "        Number of trials\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    E: float\n",
    "        Expected mean squared error of mu_h\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "\n",
    "def plot_shrinkage():\n",
    "    \"\"\"\n",
    "    Take s in a grid between 0 and 1 with 100 points.\n",
    "    Plot in function of s the mean squared error of\n",
    "    the shrinkage estimator.\n",
    "    Plot a vertical line at the maximum value\n",
    "    s =  1/(1 + n sigma^2 / (|mu|^2)) where the norm is the l2 norm\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "plot_shrinkage()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "958f7c7f",
   "metadata": {},
   "source": [
    "## Part 3\n",
    "- Compute Akaike criterion for linear estimators (y_hat = Ay)\n",
    "- Compute A for linear regression\n",
    "- Using the Akaike criterion for linear regression, select the best 3 subset of features in the diabete dataset and check that the linear regression using only the 3 best features outperform the linear regression with all features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f7f05e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def akaike(A, y, sigma):\n",
    "    \"\"\"\n",
    "    Compute the akaike estimator of the error of mu_hat = Ay\n",
    "    Recall that it is given by\n",
    "    ||mu_hat - y ||^2 + 2 sigma^2 tr(A) - sigma^2 n\n",
    "\n",
    "    Parameters\n",
    "    -----------\n",
    "    A: array of shape (n, n)\n",
    "    y: array of shape n\n",
    "    sigma: float\n",
    "\n",
    "    Return\n",
    "    -------\n",
    "    akaike: float\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "def linear_regression_coefficients(X, y):\n",
    "    \"\"\"\n",
    "    Return argmin_b ||Xb - y ||^2\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    X : array of shape (n, p)\n",
    "        Design matrix\n",
    "\n",
    "    y: array of shape n\n",
    "        Observation\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    b: array of shape p\n",
    "        regression coefficients\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "def linear_regression_operator(X):\n",
    "    \"\"\"\n",
    "    Return A such that  argmin_b ||Xb - y ||^2 = Ay\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    X : array of shape (n, p)\n",
    "        Design matrix\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    A: array of shape n, n\n",
    "        linear regression operator\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "def akaike_regression(X, y, sigma):\n",
    "    \"\"\" Akaike criterion for linear regression\n",
    "    with design matrix X\n",
    "    Find A such that Ay = b_hat = argmin_b || Xb - y ||^2\n",
    "    Return the Akaike estimate of the error of y_hat = X b_hat\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    X : array\n",
    "        Design matrix\n",
    "    y: array of shape n\n",
    "        Observation\n",
    "    sigma: float\n",
    "        standard deviation\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    akaike : float\n",
    "        The akaike estimate of the error of y_hat\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "def select_best_3_features(X, y, sigma):\n",
    "    \"\"\" Loop across features of X to find the best according to Akaike\n",
    "    Consider all possible choice of 3 features across\n",
    "    the number of features of X.\n",
    "    For each choice of 3 features:\n",
    "    - Build the matrix Xp which is of size (n, 3) and only contains\n",
    "    the chosen features.\n",
    "    - Compute the akaike regression criterion of Xp\n",
    "    Return Xp which yields the lowest akaike regression criterion\n",
    "\n",
    "    Parameters\n",
    "    ---------\n",
    "    X : array (n_samples, n_features)\n",
    "        Design matrix\n",
    "    y: array of shape n\n",
    "        Observation\n",
    "    sigma: float\n",
    "        standard deviation\n",
    "\n",
    "    Return\n",
    "    -------\n",
    "    S: boolean array (n_features,)\n",
    "        S[i] = True if feature i is selected\n",
    "        S has at exactly 3 non-False coefficients\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "def test_diabetes():\n",
    "    from sklearn.datasets import load_diabetes\n",
    "    from sklearn.model_selection import train_test_split\n",
    "    data = load_diabetes()\n",
    "    X, y = data.data, data.target\n",
    "    print(X.shape)\n",
    "    X = np.column_stack([X, np.ones(len(X)).reshape(-1, 1)])\n",
    "    print(X.shape)\n",
    "    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.9, random_state=0)\n",
    "    S = select_best_3_features(X_train, y_train, sigma)\n",
    "    b = linear_regression_coefficients(X_train, y_train)\n",
    "    bS = linear_regression_coefficients(X_train[:, S], y_train)\n",
    "    error = np.mean((X_test.dot(b) - y_test)**2)\n",
    "    errorS = np.mean((X_test[:, S].dot(bS) - y_test)**2)\n",
    "    print(\"Error without feature selection\", error)\n",
    "    print(\"Error with feature selection\", errorS)\n",
    "    assert errorS < error\n",
    "\n",
    "\n",
    "test_diabetes()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0033a2ef",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "# Bonus\n",
    "\n",
    "Take Ridge Regression min_b ||Xb - y ||^2 + lambda ||b||^2\n",
    "And find the best lambda using akaike criterion\n",
    "\n",
    "Check that it gives better results in practice (diabete dataset)"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "auto:percent,ipynb"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}