<jupyter_start><jupyter_text>Lecture 7: Bias-Variance Tradeoff, Regularization  04/10/2021 Hosted by and maintained by the [Statistics Undergraduate Students Association (SUSA)](https://susa.berkeley.edu). Authored by [<NAME>](mailto:<EMAIL>), [<NAME>](mailto:<EMAIL>), [<NAME>](mailto:<EMAIL>), [<NAME>](mailto:<EMAIL>), and [<NAME>](mailto:<EMAIL>). Table Of Contents* [Recap](recap)* [Bias-Variance Tradeoff](bv-tradeoff)    * [Bias](bias)    * [Variance](variance)    * [The Tradeoff](the-tradeoff)    * [Polynomial Regression](polynomial-regression)* [Regularization](regularization)    * [Ridge](ridge)    * [LASSO](lasso)    * [Visualizing Ridge and Lasso](visualizing-ridge-and-lasso)    * [Regularization and Bias Variance](regularization-and-bias-variance)    * [Lambda](lambda)    * [Validation on Lambda](validation-on-lambda)* [Exercises](exercises)<jupyter_code>import matplotlib.pyplot as plt
import random
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from plotting import overfittingDemo, plot_multiple_linear_regression, ridgeRegularizationDemo, lassoRegularizationDemo
from scipy.optimize import curve_fit
from sklearn.metrics import mean_squared_error
import warnings
warnings.filterwarnings('ignore')

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Recap ![alt text](fit_graphs.png "Fit Graphs") High bias corresponds to underfitting. If we look at the first model, the points seem to follow some sort of curve, but our predictor is linear and therefore, unable to capture all the points. In this case, we have chosen a model which is not complex enough to accurately capture all the information from our data set. If we look at the last model, the predictor is now overly complex because it adjusts based on every point in order to get as close to every data point as possible. In this case, the model changes too much based on small fluctuations caused by insignificant details in the data. This model is fitting more to noise than signal.   Bias-Variance Tradeoff Today we'll perform **model evaluation**, where we'll judge how our linear regression models actually perform. Last week, we talked about **loss functions**, which describes a numerical value for how far your model is from the true values.$$\text{Mean Squared Error: } \frac{1}{n}\sum_{i=1}^n \left(y_i - f(x_i)\right)^2$$In this loss function, $y_i$ is a scalar, and $x_i$ is a $p$ dimensional vector, because there are $p$ features. This loss is called **mean squared error**, or **MSE**.Now, we'll talk about other ways to evaluate a model. **First, let's define some terms.** We can say that everything in the universe can be described with the following equation:$$y = h(x) + \epsilon$$- $y$ is the quantity you are trying to model or approximate- $x$ are the features (independent variables)- $h$ is the **true model** for $y$ in terms of $x$- $\epsilon$ represents **noise**, a random variable which has mean zeroLet $f$ be your model for $y$ in terms of $x$.  BiasWhen evaluating a model, the most intuitive first step is to look at how well the model performs. For classification, this may be the percentage of data points correctly classified, or for regression it may be how close the predicted values are to actual. The **bias** of a model is *a measure of how close our prediction is to the actual value on average from an average model*. Note that bias is **not** a measure of a single model, it encapuslates the scenario in which we collect many datasets, create models for each dataset, and average the error over all of models. Bias is not a measure of error for a single model, but a more abstract concept describing the average error over all errors. A low value for the bias of a model describes that on average, our predictions are similar to the actual values.  VarianceThe **variance** of  a model relates to the variance of the distribution of all models. In the previous section about bias, we envisoned the scenario of collecting many datasets, creating models for each dataset, and averaging the error overall the datasets. Instead, the variance of a model describes the variance in prediction. While we might be able to predict a value very well on average, if the variance of predictions is very high this may not be very helpful, as when we train a model we only have one such instance, and a high model variance tells us little about the true nature of the predictions. A low variance describes that our model will not predict very different values for different datasets.**We can take a look at how bias and variance differ and can be explained in a dataset with the following diagram:** ![alt text](BiasVariance.jpg "Bias Variance Visualization")Image from http://scott.fortmann-roe.com/docs/BiasVariance.html The image describes what bias and variance are in a more simplified example. Consider that we would like to create a model that selects a point close to the center. The models on the top row have low bias, meaning the center of the cluster is close to the red dot on the target. The models on the left column have low variance, the clusters are quite tight, meaning our predictions are close together.**Question: What do the blue dots represent? What about the bullseye?****Question: What is the order of best scenarios?**  The Tradeoff We are trying to minimize **expected error**, or the average **MSE** over all datasets. It turns out (with some advanced probability gymnastics), that:$$\text{Mean Squared Error} = \text{Noise Variance} + \text{Bias}^2 + \text{Variance}$$Note that $\text{Noise Variance}$ is constant: we assume there is some noise, and $\text{moise variance}$ is simply a value that describes how noisy your dataset will be on average. This is often also called "irreducible noise", as it is literally irreducible, we cannot avoid this.  Furthermore, notice that the equation above is the sum of (squared) bias and variance. Thus there is a literal tradeoff between these two, since decreasing one increases the other. This defines what is known as the **bias variance tradeoff**. ![alt text](BiasVarianceTradeoff.png "Bias Variance Tradeoff")Image from http://scott.fortmann-roe.com/docs/BiasVariance.html **Why does this happen?**At some point, as we decrease **bias**, instead of getting closer to the **true model** $h$, we go past and try to fit to the $\epsilon$ (noise) that is part of our current dataset. This is equivalent to making our model more noisy, or **overfit** on our dataset, which means that over all datasets, it has more **variance**. **Questions for understanding**:> 1. Where does underfitting and overfitting lie in the graph above? How do they relate to bias and variance?> 2. Why can't we usually just make a bunch of models with low bias and high variance and average them?> 3. Why is low variance important in models?  Polynomial RegressionLet's revisit the polynomial problem that we have discussed.In this case, if our model has degree $d$, we have $d + 1$ features: $x = [x^0, x^1, ..., x^d]$. Now, we have a linear model with $d + 1$ features:$$\hat{f}(x) = \sum_{i=0}^{d} a_i x^i$$ Model complexity in this case is the degree of the polynomial. As we saw last week, as $d$ increases, model complexity increases. The model gets better, but then gets erratic. This directly corresponds to the bias-variance graph above.<jupyter_code>overfittingDemo()<jupyter_output><empty_output><jupyter_text>As we saw from last time, the best model was a degree 3 model.  RegularizationWe talked about **validation** as a means of combating overfitting. However, this is not the only method to combat overfitting. Another method to do so is to add *regularization* terms to our loss function. **Regularization** basically penalizes complexity in our models. This allows us to add explanatory variables to our model without worrying as much about overfitting. Here's what our ordinary least squares model looks like with a regularization term:$$\hat{\boldsymbol{\theta}} = \arg\!\min_\theta \sum_{i=1}^n (y_i - f_\boldsymbol{\theta}(x_i))^2 + \lambda R(\boldsymbol{\theta})$$We've written the model a little differently here, but the first term is the same as the ordinary least squares regression model you learned last week. This time it's just generalized to any function of $x$ where $\theta$ is a list of parameters, or weights on our explanatory variables, such as coefficients to a polynomial. We're minimizing a loss function to find the best coefficients for our model. The second term is the **regularization** term. The $\lambda$ parameter in front of it dictates how much we care about our regularization term – the higher $\lambda$ is, the more we penalize large weights, and the more the regularization makes our weights deviate from OLS. **Question**: What happens when $\lambda = 0$? So, what is $R(\theta)$ in the equation? There are a variety of different regularization functions that could take its place, but today, we'll just talk about the two most common types of functions: **ridge regression** and **LASSO regression**.$$\begin{align}\text{ Ridge (L2 Norm)}: &\  R(\boldsymbol{\theta}) = \|\theta\|_2^2 = \sum_{i=1}^p \theta_i^2\\\text{ LASSO (L1 Norm)}: &\ R(\boldsymbol{\theta}) = \|\theta\|_1=\sum_{i=1}^p \lvert \theta_i\rvert\end{align}$$  Ridge $$\hat{\boldsymbol{\theta}} = \arg\!\min_\theta \sum_{i=1}^n (y_i - f_\boldsymbol{\theta}(x_i))^2 + \lambda \|\theta\|_2^2$$In **ridge** regression, the regularization function is the sum of squared weights. One nice thing about ridge regression is that there is always a unique, mathematical solution to minimizing the loss of that term. The solution involves some linear algebra, which we won't get into in this notebook, but the existence of this formula makes this minimization computationally easy to solve!$$\hat{\boldsymbol{\theta}} = \left(\boldsymbol{X}^T \boldsymbol{X} + \lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{Y}$$If you recall, the solution to linear regression was of the form:$$\hat{\boldsymbol{\theta}} = \left(\boldsymbol{X}^T \boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{Y}$$And we said that the $\boldsymbol{X}^T\boldsymbol{X}$ isn't always invertible. **What about $\boldsymbol{X}^T \boldsymbol{X} + \lambda\boldsymbol{I}$?**Turns out, this is always invertible! If you are familiar with linear algebra, this is equivalent to adding $\lambda$ to all the eigenvalues of $X^TX$.  **To see this in practice**, we'll first create a regular linear regression model, and compare how it does against models using regularization on the `mpg` dataset we used from last week! We'll be constructing models of `displacement` vs. `mpg`, and seeing the difference from there! First, let's construct the `mpg_train` dataset!<jupyter_code>mpg = pd.read_csv("mpg.csv", index_col='name')# load mpg dataset
mpg = mpg.loc[mpg["horsepower"] != '?'].astype(float) # remove columns with missing horsepower values
mpg_train, mpg_test = train_test_split(mpg, test_size = .2, random_state = 0) # split into training set and test set
mpg_train, mpg_validation = train_test_split(mpg_train, test_size = .5, random_state = 0)
mpg_train.head()<jupyter_output><empty_output><jupyter_text>**Exercise:** Now, let's create a regular linear regression model using the same process we've learned before (fitting, predicting, finding the loss)!<jupyter_code>from sklearn.linear_model import LinearRegression

x_train = np.vander(mpg_train["displacement"], 13)
y_train = mpg_train[["mpg"]]

x_validation = np.vander(mpg_validation["displacement"], 13)
y_validation = mpg_validation[["mpg"]]

# instantiate your model
linear_model = LinearRegression()

# fit the model
linear_model.fit(x_train, y_train)
# make predictions on validation set
linear_prediction = linear_model.predict(x_validation)
# find error
linear_loss = mean_squared_error(linear_prediction, y_validation)

print("Root Mean Squared Error of linear model: {:.2f}".format(linear_loss))<jupyter_output>Root Mean Squared Error of linear model: 25.09<jupyter_text>**Exercise:** Using what you did above as reference, do the same using a Ridge regression model!<jupyter_code>from sklearn.linear_model import Ridge

ridge_model = Ridge()
ridge_model.fit(x_train, y_train)
ridge_prediction = ridge_model.predict(x_validation)
ridge_loss = mean_squared_error(ridge_prediction, y_validation) # mean squared error of ridge model

print("Root Mean Squared Error of Linear Model: {:.2f}".format(linear_loss))
print("Root Mean Squared Error of Ridge Model: {:.2f}".format(ridge_loss))<jupyter_output>Root Mean Squared Error of Linear Model: 25.09
Root Mean Squared Error of Ridge Model: 20.11<jupyter_text>LASSO$$\hat{\boldsymbol{\theta}} = \arg\!\min_\theta \sum_{i=1}^n (y_i - f_\boldsymbol{\theta}(x_i))^2 + \lambda \|\theta\|_1$$In **LASSO** regression, the regularization function is **the sum of absolute values of the weights**. One key thing to note about **LASSO** is that it is **sparsity inducing**, meaning it forces weights to be zero values rather than really small values (which can happen in **Ridge Regression**), leaving you with fewer explanatory variables in the resulting model! Unlike Ridge Regression, LASSO doesn't necessarily have a unique solution that can be solved for with linear algebra, so there's no formula that determines what the optimal weights should be.<jupyter_code>from sklearn.linear_model import Lasso

lasso_model = Lasso()
lasso_model.fit(x_train, y_train)
lasso_prediction = lasso_model.predict(x_validation)
lasso_loss = mean_squared_error(lasso_prediction, y_validation) # mean squared error of lasso model

print("Root Mean Squared Error of Linear Model: {:.2f}".format(linear_loss))
print("Root Mean Squared Error of LASSO Model: {:.2f}".format(lasso_loss))<jupyter_output>Root Mean Squared Error of Linear Model: 25.09
Root Mean Squared Error of LASSO Model: 21.44<jupyter_text>As we can see, both **Ridge Regression and LASSO Regression** minimized the loss of our linear regression models, so maybe penalizing some features allowed us to prevent overfitting on our dataset!  Visualizing Ridge and LASSOWe went through a lot about **ridge** and **LASSO**, but we didn't really get into what they look like for understanding! And so, here are some visualizations that might help build the intution behind some of the characteristics of these two regularization methods. Another way to describe the modified minimization function above is that it's the same loss function as before, with the *additional constraint* that $R(\boldsymbol{\theta}) \leq t$. Now, $t$ is related to $\lambda$ but the exact relationship between the two parameters depends on your data. Regardless, let's take a look at what this means in the two-dimensional case. For ridge,$$\theta_0^2 + \theta_1^2 = t$$Lasso is of the form $$\left|\theta_0\right| + \left|\theta_1\right| =t$$  Norm BallsLet's take at another visualization that may help build some intuition behind how both of these regularization methods work!Image from https://upload.wikimedia.org/wikipedia/commons/f/f8/L1_and_L2_balls.svg.Image from https://towardsdatascience.com/regression-analysis-lasso-ridge-and-elastic-net-9e65dc61d6d3.The rhombus and circle as a visualization of the regularization term, while the blue circles are the topological curves/level sets representing the loss function based on the weights. You want to minimize the sum of these, which means you want to minimize each of those. The point that minimizes the sum is the minimum point at which they intersect.**Question**: Based on these visualizations, could you explain why LASSO is sparsity-inducing?Turns out that the $L2-norm$ is always some sort of smooth surface, from a circle in 2D to a sphere in 3D. On the other hand, LASSO always has sharp corners. In higher dimensions, it forms an octahedron. This is exactly the feature that makes it sparsiy-inducing. As you might imagine, just as humans are more likely to bump into sharp corners than smooth surfaces, the loss term is also most likely to intersect the $L2-norm$ at one of the corners.  Regularization and Bias VarianceAs we mentioned earlier, **bias** is the average linear squares loss term across multiple models of the same family (e.g. same degree polynomial) trained on separate datasets. **Variance** is the average variance of the weight vectors (coefficients) on your features. Without the regularization term, we’re just minimizing bias; the regularization term means we won’t get the lowest possible bias, but we’re exchanging that for some lower variance so that our model does better at generalizing to data points outside of our training data.  Lambda We said that $\lambda$ is how much we care about the regularization term, but what does that look like? Let's return to the polynomial example from last week, and see what the resulting models look like with different values of $\lambda$ given a degree 8 polynomial.<jupyter_code>ridgeRegularizationDemo([0, 0.5, 1.0, 5.0, 10.0], 8)<jupyter_output><empty_output><jupyter_text>From the diagram above, it's difficult to determine which lambda value help fit our model the closest to the true data distribution. So, **how do we know what to use for $\lambda$ (or `alpha` in the `sklearn.linear_model` constructors)?**That's right, let's use the process of **validation** here! In this case, we'd be finding the value for lambda that **minimizes the loss for ridge regression, and then the one that minimizes the loss for LASSO regression**!  Validation on LambdaLet's try to find the best $\lambda$ for the degree 20 polynomial on `displacement` from above.<jupyter_code>lambdas = np.arange(0, 200) # create a list of potential lambda values

# create a list containing the corresponding mean_squared_error for each lambda usinb both ridge and lasso regression
ridge_errors = [] 
lasso_errors = []

for l in lambdas:
    ridge_model = Ridge(l)
    ridge_model.fit(x_train, y_train)
    ridge_predictions = ridge_model.predict(x_validation)
    ridge_errors.append(mean_squared_error(ridge_predictions, y_validation))
    
    
    lasso_model = Lasso(l)
    lasso_model.fit(x_train, y_train)
    lasso_predictions = lasso_model.predict(x_validation)
    lasso_errors.append(mean_squared_error(lasso_predictions, y_validation))

answer = ridge_errors.index(min(ridge_errors)), lasso_errors.index(min(lasso_errors))
answer<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08642779587404995, 'fraction_numerical': 0.0061346362649294245, 'mean_word_length': 5.272046305754171, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 43, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '98221', 'license': 'MIT', 'n_tokens_mistral': 5147, 'n_tokens_neox': 4825, 'n_words': 2399, 'path': 'Lectures/Lecture7_BiasVariance/Solutions.ipynb', 'repo_name': 'Mphan6139/Spring-2021-Career-Exploration', 'size': '18516'}