# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Videos and Exercises for Session 11: Regression and Regularization
#
# In this combined teaching module and exercise set, you will learn about linear regression models in a machine learning perspective. We will see how overfitting can arise and how we can tackle it with a modification of the linear regression model.
#
# The structure of this notebook is as follows:
# 1. Linear Regression Mechanics
# 2. Overfitting and Underfitting in Linear Regression
#     - Exploring Overfitting in Linear Regression
#     - A Cure for Overfitting in Linear Regression
# 3. Modelling Houseprices (Exercise)

# ## Packages
# First, we need to import our standard stuff. Notice that we are not interested in seeing the convergence warning in scikit-learn, so we suppress them for now.

# +
import warnings
from sklearn.exceptions import ConvergenceWarning
warnings.filterwarnings(action='ignore', category=ConvergenceWarning)

import matplotlib.pyplot as plt
import numpy as np 
import pandas as pd 
import seaborn as sns

# %matplotlib inline
# -

# # Part 1: Linear Regression Mechanics
# ## Implementing and evaluating the gradient decent 
#  
# Normally we use OLS to estimate linear regression models, but this is only way of solving the problem of minimizing the least squares problem (that minimizes the sum of squared errors). In the video below we show how to implement gradient descent below and compare it along with other approximate solutions to OLS. 
#
# You may find PML pp. 310-312, 319-324 useful as background reading.

from IPython.display import YouTubeVideo
YouTubeVideo('pl3ep6qRaZw', width=640, height=360)

# Before continuing, we want to provide some brief intuition for how taking steps in the opposite direction of the gradient of a loss function helps us find the weights that minimize the loss. In the video, encountered the following derivation of how weights are updated with gradient descent in a regression type problem: 
#
# \begin{align}\frac{\partial SSE}{\partial\hat{\textbf{w}}} & =-\textbf{X}^{T}\textbf{e}\qquad\text{(the gradient)}\\
# \Delta\hat{\textbf{w}} & =-\eta\cdot\frac{\partial SSE}{\partial\hat{\textbf{w}}}\qquad\text{(gradient descent)}\\
#  & =\eta\cdot\textbf{X}^{T}\textbf{e}\\
#  & =\eta\cdot\textbf{X}^{T}(\textbf{y}-\textbf{X}\hat{\textbf{w}})
# \end{align}
#
# You may ask: Why do we take steps in the opposite direction of the gradient? Consider the illustration of a one-dimensional problem below. Here $w$ is our weight, and $J(w)$ is the loss function that we want to minimize by choosing an appropriate weight. In the example, we started out by guessing a too _high_ value for the weight. At this point, the loss function is _increasing_ in the size of the weight. So what do we do? We take a step in the opposite direction of the gradient and _decrease_ the size of the weight. If the gradient is steep, we take a big step (we know that we are probably relatively far away from the solution), and if the gradient is relatively flat, we take a small step (we might be close to the solution).

# <center><img src='https://i.stack.imgur.com/yk1mk.png' width="600"></center>

# We continue straight to an exercise where you are to implement a new estimator that we code up from scratch. We solve the numerical optimization using the gradient decent algorithm. This will be very similar to what we just saw in the video, but we will pay a bit more attention to each step in the process.
#
# Using our algorithm, we will fit it to some data, and compare our own solution to the standard solution from `sklearn`

# > **Ex. 11.1.0**: Import the dataset `tips` from the `seaborn`.
#
#
# *Hint*: use the `load_dataset` method in seaborn

# + nbgrader={"grade": false, "grade_id": "cell-2dd56f36f76bce57", "locked": false, "schema_version": 3, "solution": true, "task": false}
### BEGIN SOLUTION
tips = sns.load_dataset("tips")
### END SOLUTION
# -

# > **Ex. 11.1.1**: Convert non-numeric variables to dummy variables for each category (remember to leave one column out for each catagorical variable, so you have a reference). Restructure the data so we get a dataset `y` containing the variable tip, and a dataset `X` containing the 
# features. 
#
# > *Hint*: You might want to use the `get_dummies` method in pandas, with the `drop_first = True` parameter. 

# + nbgrader={"grade": false, "grade_id": "cell-2e42eb4f59160bed", "locked": false, "schema_version": 3, "solution": true, "task": false}
### BEGIN SOLUTION
tips_num = pd.get_dummies(tips, drop_first=True)

X = tips_num.drop('tip', axis = 1)
y = tips_num['tip']
### END SOLUTION
# -

