{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T16:59:50.806061Z",
     "start_time": "2025-09-28T16:59:50.490698Z"
    }
   },
   "outputs": [],
   "source": [
    "# Useful starting lines\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "from test_utils import test"
   ]
  },
  {
   "cell_type": "markdown",
   "execution_count": null,
   "metadata": {},
   "source": [
    "# Cross-Validation and Bias-Variance decomposition\n",
    "## Cross-Validation\n",
    "Implementing 4-fold cross-validation below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T16:59:52.086679Z",
     "start_time": "2025-09-28T16:59:52.074045Z"
    }
   },
   "outputs": [],
   "source": [
    "from helpers import load_data\n",
    "\n",
    "# load dataset\n",
    "x, y = load_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T16:59:54.105196Z",
     "start_time": "2025-09-28T16:59:54.096977Z"
    }
   },
   "outputs": [],
   "source": [
    "def build_k_indices(y, k_fold, seed):\n",
    "    \"\"\"build k indices for k-fold.\n",
    "\n",
    "    Args:\n",
    "        y:      shape=(N,)\n",
    "        k_fold: K in K-fold, i.e. the fold num\n",
    "        seed:   the random seed\n",
    "\n",
    "    Returns:\n",
    "        A 2D array of shape=(k_fold, N/k_fold) that indicates the data indices for each fold\n",
    "\n",
    "    >>> build_k_indices(np.array([1., 2., 3., 4.]), 2, 1)\n",
    "    array([[3, 2],\n",
    "           [0, 1]])\n",
    "    \"\"\"\n",
    "    num_row = y.shape[0]\n",
    "    interval = int(num_row / k_fold)\n",
    "    np.random.seed(seed)\n",
    "    indices = np.random.permutation(num_row)\n",
    "    k_indices = [indices[k * interval : (k + 1) * interval] for k in range(k_fold)]\n",
    "    return np.array(k_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T16:59:55.631894Z",
     "start_time": "2025-09-28T16:59:55.601895Z"
    }
   },
   "outputs": [],
   "source": [
    "test(build_k_indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "execution_count": null,
   "metadata": {},
   "source": [
    "For the following cross_validation( ) function you need to implement, you can help yourselves of the build_poly( ) and ridge_regression( ) functions that you implemented in lab 3. Copy paste the code in the build_polynomial.py and ridge_regression.py files, they should pass the two following tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T16:59:58.309835Z",
     "start_time": "2025-09-28T16:59:58.291747Z"
    }
   },
   "outputs": [],
   "source": [
    "from costs import compute_mse\n",
    "from ridge_regression import ridge_regression\n",
    "from build_polynomial import build_poly\n",
    "\n",
    "test(build_poly)\n",
    "test(ridge_regression)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T17:01:43.184416Z",
     "start_time": "2025-09-28T17:01:43.152215Z"
    }
   },
   "outputs": [],
   "source": [
    "def cross_validation(y, x, k_indices, k, lambda_, degree):\n",
    "    \"\"\"return the loss of ridge regression for a fold corresponding to k_indices\n",
    "\n",
    "    Args:\n",
    "        y:          shape=(N,)\n",
    "        x:          shape=(N,)\n",
    "        k_indices:  2D array returned by build_k_indices()\n",
    "        k:          scalar, the k-th fold (N.B.: not to confused with k_fold which is the fold nums)\n",
    "        lambda_:    scalar, cf. ridge_regression()\n",
    "        degree:     scalar, cf. build_poly()\n",
    "\n",
    "    Returns:\n",
    "        train and test root mean square errors rmse = sqrt(2 mse)\n",
    "\n",
    "    >>> cross_validation(np.array([1.,2.,3.,4.]), np.array([6.,7.,8.,9.]), np.array([[3,2], [0,1]]), 1, 2, 3)\n",
    "    (0.019866645527596753, 0.33555914361295197)\n",
    "    \"\"\"\n",
    "\n",
    "    # ***************************************************\n",
    "    # INSERT YOUR CODE HERE\n",
    "    # get k'th subgroup in test, others in train: TODO\n",
    "    # ***************************************************\n",
    "    raise NotImplementedError\n",
    "    # ***************************************************\n",
    "    # INSERT YOUR CODE HERE\n",
    "    # form data with polynomial degree: TODO\n",
    "    # ***************************************************\n",
    "    raise NotImplementedError\n",
    "    # ***************************************************\n",
    "    # INSERT YOUR CODE HERE\n",
    "    # ridge regression: TODO\n",
    "    # ***************************************************\n",
    "    raise NotImplementedError\n",
    "    # ***************************************************\n",
    "    # INSERT YOUR CODE HERE\n",
    "    # calculate the loss for train and test data: TODO\n",
    "    # ***************************************************\n",
    "    raise NotImplementedError\n",
    "    return loss_tr, loss_te"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T17:01:45.284163Z",
     "start_time": "2025-09-28T17:01:45.252178Z"
    }
   },
   "outputs": [],
   "source": [
    "# can lead to a numerical error if you use an older version than Python 3.9\n",
    "test(cross_validation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T17:02:03.385114Z",
     "start_time": "2025-09-28T17:02:03.184707Z"
    }
   },
   "outputs": [],
   "source": [
    "from plots import cross_validation_visualization\n",
    "\n",
    "\n",
    "def cross_validation_demo(degree, k_fold, lambdas):\n",
    "    \"\"\"cross validation over regularisation parameter lambda.\n",
    "\n",
    "    Args:\n",
    "        degree: integer, degree of the polynomial expansion\n",
    "        k_fold: integer, the number of folds\n",
    "        lambdas: shape = (p, ) where p is the number of values of lambda to test\n",
    "    Returns:\n",
    "        best_lambda : scalar, value of the best lambda\n",
    "        best_rmse : scalar, the associated root mean squared error for the best lambda\n",
    "    \"\"\"\n",
    "\n",
    "    seed = 12\n",
    "    degree = degree\n",
    "    k_fold = k_fold\n",
    "    lambdas = lambdas\n",
    "    # split data in k fold\n",
    "    k_indices = build_k_indices(y, k_fold, seed)\n",
    "    # define lists to store the loss of training data and test data\n",
    "    rmse_tr = []\n",
    "    rmse_te = []\n",
    "    # ***************************************************\n",
    "    # INSERT YOUR CODE HERE\n",
    "    # cross validation over lambdas: TODO\n",
    "    # ***************************************************\n",
    "    raise NotImplementedError\n",
    "\n",
    "    cross_validation_visualization(lambdas, rmse_tr, rmse_te)\n",
    "    print(\n",
    "        \"For polynomial expansion up to degree %.f, the choice of lambda which leads to the best test rmse is %.5f with a test rmse of %.3f\"\n",
    "        % (degree, best_lambda, best_rmse)\n",
    "    )\n",
    "    return best_lambda, best_rmse\n",
    "\n",
    "\n",
    "best_lambda, best_rmse = cross_validation_demo(7, 4, np.logspace(-4, 0, 30))"
   ]
  },
  {
   "cell_type": "markdown",
   "execution_count": null,
   "metadata": {},
   "source": [
    "Your output should look like this for seed = 12, degree = 7 and k_fold = 4:\n",
    "\n",
    "![alt text](cross_validation2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "execution_count": null,
   "metadata": {},
   "source": [
    "You can play around with the number of folds and the degree of your polynomial expansion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T17:02:23.559133Z",
     "start_time": "2025-09-28T17:02:23.343141Z"
    }
   },
   "outputs": [],
   "source": [
    "best_lambda, best_rmse = cross_validation_demo(10, 4, np.logspace(-10, -2, 30))"
   ]
  },
  {
   "cell_type": "markdown",
   "execution_count": null,
   "metadata": {},
   "source": [
    "In the previous task we did a grid search over several values of $\\lambda$ for a fixed degree. We can also perform a grid search amongst $\\lambda$ and degrees simultaneously:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T17:03:59.889173Z",
     "start_time": "2025-09-28T17:03:59.855624Z"
    }
   },
   "outputs": [],
   "source": [
    "def best_degree_selection(degrees, k_fold, lambdas, seed=1):\n",
    "    \"\"\"cross validation over regularisation parameter lambda and degree.\n",
    "\n",
    "    Args:\n",
    "        degrees: shape = (d,), where d is the number of degrees to test\n",
    "        k_fold: integer, the number of folds\n",
    "        lambdas: shape = (p, ) where p is the number of values of lambda to test\n",
    "    Returns:\n",
    "        best_degree : integer, value of the best degree\n",
    "        best_lambda : scalar, value of the best lambda\n",
    "        best_rmse : value of the rmse for the couple (best_degree, best_lambda)\n",
    "\n",
    "    >>> best_degree_selection(np.arange(2,11), 4, np.logspace(-4, 0, 30))\n",
    "    (7, 0.004520353656360241, 0.2895728056656946)\n",
    "    \"\"\"\n",
    "\n",
    "    # split data in k fold\n",
    "    k_indices = build_k_indices(y, k_fold, seed)\n",
    "\n",
    "    # ***************************************************\n",
    "    # INSERT YOUR CODE HERE\n",
    "    # cross validation over degrees and lambdas: TODO\n",
    "    # ***************************************************\n",
    "    raise NotImplementedError\n",
    "\n",
    "    return int(best_degree), float(best_lambda), float(best_rmse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T17:04:01.609582Z",
     "start_time": "2025-09-28T17:04:01.477625Z"
    }
   },
   "outputs": [],
   "source": [
    "# can lead to a numerical error if you use an older version than Python 3.9\n",
    "test(best_degree_selection)\n",
    "\n",
    "best_degree, best_lambda, best_rmse = best_degree_selection(\n",
    "    np.arange(2, 11), 4, np.logspace(-4, 0, 30)\n",
    ")\n",
    "print(\n",
    "    \"The best rmse of %.3f is obtained for a degree of %.f and a lambda of %.5f.\"\n",
    "    % (best_rmse, best_degree, best_lambda)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "execution_count": null,
   "metadata": {},
   "source": [
    "## Bias-Variance Decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T17:04:05.923229Z",
     "start_time": "2025-09-28T17:04:05.911253Z"
    }
   },
   "outputs": [],
   "source": [
    "# true function we want to learn\n",
    "def f_star(x):\n",
    "    return x**3 - x**2 + 0.5\n",
    "\n",
    "\n",
    "# plotting function for f_star\n",
    "def plot_fstar(ax):\n",
    "    xvals = np.arange(-1, 1, 0.01)\n",
    "    ax.plot(xvals, f_star(xvals), linestyle=\"--\", color=\"k\", label=\"f_star\")\n",
    "    ax.set_ylim(-2, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T17:04:07.211306Z",
     "start_time": "2025-09-28T17:04:07.199147Z"
    }
   },
   "outputs": [],
   "source": [
    "# helper plot function\n",
    "def plot_poly(x, y, weights, degree, ax, alpha=0.3):\n",
    "    xvals = np.arange(-1, 1, 0.01)\n",
    "    tx = build_poly(xvals, degree)\n",
    "    f = tx.dot(weights)\n",
    "    ax.plot(xvals, f, color=\"orange\", alpha=alpha)\n",
    "    ax.scatter(x, y, color=\"b\", alpha=alpha, s=10)\n",
    "    ax.set_xlabel(\"x\")\n",
    "    ax.set_ylabel(\"y\")\n",
    "    ax.set_title(\"Polynomial degree \" + str(degree))\n",
    "    ax.set_ylim(-1, 2)\n",
    "\n",
    "\n",
    "# helper plot function\n",
    "def plot_f(weights, degree, ax, label=None):\n",
    "    xvals = np.arange(-1, 1, 0.01)\n",
    "    tx = build_poly(xvals, degree)\n",
    "    f = tx.dot(weights)\n",
    "    ax.plot(xvals, f, color=\"black\", alpha=1, label=label)\n",
    "    ax.set_ylim(-1, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "execution_count": null,
   "metadata": {},
   "source": [
    "Implement the following function: for 15 random datapoints, it finds the optimal fit (using the least square formula, with no regularisation λ) for a polynomial expansion of degree 1, 3 and 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T17:04:10.186222Z",
     "start_time": "2025-09-28T17:04:10.042062Z"
    }
   },
   "outputs": [],
   "source": [
    "from least_squares import least_squares\n",
    "\n",
    "\n",
    "def bias_variance_one_seed(sigma, degrees, seed):\n",
    "    \"\"\"One run of the optimal fit for 15 random points and different polynomial expansion degrees.\n",
    "\n",
    "    Args:\n",
    "        sigma: scalar, noise variance\n",
    "        degrees: shape = (3,), 3 different degrees to consider\n",
    "        seed: integer, random see\n",
    "    Returns:\n",
    "    \"\"\"\n",
    "\n",
    "    # we will generate 15 random datapoints from the [-1, 1] uniform distribuion\n",
    "    num_data = 15\n",
    "    np.random.seed(seed)  # set random seed for reproducibility\n",
    "    xs = np.random.uniform(-1, 1, num_data)\n",
    "    # the outputs will be f_star(x) + some random gaussian noise of variance sigma**2\n",
    "    ys = f_star(xs) + sigma * np.random.randn(num_data)\n",
    "\n",
    "    fig, axs = plt.subplots(1, len(degrees), figsize=(20, 5))\n",
    "    for index_degree, degree in enumerate(degrees):\n",
    "        # ***************************************************\n",
    "        # INSERT YOUR CODE HERE\n",
    "        # ***************************************************\n",
    "\n",
    "        plot_fstar(axs[index_degree])\n",
    "        axs[index_degree].legend()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "bias_variance_one_seed(0.1, [1, 3, 6], seed=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "execution_count": null,
   "metadata": {},
   "source": [
    "Your output should ressemble (for seed = 2) to this: \n",
    "![alt text](bias_variance_one_run.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "execution_count": null,
   "metadata": {},
   "source": [
    "Now to illustrate the bias variance tradeoff we will repeat many times the previous experiment but using a different random seed each time. We also plot (in plain black) the mean of all the orange functions obtained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-28T17:04:37.227604Z",
     "start_time": "2025-09-28T17:04:33.986346Z"
    }
   },
   "outputs": [],
   "source": [
    "def bias_variance_demo(sigma, degrees):\n",
    "    \"\"\"Illustration of the bias-variance tradeoff.\n",
    "\n",
    "    Args:\n",
    "        sigma: scalar, noise variance\n",
    "        degrees: shape = (3,), 3 different degrees to consider\n",
    "    Returns:\n",
    "    \"\"\"\n",
    "    # define parameters\n",
    "    seeds = range(400)  # number of runs\n",
    "    num_data = 15\n",
    "\n",
    "    fig, axs = plt.subplots(1, len(degrees), figsize=(20, 5))\n",
    "    for index_degree, degree in enumerate(degrees):\n",
    "        # ***************************************************\n",
    "        # INSERT YOUR CODE HERE\n",
    "        # ***************************************************\n",
    "\n",
    "        plot_fstar(axs[index_degree])\n",
    "        axs[index_degree].legend()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "bias_variance_demo(0.1, [1, 3, 6])"
   ]
  },
  {
   "cell_type": "markdown",
   "execution_count": null,
   "metadata": {},
   "source": [
    "Your output should ressemble to this: \n",
    "![alt text](bias_variance.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}