regularization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Biol 359A | Regularization\n",
    "### Spring 2024, Week 6\n",
    "Objectives:\n",
    "- Understand the Bias-Variance Tradeoff: Learn how increasing model complexity impacts bias and variance.\n",
    "- Gain intuition for Regularization Techniques: Implement Ridge and LASSO to prevent overfitting.\n",
    "- Implement Model Evaluation and Selection: Employ train-test splits for optimal model selection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.preprocessing import PolynomialFeatures, StandardScaler\n",
    "from sklearn.linear_model import LinearRegression, Ridge, Lasso\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "from ipywidgets import interactive, IntSlider, FloatSlider, Dropdown, fixed, Output\n",
    "import seaborn as sns\n",
    "import textwrap\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Today we will start by working with in-silico data. The code below will generate the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set random seed for reproducibility\n",
    "np.random.seed(9)\n",
    "\n",
    "# Generate x values drawn from normal distribution with mean 4.5 and std deviation of 5\n",
    "x = np.linspace(0, 10, 10)\n",
    "\n",
    "# Define a 5th order polynomial function\n",
    "# Coefficients for the polynomial: ax^5 + bx^4 + cx^3 + dx^2 + ex + f\n",
    "a, b, c, d, e, f = 0.0005, -0.001, 0.002, -0.001, 0.005, 2\n",
    "\n",
    "# Calculate y values based on the polynomial function\n",
    "y = a*x**5 + b*x**4 + c*x**3 + d*x**2 + e*x + f\n",
    "\n",
    "# Add random noise\n",
    "noise = np.random.randn(10) * 3\n",
    "y_noisy = y + noise\n",
    "\n",
    "# Create a new array that combines x and y, sort by x, then separate them\n",
    "combined = np.column_stack((x, y_noisy))\n",
    "sorted_combined = combined[np.argsort(combined[:, 0])]\n",
    "x_sorted = sorted_combined[:, 0]\n",
    "y_sorted = sorted_combined[:, 1]\n",
    "\n",
    "# Plot the sorted data\n",
    "plt.scatter(x_sorted, y_sorted, marker=\".\")\n",
    "plt.title('In-silico Data')\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to build a linear regression model that will help us learn about the system that generated these data. To make the coefficients of our model more comparable to each other, we need to start by standardizing our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_standardized_distributions(data, title, color):\n",
    "    # Calculate mean and standard deviation\n",
    "    mean = np.mean(data)\n",
    "    std = np.std(data)\n",
    "    \n",
    "    # Standardize the data\n",
    "    data_zscore = (data - mean) / std\n",
    "    \n",
    "    # Calculate z-scored mean and std deviation\n",
    "    zscored_mean = np.mean(data_zscore)\n",
    "    zscored_std = np.std(data_zscore)\n",
    "    \n",
    "    # Plotting the original and standardized data\n",
    "    plt.figure(figsize=(10, 5))\n",
    "    sns.histplot(data, kde=True, color=color, label='Original Data', element='step', stat='density')\n",
    "    sns.histplot(data_zscore, kde=True, color='orange', label='Standardized Data', element='step', stat='density')\n",
    "    plt.xlabel(f'{title} values')\n",
    "    plt.ylabel('Density')\n",
    "    plt.axvline(mean, color=color, linestyle='-', label='Mean')\n",
    "    plt.axvline(mean + std, color=color, linestyle='--', label='Mean ± Std Dev')\n",
    "    plt.axvline(mean - std, color=color, linestyle='--')\n",
    "    plt.axvline(zscored_mean, color='orange', linestyle='-')\n",
    "    plt.axvline(zscored_mean + zscored_std, color='orange', linestyle='--')\n",
    "    plt.axvline(zscored_mean - zscored_std, color='orange', linestyle='--')\n",
    "    plt.title(f'{title} Data Distribution')\n",
    "    legend = plt.legend()\n",
    "    # Set all line colors in the legend to black\n",
    "    for line in legend.get_lines():\n",
    "        line.set_color('black')\n",
    "    plt.show()\n",
    "\n",
    "# Plot for x data\n",
    "plot_standardized_distributions(x, 'X', 'blue')\n",
    "plot_standardized_distributions(y_noisy, 'Y', 'blue')\n",
    "\n",
    "# Standardizing the data\n",
    "scaler_x = StandardScaler()\n",
    "x_standardized = scaler_x.fit_transform(x_sorted[:, np.newaxis]).flatten()\n",
    "scaler_y = StandardScaler()\n",
    "y_standardized = scaler_y.fit_transform(y_sorted[:, np.newaxis]).flatten()\n",
    "\n",
    "# Plot the unstandardized and standardized data\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.scatter(x_sorted, y_sorted, color='blue', alpha=0.6, label='Unstandardized Data', marker=\".\")\n",
    "plt.scatter(x_standardized, y_standardized, color='orange', alpha=0.6, label='Standardized Data', marker=\".\")\n",
    "plt.title('Comparison of Standardized and Unstandardized Data')\n",
    "plt.xlabel('Feature X (Blue: Original, Orange: Standardized)')\n",
    "plt.ylabel('Target Y')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ASSIGNMENT QUESTIONS:\n",
    "- Please answer question 10 in the assignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Bias-Variance Tradeoff\n",
    "\n",
    "- Bias: Bias is the error introduced by approximating a real-world problem, which might be highly complex, with a simpler model. It reflects how well the model captures the underlying patterns of the data. High bias can cause the model to miss relevant relations between features and target outputs (underfitting).\n",
    "\n",
    "- Variance: Variance measures how much the model's predictions vary for a given dataset. A model with high variance pays a lot of attention to training data and learns noise as well as signals, leading to less generalizable models on new data (overfitting)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interactive plot below allows you to explore this concept. The interactive slider allows you to adjust the degree of the polynomial regression model applied to synthetic data. The polynomial features (like x, x^2) and the corresponding model coefficients are displayed in the model equation and the resulting model predictions are plotted on the synthetic data. As you increase the polynomial degree, observe how the model predictions change:\n",
    "\n",
    "- Low-degree polynomials may show high bias, where the model is too simplistic to capture the essential patterns in the data, resulting in a poor fit.\n",
    "- High-degree polynomials may exhibit high variance, fitting not only the underlying data pattern but also the noise. This is seen as the model begins to closely trace minor fluctuations in the data, indicating overfitting.\n",
    "\n",
    "As you adjust the slider, observe changes in the Mean Squared Error (MSE) and R-squared values below the plot. The MSE will help you understand the average squared difference between the observed and predicted values, while the R-squared value indicates the proportion of the variance in the dependent variable that is predictable from the independent variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_polynomial_complexity_MLR(degree):\n",
    "    polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)\n",
    "    x_poly = polynomial_features.fit_transform(x_standardized[:, np.newaxis])\n",
    "    \n",
    "    model = LinearRegression()\n",
    "    model.fit(x_poly, y_standardized)\n",
    "    y_poly_pred = model.predict(x_poly)\n",
    "    \n",
    "    # Generating polynomial feature names and equation\n",
    "    feature_names = ['x^{}'.format(i) if i > 1 else 'x' for i in range(1, degree + 1)]\n",
    "    coefficients = model.coef_\n",
    "    equation_terms = [f\"{coeff:.3f}*{name}\" for coeff, name in zip(coefficients, feature_names)]\n",
    "    equation = \" + \".join(equation_terms)\n",
    "    wrapped_equation = '\\n               '.join(textwrap.wrap(equation, width=80))  # Wrap the equation\n",
    "    \n",
    "    # Print the complete equation\n",
    "    print(\"Equation: y = \" + wrapped_equation)\n",
    "    \n",
    "    # Plotting\n",
    "    plt.figure(figsize=(10, 5))\n",
    "    plt.scatter(x_standardized, y_standardized, s=20, label=\"Noisy Data Points\")\n",
    "    plt.plot(x_standardized, y_poly_pred, color='r', label=f\"Polynomial Degree {degree} Fit\")\n",
    "    plt.xlabel(\"Feature X\")\n",
    "    plt.ylabel(\"Target Y\")\n",
    "    plt.title(\"Polynomial Regression Fit\")\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    \n",
    "    # Calculate and print errors as a measure of Bias and Variance\n",
    "    mse = mean_squared_error(y_standardized, y_poly_pred)\n",
    "    rsquared = model.score(x_poly, y_standardized)\n",
    "    print(f\"Mean Squared Error (MSE) for Degree {degree}: {mse:.2f}\")\n",
    "    print(f\"R-Squared for Degree {degree}: {rsquared:.2f}\")\n",
    "\n",
    "# Create interactive widget\n",
    "degree_slider = IntSlider(min=1, max=11, step=1, value=1, description=\"Polynomial Degree\",\n",
    "                          style={'description_width': 'initial'},\n",
    "                          layout={'width': '500px'})\n",
    "interactive_plot = interactive(plot_polynomial_complexity_MLR, degree=degree_slider)\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DISCUSSION QUESTION:\n",
    "- Are there $\\beta_0$ terms in any of the model equations? Why or why not?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ASSIGNMENT QUESTIONS:\n",
    "- Please answer questions 11, 12, and 13 in the assignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ridge regression\n",
    "\n",
    "The interactive visualization below performs Ridge regression on these data. Use the sliders to adjust the polynomial degree and the lambda value, which controls the strength of the regularization:\n",
    "\n",
    "- Polynomial Degree: Influences model complexity. Lower degrees might not capture all the data variability (underfitting), while higher degrees might model the noise as well (overfitting).\n",
    "- Lambda Value: Adjusts the Ridge regularization strength. Increasing lambda enhances the model's generalization by penalizing large coefficients, helping to prevent overfitting.\n",
    "\n",
    "Observe how adjusting the lambda values impacts the coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_polynomial_complexity_Ridge(degree, alpha):\n",
    "    polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)\n",
    "    x_poly = polynomial_features.fit_transform(x_standardized[:, np.newaxis])\n",
    "    \n",
    "    # Using Ridge regression model\n",
    "    model = Ridge(alpha=alpha)\n",
    "    model.fit(x_poly, y_standardized)\n",
    "    y_poly_pred = model.predict(x_poly)\n",
    "    \n",
    "    # Generating polynomial feature names and equation\n",
    "    feature_names = ['x^{}'.format(i) if i > 1 else 'x' for i in range(1, degree + 1)]\n",
    "    coefficients = model.coef_\n",
    "    equation_terms = [f\"{coeff:.3f}*{name}\" for coeff, name in zip(coefficients, feature_names)]\n",
    "    equation = \" + \".join(equation_terms)\n",
    "    wrapped_equation = '\\n               '.join(textwrap.wrap(equation, width=80))  # Wrap the equation\n",
    "    \n",
    "    # Print the complete equation\n",
    "    print(\"Equation: y = \" + wrapped_equation)\n",
    "    \n",
    "    # Plotting\n",
    "    plt.figure(figsize=(10, 5))\n",
    "    plt.scatter(x_standardized, y_standardized, s=20, label=\"Standardized Data Points\")\n",
    "    plt.plot(x_standardized, y_poly_pred, color='r', label=f\"Ridge Degree {degree} Fit, Alpha {alpha}\")\n",
    "    plt.xlabel(\"Standardized Feature X\")\n",
    "    plt.ylabel(\"Standardized Target Y\")\n",
    "    plt.title(\"Ridge Polynomial Regression Fit\")\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    \n",
    "    # Calculate and print errors as a measure of Bias and Variance\n",
    "    mse = mean_squared_error(y_standardized, y_poly_pred)\n",
    "    rsquared = model.score(x_poly, y_standardized)\n",
    "    print(f\"Mean Squared Error (MSE) for Degree {degree}: {mse:.2f}\")\n",
    "    print(f\"R-Squared for Degree {degree}: {rsquared:.2f}\")\n",
    "\n",
    "# Create interactive widgets\n",
    "degree_slider = IntSlider(min=1, max=11, step=1, value=1, description=\"Polynomial Degree\",\n",
    "                          style={'description_width': 'initial'},\n",
    "                          layout={'width': '500px'})\n",
    "lambda_slider = FloatSlider(min=.01, max=1, step=.01, value=.1, description=\"Lambda\",\n",
    "                         style={'description_width': 'initial'},\n",
    "                         layout={'width': '500px'})\n",
    "\n",
    "# Keen observers may notice that we rename lambda to alpha in the function (see below). This is because lambda is a reserved keyword in Python.\n",
    "interactive_plot = interactive(plot_polynomial_complexity_Ridge, degree=degree_slider, alpha=lambda_slider)\n",
    "interactive_plot\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ASSIGNMENT QUESTION:\n",
    "- Please answer question 14 in the assignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Lasso\n",
    "\n",
    "This interactive visualization showcases Lasso regression applied to polynomial models. Adjust the polynomial degree and alpha value using the provided sliders to explore their effects:\n",
    "\n",
    "- Polynomial Degree: Controls the complexity of the model. Lower degrees might result in underfitting, where the model is too simplistic, while higher degrees might lead to overfitting, capturing noise instead of just the underlying data pattern.\n",
    "- Lambda Value: Manages the Lasso regularization strength, promoting sparsity in the model coefficients. Higher lambda values can lead to more coefficients being reduced to zero, simplifying the model and potentially improving its generalizability.\n",
    "\n",
    "Observe how changing the lambda value affects the number of parameters in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_polynomial_complexity_Lasso(degree, alpha):\n",
    "    polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)\n",
    "    x_poly = polynomial_features.fit_transform(x_standardized[:, np.newaxis])\n",
    "    \n",
    "    # Using Lasso regression model\n",
    "    model = Lasso(alpha=alpha, max_iter=10000)\n",
    "    model.fit(x_poly, y_standardized)\n",
    "    y_poly_pred = model.predict(x_poly)\n",
    "    \n",
    "    # Generating polynomial feature names and equation\n",
    "    feature_names = ['x^{}'.format(i) if i > 1 else 'x' for i in range(1, degree + 1)]\n",
    "    coefficients = model.coef_\n",
    "    equation_terms = [f\"{coeff:.3f}*{name}\" if coeff != 0 else \"\" for coeff, name in zip(coefficients, feature_names)]\n",
    "    equation_terms = [term for term in equation_terms if term]  # Remove empty strings\n",
    "    equation = \" + \".join(equation_terms)\n",
    "    wrapped_equation = '\\n               '.join(textwrap.wrap(equation, width=80))  # Wrap the equation\n",
    "    \n",
    "    # Print the complete equation\n",
    "    print(\"Equation: y = \" + wrapped_equation)\n",
    "    \n",
    "    # Plotting\n",
    "    plt.figure(figsize=(10, 5))\n",
    "    plt.scatter(x_standardized, y_standardized, s=20, label=\"Standardized Data Points\")\n",
    "    plt.plot(x_standardized, y_poly_pred, color='r', label=f\"Lasso Degree {degree} Fit, Alpha {alpha}\")\n",
    "    plt.xlabel(\"Standardized Feature X\")\n",
    "    plt.ylabel(\"Standardized Target Y\")\n",
    "    plt.title(\"Lasso Polynomial Regression Fit\")\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    \n",
    "    # Calculate and print errors as a measure of Bias and Variance\n",
    "    mse = mean_squared_error(y_standardized, y_poly_pred)\n",
    "    rsquared = model.score(x_poly, y_standardized)\n",
    "    print(f\"Mean Squared Error (MSE) for Degree {degree}: {mse:.2f}\")\n",
    "    print(f\"R-Squared for Degree {degree}: {rsquared:.2f}\")\n",
    "\n",
    "# Create interactive widgets\n",
    "degree_slider = IntSlider(min=1, max=11, step=1, value=1, description=\"Polynomial Degree\",\n",
    "                          style={'description_width': 'initial'},\n",
    "                          layout={'width': '500px'})\n",
    "lambda_slider = FloatSlider(min=0.01, max=1, step=.01, value=.1, description=\"Lambda\",\n",
    "                         style={'description_width': 'initial'},\n",
    "                         layout={'width': '500px'})\n",
    "\n",
    "interactive_plot = interactive(plot_polynomial_complexity_Lasso, degree=degree_slider, alpha=lambda_slider)\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ASSIGNMENT QUESTION\n",
    "- Please answer question 15 in the assignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are not going to visualize the model fit with the data, we are just going to return a table. This interactive tool allows you to compare three different regression models: Standard Multilinear Regression (MLR), Lasso, and Ridge. Each model is evaluated using polynomial features of the data, which you can adjust using the provided sliders.\n",
    "\n",
    "Use the sliders to:\n",
    "- Adjust the Polynomial Degree: Use the degree slider to change the complexity of the polynomial features from 1 to 11. This adjustment affects how well each model can potentially fit the data.\n",
    "- Set the Regularization Strength (Lambda): The lambda slider controls the alpha value for the Lasso and Ridge models. It adjusts the intensity of the penalty applied to the models' coefficients to prevent overfitting. For Standard MLR, this setting has no effect as there is no regularization.\n",
    "\n",
    "The results table reports:\n",
    "- Sum of Squared Errors (SSE): the total error of the predictions. Lower values indicate a better fit to the data.\n",
    "- Penalties (L1 and L2): the sum of the absolute values of the coefficients (L1) and the sum of the squares of the coefficients (L2). These values are relevant for Lasso and Ridge, respectively, and help in assessing the impact of regularization.\n",
    "- Number of Parameters: how many coefficients are actively used in the model. Lasso may set some coefficients to zero, reducing the effective number of parameters.\n",
    "- Full Cost: combines the SSE with the penalties, giving a comprehensive measure of model performance considering both fit and complexity.\n",
    "\n",
    "Review the full cost and the number of parameters to gauge each model's efficiency and complexity. Consider models with lower SSE and appropriate costs, which balance fit and simplicity, to prevent overfitting. Use the R-squared values to understand how well each model's variability explains the dependent variable.\n",
    "\n",
    "Choose the model (type, polynomial degree, and lambda if relevant) that you believe best fits the data based on these metrics. Plot that model's parity plot using the next code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_models(degree, alpha_value):\n",
    "    # Prepare the polynomial features\n",
    "    polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)\n",
    "    x_poly = polynomial_features.fit_transform(x_standardized[:, np.newaxis])\n",
    "\n",
    "    # Initialize models with a dynamic alpha value for regularization\n",
    "    models = {\n",
    "        'Linear': LinearRegression(),\n",
    "        'Lasso': Lasso(alpha=alpha_value, max_iter=10000),\n",
    "        'Ridge': Ridge(alpha=alpha_value)\n",
    "    }\n",
    "\n",
    "    # Table to store results\n",
    "    results = []\n",
    "\n",
    "    # Fit models and calculate errors and penalties\n",
    "    for name, model in models.items():\n",
    "        model.fit(x_poly, y_standardized)\n",
    "        y_poly_pred = model.predict(x_poly)\n",
    "\n",
    "        # Calculating SSE\n",
    "        sse = np.sum((y_standardized - y_poly_pred) ** 2)\n",
    "\n",
    "        # Calculating L1 and L2 penalties\n",
    "        L1_penalty = np.sum(np.abs(model.coef_))\n",
    "        L2_penalty = np.sum(np.square(model.coef_))\n",
    "\n",
    "        # Number of parameters\n",
    "        num_params = len([param for param in model.coef_ if param != 0])\n",
    "        \n",
    "        results.append({\n",
    "            'Model': name,\n",
    "            'SSE': sse,\n",
    "            'Penalty': 'L1' if name == 'Lasso' else 'L2' if name == 'Ridge' else 'None',\n",
    "            'L1': L1_penalty,\n",
    "            'L2': L2_penalty,\n",
    "            'Number of Parameters' : num_params,\n",
    "        })\n",
    "\n",
    "    # Create a DataFrame to display the results\n",
    "    results_df = pd.DataFrame(results)\n",
    "    display(results_df.set_index('Model'))  # Use display to show the DataFrame in the interactive output\n",
    "\n",
    "# Create interactive widgets for the degree of the polynomial and alpha value\n",
    "degree_slider = IntSlider(min=1, max=11, step=1, value=3, description=\"Polynomial Degree\",\n",
    "                          style={'description_width': 'initial'},\n",
    "                          layout={'width': '500px'})\n",
    "lambda_slider = FloatSlider(min=0.01, max=1.0, step=0.01, value=0.1, description=\"Lambda\",\n",
    "                           style={'description_width': 'initial'},\n",
    "                           layout={'width': '500px'})\n",
    "\n",
    "interactive_plot = interactive(compare_models, degree=degree_slider, alpha_value=lambda_slider)\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_parity(model_type, degree, alpha_value):\n",
    "    # Prepare the polynomial features\n",
    "    polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)\n",
    "    x_poly = polynomial_features.fit_transform(x_standardized[:, np.newaxis])\n",
    "\n",
    "    # Choose the model based on user input\n",
    "    if model_type == 'Standard MLR':\n",
    "        model = LinearRegression()\n",
    "    elif model_type == 'Lasso':\n",
    "        model = Lasso(alpha=alpha_value, max_iter=10000)\n",
    "    elif model_type == 'Ridge':\n",
    "        model = Ridge(alpha=alpha_value)\n",
    "\n",
    "    # Fit the model\n",
    "    model.fit(x_poly, y_standardized)\n",
    "    y_poly_pred = model.predict(x_poly)\n",
    "\n",
    "    # Calculating the R^2 score\n",
    "    r_squared = r2_score(y_standardized, y_poly_pred)\n",
    "\n",
    "    # Plotting\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    plt.scatter(y_standardized, y_poly_pred, alpha=0.7, color='dodgerblue', label=f'Predicted vs. Actual')\n",
    "    plt.plot(y_standardized, y_standardized, color='red', linestyle='--', label='Ideal Fit')\n",
    "    plt.xlabel('Actual Values')\n",
    "    plt.ylabel('Predicted Values')\n",
    "    plt.title(f'Parity Plot: {model_type} (Degree {degree}, Lambda {alpha_value})\\n$R^2={r_squared:.3f}$')\n",
    "    plt.legend()\n",
    "    plt.grid(True)\n",
    "    plt.show()\n",
    "\n",
    "# Creating interactive widgets\n",
    "model_dropdown = Dropdown(\n",
    "    options=['Standard MLR', 'Lasso', 'Ridge'],\n",
    "    value='Standard MLR',\n",
    "    description='Model Type:',\n",
    "    disabled=False,\n",
    ")\n",
    "\n",
    "degree_slider = IntSlider(\n",
    "    min=1, \n",
    "    max=12, \n",
    "    step=1, \n",
    "    value=3, \n",
    "    description=\"Degree\",\n",
    "    style={'description_width': 'initial'},\n",
    "    layout={'width': '400px'}\n",
    ")\n",
    "\n",
    "lambda_slider = FloatSlider(\n",
    "    min=0.01, \n",
    "    max=1.0, \n",
    "    step=0.01, \n",
    "    value=0.1, \n",
    "    description=\"Lambda\",\n",
    "    style={'description_width': 'initial'},\n",
    "    layout={'width': '400px'}\n",
    ")\n",
    "\n",
    "# Create interactive output\n",
    "interactive_output = interactive(plot_parity, model_type=model_dropdown, degree=degree_slider, alpha_value=lambda_slider)\n",
    "\n",
    "# Display the interactive output and sliders\n",
    "display(interactive_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discussion question:\n",
    "- What degree did you choose? Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, Now we are going to pretend we gathered a lot more data of the same system. Set the variables in the first three lines to make the function plot the model you selected in the exercise above. How good is the fit? What changes can you make to the model to improve the predictions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SET HYPERPARAMETERS HERE:\n",
    "initial_model_type = 'Standard MLR' # Options are 'Standard MLR', 'Lasso', 'Ridge'\n",
    "initial_degree = 1\n",
    "initial_lambda = 0.1\n",
    "# ------------------------------------\n",
    "\n",
    "# Generate more data\n",
    "x_new = np.linspace(0, 10, 100)  # More data points for better visualization\n",
    "y_new = 0.0005*x_new**5 - 0.001*x_new**4 + 0.002*x_new**3 - 0.001*x_new**2 + 0.005*x_new + 2\n",
    "noise_new = np.random.randn(100) * 3  # Same noise level\n",
    "y_noisy_new = y_new + noise_new\n",
    "\n",
    "# Standardize the new data\n",
    "scaler_x = StandardScaler()\n",
    "x_standardized_new = scaler_x.fit_transform(x_new[:, np.newaxis]).flatten()\n",
    "scaler_y = StandardScaler()\n",
    "y_standardized_new = scaler_y.fit_transform(y_noisy_new[:, np.newaxis]).flatten()\n",
    "\n",
    "def plot_model_predictions(model_type, degree, alpha):\n",
    "    polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)\n",
    "    x_poly_new = polynomial_features.fit_transform(x_standardized_new[:, np.newaxis])\n",
    "    \n",
    "    if model_type == 'Standard MLR':\n",
    "        model = LinearRegression()\n",
    "    elif model_type == 'Lasso':\n",
    "        model = Lasso(alpha=alpha, max_iter=10000)\n",
    "    elif model_type == 'Ridge':\n",
    "        model = Ridge(alpha=alpha)\n",
    "        \n",
    "    model.fit(x_poly_new, y_standardized_new)\n",
    "    y_poly_pred_new = model.predict(x_poly_new)\n",
    "\n",
    "    plt.figure(figsize=(10, 6))\n",
    "    plt.scatter(x_standardized_new, y_standardized_new, color='blue', label='Actual Data')\n",
    "    plt.plot(x_standardized_new, y_poly_pred_new, color='red', label='Model Prediction')\n",
    "    plt.title(f'Model Type: {model_type}, Degree: {degree}, Lambda: {alpha}\\nR^2: {r2_score(y_standardized_new, y_poly_pred_new):.3f}')\n",
    "    plt.xlabel('Standardized X')\n",
    "    plt.ylabel('Standardized Y')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "# Interactive widgets setup\n",
    "model_dropdown = Dropdown(options=['Standard MLR', 'Lasso', 'Ridge'], value=initial_model_type, description='Model Type:')\n",
    "degree_slider = IntSlider(min=1, max=11, value=initial_degree, description='Polynomial Degree', continuous_update=False)\n",
    "lambda_slider = FloatSlider(min=0.01, max=1.0, value=initial_lambda, step=0.01, description='Lambda', continuous_update=False)\n",
    "\n",
    "interactive_plot = interactive(plot_model_predictions, model_type=model_dropdown, degree=degree_slider, alpha=lambda_slider)\n",
    "display(interactive_plot)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the time, we can't collect a bunch of new data to see how well our model generalizes. Instead, we have to split our existing data into \"training data\", which we use to build our model and select which model we think best describes our data, and \"testing data\" which we use to assess our model's generalizability. The code below splits the data from above into a training data set (which is a randomly selected 70% of the data) and a testing set (the remaining 30%)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Separates data into training and test sets\n",
    "def train_test_split(x, y, test_size=0.2):\n",
    "    # Shuffle the data\n",
    "    indices = np.random.permutation(len(x))\n",
    "    x_shuffled = x[indices]\n",
    "    y_shuffled = y[indices]\n",
    "    \n",
    "    # Split the data into training and test sets\n",
    "    split_index = int(len(x) * (1 - test_size))\n",
    "    x_train, x_test = x_shuffled[:split_index], x_shuffled[split_index:]\n",
    "    y_train, y_test = y_shuffled[:split_index], y_shuffled[split_index:]\n",
    "    \n",
    "    return x_train, x_test, y_train, y_test\n",
    "\n",
    "# Split the data\n",
    "x_train, x_test, y_train, y_test = train_test_split(x_standardized_new, y_standardized_new)\n",
    "\n",
    "# Plot training data and test data\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.scatter(x_train, y_train, color='blue', label='Training Data')\n",
    "plt.scatter(x_test, y_test, color='red', label='Test Data')\n",
    "plt.xlabel('Standardized X')\n",
    "plt.ylabel('Standardized Y')\n",
    "plt.title('Training and Test Data')\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets repeat our model exploration on just the training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_models(degree, alpha_value):\n",
    "    # Prepare the polynomial features\n",
    "    polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)\n",
    "    x_poly = polynomial_features.fit_transform(x_train[:, np.newaxis])\n",
    "\n",
    "    # Initialize models with a dynamic alpha value for regularization\n",
    "    models = {\n",
    "        'Linear': LinearRegression(),\n",
    "        'Lasso': Lasso(alpha=alpha_value, max_iter=10000),\n",
    "        'Ridge': Ridge(alpha=alpha_value)\n",
    "    }\n",
    "\n",
    "    # Table to store results\n",
    "    results = []\n",
    "\n",
    "    # Fit models and calculate errors and penalties\n",
    "    for name, model in models.items():\n",
    "        model.fit(x_poly, y_train)\n",
    "        y_poly_pred = model.predict(x_poly)\n",
    "\n",
    "        # Calculating SSE\n",
    "        sse = np.sum((y_train - y_poly_pred) ** 2)\n",
    "\n",
    "        # Calculating L1 and L2 penalties\n",
    "        L1_penalty = np.sum(np.abs(model.coef_))\n",
    "        L2_penalty = np.sum(np.square(model.coef_))\n",
    "\n",
    "        # Number of parameters\n",
    "        num_params = len([param for param in model.coef_ if param != 0])\n",
    "        \n",
    "        results.append({\n",
    "            'Model': name,\n",
    "            'SSE': sse,\n",
    "            'Penalty': 'L1' if name == 'Lasso' else 'L2' if name == 'Ridge' else 'None',\n",
    "            'L1': L1_penalty,\n",
    "            'L2': L2_penalty,\n",
    "            'Number of Parameters' : num_params,\n",
    "        })\n",
    "\n",
    "    # Create a DataFrame to display the results\n",
    "    results_df = pd.DataFrame(results)\n",
    "    display(results_df.set_index('Model'))  # Use display to show the DataFrame in the interactive output\n",
    "\n",
    "# Create interactive widgets for the degree of the polynomial and alpha value\n",
    "degree_slider = IntSlider(min=1, max=11, step=1, value=3, description=\"Polynomial Degree\",\n",
    "                          style={'description_width': 'initial'},\n",
    "                          layout={'width': '500px'})\n",
    "lambda_slider = FloatSlider(min=0.01, max=1.0, step=0.01, value=0.1, description=\"Lambda\",\n",
    "                           style={'description_width': 'initial'},\n",
    "                           layout={'width': '500px'})\n",
    "\n",
    "interactive_plot = interactive(compare_models, degree=degree_slider, alpha_value=lambda_slider)\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets evaluate how this model fits the test data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SET HYPERPARAMETERS HERE:\n",
    "initial_model_type = 'Standard MLR' # Options are 'Standard MLR', 'Lasso', 'Ridge'\n",
    "initial_degree = 1\n",
    "initial_lambda = 0.1\n",
    "# ------------------------------------\n",
    "\n",
    "# Standardize the new data\n",
    "scaler_x = StandardScaler()\n",
    "x_standardized_new = scaler_x.fit_transform(x_test[:, np.newaxis]).flatten()\n",
    "scaler_y = StandardScaler()\n",
    "y_standardized_new = scaler_y.fit_transform(y_test[:, np.newaxis]).flatten()\n",
    "\n",
    "def plot_model_predictions(model_type, degree, alpha):\n",
    "    polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)\n",
    "    x_poly_test = polynomial_features.fit_transform(x_test[:, np.newaxis])\n",
    "    \n",
    "    if model_type == 'Standard MLR':\n",
    "        model = LinearRegression()\n",
    "    elif model_type == 'Lasso':\n",
    "        model = Lasso(alpha=alpha, max_iter=10000)\n",
    "    elif model_type == 'Ridge':\n",
    "        model = Ridge(alpha=alpha)\n",
    "        \n",
    "    model.fit(x_poly_test, y_test)\n",
    "    y_poly_pred_new = model.predict(x_poly_test)\n",
    "\n",
    "    plt.figure(figsize=(10, 6))\n",
    "    plt.scatter(x_test, y_test, color='blue', label='Actual Data')\n",
    "    plt.plot(sorted(x_test), sorted(y_poly_pred_new), color='red', label='Model Prediction')\n",
    "    plt.title(f'Model Type: {model_type}, Degree: {degree}, Lambda: {alpha}\\nR^2: {r2_score(y_standardized_new, y_poly_pred_new):.3f}')\n",
    "    plt.xlabel('Standardized X')\n",
    "    plt.ylabel('Standardized Y')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "# Interactive widgets setup\n",
    "model_dropdown = Dropdown(options=['Standard MLR', 'Lasso', 'Ridge'], value=initial_model_type, description='Model Type:')\n",
    "degree_slider = IntSlider(min=1, max=11, value=initial_degree, description='Polynomial Degree', continuous_update=False)\n",
    "lambda_slider = FloatSlider(min=0.01, max=1.0, value=initial_lambda, step=0.01, description='Lambda', continuous_update=False)\n",
    "\n",
    "interactive_plot = interactive(plot_model_predictions, model_type=model_dropdown, degree=degree_slider, alpha=lambda_slider)\n",
    "display(interactive_plot)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ASSIGNMENT QUESTION:\n",
    "- Please complete questions 16, 17 and 18 in the assignment."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}