Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis

Further Exercises

Exercise 1

We will generate our own dataset for a function $ y(x) $ where $ x \in [0,1] $ and defined by random numbers computed with the uniform distribution. The function $ y $ is a quadratic polynomial in $ x $ with added stochastic noise according to the normal distribution $ \cal {N}(0,1) $. The following simple Python instructions define our $ x $ and $ y $ values (with 100 data points).

x = np.random.rand(100,1)
y = 5*x*x+0.1*np.random.randn(100,1)

Write your own code (following the examples above) for computing the parametrization of the data set fitting a second-order polynomial.
Use thereafter scikit-learn (see again the examples in the regression slides) and compare with your own code.
Using scikit-learn, compute also the mean square error, a risk metric corresponding to the expected value of the squared (quadratic) error defined as

$$ MSE(\hat{y},\hat{\tilde{y}}) = \frac{1}{n} \sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2, $$ and the $ R^2 $ score function. If $ \tilde{\hat{y}}_i $ is the predicted value of the $ i-th $ sample and $ y_i $ is the corresponding true value, then the score $ R^2 $ is defined as $$ R^2(\hat{y}, \tilde{\hat{y}}) = 1 - \frac{\sum_{i=0}^{n - 1} (y_i - \tilde{y}_i)^2}{\sum_{i=0}^{n - 1} (y_i - \bar{y})^2}, $$ where we have defined the mean value of $ \hat{y} $ as $$ \bar{y} = \frac{1}{n} \sum_{i=0}^{n - 1} y_i. $$

You can use the functionality included in scikit-learn. If you feel for it, you can use your own program and define functions which compute the above two functions. Discuss the meaning of these results. Try also to vary the coefficient in front of the added stochastic noise term and discuss the quality of the fits.

Exercise 2, variance of the parameters $ \beta $ in linear regression

Show that the variance of the parameters $ \beta $ in the linear regression method (chapter 3, equation (3.8) of Trevor Hastie, Robert Tibshirani, Jerome H. Friedman, The Elements of Statistical Learning, Springer) is given as $$ \mathrm{Var}(\hat{\beta}) = \left(\hat{X}^T\hat{X}\right)^{-1}\sigma^2, $$ with $$ \sigma^2 = \frac{1}{N-p-1}\sum_{i=1}^{N} (y_i-\tilde{y}_i)^2, $$ where we have assumed that we fit a function of degree $ p-1 $ (for example a polynomial in $ x $).

Exercise 3

This exercise is a continuation of exercise 1. We will use the same function to generate our data set, still staying with a simple function $ y(x) $ which we want to fit using linear regression, but now extending the analysis to include the Ridge and the Lasso regression methods. You can use the code under the Regression as an example on how to use the Ridge and the Lasso methods.

We will thus again generate our own dataset for a function $ y(x) $ where $ x \in [0,1] $ and defined by random numbers computed with the uniform distribution. The function $ y $ is a quadratic polynomial in $ x $ with added stochastic noise according to the normal distribution $ \cal{N}(0,1) $.

The following simple Python instructions define our $ x $ and $ y $ values (with 100 data points).

x = np.random.rand(100,1)
y = 5*x*x+0.1*np.random.randn(100,1)

Write your own code for the Ridge method and compute the parametrization for different values of $ \lambda $. Compare and analyze your results with those from exercise 1. Study the dependence on $ \lambda $ while also varying the strength of the noise in your expression for $ y(x) $.
Repeat the above but using the functionality of scikit-learn. Compare your code with the results from scikit-learn. Remember to run with the same random numbers for generating $ x $ and $ y $.
Our next step is to study the variance of the parameters $ \beta_1 $ and $ \beta_2 $ (assuming that we are parametrizing our function with a second-order polynomial. We will use standard linear regression and the Ridge regression. You can now opt for either writing your own function that calculates the variance of these paramaters (recall that this is equal to the diagonal elements of the matrix $ (\hat{X}^T\hat{X})+\lambda\hat{I})^{-1} $) or use the functionality of scikit-learn and compute their variances. Discuss the results of these variances as functions
Repeat the previous step but add now the Lasso method. Discuss your results and compare with standard regression and the Ridge regression results.
Try to implement the cross-validation as well.
Finally, using scikit-learn or your own code, compute also the mean square error, a risk metric corresponding to the expected value of the squared (quadratic) error defined as

Exercise 4

We will study how to fit polynomials to a specific two-dimensional function called Franke's function. This is a function which has been widely used when testing various interpolation and fitting algorithms. Furthermore, after having established the model and the method, we will employ resamling techniques such as the cross-validation and/or the bootstrap methods, in order to perform a proper assessment of our models.

