Add kfold cross validation for comparison against bootstrapping and comparison between cost functions

src/datamanip.py

@@ -2,6 +2,7 @@ import numpy as np
 from sklearn.preprocessing import StandardScaler  # type: ignore
 from sklearn.metrics import mean_squared_error, r2_score  # type: ignore
 from sklearn.utils import resample  # type: ignore
+from sklearn.model_selection import KFold  # type: ignore
 
 
 def polynomial_features(x: np.ndarray, p: int, intercept: bool = True) -> np.ndarray:
@@ -99,3 +100,23 @@ def bootstrap_resample(
         X_resample, y_resample = resample(X, y)
         resamples.append((X_resample, y_resample))
     return resamples
+
+
+def k_fold_split(
+    X: np.ndarray, y: np.ndarray, k: int
+) -> list[tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]]:
+    """Splits the dataset into k folds for cross-validation.
+    Args:
+        X: The input data matrix of shape (n_samples, n_features).
+        y: The target vector of shape (n_samples,).
+        k: The number of folds.
+    Returns:
+        A list of tuples, each containing (X_train, y_train, X_val, y_val) for each fold.
+    """
+    kf = KFold(n_splits=k, shuffle=True)
+    folds = []
+    for train_index, val_index in kf.split(X):
+        X_train, X_val = X[train_index], X[val_index]
+        y_train, y_val = y[train_index], y[val_index]
+        folds.append((X_train, y_train, X_val, y_val))
+    return folds

src/main.ipynb

@@ -650,7 +650,7 @@
     "y = datamanip.noise_data(datamanip.runge_function(x), 0.1)\n",
     "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)\n",
     "\n",
-    "polynomial_degrees = np.arange(1, 30)\n",
+    "polynomial_degrees = np.arange(1, 50)\n",
     "n_bootstraps = len(x_train)\n",
     "mses = np.zeros((len(polynomial_degrees), n_bootstraps))\n",
     "biases = np.zeros((len(polynomial_degrees), n_bootstraps))\n",
@@ -680,7 +680,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "FILL_BETWEEN = False\n",
+    "FILL_BETWEEN = True\n",
     "fig, axs = plt.subplots(1, 2, figsize=plotting.get_figsize(0.4))\n",
     "axs[0].plot(polynomial_degrees, np.mean(mses, axis=1), label=\"MSE (Test)\", color=\"C0\")\n",
     "\n",
@@ -716,6 +716,7 @@
     "\n",
     "for ax in axs:\n",
     "    ax.set_xlabel(\"Polynomial Degree\")\n",
+    "    ax.set_yscale(\"log\")\n",
     "    ax.legend()\n",
     "\n",
     "fig.tight_layout()\n",
@@ -728,6 +729,145 @@
    "id": "20",
    "metadata": {},
    "outputs": [],
