🏋️ Exercises - Statistical Learning with Python

Exercise 1: Linear Regression Recap¶

Please implement regression models using statsmodels, sklearn, and numpy. Which degree of polynomial is most suited for the synthetic data provided in the code cell below?

Plot the resulting regression models.

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

# Generate synthetic data
n_samples = 100
X = np.linspace(-2, 2, n_samples).reshape(-1, 1)
y = X**2 + np.random.normal(scale=0.5, size=X.shape)

# TODO: Statsmodels
...

# TODO: Numpy
...

# TODO: Scikit-learn
...

# TODO: Plot the predictions
...

Exercise 2: Bias-variance Tradeoff¶

To get a better understanding of the bias-variance tradeoff, we will fit polynomial regression models to synthetic data from a known function, $y = \sin(x) + \varepsilon$ .

Please perform the following tasks:

Visualize the data. Which model do you think would be optimal?
Split the data into a training set (70%) and a test set (30%). Use random_state=42 so that your results are reproducible.
Fit polynomial regression models for degrees 1 to 12 and calculate the training and test MSE for each of them
Plot the training and test MSE against the polynomial degree
Which degree has the lowest test MSE? Does it match your guess from task 1? Now repeat the split with a few different values for random_state. Does the best degree change? How much does the test MSE vary?

Bonus: Because we know the true function, we can also measure bias and variance directly, just like in the bias-variance session. Simulate many datasets with the same X but new noise, fit a model of each degree to each dataset, and calculate bias² and variance of the predictions.

Hint: You can split the data with the train_test_split() function from sklearn.model_selection.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

# Generate synthetic data
np.random.seed(55)
n_samples = 100
X = np.linspace(0, 2*np.pi, n_samples).reshape(-1, 1)  # Reshape for sklearn
y = np.sin(X) + np.random.normal(scale=0.5, size=X.shape)

# 1. TODO: Visualize the data
...

# 2. TODO: Split the data into training and test sets (use random_state=42)
...

# 3. TODO: Fit polynomial regression models for degrees 1 to 12
degrees = range(1, 13)
...

# 4. TODO: Plot the training and test MSE against the polynomial degree
...

# 5. TODO: Repeat the split with different values for random_state
...

# Bonus: TODO: Estimate bias² and variance by simulating many datasets
...

Exercise 3: Resampling Methods¶

The dataset we are using for the exercise is the California Housing Dataset. It contains 20640 samples and 8 features. In this dataset, we have information regarding the demography (income, population, house occupancy) in the districts, the location of the districts (latitude, longitude), and general information regarding the house in the districts (number of rooms, number of bedrooms, age of the house). Since these statistics are at the granularity of the district, they correspond to averages or medians.

Familiarize yourself with the dataset by exploring the documentation and looking at the data. What are the features and target?

Hints:

print(data.DESCR) prints the documentation of the dataset
df.head() shows the first rows of the DataFrame, and df.describe() gives summary statistics for every column
The names of the features and the target are stored in data.feature_names and data.target_names

from sklearn.datasets import fetch_california_housing

data = fetch_california_housing(as_frame=True)
df = data.frame
...

Let’s have a quick look at the distribution of these features by plotting their histograms.

Hint: DataFrames have a .hist() method that plots a histogram for every column. Use its figsize and bins arguments to make the plot readable.

# TODO: Plot histograms
...

Prepare the data for cross validation. This will require you to have a variable (e.g. X) for the features and a variable for the target (e.g. y).
Set up a k-fold cross validation for a linear regression
- Choose an appropriate k
- Define the model
- Perform cross validation
- Use the mean squared error (MSE) to assess model performance
Visualize which rows end up in the validation fold in each of the 5 rounds, once with KFold(n_splits=5) and once with KFold(n_splits=5, shuffle=True). Then run the cross-validation with both splitters and compare the MSE of the individual folds. What do you notice?

Hints:

For 1: You can achieve this by e.g. creating them from the DataFrame, or by using the return_X_y parameter on the fetch_california_housing() function
For 2: You can use the LinearRegression() model from sklearn. You can further evaluate the model and specify the (negative) MSE as a performance measure in cross_val_score()
For 3: Have a look at the plot_cv_indices() function in the scikit-learn example on visualizing cross-validation behavior as an example. A simpler version: create an array of zeros with one row per fold and one column per observation, set the validation rows of each fold to 1, and plot it with plt.imshow(..., aspect='auto').
The scikit-learn user guide on cross-validation explains all splitters and functions used in this exercise, with examples.

# TODO: Prepare data
...

from sklearn.model_selection import KFold
from sklearn.linear_model import LinearRegression 
from sklearn.model_selection import cross_val_score 

# TODO: Implement CV
...

# 3. TODO: Visualize the folds with and without shuffling and compare the fold MSEs
...

Exercise 4: LOOCV¶

Draw a random subset of 5000 observations (LOOCV on all 20640 observations would take a while). Run LOOCV on this subset and compare the average MSE with a 5-fold CV on the same subset.
Get the minimum and maximum MSE value. Discuss the range!
Plot the MSE values in a histogram (the x-range should be from 0 to 6)
Calculate the median MSE and discuss if it might be a more appropriate measure than the mean

Hint: You can draw a random subset of row indices with np.random.default_rng(42).choice(len(y), size=5000, replace=False).

# TODO: Implement LOOCV
...