The Franke function, which is a weighted sum of four exponentials reads as follows $$ \begin{align*} f(x,y) &= \frac{3}{4}\exp{\left(-\frac{(9x-2)^2}{4} - \frac{(9y-2)^2}{4}\right)}+\frac{3}{4}\exp{\left(-\frac{(9x+1)^2}{49}- \frac{(9y+1)}{10}\right)} \\ &+\frac{1}{2}\exp{\left(-\frac{(9x-7)^2}{4} - \frac{(9y-3)^2}{4}\right)} -\frac{1}{5}\exp{\left(-(9x-4)^2 - (9y-7)^2\right) }. \end{align*} $$

The function will be defined for $ x,y\in [0,1] $. Our first step will be to perform an OLS regression analysis of this function, trying out a polynomial fit with an $ x $ and $ y $ dependence of the form $ [x, y, x^2, y^2, xy, \dots] $. We will also include cross-validation and bootstrap as resampling techniques. As in homeworks 1 and 2, we can use a uniform distribution to set up the arrays of values for $ x $ and $ y $, or as in the example below just a fix values for $ x $ and $ y $ with a given step size. In this case we will have two predictors and need to fit a function (for example a polynomial) of $ x $ and $ y $. Thereafter we will repeat much of the same procedure using the the Ridge and Lasso regression methods, introducing thus a dependence on the bias (penalty) $ \lambda $.

The Python fucntion for the Franke function is included here (it performs also a three-dimensional plot of it)

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
from random import random, seed

fig = plt.figure()
ax = fig.gca(projection='3d')

# Make data.
x = np.arange(0, 1, 0.05)
y = np.arange(0, 1, 0.05)
x, y = np.meshgrid(x,y)


def FrankeFunction(x,y):
    term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))
    term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))
    term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))
    term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)
    return term1 + term2 + term3 + term4


z = FrankeFunction(x, y)

# Plot the surface.
surf = ax.plot_surface(x, y, z, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

# Customize the z axis.
ax.set_zlim(-0.10, 1.40)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

# Add a color bar which maps values to colors.
fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()

We will thus again generate our own dataset for a function $ \mathrm{FrankeFunction}(x,y) $ where $ x,y \in [0,1] $ could be defined by random numbers computed with the uniform distribution. The function $ f(x,y) $ is the Franke function. You should explore also the addition an added stochastic noise to this function using the normal distribution $ \cal{N}(0,1) $.

Write your own code (using either a matrix inversion or a singular value decomposition from e.g., numpy ) or use your code from exercises 1 and 3 and perform a standard least square regression analysis using polynomials in $ x $ and $ y $ up to fifth order. Find the confidence intervals of the parameters $ \beta $ by computing their variances, evaluate the Mean Squared error (MSE) $$ MSE(\hat{y},\hat{\tilde{y}}) = \frac{1}{n} \sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2, $$ and the $ R^2 $ score function. If $ \tilde{\hat{y}}_i $ is the predicted value of the $ i-th $ sample and $ y_i $ is the corresponding true value, then the score $ R^2 $ is defined as $$ R^2(\hat{y}, \tilde{\hat{y}}) = 1 - \frac{\sum_{i=0}^{n - 1} (y_i - \tilde{y}_i)^2}{\sum_{i=0}^{n - 1} (y_i - \bar{y})^2}, $$ where we have defined the mean value of $ \hat{y} $ as $$ \bar{y} = \frac{1}{n} \sum_{i=0}^{n - 1} y_i. $$

Perform a resampling of the data where you split the data in training data and test data. Implement the $ k $-fold cross-validation algorithm and/or the bootstrap algorithm and evaluate again the MSE and the $ R^2 $ functions resulting from the test data. Evaluate also the bias and variance of the final models.

Write then your own code for the Ridge method, either using matrix inversion or the singular value decomposition as done for standard OLS. Perform the same analysis as in the previous exercise (for the same polynomials and include resampling techniques) but now for different values of $ \lambda $. Compare and analyze your results with those obtained with standard OLS. Study the dependence on $ \lambda $ while also varying eventually the strength of the noise in your expression for $ \mathrm{FrankeFunction}(x,y) $.

Then perform the same studies but now with Lasso regression. Use the functionalities of scikit-learn. Give a critical discussion of the three methods and a judgement of which model fits the data best.

Further Exercises

Exercise 1

Exercise 2, variance of the parameters \( \beta \) in linear regression

Exercise 3

Exercise 4