# > **Ex. 11.1.2**: Divide the features and target into test and train data. Make the split 50 pct. of each. The split data should be called `X_train`, `X_test`, `y_train`, `y_test`.
#
# > *Hint*: You may use `train_test_split` in `sklearn.model_selection`.

# + nbgrader={"grade": false, "grade_id": "cell-ba197171f1b2bfef", "locked": false, "schema_version": 3, "solution": true, "task": false}
### BEGIN SOLUTION
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=.5, random_state=161193)
### END SOLUTION
# -

# > **Ex. 11.1.3**: Normalize your features by converting to zero mean and one std. deviation.
#
# > *Hint*: Take a look at `StandardScaler` in `sklearn.preprocessing`. If in doubt about which distribution to scale, you may read [this post](https://stats.stackexchange.com/questions/174823/how-to-apply-standardization-normalization-to-train-and-testset-if-prediction-i).

# + nbgrader={"grade": false, "grade_id": "cell-8ab591d5927be1d5", "locked": false, "schema_version": 3, "solution": true, "task": false}
### BEGIN SOLUTION
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
X_train = sc.fit_transform(X_train)
X_test = sc.transform(X_test)
### END SOLUTION
# -

# > **Ex. 11.1.4**: Make a function called `compute_error` to compute the prediction errors given input target `y_`, input features `X_` and input weights `w_`. You should use matrix multiplication.
# >
# > *Hint:* You can use the net-input fct. from yesterday.
#
#

# + nbgrader={"grade": false, "grade_id": "cell-a70101715bbbb443", "locked": false, "schema_version": 3, "solution": true, "task": false}
### BEGIN SOLUTION
def net_input(X_, w_):    
    ''' Computes the matrix product between X and w. Note that
    X is assumed not to contain a bias/intercept column.'''
    return np.dot(X_, w_[1:]) + w_[0]   # We have to add w_[0] separately because this is the constant term. We could also have added a constant term (columns of 1's to X_ and multipliced it to all of w_)

def compute_error(y_, X_, w_):
    return y_ - net_input(X_, w_)
### END SOLUTION


# -

# > **Ex. 11.1.5**: Make a function to update the weights given input target `y_`, input features `X_` and input weights `w_` as well as learning rate, $\eta$, i.e. greek `eta`. You should use matrix multiplication.

# +
# INCLUDED IN ASSIGNMENT 2
# -

# > **Ex. 11.1.6**: Use the code below to initialize weights `w` at zero given feature set `X`. Notice how we include an extra weight that includes the bias term. Set the learning rate `eta` to 0.001. Make a loop with 50 iterations where you iteratively apply your weight updating function. 
#
# >```python
# w = np.zeros(1+X_train.shape[1])
# ```

# +
# INCLUDED IN ASSIGNMENT 2
# -

# > **Ex. 11.1.7**: Make a function to compute the mean squared error. Alter the loop so it makes 100 iterations and computes the MSE for test and train after each iteration, plot these in one figure. 
#
# > Hint: You can use the following code to check that your model works:
# >```python
# from sklearn.linear_model import LinearRegression
# reg = LinearRegression()
# reg.fit(X_train, y_train)
# assert((w[1:] - reg.coef_).sum() < 0.01)
# ```

# + nbgrader={"grade": false, "grade_id": "cell-5cabc75ac6152434", "locked": false, "schema_version": 3, "solution": true, "task": false}
### BEGIN SOLUTION
def MSE(y_,X_,w_):
    return (compute_error(y_, X_, w_)**2).mean()
    
w = np.zeros(1+X_train.shape[1])

MSE_train = [MSE(y_train, X_train, w)]
MSE_test = [MSE(y_test, X_test, w)]

for i in range(100):
    w = update_weights(y_train, X_train, w)
    MSE_train.append(MSE(y_train, X_train, w))
    MSE_test.append(MSE(y_test, X_test, w))
### END SOLUTION

pd.Series(MSE_train).plot()
pd.Series(MSE_test).plot()
# -

# The following bonus exercises are for those who have completed all other exercises until now and have a deep motivation for learning more.

# > **Ex. 11.1.8 (BONUS)**: Implement your linear regression model as a class.

# > A solution is found on p. 320 in Python for Machine Learning.

# # Part 2: Overfitting and Underfitting in Linear Regression 
#
# ## Exploring Overfitting in Linear Regression
# How does overfitting manifest itself in linear regression? In the video below we simulate what happens as make a better and better taylor approximation, i.e. we estimate a polynomial of higher and higher order. Two issues arise simultaneously - one is related to the number of parameters and the to the size of the parameters. 
#
# You may find PML pp. 334-339 useful as background reading.

