bias_variance_ridge_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e093828d",
   "metadata": {},
   "source": [
    "Sascha Spors,\n",
    "Professorship Signal Theory and Digital Signal Processing,\n",
    "Institute of Communications Engineering (INT),\n",
    "Faculty of Computer Science and Electrical Engineering (IEF),\n",
    "University of Rostock,\n",
    "Germany\n",
    "\n",
    "# Data Driven Audio Signal Processing - A Tutorial with Computational Examples\n",
    "\n",
    "Winter Semester 2023/24 (Master Course #24512)\n",
    "\n",
    "- lecture: https://github.com/spatialaudio/data-driven-audio-signal-processing-lecture\n",
    "- tutorial: https://github.com/spatialaudio/data-driven-audio-signal-processing-exercise\n",
    "\n",
    "Feel free to contact lecturer frank.schultz@uni-rostock.de"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c86e15f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from statsmodels.api import OLS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97708e79",
   "metadata": {},
   "source": [
    "# Trade-Off Between Bias^2 / Variance; Regularization by Ridge Regression\n",
    "\n",
    "- we use plain **regularized least squares** based **ridge regression** to discuss a very fundamental aspect when we learn from data, i.e. we create prediction models\n",
    "- this aspect is known as bias-variance trade-off\n",
    "- in general we can split the **sum of (true model data - predicted model data)^2** into three components\n",
    "\n",
    "$$\\text{model bias}^2 + \\text{model variance} + \\text{noise variance}$$\n",
    "\n",
    "- a model will never explain all variance (which is actually not wanted for a useful robust model), so certain noise variance remains\n",
    "- we can influence the model bias and model variance obviously by the choice of the model, see [bias_variance_linear_regression.ipynb](bias_variance_linear_regression.ipynb) for usage of different design/feature matrices that set up models with different complexity\n",
    "- however, we cannot at the same time have lowest model bias *and* lowest model variance to reduce the overall error for predictions\n",
    "- we therefore need to find a good compromise between bias and variance and especially we need to avoid two extremes\n",
    "    - underfit case, with typically too low model complexity yielding high bias and low variance\n",
    "    - overfit case, with typically too high model complexity yielding low bias and high variance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e419572",
   "metadata": {},
   "source": [
    "In this notebook we use a model that is in principle capable of overfitting, i.e. it can fit some amount of the noise due to its comparably too high model complexity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e1db974",
   "metadata": {},
   "source": [
    "One way to avoid overfitting to a certain degree is the regularization of the inverse problem. A most simple form of regularization can be used for our well known linear model (i.e. the ordinary least squares (OLS))\n",
    "\n",
    "$$\\min_{\\text{wrt }\\beta} (||\\mathbf{y} - \\mathbf{X} \\beta||_2^2)$$\n",
    "\n",
    "by adding a penalty onto $||\\beta||^2_2$, which leads to the optimization problem\n",
    "\n",
    "$$\\min_{\\text{wrt }\\beta} (||\\mathbf{y} - \\mathbf{X} \\beta||_2^2 + \\alpha^2 ||\\beta||^2_2).$$\n",
    "\n",
    "This is known as **ridge regression** or **Tikhonov regularization**.\n",
    "The approach has the model parameters $\\beta$, but compared to OLS it has **one additional hyper parameter**, which is the regularization value / penalty weight $\\alpha^2>0$ (it is a real-valued scalar).\n",
    "When $\\alpha^2=0$ we obtain ordinary least squares case.\n",
    "\n",
    "Hyper parameters are parameters that are linked to the actual chosen optimization problem, but not to the model architecture, rather to the optimization 'algorithm'/approach.\n",
    "Optimum hyper parameters should be learned within the model's training/fitting stage.\n",
    "We do this later on in the course, cf. [exercise12_HyperParameterTuning.ipynb](exercise12_HyperParameterTuning.ipynb).\n",
    "\n",
    "Note that in textbooks $\\alpha^2$ is very often denoted as variable $\\lambda$.\n",
    "We here rather use a squared variable $\\alpha^2$ to have consistent quantities in the denominator of the \n",
    "singular value inversion\n",
    "$$\\frac{\\sigma_i}{\\sigma_i^2 + \\alpha^2}$$\n",
    "that holds for the ridge regression, cf. [exercise07_left_inverse_SVD_QR.ipynb](exercise07_left_inverse_SVD_QR.ipynb).\n",
    "Just for recap, see the plot we already had below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b562afa",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 1 / 10\n",
    "lmb = alpha**2\n",
    "\n",
    "singval = np.logspace(-4, 4, 2**6)\n",
    "# ridge regression\n",
    "inv_singval = singval / (singval**2 + alpha**2)\n",
    "\n",
    "plt.plot(singval, 1 / singval, label=\"no penalty\")\n",
    "plt.plot(singval, inv_singval, label=\"penalty\")\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xticks(10.0 ** np.arange(-4, 5))\n",
    "plt.yticks(10.0 ** np.arange(-4, 5))\n",
    "plt.axis(\"equal\")\n",
    "plt.xlabel(r\"singular value $\\sigma_i$\")\n",
    "plt.ylabel(\n",
    "    r\"ridge inverted singular value $\\sigma_i \\,\\,\\,/\\,\\,\\, (\\sigma_i^2 + \\alpha^2)$\"\n",
    ")\n",
    "plt.title(\n",
    "    r\"ridge penalty $\\alpha =$\"\n",
    "    + str(alpha)\n",
    "    + r\", $\\alpha^2 =$\"\n",
    "    + str(alpha**2)\n",
    ")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "print(\"alpha =\", alpha, \"alpha^2 = lambda =\", lmb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "967d2b7a",
   "metadata": {},
   "source": [
    "The solution for the ridge regression optimization problem (in machine learning wording: **we trained the model**) can be analytically given and is well known as\n",
    "\n",
    "$$\\hat{\\beta} = (\\mathbf{X}^\\mathrm{H} \\mathbf{X} + \\alpha^2 \\mathbf{I})^{-1} \\mathbf{X}^\\mathrm{H} \\mathbf{y}$$\n",
    "\n",
    "and consistently results in the left inverse solution for $\\alpha^2=0$, i.e. again linear regression with ordinary least squares (OLS).\n",
    "\n",
    "We see, that the estimated model parameters $\\hat{\\beta}$ are influenced by the hyper parameter $\\alpha^2$.\n",
    "We precisely know how exactly, as we've discussed the regularized inverted singular values $\\frac{\\sigma_i}{\\sigma_i^2 + \\alpha^2}$ in [exercise07_left_inverse_SVD_QR.ipynb](exercise07_left_inverse_SVD_QR.ipynb).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d4b7e45",
   "metadata": {},
   "source": [
    "In the examples below we will see that\n",
    "- very **small** $\\alpha^2$ produces **high var**iance, but **low** squared **bias**; hence, we potentially **over**fit the model\n",
    "- very **large** $\\alpha^2$ produces **low** **var**iance, but **high** squared **bias**; hence, we potentially **under**fit the model\n",
    "\n",
    "We could consider an **optimum** regularization amount $\\alpha^2$, where the sum $\\text{model bias}^2 + \\text{model variance}$ is **minimum**. This mean that we found an optimum hyper parameter set, here only the one value $\\alpha^2$.\n",
    "\n",
    "We should realize that regularization does not solve the general problem of choosing an appropriate design/matrix, i.e. an appropriate model and features. For example, if the true model has $f(x^2)$ and the prediction model is set up for $f(x^3)$, it will be hard to train/predict for negative $x$-values, simply because the true model and the prediction model have not too much in common for these independent variables. So, regularization can only help a little in this example. We could try this as a toy example on our own..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "562808a4",
   "metadata": {},
   "source": [
    "Useful chapters in textbooks on bias-variance-tradeoff and ridge regression:\n",
    "- [Bishop 2006] Christopher M. Bishop, *Pattern Recognition and Machine Learning*, Springer, 2006, Chapter 3.2, 3.1.4\n",
    "- Sergios Theodoridis, *Machine Learning*, Academic Press, 2020, 2nd ed., Chapter 3.9, 3.8\n",
    "- Kevin P. Murphy, *Machine Learning-A Probabilistic Perspective*, MIT Press, 2012, 1st ed., Chapter 6.4.4, 7.5\n",
    "- Kevin P. Murphy, *Probabilistic Machine Learning-An Introduction*, MIT Press, 2022, Chapter 4.7.6.3, 11.3\n",
    "- Trevor Hastie, Robert Tibshirani, Jerome Friedman, *The Elements of  Statistical Learning: Data Mining, Inference, and Prediction*, Springer, 2009, 2nd ed., Chapter 2.9, 3.4\n",
    "- Gareth James, Daniela Witten, Trevor Hastie, Rob Tibshirani, *An Introduction to Statistical Learning with Applications in R*, Springer, 2021, 2nd ed., Chapter 2.2.2, 6.2.1\n",
    "- Richard O. Duda, Peter E. Hart, David G. Stork, *Pattern Classification*, Wiley, 2000, 2nd ed., Chapter 9.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6764d6da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for reproducible outputs\n",
    "rng = np.random.default_rng(12345)  # used for data creation and shuffling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "d8ef032a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_dataset(M, split=0.8, noise_scale=2, shuffled=True):\n",
    "    x = np.linspace(0, 2 * np.pi, M)  # lin increase\n",
    "    # shuffle data for simple train/test split handling using [:Ns], [Ns:]\n",
    "    if shuffled:\n",
    "        rng.shuffle(x)\n",
    "\n",
    "    # design/feature matrix of the true model\n",
    "    X = np.column_stack(\n",
    "        (np.cos(1 * x), np.sin(2 * x), np.cos(5 * x), np.cos(6 * x))\n",
    "    )\n",
    "    # add bias/intercept column to the design/feature matrix\n",
    "    X = np.hstack((np.ones((M, 1)), X))\n",
    "    # some nice numbers for the true model parameters beta\n",
    "    beta = np.array([3, 2, 1, 1 / 2, 1 / 4])\n",
    "    # outcome of true model, i.e. linear combination\n",
    "    y_true = (X @ beta)[:, None]\n",
    "    # add measurement noise\n",
    "    noise = rng.normal(loc=0, scale=noise_scale, size=(M, 1))\n",
    "    y = y_true + noise\n",
    "\n",
    "    # design/feature matrix of the prediction model\n",
    "    # we create a model that can overfit the noisy data\n",
    "    # as the feature/design matrix contains also non-matching Fourier series\n",
    "    # components, thus we firstly go for:\n",
    "    # true model Fourier components, same as above\n",
    "    X = np.column_stack(\n",
    "        (np.cos(1 * x), np.sin(2 * x), np.cos(5 * x), np.cos(6 * x))\n",
    "    )\n",
    "    # and then add additional Fourier components, that do not explain our\n",
    "    # y_true but will be sensible to the measurement noise contained in y\n",
    "    if True:  # if False can be used for debug\n",
    "        X = np.column_stack(\n",
    "            (\n",
    "                X,\n",
    "                np.cos(2 * x),\n",
    "                np.cos(3 * x),\n",
    "                np.cos(4 * x),\n",
    "                np.cos(7 * x),\n",
    "                np.cos(8 * x),\n",
    "                np.cos(9 * x),\n",
    "                np.cos(10 * x),\n",
    "                np.sin(1 * x),\n",
    "                np.sin(3 * x),\n",
    "                np.sin(4 * x),\n",
    "                np.sin(5 * x),\n",
    "                np.sin(6 * x),\n",
    "                np.sin(7 * x),\n",
    "                np.sin(8 * x),\n",
    "                np.sin(9 * x),\n",
    "                np.sin(10 * x),\n",
    "            )\n",
    "        )\n",
    "    X = np.hstack((np.ones((M, 1)), X))  # add bias\n",
    "\n",
    "    # split data set into a training data set and a test data set\n",
    "    Ns = int(split * M)\n",
    "    X_trn, X_tst = X[:Ns, :], X[Ns:, :]\n",
    "    y_true_trn, y_true_tst = y_true[:Ns], y_true[Ns:]  # without noise\n",
    "    y_trn, y_tst = y[:Ns], y[Ns:]  # with measurement noise\n",
    "\n",
    "    return x, X_trn, X_tst, y_true_trn, y_true_tst, y_trn, y_tst"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9320e8f5",
   "metadata": {},
   "source": [
    "## Empirical Estimators for Bias and Variance of Ridge Regression Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "e21f304f",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_scale = 5  # standard deviation for added measurement noise\n",
    "\n",
    "M = 2**10  # no of rows in X = no of samples in one full data set\n",
    "\n",
    "Ndatasampling = 2**6  # number of data sets\n",
    "\n",
    "alpha2_min, alpha2_max = -2, +2\n",
    "Nalpha = 2**6 + 1\n",
    "\n",
    "alpha2_vec = np.logspace(alpha2_min, alpha2_max, Nalpha - 1)\n",
    "# insert 0 at end for the case 'no regularization'\n",
    "alpha2_vec = np.insert(alpha2_vec, 0, 0)\n",
    "\n",
    "# create only one data set to get value for Mtest\n",
    "x, X_trn, X_tst, y_true_trn, y_true_tst, y_trn, y_tst = create_dataset(M)\n",
    "Mtest = y_tst.shape[0]\n",
    "\n",
    "# we use capital Y variable to refer to output data\n",
    "# of many models stored into columns of matrices Yxxxx\n",
    "Yh_tst = np.zeros((Mtest, Ndatasampling, Nalpha))\n",
    "Y_true_tst = np.zeros((Mtest, Ndatasampling, Nalpha))\n",
    "\n",
    "# in total we train (Ndatasampling x Nalpha) models in the next for loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "9afb4067",
   "metadata": {},
   "outputs": [],
   "source": [
    "for data_idx in range(Ndatasampling):\n",
    "    # get new data set, split to train/test\n",
    "    x, X_trn, X_tst, y_true_trn, y_true_tst, y_trn, y_tst = create_dataset(\n",
    "        M, noise_scale=noise_scale\n",
    "    )\n",
    "    # set up model with training data\n",
    "    model = OLS(y_trn, X_trn, hasconst=True)\n",
    "    # vary the ridge regression hyper parameter alpha^2\n",
    "    for alpha2_idx, alpha2 in enumerate(alpha2_vec):\n",
    "        # train a model\n",
    "        results = model.fit_regularized(\n",
    "            alpha=alpha2, L1_wt=0, profile_scale=False\n",
    "        )\n",
    "        # predict data using the actual model\n",
    "        Yh_tst[:, data_idx, alpha2_idx] = results.predict(X_tst)\n",
    "        # either: true model data without noise:\n",
    "        Y_true_tst[:, data_idx, alpha2_idx] = np.squeeze(y_true_tst)\n",
    "        # or: take the noisy data\n",
    "        # Y_true_tst[:, data_idx, alpha2_idx] = np.squeeze(y_tst)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c8d2785",
   "metadata": {},
   "source": [
    "Next, we calculate the bias and variance using the data we obtained from many training procedures / many models.\n",
    "\n",
    "Bias and variance is calculated for models that were trained on different train/test data, but for the same hyper parameter $\\alpha^2$.\n",
    "So, for each $\\alpha^2$ we obtain a bias and a variance quantity.\n",
    "We then can estimate that one **optimum hyper parameter** $\\alpha^2_\\mathrm{opt}$, where the sum bias$^2$ + variance is **minimum**.\n",
    "\n",
    "See [bias_variance_linear_regression.ipynb](bias_variance_linear_regression.ipynb) for details on the equations for bias and variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9c3ccf25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get average prediction, i.e. mean over the L models\n",
    "# which is a numerical eval of the expectation:\n",
    "ym = np.mean(Yh_tst, axis=1)  # (3.45) in [Bishop 2006]\n",
    "ym_true = np.mean(Y_true_tst, axis=1)\n",
    "\n",
    "# get integrated squared bias (numerical eval of the expectation):\n",
    "# (3.42), (3.46) in [Bishop 2006]\n",
    "bias_squared = np.mean((ym - ym_true) ** 2, axis=0)\n",
    "\n",
    "# get integrated variance (numerical eval of the expectation):\n",
    "# (3.43), (3.47) in [Bishop 2006]\n",
    "variance = np.mean(\n",
    "    np.mean((Yh_tst - np.expand_dims(ym, axis=1)) ** 2, axis=1), axis=0\n",
    ")\n",
    "\n",
    "# find min for bias_squared+variance\n",
    "idx = np.argmin(bias_squared + variance)\n",
    "# get specific alpha^2 for this min\n",
    "alpha2_opt = alpha2_vec[idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49248a49",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 1, figsize=(8, 4))\n",
    "axs.plot(alpha2_vec, bias_squared, \"C0\", label=r\"bias$^2$\", lw=3)\n",
    "axs.plot(alpha2_vec, variance, \"C1\", label=r\"var\")\n",
    "axs.plot(alpha2_vec, bias_squared + variance, \"C2\", label=r\"bias$^2$+var\")\n",
    "\n",
    "axs.plot(alpha2_opt, bias_squared[idx], \"C0o\")\n",
    "axs.plot(alpha2_opt, variance[idx], \"C1o\")\n",
    "axs.plot(alpha2_opt, bias_squared[idx] + variance[idx], \"C2o\")\n",
    "\n",
    "axs.set_xscale(\"log\")\n",
    "axs.set_yscale(\"log\")\n",
    "axs.set_xlabel(r\"underfit region                                   regularization value $\\alpha^2$                                   overfit region\")\n",
    "axs.set_title(r\"$\\alpha^2_\\mathrm{opt}$=\" + \"{:4.3f}\".format(alpha2_vec[idx]))\n",
    "axs.legend()\n",
    "axs.set_xlim(10**alpha2_min, 10**alpha2_max)\n",
    "axs.set_ylim(1e-2, 1e1)\n",
    "axs.grid(True)\n",
    "axs.xaxis.set_inverted(True)\n",
    "plt.tight_layout()\n",
    "plt.savefig('bias_variance_plots/bias_var_l2_regularisation.png', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d8ee910",
   "metadata": {},
   "source": [
    "## Compare prediction data to true model data and to noisy data for different regularization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22f3db48",
   "metadata": {},
   "source": [
    "Let us next check and visualize what different hyper parameters $\\alpha^2$ do in terms of regularization, i.e. in terms of model complexity, when we predict data and compare this to the true model data and the noisy data.\n",
    "As we draw data from the same sampling distribution as above, we can assume that the estimated $\\alpha^2_\\mathrm{opt}$\n",
    "still holds.\n",
    "So we present four cases for $\\alpha^2$:\n",
    "- $\\alpha^2 = 0$\n",
    "- $\\alpha^2 = 0.01$\n",
    "- $\\alpha^2 = \\alpha^2_\\mathrm{opt}$\n",
    "- $\\alpha^2 = 100$\n",
    "\n",
    "Recall that typically\n",
    "- very **small** $\\alpha^2$ produces **high var**iance, but **low** squared **bias**; hence, we potentially **over**fit the model\n",
    "- very **large** $\\alpha^2$ produces **low** **var**iance, but **high** squared **bias**; hence, we potentially **under**fit the model\n",
    "\n",
    "For convenience (i.e. easier data handling), we use the same data for training/fitting and testing/predicting.\n",
    "We should never do this in real applications, but here we can go for it, as we are interested in the essence of what's going on with simple regularization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4116e08a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# split=1 here, i.e. we realize handling: train data == test data\n",
    "# we do this to conviently show y(x) for shuffled=False data\n",
    "x, X_trn, X_tst, y_true_trn, y_true_tst, y_trn, y_tst = create_dataset(\n",
    "    M, split=1, noise_scale=noise_scale, shuffled=False\n",
    ")\n",
    "model = OLS(y_trn, X_trn, hasconst=True)\n",
    "\n",
    "fig, axs = plt.subplots(2, 2, figsize=(10, 5))\n",
    "fig2, axs2 = plt.subplots(1, 1, figsize=(10, 3))\n",
    "\n",
    "# model with alpha=0\n",
    "results = model.fit_regularized(alpha=0, L1_wt=0, profile_scale=False)\n",
    "yh = results.predict(X_trn)\n",
    "axs2.plot(\n",
    "    np.arange(X_trn.shape[1]),\n",
    "    results.params,\n",
    "    \"C0o-\",\n",
    "    label=r\"$\\alpha^2$=\" + \"{:4.3f}\".format(alpha2_vec[0]),\n",
    ")\n",
    "axs[0, 0].plot(x, y_trn, \"C0\")\n",
    "axs[0, 0].plot(x, y_true_trn, \"C1\")\n",
    "axs[0, 0].plot(x, yh, \"C3\")\n",
    "axs[0, 0].set_title(r\"$\\alpha^2$=\" + \"{:4.3f}\".format(alpha2_vec[0]))\n",
    "\n",
    "# model with alpha^2=0.01 (stored in alpha2_vec[0])\n",
    "results = model.fit_regularized(\n",
    "    alpha=alpha2_vec[0], L1_wt=0, profile_scale=False\n",
    ")\n",
    "yh = results.predict(X_trn)\n",
    "axs2.plot(\n",
    "    np.arange(X_trn.shape[1]),\n",
    "    results.params,\n",
    "    \"C1o:\",\n",
    "    label=r\"$\\alpha^2$=\" + \"{:4.3f}\".format(alpha2_vec[1]),\n",
    ")\n",
    "axs[0, 1].plot(x, y_trn, \"C0\")\n",
    "axs[0, 1].plot(x, y_true_trn, \"C1\")\n",
    "axs[0, 1].plot(x, yh, \"C3\")\n",
    "axs[0, 1].set_title(r\"$\\alpha^2$=\" + \"{:4.3f}\".format(alpha2_vec[1]))\n",
    "\n",
    "# model with optimum alpha (stored in alpha_opt)\n",
    "results = model.fit_regularized(alpha=alpha2_opt, L1_wt=0, profile_scale=False)\n",
    "yh = results.predict(X_trn)\n",
    "axs2.plot(\n",
    "    np.arange(X_trn.shape[1]),\n",
    "    results.params,\n",
    "    \"C2o-\",\n",
    "    label=r\"$\\alpha^2_\\mathrm{opt}$=\" + \"{:4.3f}\".format(alpha2_opt),\n",
    ")\n",
    "axs[1, 0].plot(x, y_trn, \"C0\")\n",
    "axs[1, 0].plot(x, y_true_trn, \"C1\")\n",
    "axs[1, 0].plot(x, yh, \"C3\")\n",
    "axs[1, 0].set_title(r\"$\\alpha^2_\\mathrm{opt}$=\" + \"{:4.3f}\".format(alpha2_opt))\n",
    "\n",
    "# model with alpha^2=100 (stored in alpha2_vec[-1])\n",
    "results = model.fit_regularized(\n",
    "    alpha=alpha2_vec[-1], L1_wt=0, profile_scale=False\n",
    ")\n",
    "yh = results.predict(X_trn)\n",
    "axs2.plot(\n",
    "    np.arange(X_trn.shape[1]),\n",
    "    results.params,\n",
    "    \"C3o-\",\n",
    "    label=r\"$\\alpha^2$=\" + \"{:4.3f}\".format(alpha2_vec[-1]),\n",
    ")\n",
    "axs[1, 1].plot(x, y_trn, \"C0\", label=\"y_train (with noise)\")\n",
    "axs[1, 1].plot(x, y_true_trn, \"C1\", label=\"y_true_train (w/o noise)\")\n",
    "axs[1, 1].plot(x, yh, \"C3\", label=\"predicted y from X_trn\")\n",
    "axs[1, 1].set_title(r\"$\\alpha^2$=\" + \"{:4.3f}\".format(alpha2_vec[-1]))\n",
    "\n",
    "\n",
    "for i in range(2):\n",
    "    for j in range(2):\n",
    "        axs[i, j].grid(True)\n",
    "        axs[i, j].set_xlabel(\"x\")\n",
    "        axs[i, j].set_ylabel(\"y\")\n",
    "        axs[i, j].set_xlim(x[0], x[-1])\n",
    "        axs[i, j].set_ylim(-7, 13)\n",
    "axs[1, 1].legend()\n",
    "fig.tight_layout()\n",
    "\n",
    "axs2.legend()\n",
    "axs2.set_xlabel(\n",
    "    r\"$\\beta$ coefficient index, true model features 0...4, features >4 contribute to overfit\"\n",
    ")\n",
    "axs2.set_ylabel(r\"$\\beta$ value\")\n",
    "axs2.set_title(\"prediction model parameters\")\n",
    "axs2.set_xticks(np.arange(X_trn.shape[1]))\n",
    "axs2.grid(True)\n",
    "fig2.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09838a13",
   "metadata": {},
   "source": [
    "## Estimate optimum regularization value for one specific data set\n",
    "\n",
    "The approach discussed below is known as **hyper parameter tuning**.\n",
    "The essence is as follows:\n",
    "\n",
    "We have certain finite amount of data and want to train a regularized model.\n",
    "We here go for ridge regression, so the model parameters are $\\beta$ and there is one hyper parameter $\\alpha^2$.\n",
    "We split data to training set and test set.\n",
    "We train and test many models where we vary $\\alpha^2$.\n",
    "The model with best performance exhibits $\\alpha^2 = \\alpha^2_\\mathrm{opt}$.\n",
    "Finding this value is known as hyper parameter tuning.\n",
    "The model is afterwards trained for the optimum $\\beta$ weights.\n",
    "\n",
    "**Note that for real applications**, the final performance check should use a test data set, which was never used for the hyper parameter training.\n",
    "This implies that the full data set is actually split into at least\n",
    "\n",
    "- **train** (train model in hyper parameter tuning stage)\n",
    "- **dev** (predict data in hyper parameter tuning stage)\n",
    "- **test** (use only once to predict data for final performance estimation)\n",
    "\n",
    "sets.\n",
    "For final training (i.e. using the optimum hyper parameter $\\alpha^2$ to estimate optimum model weights $\\beta$) often either the train data or the train+dev data is used for training.\n",
    "\n",
    "Below, for convenience, we only use a split into two data sets. Splitting into train, dev, test might be a nice homework..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "777d78d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we use very noisy y data\n",
    "# we have split=1 and shuffled=False here, as we split and shuffle manually\n",
    "# below\n",
    "# we do this because we want to concatenate and re-sort the data after\n",
    "# training / testing in order to plot the data nicely\n",
    "x, X_trn, _, y_true_trn, _, y_trn, _ = create_dataset(\n",
    "    M, split=1, noise_scale=10, shuffled=False\n",
    ")\n",
    "\n",
    "if False:  # we could use just the features that correspond to the true model\n",
    "    X_trn = X_trn[:, 0:5]\n",
    "\n",
    "# for shuffling data\n",
    "idx = np.arange(M)\n",
    "rng.shuffle(idx)\n",
    "# shuffle data\n",
    "x = x[idx]\n",
    "X = X_trn[idx, :]\n",
    "y_true = y_true_trn[idx]\n",
    "y = y_trn[idx]\n",
    "# split data\n",
    "split = 0.8  # 80 % go into training data, 20% into test data\n",
    "Ns = int(split * M)\n",
    "X_trn, X_tst = X[:Ns, :], X[Ns:, :]\n",
    "y_true_trn, y_true_tst = y_true[:Ns], y_true[Ns:]  # without noise\n",
    "y_trn, y_tst = y[:Ns], y[Ns:]  # with measurement noise\n",
    "# set up OLS model\n",
    "model = OLS(y_trn, X_trn, hasconst=True)\n",
    "# we use capital Y to refer to output data of many models stored into matrices\n",
    "Yh_trn = np.zeros((Ns, Nalpha))\n",
    "Yh_tst = np.zeros((M - Ns, Nalpha))\n",
    "# train/predict for different regularization\n",
    "for alpha2_idx, alpha2 in enumerate(alpha2_vec):\n",
    "    results = model.fit_regularized(alpha=alpha2, L1_wt=0, profile_scale=False)\n",
    "    Yh_trn[:, alpha2_idx] = results.predict(X_trn)\n",
    "    Yh_tst[:, alpha2_idx] = results.predict(X_tst)\n",
    "# residual check for train / test data compared with\n",
    "# true! data (which we don't have in pratice)\n",
    "SSE_trn = np.sum((Yh_trn - y_true_trn) ** 2, axis=0)\n",
    "SSE_tst = np.sum((Yh_tst - y_true_tst) ** 2, axis=0)\n",
    "# get optimum alpha2 where smallest SSE_trn+SSE_tst is obtained\n",
    "alpha2_opt_idx = np.argmin(SSE_trn + SSE_tst)\n",
    "alpha2_opt = alpha2_vec[alpha2_opt_idx]\n",
    "print(alpha2_opt)\n",
    "\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.plot(alpha2_vec, SSE_trn, \"C0\", label=\"sum squared errors train data\")\n",
    "plt.plot(alpha2_vec, SSE_tst, \"C2\", label=\"sum squared errors test data\")\n",
    "plt.plot(alpha2_vec, SSE_trn + SSE_tst, \"C3\", label=\"SSE train + SSE test\")\n",
    "plt.plot(\n",
    "    alpha2_vec[alpha2_opt_idx],\n",
    "    SSE_trn[alpha2_opt_idx] + SSE_tst[alpha2_opt_idx],\n",
    "    \"C3o\",\n",
    ")\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlim(10**alpha2_min, 10**alpha2_max)\n",
    "plt.xlabel(r\"regularization value $\\alpha^2$\")\n",
    "plt.title(r\"$\\alpha^2_\\mathrm{opt}$=\" + \"{:4.3f}\".format(alpha2_opt))\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffedd52b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# concatenate train and test data\n",
    "# re-shuffle them, so bring it into original order, such that x is increasing\n",
    "y_true = np.concatenate((y_true_trn, y_true_tst))\n",
    "y_true = y_true[np.argsort(idx)]\n",
    "y = np.concatenate((y_trn, y_tst))\n",
    "y = y[np.argsort(idx)]\n",
    "yh = np.concatenate((Yh_trn[:, alpha2_opt_idx], Yh_tst[:, alpha2_opt_idx]))\n",
    "yh = yh[np.argsort(idx)]\n",
    "# plot y(x)\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.plot(x[np.argsort(idx)], y, \"C0\", label=\"measured y (with noise)\")\n",
    "plt.plot(x[np.argsort(idx)], y_true, \"C1\", label=\"true y (w/o noise)\")\n",
    "plt.plot(x[np.argsort(idx)], yh, \"C3\", label=\"predicted y\")\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.xlim(x[np.argsort(idx)][0], x[np.argsort(idx)][-1])\n",
    "plt.ylim(-5, 10)\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e301efb9",
   "metadata": {},
   "source": [
    "## Copyright\n",
    "\n",
    "- the notebooks are provided as [Open Educational Resources](https://en.wikipedia.org/wiki/Open_educational_resources)\n",
    "- the text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/)\n",
    "- the code of the IPython examples is licensed under the [MIT license](https://opensource.org/licenses/MIT)\n",
    "- feel free to use the notebooks for your own purposes\n",
    "- please attribute the work as follows: *Frank Schultz, Data Driven Audio Signal Processing - A Tutorial Featuring Computational Examples, University of Rostock* ideally with relevant file(s), github URL https://github.com/spatialaudio/data-driven-audio-signal-processing-exercise, commit number and/or version tag, year."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "myddasp",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}