diff --git a/doc/Projects/2020/hw2/pdf/hw2.p.tex b/doc/Projects/2020/hw2/pdf/hw2.p.tex new file mode 100644 index 000000000..85004d1d1 --- /dev/null +++ b/doc/Projects/2020/hw2/pdf/hw2.p.tex @@ -0,0 +1,742 @@ +%% +%% Automatically generated file from DocOnce source +%% (https://github.com/hplgit/doconce/) +%% +%% +% #ifdef PTEX2TEX_EXPLANATION +%% +%% The file follows the ptex2tex extended LaTeX format, see +%% ptex2tex: http://code.google.com/p/ptex2tex/ +%% +%% Run +%% ptex2tex myfile +%% or +%% doconce ptex2tex myfile +%% +%% to turn myfile.p.tex into an ordinary LaTeX file myfile.tex. +%% (The ptex2tex program: http://code.google.com/p/ptex2tex) +%% Many preprocess options can be added to ptex2tex or doconce ptex2tex +%% +%% ptex2tex -DMINTED myfile +%% doconce ptex2tex myfile envir=minted +%% +%% ptex2tex will typeset code environments according to a global or local +%% .ptex2tex.cfg configure file. doconce ptex2tex will typeset code +%% according to options on the command line (just type doconce ptex2tex to +%% see examples). If doconce ptex2tex has envir=minted, it enables the +%% minted style without needing -DMINTED. +% #endif + +% #define PREAMBLE + +% #ifdef PREAMBLE +%-------------------- begin preamble ---------------------- + +\documentclass[% +oneside, % oneside: electronic viewing, twoside: printing +final, % draft: marks overfull hboxes, figures with paths +10pt]{article} + +\listfiles % print all files needed to compile this document + +\usepackage{relsize,makeidx,color,setspace,amsmath,amsfonts,amssymb} +\usepackage[table]{xcolor} +\usepackage{bm,ltablex,microtype} + +\usepackage[pdftex]{graphicx} + +\usepackage{ptex2tex} +% #ifdef MINTED +\usepackage{minted} +\usemintedstyle{default} +% #endif + +\usepackage[T1]{fontenc} +%\usepackage[latin1]{inputenc} +\usepackage{ucs} +\usepackage[utf8x]{inputenc} + +\usepackage{lmodern} % Latin Modern fonts derived from Computer Modern + +% Hyperlinks in PDF: +\definecolor{linkcolor}{rgb}{0,0,0.4} +\usepackage{hyperref} +\hypersetup{ + breaklinks=true, + colorlinks=true, + linkcolor=linkcolor, + urlcolor=linkcolor, + citecolor=black, + filecolor=black, + %filecolor=blue, + pdfmenubar=true, + pdftoolbar=true, + bookmarksdepth=3 % Uncomment (and tweak) for PDF bookmarks with more levels than the TOC + } +%\hyperbaseurl{} % hyperlinks are relative to this root + +\setcounter{tocdepth}{2} % levels in table of contents + +% --- fancyhdr package for fancy headers --- +\usepackage{fancyhdr} +\fancyhf{} % sets both header and footer to nothing +\renewcommand{\headrulewidth}{0pt} +\fancyfoot[LE,RO]{\thepage} +% Ensure copyright on titlepage (article style) and chapter pages (book style) +\fancypagestyle{plain}{ + \fancyhf{} + \fancyfoot[C]{{\footnotesize \copyright\ 1999-2020, "Data Analysis and Machine Learning FYS-STK3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html". Released under CC Attribution-NonCommercial 4.0 license}} +% \renewcommand{\footrulewidth}{0mm} + \renewcommand{\headrulewidth}{0mm} +} +% Ensure copyright on titlepages with \thispagestyle{empty} +\fancypagestyle{empty}{ + \fancyhf{} + \fancyfoot[C]{{\footnotesize \copyright\ 1999-2020, "Data Analysis and Machine Learning FYS-STK3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html". Released under CC Attribution-NonCommercial 4.0 license}} + \renewcommand{\footrulewidth}{0mm} + \renewcommand{\headrulewidth}{0mm} +} + +\pagestyle{fancy} + + +% prevent orhpans and widows +\clubpenalty = 10000 +\widowpenalty = 10000 + +\newenvironment{doconceexercise}{}{} +\newcounter{doconceexercisecounter} + + +% ------ header in subexercises ------ +%\newcommand{\subex}[1]{\paragraph{#1}} +%\newcommand{\subex}[1]{\par\vspace{1.7mm}\noindent{\bf #1}\ \ } +\makeatletter +% 1.5ex is the spacing above the header, 0.5em the spacing after subex title +\newcommand\subex{\@startsection{paragraph}{4}{\z@}% + {1.5ex\@plus1ex \@minus.2ex}% + {-0.5em}% + {\normalfont\normalsize\bfseries}} +\makeatother + + +% --- end of standard preamble for documents --- + + +% insert custom LaTeX commands... + +\raggedbottom +\makeindex +\usepackage[totoc]{idxlayout} % for index in the toc +\usepackage[nottoc]{tocbibind} % for references/bibliography in the toc + +%-------------------- end preamble ---------------------- + +\begin{document} + +% matching end for #ifdef PREAMBLE +% #endif + +\newcommand{\exercisesection}[1]{\subsection*{#1}} + + +% ------------------- main content ---------------------- + + + +% ----------------- title ------------------------- + +\thispagestyle{empty} + +\begin{center} +{\LARGE\bf +\begin{spacing}{1.25} +Homework 2, weeks 36 and 37 +\end{spacing} +} +\end{center} + +% ----------------- author(s) ------------------------- + +\begin{center} +{\bf \href{{http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html}}{Data Analysis and Machine Learning FYS-STK3155/FYS4155}} +\end{center} + + \begin{center} +% List of all institutions: +\centerline{{\small Department of Physics, University of Oslo, Norway}} +\end{center} + +% ----------------- end author(s) ------------------------- + +% --- begin date --- +\begin{center} +Sep 1, 2020 +\end{center} +% --- end date --- + +\vspace{1cm} + + +% --- begin exercise --- +\begin{doconceexercise} +\refstepcounter{doconceexercisecounter} + +\exercisesection{Exercise \thedoconceexercisecounter: Adding Ridge and Lasso Regression} + + +This exercise is a continuation of exercise 3 from exercise set 1 (week 35). 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, see the \href{{https://compphysics.github.io/MachineLearning/doc/pub/Regression/html/Regression-bs.html}}{regression slides}). + +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). +\bpycod +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) +\epycod + + +\subex{a)} +Write your own code for the Ridge method (see chapter 3.4 of Hastie \emph{et al.}, equations (3.43) and (3.44)) and compute the parametrization for different values of $\lambda$. Compare and analyze your results with those from exercise 3. Study the dependence on $\lambda$ while also varying the strength of the noise in your expression for $y(x)$. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +The code here allows you to perform your own Ridge calculation and perform calculations for various values of the regularization parameter $\lambda$. This program can easily be extended upon. +\bpycod +import os +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# and then make the prediction +ytildeOLS = X_train @ OLSbeta +print("Training R2 for OLS") +print(R2(y_train,ytildeOLS)) +print("Training MSE for OLS") +print(MSE(y_train,ytildeOLS)) +ypredictOLS = X_test @ OLSbeta +print("Test R2 for OLS") +print(R2(y_test,ypredictOLS)) +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 20 +MSEPredict = np.zeros(nlambdas) +MSETrain = np.zeros(nlambdas) +lambdas = np.logspace(-4, 1, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train + # and then make the prediction + ytildeRidge = X_train @ Ridgebeta + ypredictRidge = X_test @ Ridgebeta + MSEPredict[i] = MSE(y_test,ypredictRidge) + MSETrain[i] = MSE(y_train,ytildeRidge) +# Now plot the resulys +plt.figure() +plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train') +plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() +\epycod + +% --- end solution of exercise --- + +\subex{b)} +Repeat the above but using the functionality of \textbf{Scikit-Learn}. Compare your code with the results from \textbf{Scikit-Learn}. Remember to run with the same random numbers for generating $x$ and $y$. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +To use \textbf{scikit-learn} with Ridge, we simply need to add the relevant function \textbf{Ridge()}, as done in the code here. +\bpycod +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler +import sklearn.linear_model as skl + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# and then make the prediction +ytildeOLS = X_train @ OLSbeta +print("Training R2 for OLS") +print(R2(y_train,ytildeOLS)) +print("Training MSE for OLS") +print(MSE(y_train,ytildeOLS)) +ypredictOLS = X_test @ OLSbeta +print("Test R2 for OLS") +print(R2(y_test,ypredictOLS)) +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 100 +MSEPredict = np.zeros(nlambdas) +MSEPredictSKL = np.zeros(nlambdas) +MSETrain = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + # add ridge + clf_ridge = skl.Ridge(alpha=lmb).fit(X_train, y_train) + yridge = clf_ridge.predict(X_test) + Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train + # and then make the prediction + ytildeRidge = X_train @ Ridgebeta + ypredictRidge = X_test @ Ridgebeta + MSEPredict[i] = MSE(y_test,ypredictRidge) + MSEPredictSKL[i] = MSE(y_test,yridge) + MSETrain[i] = MSE(y_train,ytildeRidge) +#then plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train') +plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test') +plt.plot(np.log10(lambdas), MSEPredictSKL, 'g--', label = 'MSE Ridge sickit-learn Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() +\epycod + +% --- end solution of exercise --- + +\subex{c)} +Our next step is to study the variance of the parameters $\beta_1$ and $\beta_2$ (assuming that we are parameterizing 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 or using \textbf{Scikit-Learn} to find the parameters $\beta$. From your results calculate the variance of these parameters (recall that this is equal to the diagonal elements of the matrix $(\hat{X}^T\hat{X})+\lambda\hat{I})^{-1}$). Discuss the results of these variances as functions of $\lambda$. In particular, try to link your discussion with the discussion in Hastie \emph{et al.} and their figures 3.10 and 3.11. \textbf{Scikit-Learn} may not provide the variance of the parameters $\beta$. This needs to be checked. With your own code you can however do so. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +\bpycod +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler +import sklearn.linear_model as skl + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# The variance is given by the inverse of the matrix X^TX +print(np.linalg.inv(X_train.T @ X_train)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 10 +MSEPredict = np.zeros(nlambdas) +MSEPredictSKL = np.zeros(nlambdas) +MSETrain = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train + print(np.linalg.inv(X_train.T @ X_train+lmb*I)) + + + +\epycod + +% --- end solution of exercise --- + +\subex{d)} +Repeat the previous step but add now the Lasso method, see equation (3.53) of Hastie \emph{et al.}. Discuss your results and compare with standard regression and the Ridge regression results. You can write your own code or use the functionality of \textbf{scikit-learn}. We recommend the latter since we have not yet discussed how to solve the Lasso equations numerically. Also, you do not need to compute the variance of the parameters $\beta$ but you can extract their values and study their behavior as functions of the regularization parameter $\lambda$. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +\bpycod +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler +import sklearn.linear_model as skl +#from sklearn.linear_model import LinearRegression, Ridge, Lasso +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# and then make the prediction +ytildeOLS = X_train @ OLSbeta +print("Training R2 for OLS") +print(R2(y_train,ytildeOLS)) +print("Training MSE for OLS") +print(MSE(y_train,ytildeOLS)) +ypredictOLS = X_test @ OLSbeta +print("Test R2 for OLS") +print(R2(y_test,ypredictOLS)) +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 100 +MSEPredictLasso = np.zeros(nlambdas) +MSEPredictRidge = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + # add ridge + clf_ridge = skl.Ridge(alpha=lmb).fit(X_train, y_train) + clf_lasso = skl.Lasso(alpha=lmb).fit(X_train, y_train) + yridge = clf_ridge.predict(X_test) + ylasso = clf_lasso.predict(X_test) + MSEPredictLasso[i] = MSE(y_test,ylasso) + MSEPredictRidge[i] = MSE(y_test,yridge) +#then plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSEPredictRidge, 'r--', label = 'MSE Ridge Test') +plt.plot(np.log10(lambdas), MSEPredictLasso, 'g--', label = 'MSE Lasso Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() +\epycod + +% --- end solution of exercise --- + +\subex{e)} +Finally, using \textbf{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 +\[ 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. +\] +Discuss these quantities as functions of the variable $\lambda$ in the Ridge and Lasso regression methods. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +These results can all be studied with the codes we have above. These scores are included in the codes above. + +% --- end solution of exercise --- + + +\end{doconceexercise} +% --- end exercise --- + + + + +% --- begin exercise --- +\begin{doconceexercise} +\refstepcounter{doconceexercisecounter} + +\exercisesection{Exercise \thedoconceexercisecounter: Normalizing our data} + + +A much used approach before starting to train the data is to preprocess our +data. Normally the data may need a rescaling and/or may be sensitive +to extreme values. Scaling the data renders our inputs much more +suitable for the algorithms we want to employ. + +\textbf{Scikit-Learn} has several functions which allow us to rescale the +data, normally resulting in much better results in terms of various +accuracy scores. The \textbf{StandardScaler} function in \textbf{Scikit-Learn} +ensures that for each feature/predictor we study the mean value is +zero and the variance is one (every column in the design/feature +matrix). This scaling has the drawback that it does not ensure that +we have a particular maximum or minimum in our data set. Another +function included in \textbf{Scikit-Learn} is the \textbf{MinMaxScaler} which +ensures that all features are exactly between $0$ and $1$. The + + +The \textbf{Normalizer} scales each data +point such that the feature vector has a euclidean length of one. In other words, it +projects a data point on the circle (or sphere in the case of higher dimensions) with a +radius of 1. This means every data point is scaled by a different number (by the +inverse of it’s length). +This normalization is often used when only the direction (or angle) of the data matters, +not the length of the feature vector. + +The \textbf{RobustScaler} works similarly to the StandardScaler in that it +ensures statistical properties for each feature that guarantee that +they are on the same scale. However, the RobustScaler uses the median +and quartiles, instead of mean and variance. This makes the +RobustScaler ignore data points that are very different from the rest +(like measurement errors). These odd data points are also called +outliers, and might often lead to trouble for other scaling +techniques. + + +It also common to split the data in a \textbf{training} set and a \textbf{testing} set. A typical split is to use $80\%$ of the data for training and the rest +for testing. This can be done as follows with our design matrix $\bm{X}$ and data $\bm{y}$ (remember to import \textbf{scikit-learn}) +\bpycod +# split in training and test data +X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.2) +\epycod +Then we can use the standard scaler to scale our data as +\bpycod +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) +\epycod + + +In this exercise we want you to to compute the MSE for the training +data and the test data as function of the complexity of a polynomial, +that is the degree of a given polynomial. We want you also to compute the $R2$ score as function of the complexity of the model for both training data and test data. You should also run the calculation with and without scaling. + +One of +the aims is to reproduce Figure 2.11 of \href{{https://github.com/CompPhysics/MLErasmus/blob/master/doc/Textbooks/elementsstat.pdf}}{Hastie et al}. +We will also use Ridge and Lasso regression. + + +Our data is defined by $x\in [-3,3]$ with a total of for example $100$ data points. +\bpycod +np.random.seed() +n = 100 +maxdegree = 14 +# Make data set. +x = np.linspace(-3, 3, n).reshape(-1, 1) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape) +\epycod +where $y$ is the function we want to fit with a given polynomial. + + +\subex{a)} +Write a first code which sets up a design matrix $X$ defined by a fifth-order polynomial. Scale your data and split it in training and test data. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +\bpycod +import matplotlib.pyplot as plt +import numpy as np +from sklearn.linear_model import LinearRegression, Ridge, Lasso +from sklearn.preprocessing import PolynomialFeatures +from sklearn.model_selection import train_test_split +from sklearn.pipeline import make_pipeline + + +np.random.seed(2018) +n = 50 +maxdegree = 5 +# Make data set. +x = np.linspace(-3, 3, n).reshape(-1, 1) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape) +TestError = np.zeros(maxdegree) +TrainError = np.zeros(maxdegree) +polydegree = np.zeros(maxdegree) +x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +x_train_scaled = scaler.transform(x_train) +x_test_scaled = scaler.transform(x_test) + +for degree in range(maxdegree): + model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False)) + clf = model.fit(x_train_scale,y_train) + y_fit = clf.predict(x_train_scaled) + y_pred = clf.predict(x_test_scaled) + polydegree[degree] = degree + TestError[degree] = np.mean( np.mean((y_test - y_pred)**2) ) + TrainError[degree] = np.mean( np.mean((y_train - y_fit)**2) ) + +plt.plot(polydegree, TestError, label='Test Error') +plt.plot(polydegree, TrainError, label='Train Error') +plt.legend() +plt.show() +\epycod + +% --- end solution of exercise --- + +\subex{b)} +Perform an ordinary least squares and compute the means squared error and the $R2$ factor for the training data and the test data, with and without scaling. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +This requires a simple extension to the above code where you simply add a statement calling the $R2$ function included in the same code. + +% --- end solution of exercise --- + +\subex{c)} +Add now a model which allows you to make polynomials up to degree $15$. Perform a standard OLS fitting of the training data and compute the MSE and $R2$ for the training and test data and plot both test and training data MSE and $R2$ as functions of the polynomial degree. Compare what you see with Figure 2.11 of Hastie et al.~Comment your results. For which polynomial degree do you find an optimal MSE (smallest value)? + + +% --- begin solution of exercise --- +\paragraph{Solution.} +Here you simply need to change the degree of the polynomial in the above code to $n=15$. + +% --- end solution of exercise --- + +\subex{d)} +Repeat part (2c) but now using Ridge regressions with various hyperparameters $\lambda$. Make the same plots for the optimal $\lambda$ value for each polynomial degree. Compare these results with those from the standard OLS approach. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +Here you need to add for example the same loop over the parameters $\lambda$ as you did in the first exercise, that is add +\bpycod +nlambdas = 100 +MSEPredictRidge = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + # add ridge + clf_ridge = skl.Ridge(alpha=lmb).fit(X_train_scaled, y_train) + +\epycod +The plotting functionality of the first exercise can be reused here as well. + +% --- end solution of exercise --- + + + + + +\end{doconceexercise} +% --- end exercise --- + + +% ------------------- end of main content --------------- + +% #ifdef PREAMBLE +\end{document} +% #endif + diff --git a/doc/Projects/2020/hw2/pdf/hw2.pdf b/doc/Projects/2020/hw2/pdf/hw2.pdf new file mode 100644 index 000000000..de91d0a65 Binary files /dev/null and b/doc/Projects/2020/hw2/pdf/hw2.pdf differ diff --git a/doc/Projects/2020/hw2/pdf/hw2.tex b/doc/Projects/2020/hw2/pdf/hw2.tex new file mode 100644 index 000000000..cc3c2b843 --- /dev/null +++ b/doc/Projects/2020/hw2/pdf/hw2.tex @@ -0,0 +1,710 @@ +%% +%% Automatically generated file from DocOnce source +%% (https://github.com/hplgit/doconce/) +%% +%% + + +%-------------------- begin preamble ---------------------- + +\documentclass[% +oneside, % oneside: electronic viewing, twoside: printing +final, % draft: marks overfull hboxes, figures with paths +10pt]{article} + +\listfiles % print all files needed to compile this document + +\usepackage{relsize,makeidx,color,setspace,amsmath,amsfonts,amssymb} +\usepackage[table]{xcolor} +\usepackage{bm,ltablex,microtype} + +\usepackage[pdftex]{graphicx} + +\usepackage{fancyvrb} % packages needed for verbatim environments + +\usepackage[T1]{fontenc} +%\usepackage[latin1]{inputenc} +\usepackage{ucs} +\usepackage[utf8x]{inputenc} + +\usepackage{lmodern} % Latin Modern fonts derived from Computer Modern + +% Hyperlinks in PDF: +\definecolor{linkcolor}{rgb}{0,0,0.4} +\usepackage{hyperref} +\hypersetup{ + breaklinks=true, + colorlinks=true, + linkcolor=linkcolor, + urlcolor=linkcolor, + citecolor=black, + filecolor=black, + %filecolor=blue, + pdfmenubar=true, + pdftoolbar=true, + bookmarksdepth=3 % Uncomment (and tweak) for PDF bookmarks with more levels than the TOC + } +%\hyperbaseurl{} % hyperlinks are relative to this root + +\setcounter{tocdepth}{2} % levels in table of contents + +% --- fancyhdr package for fancy headers --- +\usepackage{fancyhdr} +\fancyhf{} % sets both header and footer to nothing +\renewcommand{\headrulewidth}{0pt} +\fancyfoot[LE,RO]{\thepage} +% Ensure copyright on titlepage (article style) and chapter pages (book style) +\fancypagestyle{plain}{ + \fancyhf{} + \fancyfoot[C]{{\footnotesize \copyright\ 1999-2020, "Data Analysis and Machine Learning FYS-STK3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html". Released under CC Attribution-NonCommercial 4.0 license}} +% \renewcommand{\footrulewidth}{0mm} + \renewcommand{\headrulewidth}{0mm} +} +% Ensure copyright on titlepages with \thispagestyle{empty} +\fancypagestyle{empty}{ + \fancyhf{} + \fancyfoot[C]{{\footnotesize \copyright\ 1999-2020, "Data Analysis and Machine Learning FYS-STK3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html". Released under CC Attribution-NonCommercial 4.0 license}} + \renewcommand{\footrulewidth}{0mm} + \renewcommand{\headrulewidth}{0mm} +} + +\pagestyle{fancy} + + +% prevent orhpans and widows +\clubpenalty = 10000 +\widowpenalty = 10000 + +\newenvironment{doconceexercise}{}{} +\newcounter{doconceexercisecounter} + + +% ------ header in subexercises ------ +%\newcommand{\subex}[1]{\paragraph{#1}} +%\newcommand{\subex}[1]{\par\vspace{1.7mm}\noindent{\bf #1}\ \ } +\makeatletter +% 1.5ex is the spacing above the header, 0.5em the spacing after subex title +\newcommand\subex{\@startsection*{paragraph}{4}{\z@}% + {1.5ex\@plus1ex \@minus.2ex}% + {-0.5em}% + {\normalfont\normalsize\bfseries}} +\makeatother + + +% --- end of standard preamble for documents --- + + +% insert custom LaTeX commands... + +\raggedbottom +\makeindex +\usepackage[totoc]{idxlayout} % for index in the toc +\usepackage[nottoc]{tocbibind} % for references/bibliography in the toc + +%-------------------- end preamble ---------------------- + +\begin{document} + +% matching end for #ifdef PREAMBLE + +\newcommand{\exercisesection}[1]{\subsection*{#1}} + + +% ------------------- main content ---------------------- + + + +% ----------------- title ------------------------- + +\thispagestyle{empty} + +\begin{center} +{\LARGE\bf +\begin{spacing}{1.25} +Homework 2, weeks 36 and 37 +\end{spacing} +} +\end{center} + +% ----------------- author(s) ------------------------- + +\begin{center} +{\bf \href{{http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html}}{Data Analysis and Machine Learning FYS-STK3155/FYS4155}} +\end{center} + + \begin{center} +% List of all institutions: +\centerline{{\small Department of Physics, University of Oslo, Norway}} +\end{center} + +% ----------------- end author(s) ------------------------- + +% --- begin date --- +\begin{center} +Sep 1, 2020 +\end{center} +% --- end date --- + +\vspace{1cm} + + +% --- begin exercise --- +\begin{doconceexercise} +\refstepcounter{doconceexercisecounter} + +\exercisesection*{Exercise \thedoconceexercisecounter: Adding Ridge and Lasso Regression} + + +This exercise is a continuation of exercise 3 from exercise set 1 (week 35). 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, see the \href{{https://compphysics.github.io/MachineLearning/doc/pub/Regression/html/Regression-bs.html}}{regression slides}). + +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). +\begin{verbatim} +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) +\end{verbatim} + + +\subex{a)} +Write your own code for the Ridge method (see chapter 3.4 of Hastie \emph{et al.}, equations (3.43) and (3.44)) and compute the parametrization for different values of $\lambda$. Compare and analyze your results with those from exercise 3. Study the dependence on $\lambda$ while also varying the strength of the noise in your expression for $y(x)$. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +The code here allows you to perform your own Ridge calculation and perform calculations for various values of the regularization parameter $\lambda$. This program can easily be extended upon. +\begin{verbatim} +import os +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# and then make the prediction +ytildeOLS = X_train @ OLSbeta +print("Training R2 for OLS") +print(R2(y_train,ytildeOLS)) +print("Training MSE for OLS") +print(MSE(y_train,ytildeOLS)) +ypredictOLS = X_test @ OLSbeta +print("Test R2 for OLS") +print(R2(y_test,ypredictOLS)) +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 20 +MSEPredict = np.zeros(nlambdas) +MSETrain = np.zeros(nlambdas) +lambdas = np.logspace(-4, 1, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train + # and then make the prediction + ytildeRidge = X_train @ Ridgebeta + ypredictRidge = X_test @ Ridgebeta + MSEPredict[i] = MSE(y_test,ypredictRidge) + MSETrain[i] = MSE(y_train,ytildeRidge) +# Now plot the resulys +plt.figure() +plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train') +plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() +\end{verbatim} + +% --- end solution of exercise --- + +\subex{b)} +Repeat the above but using the functionality of \textbf{Scikit-Learn}. Compare your code with the results from \textbf{Scikit-Learn}. Remember to run with the same random numbers for generating $x$ and $y$. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +To use \textbf{scikit-learn} with Ridge, we simply need to add the relevant function \textbf{Ridge()}, as done in the code here. +\begin{verbatim} +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler +import sklearn.linear_model as skl + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# and then make the prediction +ytildeOLS = X_train @ OLSbeta +print("Training R2 for OLS") +print(R2(y_train,ytildeOLS)) +print("Training MSE for OLS") +print(MSE(y_train,ytildeOLS)) +ypredictOLS = X_test @ OLSbeta +print("Test R2 for OLS") +print(R2(y_test,ypredictOLS)) +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 100 +MSEPredict = np.zeros(nlambdas) +MSEPredictSKL = np.zeros(nlambdas) +MSETrain = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + # add ridge + clf_ridge = skl.Ridge(alpha=lmb).fit(X_train, y_train) + yridge = clf_ridge.predict(X_test) + Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train + # and then make the prediction + ytildeRidge = X_train @ Ridgebeta + ypredictRidge = X_test @ Ridgebeta + MSEPredict[i] = MSE(y_test,ypredictRidge) + MSEPredictSKL[i] = MSE(y_test,yridge) + MSETrain[i] = MSE(y_train,ytildeRidge) +#then plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train') +plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test') +plt.plot(np.log10(lambdas), MSEPredictSKL, 'g--', label = 'MSE Ridge sickit-learn Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() +\end{verbatim} + +% --- end solution of exercise --- + +\subex{c)} +Our next step is to study the variance of the parameters $\beta_1$ and $\beta_2$ (assuming that we are parameterizing 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 or using \textbf{Scikit-Learn} to find the parameters $\beta$. From your results calculate the variance of these parameters (recall that this is equal to the diagonal elements of the matrix $(\hat{X}^T\hat{X})+\lambda\hat{I})^{-1}$). Discuss the results of these variances as functions of $\lambda$. In particular, try to link your discussion with the discussion in Hastie \emph{et al.} and their figures 3.10 and 3.11. \textbf{Scikit-Learn} may not provide the variance of the parameters $\beta$. This needs to be checked. With your own code you can however do so. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +\begin{verbatim} +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler +import sklearn.linear_model as skl + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# The variance is given by the inverse of the matrix X^TX +print(np.linalg.inv(X_train.T @ X_train)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 10 +MSEPredict = np.zeros(nlambdas) +MSEPredictSKL = np.zeros(nlambdas) +MSETrain = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train + print(np.linalg.inv(X_train.T @ X_train+lmb*I)) + + + +\end{verbatim} + +% --- end solution of exercise --- + +\subex{d)} +Repeat the previous step but add now the Lasso method, see equation (3.53) of Hastie \emph{et al.}. Discuss your results and compare with standard regression and the Ridge regression results. You can write your own code or use the functionality of \textbf{scikit-learn}. We recommend the latter since we have not yet discussed how to solve the Lasso equations numerically. Also, you do not need to compute the variance of the parameters $\beta$ but you can extract their values and study their behavior as functions of the regularization parameter $\lambda$. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +\begin{verbatim} +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler +import sklearn.linear_model as skl +#from sklearn.linear_model import LinearRegression, Ridge, Lasso +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# and then make the prediction +ytildeOLS = X_train @ OLSbeta +print("Training R2 for OLS") +print(R2(y_train,ytildeOLS)) +print("Training MSE for OLS") +print(MSE(y_train,ytildeOLS)) +ypredictOLS = X_test @ OLSbeta +print("Test R2 for OLS") +print(R2(y_test,ypredictOLS)) +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 100 +MSEPredictLasso = np.zeros(nlambdas) +MSEPredictRidge = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + # add ridge + clf_ridge = skl.Ridge(alpha=lmb).fit(X_train, y_train) + clf_lasso = skl.Lasso(alpha=lmb).fit(X_train, y_train) + yridge = clf_ridge.predict(X_test) + ylasso = clf_lasso.predict(X_test) + MSEPredictLasso[i] = MSE(y_test,ylasso) + MSEPredictRidge[i] = MSE(y_test,yridge) +#then plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSEPredictRidge, 'r--', label = 'MSE Ridge Test') +plt.plot(np.log10(lambdas), MSEPredictLasso, 'g--', label = 'MSE Lasso Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() +\end{verbatim} + +% --- end solution of exercise --- + +\subex{e)} +Finally, using \textbf{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 +\[ 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. +\] +Discuss these quantities as functions of the variable $\lambda$ in the Ridge and Lasso regression methods. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +These results can all be studied with the codes we have above. These scores are included in the codes above. + +% --- end solution of exercise --- + + +\end{doconceexercise} +% --- end exercise --- + + + + +% --- begin exercise --- +\begin{doconceexercise} +\refstepcounter{doconceexercisecounter} + +\exercisesection*{Exercise \thedoconceexercisecounter: Normalizing our data} + + +A much used approach before starting to train the data is to preprocess our +data. Normally the data may need a rescaling and/or may be sensitive +to extreme values. Scaling the data renders our inputs much more +suitable for the algorithms we want to employ. + +\textbf{Scikit-Learn} has several functions which allow us to rescale the +data, normally resulting in much better results in terms of various +accuracy scores. The \textbf{StandardScaler} function in \textbf{Scikit-Learn} +ensures that for each feature/predictor we study the mean value is +zero and the variance is one (every column in the design/feature +matrix). This scaling has the drawback that it does not ensure that +we have a particular maximum or minimum in our data set. Another +function included in \textbf{Scikit-Learn} is the \textbf{MinMaxScaler} which +ensures that all features are exactly between $0$ and $1$. The + + +The \textbf{Normalizer} scales each data +point such that the feature vector has a euclidean length of one. In other words, it +projects a data point on the circle (or sphere in the case of higher dimensions) with a +radius of 1. This means every data point is scaled by a different number (by the +inverse of it’s length). +This normalization is often used when only the direction (or angle) of the data matters, +not the length of the feature vector. + +The \textbf{RobustScaler} works similarly to the StandardScaler in that it +ensures statistical properties for each feature that guarantee that +they are on the same scale. However, the RobustScaler uses the median +and quartiles, instead of mean and variance. This makes the +RobustScaler ignore data points that are very different from the rest +(like measurement errors). These odd data points are also called +outliers, and might often lead to trouble for other scaling +techniques. + + +It also common to split the data in a \textbf{training} set and a \textbf{testing} set. A typical split is to use $80\%$ of the data for training and the rest +for testing. This can be done as follows with our design matrix $\bm{X}$ and data $\bm{y}$ (remember to import \textbf{scikit-learn}) +\begin{verbatim} +# split in training and test data +X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.2) +\end{verbatim} +Then we can use the standard scaler to scale our data as +\begin{verbatim} +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) +\end{verbatim} + + +In this exercise we want you to to compute the MSE for the training +data and the test data as function of the complexity of a polynomial, +that is the degree of a given polynomial. We want you also to compute the $R2$ score as function of the complexity of the model for both training data and test data. You should also run the calculation with and without scaling. + +One of +the aims is to reproduce Figure 2.11 of \href{{https://github.com/CompPhysics/MLErasmus/blob/master/doc/Textbooks/elementsstat.pdf}}{Hastie et al}. +We will also use Ridge and Lasso regression. + + +Our data is defined by $x\in [-3,3]$ with a total of for example $100$ data points. +\begin{verbatim} +np.random.seed() +n = 100 +maxdegree = 14 +# Make data set. +x = np.linspace(-3, 3, n).reshape(-1, 1) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape) +\end{verbatim} +where $y$ is the function we want to fit with a given polynomial. + + +\subex{a)} +Write a first code which sets up a design matrix $X$ defined by a fifth-order polynomial. Scale your data and split it in training and test data. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +\begin{verbatim} +import matplotlib.pyplot as plt +import numpy as np +from sklearn.linear_model import LinearRegression, Ridge, Lasso +from sklearn.preprocessing import PolynomialFeatures +from sklearn.model_selection import train_test_split +from sklearn.pipeline import make_pipeline + + +np.random.seed(2018) +n = 50 +maxdegree = 5 +# Make data set. +x = np.linspace(-3, 3, n).reshape(-1, 1) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape) +TestError = np.zeros(maxdegree) +TrainError = np.zeros(maxdegree) +polydegree = np.zeros(maxdegree) +x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +x_train_scaled = scaler.transform(x_train) +x_test_scaled = scaler.transform(x_test) + +for degree in range(maxdegree): + model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False)) + clf = model.fit(x_train_scale,y_train) + y_fit = clf.predict(x_train_scaled) + y_pred = clf.predict(x_test_scaled) + polydegree[degree] = degree + TestError[degree] = np.mean( np.mean((y_test - y_pred)**2) ) + TrainError[degree] = np.mean( np.mean((y_train - y_fit)**2) ) + +plt.plot(polydegree, TestError, label='Test Error') +plt.plot(polydegree, TrainError, label='Train Error') +plt.legend() +plt.show() +\end{verbatim} + +% --- end solution of exercise --- + +\subex{b)} +Perform an ordinary least squares and compute the means squared error and the $R2$ factor for the training data and the test data, with and without scaling. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +This requires a simple extension to the above code where you simply add a statement calling the $R2$ function included in the same code. + +% --- end solution of exercise --- + +\subex{c)} +Add now a model which allows you to make polynomials up to degree $15$. Perform a standard OLS fitting of the training data and compute the MSE and $R2$ for the training and test data and plot both test and training data MSE and $R2$ as functions of the polynomial degree. Compare what you see with Figure 2.11 of Hastie et al.~Comment your results. For which polynomial degree do you find an optimal MSE (smallest value)? + + +% --- begin solution of exercise --- +\paragraph{Solution.} +Here you simply need to change the degree of the polynomial in the above code to $n=15$. + +% --- end solution of exercise --- + +\subex{d)} +Repeat part (2c) but now using Ridge regressions with various hyperparameters $\lambda$. Make the same plots for the optimal $\lambda$ value for each polynomial degree. Compare these results with those from the standard OLS approach. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +Here you need to add for example the same loop over the parameters $\lambda$ as you did in the first exercise, that is add +\begin{verbatim} +nlambdas = 100 +MSEPredictRidge = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + # add ridge + clf_ridge = skl.Ridge(alpha=lmb).fit(X_train_scaled, y_train) + +\end{verbatim} +The plotting functionality of the first exercise can be reused here as well. + +% --- end solution of exercise --- + + + + + +\end{doconceexercise} +% --- end exercise --- + + +% ------------------- end of main content --------------- + +\end{document} + diff --git a/doc/Projects/2020/hw2/pdf/hw2.tex~ b/doc/Projects/2020/hw2/pdf/hw2.tex~ new file mode 100644 index 000000000..41030d5f5 --- /dev/null +++ b/doc/Projects/2020/hw2/pdf/hw2.tex~ @@ -0,0 +1,710 @@ +%% +%% Automatically generated file from DocOnce source +%% (https://github.com/hplgit/doconce/) +%% +%% + + +%-------------------- begin preamble ---------------------- + +\documentclass[% +oneside, % oneside: electronic viewing, twoside: printing +final, % draft: marks overfull hboxes, figures with paths +10pt]{article} + +\listfiles % print all files needed to compile this document + +\usepackage{relsize,makeidx,color,setspace,amsmath,amsfonts,amssymb} +\usepackage[table]{xcolor} +\usepackage{bm,ltablex,microtype} + +\usepackage[pdftex]{graphicx} + +\usepackage{fancyvrb} % packages needed for verbatim environments + +\usepackage[T1]{fontenc} +%\usepackage[latin1]{inputenc} +\usepackage{ucs} +\usepackage[utf8x]{inputenc} + +\usepackage{lmodern} % Latin Modern fonts derived from Computer Modern + +% Hyperlinks in PDF: +\definecolor{linkcolor}{rgb}{0,0,0.4} +\usepackage{hyperref} +\hypersetup{ + breaklinks=true, + colorlinks=true, + linkcolor=linkcolor, + urlcolor=linkcolor, + citecolor=black, + filecolor=black, + %filecolor=blue, + pdfmenubar=true, + pdftoolbar=true, + bookmarksdepth=3 % Uncomment (and tweak) for PDF bookmarks with more levels than the TOC + } +%\hyperbaseurl{} % hyperlinks are relative to this root + +\setcounter{tocdepth}{2} % levels in table of contents + +% --- fancyhdr package for fancy headers --- +\usepackage{fancyhdr} +\fancyhf{} % sets both header and footer to nothing +\renewcommand{\headrulewidth}{0pt} +\fancyfoot[LE,RO]{\thepage} +% Ensure copyright on titlepage (article style) and chapter pages (book style) +\fancypagestyle{plain}{ + \fancyhf{} + \fancyfoot[C]{{\footnotesize \copyright\ 1999-2020, "Data Analysis and Machine Learning FYS-STK3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html". Released under CC Attribution-NonCommercial 4.0 license}} +% \renewcommand{\footrulewidth}{0mm} + \renewcommand{\headrulewidth}{0mm} +} +% Ensure copyright on titlepages with \thispagestyle{empty} +\fancypagestyle{empty}{ + \fancyhf{} + \fancyfoot[C]{{\footnotesize \copyright\ 1999-2020, "Data Analysis and Machine Learning FYS-STK3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html". Released under CC Attribution-NonCommercial 4.0 license}} + \renewcommand{\footrulewidth}{0mm} + \renewcommand{\headrulewidth}{0mm} +} + +\pagestyle{fancy} + + +% prevent orhpans and widows +\clubpenalty = 10000 +\widowpenalty = 10000 + +\newenvironment{doconceexercise}{}{} +\newcounter{doconceexercisecounter} + + +% ------ header in subexercises ------ +%\newcommand{\subex}[1]{\paragraph{#1}} +%\newcommand{\subex}[1]{\par\vspace{1.7mm}\noindent{\bf #1}\ \ } +\makeatletter +% 1.5ex is the spacing above the header, 0.5em the spacing after subex title +\newcommand\subex{\@startsection*{paragraph}{4}{\z@}% + {1.5ex\@plus1ex \@minus.2ex}% + {-0.5em}% + {\normalfont\normalsize\bfseries}} +\makeatother + + +% --- end of standard preamble for documents --- + + +% insert custom LaTeX commands... + +\raggedbottom +\makeindex +\usepackage[totoc]{idxlayout} % for index in the toc +\usepackage[nottoc]{tocbibind} % for references/bibliography in the toc + +%-------------------- end preamble ---------------------- + +\begin{document} + +% matching end for #ifdef PREAMBLE + +\newcommand{\exercisesection}[1]{\subsection*{#1}} + + +% ------------------- main content ---------------------- + + + +% ----------------- title ------------------------- + +\thispagestyle{empty} + +\begin{center} +{\LARGE\bf +\begin{spacing}{1.25} +Homework 2, weeks 36 and 37 +\end{spacing} +} +\end{center} + +% ----------------- author(s) ------------------------- + +\begin{center} +{\bf \href{{http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html}}{Data Analysis and Machine Learning FYS-STK3155/FYS4155}} +\end{center} + + \begin{center} +% List of all institutions: +\centerline{{\small Department of Physics, University of Oslo, Norway}} +\end{center} + +% ----------------- end author(s) ------------------------- + +% --- begin date --- +\begin{center} +Sep 1, 2020 +\end{center} +% --- end date --- + +\vspace{1cm} + + +% --- begin exercise --- +\begin{doconceexercise} +\refstepcounter{doconceexercisecounter} + +\exercisesection*{Exercise \thedoconceexercisecounter: Adding Ridge and Lasso Regression} + + +This exercise is a continuation of exercise 3 from exercise set 1 (week 35). 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, see the \href{{https://compphysics.github.io/MachineLearning/doc/pub/Regression/html/Regression-bs.html}}{regression slides}). + +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). +\begin{print} +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) +\end{print} + + +\subex{a)} +Write your own code for the Ridge method (see chapter 3.4 of Hastie \emph{et al.}, equations (3.43) and (3.44)) and compute the parametrization for different values of $\lambda$. Compare and analyze your results with those from exercise 3. Study the dependence on $\lambda$ while also varying the strength of the noise in your expression for $y(x)$. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +The code here allows you to perform your own Ridge calculation and perform calculations for various values of the regularization parameter $\lambda$. This program can easily be extended upon. +\begin{print} +import os +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# and then make the prediction +ytildeOLS = X_train @ OLSbeta +print("Training R2 for OLS") +print(R2(y_train,ytildeOLS)) +print("Training MSE for OLS") +print(MSE(y_train,ytildeOLS)) +ypredictOLS = X_test @ OLSbeta +print("Test R2 for OLS") +print(R2(y_test,ypredictOLS)) +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 20 +MSEPredict = np.zeros(nlambdas) +MSETrain = np.zeros(nlambdas) +lambdas = np.logspace(-4, 1, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train + # and then make the prediction + ytildeRidge = X_train @ Ridgebeta + ypredictRidge = X_test @ Ridgebeta + MSEPredict[i] = MSE(y_test,ypredictRidge) + MSETrain[i] = MSE(y_train,ytildeRidge) +# Now plot the resulys +plt.figure() +plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train') +plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() +\end{print} + +% --- end solution of exercise --- + +\subex{b)} +Repeat the above but using the functionality of \textbf{Scikit-Learn}. Compare your code with the results from \textbf{Scikit-Learn}. Remember to run with the same random numbers for generating $x$ and $y$. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +To use \textbf{scikit-learn} with Ridge, we simply need to add the relevant function \textbf{Ridge()}, as done in the code here. +\begin{print} +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler +import sklearn.linear_model as skl + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# and then make the prediction +ytildeOLS = X_train @ OLSbeta +print("Training R2 for OLS") +print(R2(y_train,ytildeOLS)) +print("Training MSE for OLS") +print(MSE(y_train,ytildeOLS)) +ypredictOLS = X_test @ OLSbeta +print("Test R2 for OLS") +print(R2(y_test,ypredictOLS)) +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 100 +MSEPredict = np.zeros(nlambdas) +MSEPredictSKL = np.zeros(nlambdas) +MSETrain = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + # add ridge + clf_ridge = skl.Ridge(alpha=lmb).fit(X_train, y_train) + yridge = clf_ridge.predict(X_test) + Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train + # and then make the prediction + ytildeRidge = X_train @ Ridgebeta + ypredictRidge = X_test @ Ridgebeta + MSEPredict[i] = MSE(y_test,ypredictRidge) + MSEPredictSKL[i] = MSE(y_test,yridge) + MSETrain[i] = MSE(y_train,ytildeRidge) +#then plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train') +plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test') +plt.plot(np.log10(lambdas), MSEPredictSKL, 'g--', label = 'MSE Ridge sickit-learn Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() +\end{print} + +% --- end solution of exercise --- + +\subex{c)} +Our next step is to study the variance of the parameters $\beta_1$ and $\beta_2$ (assuming that we are parameterizing 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 or using \textbf{Scikit-Learn} to find the parameters $\beta$. From your results calculate the variance of these parameters (recall that this is equal to the diagonal elements of the matrix $(\hat{X}^T\hat{X})+\lambda\hat{I})^{-1}$). Discuss the results of these variances as functions of $\lambda$. In particular, try to link your discussion with the discussion in Hastie \emph{et al.} and their figures 3.10 and 3.11. \textbf{Scikit-Learn} may not provide the variance of the parameters $\beta$. This needs to be checked. With your own code you can however do so. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +\begin{print} +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler +import sklearn.linear_model as skl + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# The variance is given by the inverse of the matrix X^TX +print(np.linalg.inv(X_train.T @ X_train)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 10 +MSEPredict = np.zeros(nlambdas) +MSEPredictSKL = np.zeros(nlambdas) +MSETrain = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train + print(np.linalg.inv(X_train.T @ X_train+lmb*I)) + + + +\end{print} + +% --- end solution of exercise --- + +\subex{d)} +Repeat the previous step but add now the Lasso method, see equation (3.53) of Hastie \emph{et al.}. Discuss your results and compare with standard regression and the Ridge regression results. You can write your own code or use the functionality of \textbf{scikit-learn}. We recommend the latter since we have not yet discussed how to solve the Lasso equations numerically. Also, you do not need to compute the variance of the parameters $\beta$ but you can extract their values and study their behavior as functions of the regularization parameter $\lambda$. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +\begin{print} +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import StandardScaler +import sklearn.linear_model as skl +#from sklearn.linear_model import LinearRegression, Ridge, Lasso +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# and then make the prediction +ytildeOLS = X_train @ OLSbeta +print("Training R2 for OLS") +print(R2(y_train,ytildeOLS)) +print("Training MSE for OLS") +print(MSE(y_train,ytildeOLS)) +ypredictOLS = X_test @ OLSbeta +print("Test R2 for OLS") +print(R2(y_test,ypredictOLS)) +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 100 +MSEPredictLasso = np.zeros(nlambdas) +MSEPredictRidge = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + # add ridge + clf_ridge = skl.Ridge(alpha=lmb).fit(X_train, y_train) + clf_lasso = skl.Lasso(alpha=lmb).fit(X_train, y_train) + yridge = clf_ridge.predict(X_test) + ylasso = clf_lasso.predict(X_test) + MSEPredictLasso[i] = MSE(y_test,ylasso) + MSEPredictRidge[i] = MSE(y_test,yridge) +#then plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSEPredictRidge, 'r--', label = 'MSE Ridge Test') +plt.plot(np.log10(lambdas), MSEPredictLasso, 'g--', label = 'MSE Lasso Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() +\end{print} + +% --- end solution of exercise --- + +\subex{e)} +Finally, using \textbf{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 +\[ 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. +\] +Discuss these quantities as functions of the variable $\lambda$ in the Ridge and Lasso regression methods. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +These results can all be studied with the codes we have above. These scores are included in the codes above. + +% --- end solution of exercise --- + + +\end{doconceexercise} +% --- end exercise --- + + + + +% --- begin exercise --- +\begin{doconceexercise} +\refstepcounter{doconceexercisecounter} + +\exercisesection*{Exercise \thedoconceexercisecounter: Normalizing our data} + + +A much used approach before starting to train the data is to preprocess our +data. Normally the data may need a rescaling and/or may be sensitive +to extreme values. Scaling the data renders our inputs much more +suitable for the algorithms we want to employ. + +\textbf{Scikit-Learn} has several functions which allow us to rescale the +data, normally resulting in much better results in terms of various +accuracy scores. The \textbf{StandardScaler} function in \textbf{Scikit-Learn} +ensures that for each feature/predictor we study the mean value is +zero and the variance is one (every column in the design/feature +matrix). This scaling has the drawback that it does not ensure that +we have a particular maximum or minimum in our data set. Another +function included in \textbf{Scikit-Learn} is the \textbf{MinMaxScaler} which +ensures that all features are exactly between $0$ and $1$. The + + +The \textbf{Normalizer} scales each data +point such that the feature vector has a euclidean length of one. In other words, it +projects a data point on the circle (or sphere in the case of higher dimensions) with a +radius of 1. This means every data point is scaled by a different number (by the +inverse of it’s length). +This normalization is often used when only the direction (or angle) of the data matters, +not the length of the feature vector. + +The \textbf{RobustScaler} works similarly to the StandardScaler in that it +ensures statistical properties for each feature that guarantee that +they are on the same scale. However, the RobustScaler uses the median +and quartiles, instead of mean and variance. This makes the +RobustScaler ignore data points that are very different from the rest +(like measurement errors). These odd data points are also called +outliers, and might often lead to trouble for other scaling +techniques. + + +It also common to split the data in a \textbf{training} set and a \textbf{testing} set. A typical split is to use $80\%$ of the data for training and the rest +for testing. This can be done as follows with our design matrix $\bm{X}$ and data $\bm{y}$ (remember to import \textbf{scikit-learn}) +\begin{print} +# split in training and test data +X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.2) +\end{print} +Then we can use the standard scaler to scale our data as +\begin{print} +scaler = StandardScaler() +scaler.fit(X_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) +\end{print} + + +In this exercise we want you to to compute the MSE for the training +data and the test data as function of the complexity of a polynomial, +that is the degree of a given polynomial. We want you also to compute the $R2$ score as function of the complexity of the model for both training data and test data. You should also run the calculation with and without scaling. + +One of +the aims is to reproduce Figure 2.11 of \href{{https://github.com/CompPhysics/MLErasmus/blob/master/doc/Textbooks/elementsstat.pdf}}{Hastie et al}. +We will also use Ridge and Lasso regression. + + +Our data is defined by $x\in [-3,3]$ with a total of for example $100$ data points. +\begin{print} +np.random.seed() +n = 100 +maxdegree = 14 +# Make data set. +x = np.linspace(-3, 3, n).reshape(-1, 1) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape) +\end{print} +where $y$ is the function we want to fit with a given polynomial. + + +\subex{a)} +Write a first code which sets up a design matrix $X$ defined by a fifth-order polynomial. Scale your data and split it in training and test data. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +\begin{print} +import matplotlib.pyplot as plt +import numpy as np +from sklearn.linear_model import LinearRegression, Ridge, Lasso +from sklearn.preprocessing import PolynomialFeatures +from sklearn.model_selection import train_test_split +from sklearn.pipeline import make_pipeline + + +np.random.seed(2018) +n = 50 +maxdegree = 5 +# Make data set. +x = np.linspace(-3, 3, n).reshape(-1, 1) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape) +TestError = np.zeros(maxdegree) +TrainError = np.zeros(maxdegree) +polydegree = np.zeros(maxdegree) +x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2) +scaler = StandardScaler() +scaler.fit(X_train) +x_train_scaled = scaler.transform(x_train) +x_test_scaled = scaler.transform(x_test) + +for degree in range(maxdegree): + model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False)) + clf = model.fit(x_train_scale,y_train) + y_fit = clf.predict(x_train_scaled) + y_pred = clf.predict(x_test_scaled) + polydegree[degree] = degree + TestError[degree] = np.mean( np.mean((y_test - y_pred)**2) ) + TrainError[degree] = np.mean( np.mean((y_train - y_fit)**2) ) + +plt.plot(polydegree, TestError, label='Test Error') +plt.plot(polydegree, TrainError, label='Train Error') +plt.legend() +plt.show() +\end{print} + +% --- end solution of exercise --- + +\subex{b)} +Perform an ordinary least squares and compute the means squared error and the $R2$ factor for the training data and the test data, with and without scaling. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +This requires a simple extension to the above code where you simply add a statement calling the $R2$ function included in the same code. + +% --- end solution of exercise --- + +\subex{c)} +Add now a model which allows you to make polynomials up to degree $15$. Perform a standard OLS fitting of the training data and compute the MSE and $R2$ for the training and test data and plot both test and training data MSE and $R2$ as functions of the polynomial degree. Compare what you see with Figure 2.11 of Hastie et al.~Comment your results. For which polynomial degree do you find an optimal MSE (smallest value)? + + +% --- begin solution of exercise --- +\paragraph{Solution.} +Here you simply need to change the degree of the polynomial in the above code to $n=15$. + +% --- end solution of exercise --- + +\subex{d)} +Repeat part (2c) but now using Ridge regressions with various hyperparameters $\lambda$. Make the same plots for the optimal $\lambda$ value for each polynomial degree. Compare these results with those from the standard OLS approach. + + +% --- begin solution of exercise --- +\paragraph{Solution.} +Here you need to add for example the same loop over the parameters $\lambda$ as you did in the first exercise, that is add +\begin{print} +nlambdas = 100 +MSEPredictRidge = np.zeros(nlambdas) +lambdas = np.logspace(-4, 0, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + # add ridge + clf_ridge = skl.Ridge(alpha=lmb).fit(X_train_scaled, y_train) + +\end{print} +The plotting functionality of the first exercise can be reused here as well. + +% --- end solution of exercise --- + + + + + +\end{doconceexercise} +% --- end exercise --- + + +% ------------------- end of main content --------------- + +\end{document} +