+   "source": [
+    "k_folds = 5\n",
+    "k_fold_mses = np.zeros((len(polynomial_degrees), k_folds))\n",
+    "splits = datamanip.k_fold_split(x, y, k_folds)\n",
+    "\n",
+    "for i, polynomial_degree in enumerate(polynomial_degrees):\n",
+    "    for k, (x_tr, y_tr, x_val, y_val) in enumerate(splits):\n",
+    "        X_train = datamanip.polynomial_features(x_tr, polynomial_degree, False)\n",
+    "        X_val = datamanip.polynomial_features(x_val, polynomial_degree, False)\n",
+    "        X_train_scaled, X_val_scaled = datamanip.scale_data(X_train, X_val)\n",
+    "        y_train_scaled, y_val_scaled = datamanip.scale_data(y_tr, y_val)\n",
+    "\n",
+    "        beta = optimizers.Ridge_parameters(X_train_scaled, y_train_scaled, lam=1e-10)\n",
+    "        y_pred = X_val_scaled @ beta\n",
+    "        mse, _ = datamanip.evaluate_model(y_val_scaled, y_pred)\n",
+    "        k_fold_mses[i, k] = mse"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "21",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "FILL_BETWEEN = False\n",
+    "fig, ax = plt.subplots(figsize=plotting.get_figsize(0.5))\n",
+    "\n",
+    "ax.plot(\n",
+    "    polynomial_degrees, np.mean(mses, axis=1), label=\"Bootstrapping (Test)\", color=\"C0\"\n",
+    ")\n",
+    "ax.plot(\n",
+    "    polynomial_degrees,\n",
+    "    np.mean(k_fold_mses, axis=1),\n",
+    "    label=f\"{k_folds}-Fold (Test)\",\n",
+    "    color=\"C1\",\n",
+    ")\n",
+    "if FILL_BETWEEN:\n",
+    "    ax.fill_between(\n",
+    "        polynomial_degrees,\n",
+    "        np.mean(mses, axis=1) - np.std(mses, axis=1),\n",
+    "        np.mean(mses, axis=1) + np.std(mses, axis=1),\n",
+    "        color=\"C0\",\n",
+    "        alpha=0.3,\n",
+    "    )\n",
+    "    ax.fill_between(\n",
+    "        polynomial_degrees,\n",
+    "        np.mean(k_fold_mses, axis=1) - np.std(k_fold_mses, axis=1),\n",
+    "        np.mean(k_fold_mses, axis=1) + np.std(k_fold_mses, axis=1),\n",
+    "        color=\"C1\",\n",
+    "        alpha=0.3,\n",
+    "    )\n",
+    "\n",
+    "ax.set_xlabel(\"Polynomial Degree\")\n",
+    "ax.set_ylabel(\"Mean Squared Error\")\n",
+    "ax.set_yscale(\"log\")\n",
+    "ax.legend()\n",
+    "fig.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "22",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "polynomial_degrees = np.arange(1, 30)\n",
+    "k_folds = 5\n",
+    "\n",
+    "k_fold_mses_ols = np.zeros((len(polynomial_degrees), k_folds))\n",
+    "k_fold_mses_ridge = np.zeros_like(k_fold_mses_ols)\n",
+    "k_fold_mses_lasso = np.zeros_like(k_fold_mses_ols)\n",
+    "\n",
+    "splits = datamanip.k_fold_split(x, y, k_folds)\n",
+    "\n",
+    "for i, polynomial_degree in enumerate(polynomial_degrees):\n",
+    "    for k, (x_tr, y_tr, x_val, y_val) in enumerate(splits):\n",
+    "        X_train = datamanip.polynomial_features(x_tr, polynomial_degree, False)\n",
+    "        X_val = datamanip.polynomial_features(x_val, polynomial_degree, False)\n",
+    "        X_train_scaled, X_val_scaled = datamanip.scale_data(X_train, X_val)\n",
+    "        y_train_scaled, y_val_scaled = datamanip.scale_data(y_tr, y_val)\n",
+    "\n",
+    "        for optim, results_array in zip(\n",
+    "            [\n",
+    "                optimizers.OLSGradientDescent,\n",
+    "                optimizers.RidgeGradientDescent,\n",
+    "                optimizers.LASSOGradientDescent,\n",
+    "            ],\n",
+    "            [k_fold_mses_ols, k_fold_mses_ridge, k_fold_mses_lasso],\n",
+    "        ):\n",
+    "            Optimizer = optim(num_iterations=1000, learning_rate=0.1, lam=0.1)\n",
+    "            theta = Optimizer.fit(X_train_scaled, y_train_scaled)\n",
+    "            y_pred = X_val_scaled @ theta\n",
+    "            mse, _ = datamanip.evaluate_model(y_val_scaled, y_pred)\n",
+    "            results_array[i, k] = mse"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "23",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "FILL_BETWEEN = False\n",
+    "fig, ax = plt.subplots(figsize=plotting.get_figsize(0.5))\n",
+    "\n",
+    "for i, (mse_array, label) in enumerate(\n",
+    "    zip(\n",
+    "        [k_fold_mses_ols, k_fold_mses_ridge, k_fold_mses_lasso],\n",
+    "        [\"OLS\", \"Ridge Regression\", \"Lasso Regression\"],\n",
+    "    )\n",
+    "):\n",
+    "    c = f\"C{i}\"\n",
+    "    ax.plot(polynomial_degrees, np.mean(mse_array, axis=1), label=label, c=c)\n",
+    "    if FILL_BETWEEN:\n",
+    "        ax.fill_between(\n",
+    "            polynomial_degrees,\n",
+    "            np.mean(mse_array, axis=1) - np.std(mse_array, axis=1),\n",
+    "            np.mean(mse_array, axis=1) + np.std(mse_array, axis=1),\n",
+    "            color=c,\n",
+    "            alpha=0.3,\n",
+    "        )\n",
+    "\n",
+    "ax.set_xlabel(\"Polynomial Degree\")\n",
+    "ax.set_ylabel(f\"MSE ({k_folds}-fold validation)\")\n",
+    "# ax.set_yscale(\"log\")\n",
+    "# ax.set_ylim(0, 2)\n",
+    "ax.legend()\n",
+    "fig.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "24",
+   "metadata": {},
+   "outputs": [],
    "source": []
   }
  ],

src/optimizers.py

@@ -279,6 +279,9 @@ class OLSStochasticGradientDescent(OLSGradientDescent):
         self.N = len(self.y)
         self.indices = np.arange(self.N)
         self.n = self.batch_size
+        np.random.shuffle(self.indices)
+        self.X = self.X[self.indices]
+        self.y = self.y[self.indices]
 
     def _comp_step(self):
         index = np.random.randint(0, self.N)