YouTubeVideo('HbeTpK-2oeU', width=640, height=360)

# ## A Cure for Overfitting in Linear Regression
#
# How do we fix the two issues of excessively large weights/coefficients and too many spurious solutions? The video below provides a solution by directly incorporating these issues into the optimization problem.
#
# You may find PML pp. 73-76, 123-136, 332-334 useful as background reading.

YouTubeVideo('r6a8WFm9jAI', width=640, height=360)

# Above we tackled overfitting, but what about ***underfitting***? The video below shows how to address underfitting and also zooms in on some important details about regularization.
#
# You may find PML pp. 73-76, 123-136, 332-334 useful as background reading.

YouTubeVideo('IWBtYT1KI_Q', width=640, height=360)

# > **Ex. 11.2.1 (BONUS)**: Is it possible to add a penalty to our linear model above and solve this Lasso model with gradient descent? Is there a simple fix?
# >
# > *Hint:* Gradient descent essentially relies on a differentiable loss function (read more [here](https://stats.stackexchange.com/questions/177800/why-proximal-gradient-descent-instead-of-plain-subgradient-methods-for-lasso))

# + nbgrader={"grade": false, "grade_id": "cell-6aa146e9c530b2b4", "locked": false, "schema_version": 3, "solution": true, "task": false}
### BEGIN SOLUTION

ANSWER: No, we cannot exactly solve for the Lasso with gradient descent. However, we can make an approximate solution which is pretty close and quite intuitive - see good explanation here: https://stats.stackexchange.com/questions/177800/why-proximal-gradient-descent-instead-of-plain-subgradient-methods-for-lasso.

### END SOLUTION
# -

# # Part 3: Modelling Houseprices
# In this example, we will try to predict houseprices using a lot of variable (or features as they are called in Machine Learning). We are going to work with Kaggle's dataset on house prices, see information [here](https://www.kaggle.com/c/house-prices-advanced-regression-techniques). Kaggle is an organization that hosts competitions in building predictive models.

# > **Ex. 11.3.0:** Load the california housing data with scikit-learn using the code below. Now:
# > 1. Inspect *cal_house*. How are the data stored?
# > 2. Create a pandas DataFrame called *X*, using `data`. Name the columns using `feature_names`.
# > 3. Crate a pandas Series called *y* using `target`.
# > 4. Make a train test split of equal size.

# + nbgrader={"grade": false, "grade_id": "cell-5f14e576643ac94c", "locked": false, "schema_version": 3, "solution": true, "task": false}
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split

cal_house = fetch_california_housing()    

### BEGIN SOLUTION
X = pd.DataFrame(data=cal_house['data'], 
                 columns=cal_house['feature_names'])\
             .iloc[:,:-2]
y = cal_house['target']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.5, random_state=1)

X_train.describe()
### END SOLUTION
# -

#
#
# > **Ex.11.3.1**: Generate interactions between all features to third degree (make sure you **exclude** the bias/intercept term). How many variables are there? Will OLS fail? After making interactions, rescale the features to have zero mean, unit std. deviation. Should you use the distribution of the training data to rescale the test data?  
#
# > *Hint 1*: Try importing `PolynomialFeatures` from `sklearn.preprocessing`
#
# > *Hint 2*: If in doubt about which distribution to scale, you may read [this post](https://stats.stackexchange.com/questions/174823/how-to-apply-standardization-normalization-to-train-and-testset-if-prediction-i).

# +
# INCLUDED IN ASSIGNMENT 2
# -

# > **Ex.11.3.2**: Estimate the Lasso model on the rescaled train data set, using values of $\lambda$ in the range from $10^{-4}$ to $10^4$. For each $\lambda$  calculate and save the Root Mean Squared Error (RMSE) for the rescaled test and train data. Take a look at the fitted coefficients for different sizes of $\lambda$. What happens when $\lambda$ increases? Why?
#
# > *Hint 1*: use `logspace` in numpy to create the range.
#
# > *Hint 2*: read about the `coef_` feature [here](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html#sklearn.linear_model.Lasso).
#

# +
# INCLUDED IN ASSIGNMENT 2
# -

# > **Ex.11.3.3**: Make a plot with on the x-axis and the RMSE measures on the y-axis. What happens to RMSE for train and test data as $\lambda$ increases? The x-axis should be log scaled. Which one are we interested in minimizing? 
#
# > Bonus: Can you find the lambda that gives the lowest MSE-test score?

# +
# INCLUDED IN ASSIGNMENT 2

StarcoderdataJupyterScriptsDedupFiltered