+Small exercise: perform the standars scaling by including the standard deviation.
@@ -1347,7 +1350,73 @@ where \( \min(x_j) \) and \( \max(x_j) \) return the minimum and maximum value o
+Our data is defined by \( x\in [-3,3] \) with a total of for example \( 100 \) data points.
+
+
+
+
+where \( y \) is the function we want to fit with a given polynomial.
+
+
+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.
+
+
+
+
@@ -1449,10 +1518,10 @@ clf = skl.LinearRegression().fit(X_train_scaled, y_train)
-Lasso and Ridge regression
+What is presented here is a mathematical analysis of various regression algorithms (ordinary least squares, Ridge and Lasso Regression). The analysis is based on an important algorithm in linear algebra, the so-called Singular Value Decomposition (SVD).
The examples we have looked at so far are cases where we normally can
invert the matrix \( \boldsymbol{X}^T\boldsymbol{X} \). Using a polynomial expansion as we
-did both for the masses and the fitting of the equation of state,
-leads to row vectors of the design matrix which are essentially
-orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition.
+did both for the masses and the fitting of various functions leads to
+row vectors of the design matrix which are essentially orthogonal due
+to the polynomial character of our model. Obtaining the inverse of the
+design matrix is then often done via a so-called LU, QR or Cholesky
+decomposition.
-This may
+As we will also see in the first project,
+this may
however not the be case in general and a standard matrix inversion
algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.
@@ -1849,10 +1921,10 @@ $$ n > p$$
-Video of Lecture and handwritten notes
+Video of Lecture from 2020 and handwritten notes
More material will be added here, see handwritten notes also.
@@ -2492,6 +2564,334 @@ It is easy to generalize this to a matrix \( \boldsymbol{X}\in {\mathbb{R}}^{n\t
+More material will be added here.
+
+The exercises here are meant to prepare you for work with project 1. The first exercise is a follow-up of exercise 2 from week 36 August 30-September 3).
+
+
+
+
+
+This exercise is a continuation of exercise 2 from exercise set 1 (week 36, August 30-September 3)). 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.
+
+
+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).
+
+
+
+
+Write your own code for the Ridge method (see chapter 3.4 of Hastie 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) \).
+
+
+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.
+
+
+
+
+
+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 \).
+
+
+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
+
+This is a longer exercise and the aim is to study in more detail various
+regression methods, including the Ordinary Least Squares (OLS) method,
+Ridge regression and finally Lasso regression.
+This exercise forms a part of project 1.
+
+
+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.
+
+
+The Franke function, which is a weighted sum of four exponentials reads as follows
+
+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 fit a
+function (for example a polynomial) of \( x \) and \( y \). Thereafter we
+will repeat much of the same procedure using 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)
+
+
+
+
+We will generate our own dataset for a function
+\( \mathrm{FrankeFunction}(x,y) \) with \( x,y \in [0,1] \). 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 and perform a standard least square regression
+analysis using polynomials in \( x \) and \( y \) up to fifth order. You can use scikit-learn as well.
+
+
+Evaluate the Mean Squared error (MSE)
+
+
+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
+
+
+where we have defined the mean value of \( \hat{y} \) as
+
+
+You should split your data in train and test and also consider scaling the data.
+
+
+To set up the design matrix, the following code can be used
+
+
+
+
+Write then your own code for the Ridge method or use Scikit-Learn.
+Perform the same analysis as you did for ordinary Least Squares (for the same polynomials) but now for different values of \( \lambda \). Compare and
+analyze your results with those obtained with ordinary Least Squares. Study the
+dependence on \( \lambda \).
+
+
+This part is essentially a repeat of the previous ones, but now
+with Lasso regression. Write either your own code or
+use the functionalities of Scikit-Learn (recommended).
+Give a
+critical discussion of the three methods and a judgement of which
+model fits the data best.
+
+
+
+Small exercise: perform the standars scaling by including the standard deviation.
+
+Our data is defined by \( x\in [-3,3] \) with a total of for example \( 100 \) data points.
+
+
+
+
+where \( y \) is the function we want to fit with a given polynomial.
+
+
+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.
+
+
+
+
+
@@ -1521,10 +1609,10 @@ clf = skl.LinearRegression().fit(X_train_scaled, y_train)
-Lasso and Ridge regression
+What is presented here is a mathematical analysis of various regression algorithms (ordinary least squares, Ridge and Lasso Regression). The analysis is based on an important algorithm in linear algebra, the so-called Singular Value Decomposition (SVD).
The examples we have looked at so far are cases where we normally can
invert the matrix \( \boldsymbol{X}^T\boldsymbol{X} \). Using a polynomial expansion as we
-did both for the masses and the fitting of the equation of state,
-leads to row vectors of the design matrix which are essentially
-orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition.
+did both for the masses and the fitting of various functions leads to
+row vectors of the design matrix which are essentially orthogonal due
+to the polynomial character of our model. Obtaining the inverse of the
+design matrix is then often done via a so-called LU, QR or Cholesky
+decomposition.
-This may
+As we will also see in the first project,
+this may
however not the be case in general and a standard matrix inversion
algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.
@@ -1872,10 +1963,10 @@ Here we have that $${\bf X} = {\bf U}{\bf \Sigma}{\bf V}^T$$, with $$\Sigma$$ be
-Video of Lecture and handwritten notes
+Video of Lecture from 2020 and handwritten notes
More material will be added here, see handwritten notes also.
@@ -2446,6 +2537,318 @@ It is easy to generalize this to a matrix \( \boldsymbol{X}\in {\mathbb{R}}^{n\t
+More material will be added here.
+
+
+
+The exercises here are meant to prepare you for work with project 1. The first exercise is a follow-up of exercise 2 from week 36 August 30-September 3).
+
+
+
+
+
+This exercise is a continuation of exercise 2 from exercise set 1 (week 36, August 30-September 3)). 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.
+
+
+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).
+
+
+
+
+Write your own code for the Ridge method (see chapter 3.4 of Hastie 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) \).
+
+
+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.
+
+
+
+
+
+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 \).
+
+
+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
+$$ 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.
+
+
+This is a longer exercise and the aim is to study in more detail various
+regression methods, including the Ordinary Least Squares (OLS) method,
+Ridge regression and finally Lasso regression.
+This exercise forms a part of project 1.
+
+
+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.
+
+
+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 fit a
+function (for example a polynomial) of \( x \) and \( y \). Thereafter we
+will repeat much of the same procedure using 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)
+
+
+
+
+We will generate our own dataset for a function
+\( \mathrm{FrankeFunction}(x,y) \) with \( x,y \in [0,1] \). 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 and perform a standard least square regression
+analysis using polynomials in \( x \) and \( y \) up to fifth order. You can use scikit-learn as well.
+
+
+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.
+$$
+
+
+You should split your data in train and test and also consider scaling the data.
+
+
+To set up the design matrix, the following code can be used
+
+
+
+
+Write then your own code for the Ridge method or use Scikit-Learn.
+Perform the same analysis as you did for ordinary Least Squares (for the same polynomials) but now for different values of \( \lambda \). Compare and
+analyze your results with those obtained with ordinary Least Squares. Study the
+dependence on \( \lambda \).
+
+
+This part is essentially a repeat of the previous ones, but now
+with Lasso regression. Write either your own code or
+use the functionalities of Scikit-Learn (recommended).
+Give a
+critical discussion of the three methods and a judgement of which
+model fits the data best.
+
+
+
diff --git a/doc/pub/week35/html/week35.html b/doc/pub/week35/html/week35.html
index 6d38c7927..807f09764 100644
--- a/doc/pub/week35/html/week35.html
+++ b/doc/pub/week35/html/week35.html
@@ -172,11 +172,18 @@ div { text-align: justify; text-justify: inter-word; }
None,
'example-of-own-standard-scaling'),
('Min-Max Scaling', 2, None, 'min-max-scaling'),
- ('Simple preprocessing examples, Franke function and regression',
+ ('Testing the Means Squared Error as function of Complexity',
2,
None,
- 'simple-preprocessing-examples-franke-function-and-regression'),
- ('Friday September 3', 2, None, 'friday-september-3'),
+ 'testing-the-means-squared-error-as-function-of-complexity'),
+ ('More preprocessing examples, Franke function and regression',
+ 2,
+ None,
+ 'more-preprocessing-examples-franke-function-and-regression'),
+ ('Mathematical Interpretation of Ordinary Least Squares',
+ 2,
+ None,
+ 'mathematical-interpretation-of-ordinary-least-squares'),
('The singular value decomposition',
2,
None,
@@ -194,7 +201,7 @@ div { text-align: justify; text-justify: inter-word; }
('Economy-size SVD', 2, None, 'economy-size-svd'),
('Codes for the SVD', 2, None, 'codes-for-the-svd'),
('Mathematical Properties', 2, None, 'mathematical-properties'),
- ('Friday September 12', 2, None, 'friday-september-12'),
+ ('Friday September 3', 2, None, 'friday-september-3'),
('Ridge and LASSO Regression',
2,
None,
@@ -238,7 +245,19 @@ div { text-align: justify; text-justify: inter-word; }
2,
None,
'rewriting-the-covariance-and-or-correlation-matrix'),
- ('Linking with SVD', 2, None, 'linking-with-svd')]}
+ ('Linking with SVD', 2, None, 'linking-with-svd'),
+ ('Exercises for week 37, September 6-10',
+ 2,
+ None,
+ 'exercises-for-week-37-september-6-10'),
+ ('Exercise 1: Adding Ridge and Lasso Regression',
+ 2,
+ None,
+ 'exercise-1-adding-ridge-and-lasso-regression'),
+ ('Exercise: Linear Regression for a two-dimensional function',
+ 3,
+ None,
+ 'exercise-linear-regression-for-a-two-dimensional-function')]}
end of tocinfo -->
+Small exercise: perform the standars scaling by including the standard deviation.
+
+Our data is defined by \( x\in [-3,3] \) with a total of for example \( 100 \) data points.
+
+
+
+
+where \( y \) is the function we want to fit with a given polynomial.
+
+
+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.
+
+
+
+
+
@@ -1526,10 +1614,10 @@ clf = skl.
-Lasso and Ridge regression
+What is presented here is a mathematical analysis of various regression algorithms (ordinary least squares, Ridge and Lasso Regression). The analysis is based on an important algorithm in linear algebra, the so-called Singular Value Decomposition (SVD).
The examples we have looked at so far are cases where we normally can
invert the matrix \( \boldsymbol{X}^T\boldsymbol{X} \). Using a polynomial expansion as we
-did both for the masses and the fitting of the equation of state,
-leads to row vectors of the design matrix which are essentially
-orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition.
+did both for the masses and the fitting of various functions leads to
+row vectors of the design matrix which are essentially orthogonal due
+to the polynomial character of our model. Obtaining the inverse of the
+design matrix is then often done via a so-called LU, QR or Cholesky
+decomposition.
-This may
+As we will also see in the first project,
+this may
however not the be case in general and a standard matrix inversion
algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.
@@ -1877,10 +1968,10 @@ Here we have that $${\bf X} = {\bf U}{\bf \Sigma}{\bf V}^T$$, with $$\Sigma$$ be
-Video of Lecture and handwritten notes
+Video of Lecture from 2020 and handwritten notes
More material will be added here, see handwritten notes also.
@@ -2451,6 +2542,318 @@ It is easy to generalize this to a matrix \( \boldsymbol{X}\in {\mathbb{R}}^{n\t
+More material will be added here.
+
+
+
+The exercises here are meant to prepare you for work with project 1. The first exercise is a follow-up of exercise 2 from week 36 August 30-September 3).
+
+
+
+
+
+This exercise is a continuation of exercise 2 from exercise set 1 (week 36, August 30-September 3)). 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.
+
+
+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).
+
+
+
+
+Write your own code for the Ridge method (see chapter 3.4 of Hastie 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) \).
+
+
+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.
+
+
+
+
+
+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 \).
+
+
+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
+$$ 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.
+
+
+This is a longer exercise and the aim is to study in more detail various
+regression methods, including the Ordinary Least Squares (OLS) method,
+Ridge regression and finally Lasso regression.
+This exercise forms a part of project 1.
+
+
+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.
+
+
+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 fit a
+function (for example a polynomial) of \( x \) and \( y \). Thereafter we
+will repeat much of the same procedure using 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)
+
+
+
+
+We will generate our own dataset for a function
+\( \mathrm{FrankeFunction}(x,y) \) with \( x,y \in [0,1] \). 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 and perform a standard least square regression
+analysis using polynomials in \( x \) and \( y \) up to fifth order. You can use scikit-learn as well.
+
+
+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.
+$$
+
+
+You should split your data in train and test and also consider scaling the data.
+
+
+To set up the design matrix, the following code can be used
+
+
+
+
+Write then your own code for the Ridge method or use Scikit-Learn.
+Perform the same analysis as you did for ordinary Least Squares (for the same polynomials) but now for different values of \( \lambda \). Compare and
+analyze your results with those obtained with ordinary Least Squares. Study the
+dependence on \( \lambda \).
+
+
+This part is essentially a repeat of the previous ones, but now
+with Lasso regression. Write either your own code or
+use the functionalities of Scikit-Learn (recommended).
+Give a
+critical discussion of the three methods and a judgement of which
+model fits the data best.
+
+
+
diff --git a/doc/pub/week35/ipynb/ipynb-week35-src.tar.gz b/doc/pub/week35/ipynb/ipynb-week35-src.tar.gz
index 0f8bdcb3b..ae4e62b79 100644
Binary files a/doc/pub/week35/ipynb/ipynb-week35-src.tar.gz and b/doc/pub/week35/ipynb/ipynb-week35-src.tar.gz differ
diff --git a/doc/pub/week35/ipynb/week35.ipynb b/doc/pub/week35/ipynb/week35.ipynb
index 8dc9e25b3..605aaf26d 100644
--- a/doc/pub/week35/ipynb/week35.ipynb
+++ b/doc/pub/week35/ipynb/week35.ipynb
@@ -1619,17 +1619,20 @@
"display(XPandas)\n",
"print(XPandas.mean())\n",
"print(XPandas.std())\n",
- "XPandas = XPandas -XPandas,mean()\n",
+ "XPandas = (XPandas -XPandas.mean())\n",
"display(XPandas)\n",
- "scaler = StandardScaler()\n",
- "Xscaled = scaler.transform(a)\n",
- "print(Xscaled)"
+ "scaler = StandardScaler(with_std=False)\n",
+ "scaler.fit(X)\n",
+ "Xscaled = scaler.transform(X)\n",
+ "display(XPandas-Xscaled)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
+ "Small exercise: perform the standars scaling by including the standard deviation.\n",
+ "\n",
"## Min-Max Scaling\n",
"\n",
"Another commonly used scaling method is min-max scaling. This is very\n",
@@ -1654,7 +1657,93 @@
"where $\\min(x_j)$ and $\\max(x_j)$ return the minimum and maximum value of $x_j$ over the data set, respectively.\n",
"\n",
"\n",
- "## Simple preprocessing examples, Franke function and regression"
+ "## Testing the Means Squared Error as function of Complexity\n",
+ "One of \n",
+ "the aims is to reproduce Figure 2.11 of [Hastie et al](https://github.com/CompPhysics/MLErasmus/blob/master/doc/Textbooks/elementsstat.pdf).\n",
+ "We will also use Ridge and Lasso regression. \n",
+ "\n",
+ "\n",
+ "Our data is defined by $x\\in [-3,3]$ with a total of for example $100$ data points."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "editable": true
+ },
+ "outputs": [],
+ "source": [
+ "np.random.seed()\n",
+ "n = 100\n",
+ "maxdegree = 14\n",
+ "# Make data set.\n",
+ "x = np.linspace(-3, 3, n).reshape(-1, 1)\n",
+ "y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where $y$ is the function we want to fit with a given polynomial.\n",
+ "\n",
+ "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."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "editable": true
+ },
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np\n",
+ "from sklearn.linear_model import LinearRegression, Ridge, Lasso\n",
+ "from sklearn.preprocessing import PolynomialFeatures\n",
+ "from sklearn.model_selection import train_test_split\n",
+ "from sklearn.pipeline import make_pipeline\n",
+ "\n",
+ "\n",
+ "np.random.seed(2018)\n",
+ "n = 50\n",
+ "maxdegree = 5\n",
+ "# Make data set.\n",
+ "x = np.linspace(-3, 3, n).reshape(-1, 1)\n",
+ "y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)\n",
+ "TestError = np.zeros(maxdegree)\n",
+ "TrainError = np.zeros(maxdegree)\n",
+ "polydegree = np.zeros(maxdegree)\n",
+ "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)\n",
+ "scaler = StandardScaler()\n",
+ "scaler.fit(X_train)\n",
+ "x_train_scaled = scaler.transform(x_train)\n",
+ "x_test_scaled = scaler.transform(x_test)\n",
+ "\n",
+ "for degree in range(maxdegree):\n",
+ " model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False))\n",
+ " clf = model.fit(x_train_scale,y_train)\n",
+ " y_fit = clf.predict(x_train_scaled)\n",
+ " y_pred = clf.predict(x_test_scaled) \n",
+ " polydegree[degree] = degree\n",
+ " TestError[degree] = np.mean( np.mean((y_test - y_pred)**2) )\n",
+ " TrainError[degree] = np.mean( np.mean((y_train - y_fit)**2) )\n",
+ "\n",
+ "plt.plot(polydegree, TestError, label='Test Error')\n",
+ "plt.plot(polydegree, TrainError, label='Train Error')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## More preprocessing examples, Franke function and regression"
]
},
{
@@ -1764,10 +1853,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Friday September 3\n",
- "\n",
- "Lasso and Ridge regression\n",
+ "## Mathematical Interpretation of Ordinary Least Squares\n",
"\n",
+ "What is presented here is a mathematical analysis of various regression algorithms (ordinary least squares, Ridge and Lasso Regression). The analysis is based on an important algorithm in linear algebra, the so-called Singular Value Decomposition (SVD). \n",
"\n",
"\n",
"## The singular value decomposition\n",
@@ -1775,13 +1863,15 @@
"\n",
"The examples we have looked at so far are cases where we normally can\n",
"invert the matrix $\\boldsymbol{X}^T\\boldsymbol{X}$. Using a polynomial expansion as we\n",
- "did both for the masses and the fitting of the equation of state,\n",
- "leads to row vectors of the design matrix which are essentially\n",
- "orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition. \n",
+ "did both for the masses and the fitting of various functions leads to\n",
+ "row vectors of the design matrix which are essentially orthogonal due\n",
+ "to the polynomial character of our model. Obtaining the inverse of the\n",
+ "design matrix is then often done via a so-called LU, QR or Cholesky\n",
+ "decomposition.\n",
"\n",
"\n",
- "\n",
- "This may\n",
+ "As we will also see in the first project, \n",
+ "this may\n",
"however not the be case in general and a standard matrix inversion\n",
"algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.\n",
"\n",
@@ -2256,9 +2346,9 @@
"\n",
"\n",
"\n",
- "## Friday September 12\n",
+ "## Friday September 3\n",
"\n",
- "[Video of Lecture](https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h20/forelesningsvideoer/LectureSeptember11.mp4?vrtx=view-as-webpage) and [handwritten notes](https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/NotesSeptember11.pdf)\n",
+ "[Video of Lecture from 2020](https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h20/forelesningsvideoer/LectureSeptember11.mp4?vrtx=view-as-webpage) and [handwritten notes](https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/NotesSeptember11.pdf)\n",
"\n",
"More material will be added here, see handwritten notes also.\n",
"\n",
@@ -3152,7 +3242,436 @@
"It is easy to generalize this to a matrix $\\boldsymbol{X}\\in {\\mathbb{R}}^{n\\times p}$.\n",
"\n",
"\n",
- "## Linking with SVD"
+ "## Linking with SVD\n",
+ "\n",
+ "More material will be added here.\n",
+ "\n",
+ "\n",
+ "## Exercises for week 37, September 6-10\n",
+ "\n",
+ "The exercises here are meant to prepare you for work with project 1. The first exercise is a follow-up of exercise 2 from week 36 August 30-September 3).\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Exercise 1: Adding Ridge and Lasso Regression\n",
+ "\n",
+ "This exercise is a continuation of exercise 2 from exercise set 1 (week 36, August 30-September 3)). We will\n",
+ "use the same function to generate our data set, still staying with a\n",
+ "simple function $y(x)$ which we want to fit using linear regression,\n",
+ "but now extending the analysis to include the Ridge and the Lasso\n",
+ "regression methods. \n",
+ "\n",
+ "We will thus again generate our own dataset for a function $y(x)$ where \n",
+ "$x \\in [0,1]$ and defined by random numbers computed with the uniform\n",
+ "distribution. The function $y$ is a quadratic polynomial in $x$ with\n",
+ "added stochastic noise according to the normal distribution $\\cal{N}(0,1)$.\n",
+ "\n",
+ "The following simple Python instructions define our $x$ and $y$ values (with 100 data points)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "editable": true
+ },
+ "outputs": [],
+ "source": [
+ "x = np.random.rand(100)\n",
+ "y = 2.0+5*x*x+0.1*np.random.randn(100)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Write your own code for the Ridge method (see chapter 3.4 of Hastie *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)$. \n",
+ "\n",
+ "\n",
+ "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."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "editable": true
+ },
+ "outputs": [],
+ "source": [
+ "import os\n",
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "import matplotlib.pyplot as plt\n",
+ "from sklearn.model_selection import train_test_split\n",
+ "from sklearn.preprocessing import StandardScaler\n",
+ "\n",
+ "def R2(y_data, y_model):\n",
+ " return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n",
+ "def MSE(y_data,y_model):\n",
+ " n = np.size(y_model)\n",
+ " return np.sum((y_data-y_model)**2)/n\n",
+ "\n",
+ "\n",
+ "# A seed just to ensure that the random numbers are the same for every run.\n",
+ "# Useful for eventual debugging.\n",
+ "np.random.seed(3155)\n",
+ "\n",
+ "x = np.random.rand(100)\n",
+ "y = 2.0+5*x*x+0.1*np.random.randn(100)\n",
+ "\n",
+ "# number of features p (here degree of polynomial\n",
+ "p = 3\n",
+ "# The design matrix now as function of a given polynomial\n",
+ "X = np.zeros((len(x),p))\n",
+ "X[:,0] = 1.0\n",
+ "X[:,1] = x\n",
+ "X[:,2] = x*x\n",
+ "# We split the data in test and training data\n",
+ "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
+ "scaler = StandardScaler()\n",
+ "scaler.fit(X_train)\n",
+ "X_train_scaled = scaler.transform(X_train)\n",
+ "X_test_scaled = scaler.transform(X_test)\n",
+ "\n",
+ "# matrix inversion to find beta\n",
+ "OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train\n",
+ "print(OLSbeta)\n",
+ "# and then make the prediction\n",
+ "ytildeOLS = X_train @ OLSbeta\n",
+ "print(\"Training R2 for OLS\")\n",
+ "print(R2(y_train,ytildeOLS))\n",
+ "print(\"Training MSE for OLS\")\n",
+ "print(MSE(y_train,ytildeOLS))\n",
+ "ypredictOLS = X_test @ OLSbeta\n",
+ "print(\"Test R2 for OLS\")\n",
+ "print(R2(y_test,ypredictOLS))\n",
+ "print(\"Test MSE OLS\")\n",
+ "print(MSE(y_test,ypredictOLS))\n",
+ "\n",
+ "# Repeat now for Ridge regression and various values of the regularization parameter\n",
+ "I = np.eye(p,p)\n",
+ "# Decide which values of lambda to use\n",
+ "nlambdas = 20\n",
+ "MSEPredict = np.zeros(nlambdas)\n",
+ "MSETrain = np.zeros(nlambdas)\n",
+ "lambdas = np.logspace(-4, 1, nlambdas)\n",
+ "for i in range(nlambdas):\n",
+ " lmb = lambdas[i]\n",
+ " Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train\n",
+ " # and then make the prediction\n",
+ " ytildeRidge = X_train @ Ridgebeta\n",
+ " ypredictRidge = X_test @ Ridgebeta\n",
+ " MSEPredict[i] = MSE(y_test,ypredictRidge)\n",
+ " MSETrain[i] = MSE(y_train,ytildeRidge)\n",
+ "# Now plot the results\n",
+ "plt.figure()\n",
+ "plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train')\n",
+ "plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')\n",
+ "plt.xlabel('log10(lambda)')\n",
+ "plt.ylabel('MSE')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "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$. \n",
+ "\n",
+ "\n",
+ "\n",
+ "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"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "MSE(\\hat{y},\\hat{\\tilde{y}}) = \\frac{1}{n}\n",
+ "\\sum_{i=0}^{n-1}(y_i-\\tilde{y}_i)^2,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and the $R^2$ score function.\n",
+ "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"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "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},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where we have defined the mean value of $\\hat{y}$ as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\bar{y} = \\frac{1}{n} \\sum_{i=0}^{n - 1} y_i.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Discuss these quantities as functions of the variable $\\lambda$ in the Ridge and Lasso regression methods.\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "### Exercise: Linear Regression for a two-dimensional function\n",
+ "\n",
+ "This is a longer exercise and the aim is to study in more detail various\n",
+ "regression methods, including the Ordinary Least Squares (OLS) method,\n",
+ "Ridge regression and finally Lasso regression.\n",
+ "This exercise forms a part of project 1.\n",
+ "\n",
+ "We will study how to fit polynomials to a specific\n",
+ "two-dimensional function called [Franke's\n",
+ "function](http://www.dtic.mil/dtic/tr/fulltext/u2/a081688.pdf). This\n",
+ "is a function which has been widely used when testing various\n",
+ "interpolation and fitting algorithms. \n",
+ "\n",
+ "The Franke function, which is a weighted sum of four exponentials reads as follows"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\begin{align*}\n",
+ "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)} \\\\\n",
+ "&+\\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) }.\n",
+ "\\end{align*}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The function will be defined for $x,y\\in [0,1]$. Our first step will\n",
+ "be to perform an OLS regression analysis of this function, trying out\n",
+ "a polynomial fit with an $x$ and $y$ dependence of the form $[x, y,\n",
+ "x^2, y^2, xy, \\dots]$. We will fit a\n",
+ "function (for example a polynomial) of $x$ and $y$. Thereafter we\n",
+ "will repeat much of the same procedure using the Ridge and Lasso\n",
+ "regression methods, introducing thus a dependence on the bias\n",
+ "(penalty) $\\lambda$.\n",
+ "\n",
+ "\n",
+ "The Python fucntion for the Franke function is included here (it performs also a three-dimensional plot of it)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "editable": true
+ },
+ "outputs": [],
+ "source": [
+ "from mpl_toolkits.mplot3d import Axes3D\n",
+ "import matplotlib.pyplot as plt\n",
+ "from matplotlib import cm\n",
+ "from matplotlib.ticker import LinearLocator, FormatStrFormatter\n",
+ "import numpy as np\n",
+ "from random import random, seed\n",
+ "\n",
+ "fig = plt.figure()\n",
+ "ax = fig.gca(projection='3d')\n",
+ "\n",
+ "# Make data.\n",
+ "x = np.arange(0, 1, 0.05)\n",
+ "y = np.arange(0, 1, 0.05)\n",
+ "x, y = np.meshgrid(x,y)\n",
+ "\n",
+ "\n",
+ "def FrankeFunction(x,y):\n",
+ " term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))\n",
+ " term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))\n",
+ " term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))\n",
+ " term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)\n",
+ " return term1 + term2 + term3 + term4\n",
+ "\n",
+ "\n",
+ "z = FrankeFunction(x, y)\n",
+ "\n",
+ "# Plot the surface.\n",
+ "surf = ax.plot_surface(x, y, z, cmap=cm.coolwarm,\n",
+ " linewidth=0, antialiased=False)\n",
+ "\n",
+ "# Customize the z axis.\n",
+ "ax.set_zlim(-0.10, 1.40)\n",
+ "ax.zaxis.set_major_locator(LinearLocator(10))\n",
+ "ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))\n",
+ "\n",
+ "# Add a color bar which maps values to colors.\n",
+ "fig.colorbar(surf, shrink=0.5, aspect=5)\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We will generate our own dataset for a function\n",
+ "$\\mathrm{FrankeFunction}(x,y)$ with $x,y \\in [0,1]$. The function\n",
+ "$f(x,y)$ is the Franke function. You should explore also the addition\n",
+ "an added stochastic noise to this function using the normal\n",
+ "distribution $\\cal{N}(0,1)$.\n",
+ "\n",
+ "Write your own code (using either a matrix inversion or a singular\n",
+ "value decomposition from e.g., **numpy** ) or use your code and perform a standard least square regression\n",
+ "analysis using polynomials in $x$ and $y$ up to fifth order. You can use **scikit-learn** as well.\n",
+ "\n",
+ "\n",
+ "Evaluate the Mean Squared error (MSE)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "MSE(\\hat{y},\\hat{\\tilde{y}}) = \\frac{1}{n}\n",
+ "\\sum_{i=0}^{n-1}(y_i-\\tilde{y}_i)^2,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and the $R^2$ score function. If $\\tilde{\\hat{y}}_i$ is the predicted\n",
+ "value of the $i-th$ sample and $y_i$ is the corresponding true value,\n",
+ "then the score $R^2$ is defined as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "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},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where we have defined the mean value of $\\hat{y}$ as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\bar{y} = \\frac{1}{n} \\sum_{i=0}^{n - 1} y_i.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "You should split your data in train and test and also consider scaling the data.\n",
+ "\n",
+ "To set up the design matrix, the following code can be used"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "editable": true
+ },
+ "outputs": [],
+ "source": [
+ "def FrankeFunction(x,y):\n",
+ "\tterm1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))\n",
+ "\tterm2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))\n",
+ "\tterm3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))\n",
+ "\tterm4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)\n",
+ "\treturn term1 + term2 + term3 + term4\n",
+ "\n",
+ "\n",
+ "def create_X(x, y, n ):\n",
+ "\tif len(x.shape) > 1:\n",
+ "\t\tx = np.ravel(x)\n",
+ "\t\ty = np.ravel(y)\n",
+ "\n",
+ "\tN = len(x)\n",
+ "\tl = int((n+1)*(n+2)/2)\t\t# Number of elements in beta\n",
+ "\tX = np.ones((N,l))\n",
+ "\n",
+ "\tfor i in range(1,n+1):\n",
+ "\t\tq = int((i)*(i+1)/2)\n",
+ "\t\tfor k in range(i+1):\n",
+ "\t\t\tX[:,q+k] = (x**(i-k))*(y**k)\n",
+ "\n",
+ "\treturn X\n",
+ "\n",
+ "\n",
+ "# Making meshgrid of datapoints and compute Franke's function\n",
+ "n = 5\n",
+ "N = 1000\n",
+ "x = np.sort(np.random.uniform(0, 1, N))\n",
+ "y = np.sort(np.random.uniform(0, 1, N))\n",
+ "z = FrankeFunction(x, y)\n",
+ "X = create_X(x, y, n=n)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Write then your own code for the Ridge method or use **Scikit-Learn**.\n",
+ "Perform the same analysis as you did for ordinary Least Squares (for the same polynomials) but now for different values of $\\lambda$. Compare and\n",
+ "analyze your results with those obtained with ordinary Least Squares. Study the\n",
+ "dependence on $\\lambda$.\n",
+ "\n",
+ "\n",
+ "This part is essentially a repeat of the previous ones, but now\n",
+ "with Lasso regression. Write either your own code or\n",
+ "use the functionalities of **Scikit-Learn** (recommended). \n",
+ "Give a\n",
+ "critical discussion of the three methods and a judgement of which\n",
+ "model fits the data best.\n",
+ "\n",
+ ""
]
}
],
diff --git a/doc/src/week35/week35.do.txt b/doc/src/week35/week35.do.txt
index 80ce1755a..92f3f9690 100644
--- a/doc/src/week35/week35.do.txt
+++ b/doc/src/week35/week35.do.txt
@@ -965,14 +965,15 @@ XPandas = pd.DataFrame(X)
display(XPandas)
print(XPandas.mean())
print(XPandas.std())
-XPandas = XPandas -XPandas,mean()
+XPandas = (XPandas -XPandas.mean())
display(XPandas)
-scaler = StandardScaler()
-Xscaled = scaler.transform(a)
-print(Xscaled)
+scaler = StandardScaler(with_std=False)
+scaler.fit(X)
+Xscaled = scaler.transform(X)
+display(XPandas-Xscaled)
!ec
-
+Small exercise: perform the standars scaling by including the standard deviation.
!split
===== Min-Max Scaling =====
@@ -991,7 +992,67 @@ where $\min(x_j)$ and $\max(x_j)$ return the minimum and maximum value of $x_j$
!split
-===== Simple preprocessing examples, Franke function and regression =====
+===== Testing the Means Squared Error as function of Complexity =====
+One of
+the aims is to reproduce Figure 2.11 of "Hastie et al":"https://github.com/CompPhysics/MLErasmus/blob/master/doc/Textbooks/elementsstat.pdf".
+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.
+!bc pycod
+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)
+!ec
+where $y$ is the function we want to fit with a given polynomial.
+
+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.
+!bc pycod
+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()
+!ec
+
+
+
+!split
+===== More preprocessing examples, Franke function and regression =====
!bc pycod
# Common imports
@@ -1091,10 +1152,9 @@ print("R2 score for scaled data: {:.2f}".format(clf.score(X_test_scaled,y_test)
!split
-===== Friday September 3 =====
-
-Lasso and Ridge regression
+===== Mathematical Interpretation of Ordinary Least Squares =====
+What is presented here is a mathematical analysis of various regression algorithms (ordinary least squares, Ridge and Lasso Regression). The analysis is based on an important algorithm in linear algebra, the so-called Singular Value Decomposition (SVD).
!split
@@ -1104,13 +1164,15 @@ Lasso and Ridge regression
The examples we have looked at so far are cases where we normally can
invert the matrix $\bm{X}^T\bm{X}$. Using a polynomial expansion as we
-did both for the masses and the fitting of the equation of state,
-leads to row vectors of the design matrix which are essentially
-orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition.
+did both for the masses and the fitting of various functions leads to
+row vectors of the design matrix which are essentially orthogonal due
+to the polynomial character of our model. Obtaining the inverse of the
+design matrix is then often done via a so-called LU, QR or Cholesky
+decomposition.
-
-This may
+As we will also see in the first project,
+this may
however not the be case in general and a standard matrix inversion
algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.
@@ -1419,9 +1481,9 @@ Here we have that $${\bf X} = {\bf U}{\bf \Sigma}{\bf V}^T$$, with $$\Sigma$$ be
!split
-===== Friday September 12 =====
+===== Friday September 3 =====
-"Video of Lecture":"https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h20/forelesningsvideoer/LectureSeptember11.mp4?vrtx=view-as-webpage" and "handwritten notes":"https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/NotesSeptember11.pdf"
+"Video of Lecture from 2020":"https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h20/forelesningsvideoer/LectureSeptember11.mp4?vrtx=view-as-webpage" and "handwritten notes":"https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/NotesSeptember11.pdf"
More material will be added here, see handwritten notes also.
@@ -1980,4 +2042,327 @@ It is easy to generalize this to a matrix $\bm{X}\in {\mathbb{R}}^{n\times p}$.
!split
===== Linking with SVD =====
+More material will be added here.
+
+
+!split
+===== Exercises for week 37, September 6-10 =====
+
+The exercises here are meant to prepare you for work with project 1. The first exercise is a follow-up of exercise 2 from week 36 August 30-September 3).
+
+
+===== Exercise: Adding Ridge and Lasso Regression =====
+
+
+This exercise is a continuation of exercise 2 from exercise set 1 (week 36, August 30-September 3)). 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.
+
+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).
+!bc pycod
+x = np.random.rand(100)
+y = 2.0+5*x*x+0.1*np.random.randn(100)
+!ec
+
+
+Write your own code for the Ridge method (see chapter 3.4 of Hastie *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)$.
+
+
+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.
+
+!bc pycod
+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 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.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+!ec
+
+
+
+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$.
+
+
+
+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
+!bt
+\[ MSE(\hat{y},\hat{\tilde{y}}) = \frac{1}{n}
+\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2,
+\]
+!et
+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
+!bt
+\[
+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},
+\]
+!et
+where we have defined the mean value of $\hat{y}$ as
+!bt
+\[
+\bar{y} = \frac{1}{n} \sum_{i=0}^{n - 1} y_i.
+\]
+!et
+Discuss these quantities as functions of the variable $\lambda$ in the Ridge and Lasso regression methods.
+
+
+
+
+
+
+
+=== Exercise: Linear Regression for a two-dimensional function ===
+
+This is a longer exercise and the aim is to study in more detail various
+regression methods, including the Ordinary Least Squares (OLS) method,
+Ridge regression and finally Lasso regression.
+This exercise forms a part of project 1.
+
+We will study how to fit polynomials to a specific
+two-dimensional function called "Franke's
+function":"http://www.dtic.mil/dtic/tr/fulltext/u2/a081688.pdf". This
+is a function which has been widely used when testing various
+interpolation and fitting algorithms.
+
+The Franke function, which is a weighted sum of four exponentials reads as follows
+!bt
+\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*}
+!et
+
+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 fit a
+function (for example a polynomial) of $x$ and $y$. Thereafter we
+will repeat much of the same procedure using 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)
+!bc pycod
+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()
+
+!ec
+
+
+
+We will generate our own dataset for a function
+$\mathrm{FrankeFunction}(x,y)$ with $x,y \in [0,1]$. 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 and perform a standard least square regression
+analysis using polynomials in $x$ and $y$ up to fifth order. You can use _scikit-learn_ as well.
+
+
+Evaluate the Mean Squared error (MSE)
+
+!bt
+\[ MSE(\hat{y},\hat{\tilde{y}}) = \frac{1}{n}
+\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2,
+\]
+!et
+
+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
+
+!bt
+\[
+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},
+\]
+!et
+
+where we have defined the mean value of $\hat{y}$ as
+
+!bt
+\[
+\bar{y} = \frac{1}{n} \sum_{i=0}^{n - 1} y_i.
+\]
+!et
+
+
+You should split your data in train and test and also consider scaling the data.
+
+To set up the design matrix, the following code can be used
+!bc pycod
+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
+
+
+def create_X(x, y, n ):
+ if len(x.shape) > 1:
+ x = np.ravel(x)
+ y = np.ravel(y)
+
+ N = len(x)
+ l = int((n+1)*(n+2)/2) # Number of elements in beta
+ X = np.ones((N,l))
+
+ for i in range(1,n+1):
+ q = int((i)*(i+1)/2)
+ for k in range(i+1):
+ X[:,q+k] = (x**(i-k))*(y**k)
+
+ return X
+
+
+# Making meshgrid of datapoints and compute Franke's function
+n = 5
+N = 1000
+x = np.sort(np.random.uniform(0, 1, N))
+y = np.sort(np.random.uniform(0, 1, N))
+z = FrankeFunction(x, y)
+X = create_X(x, y, n=n)
+!ec
+
+
+
+Write then your own code for the Ridge method or use _Scikit-Learn_.
+Perform the same analysis as you did for ordinary Least Squares (for the same polynomials) but now for different values of $\lambda$. Compare and
+analyze your results with those obtained with ordinary Least Squares. Study the
+dependence on $\lambda$.
+
+
+This part is essentially a repeat of the previous ones, but now
+with Lasso regression. Write either your own code or
+use the functionalities of _Scikit-Learn_ (recommended).
+Give a
+critical discussion of the three methods and a judgement of which
+model fits the data best.
+
+
+
+
+
+
+
+
+
+
+
+
Simple preprocessing examples, Franke function and regression
+Testing the Means Squared Error as function of Complexity
+One of
+the aims is to reproduce Figure 2.11 of Hastie et al.
+We will also use Ridge and Lasso regression.
+
+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)
+
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()
+
More preprocessing examples, Franke function and regression
Friday September 3
+Mathematical Interpretation of Ordinary Least Squares
Friday September 12
+Friday September 3
Linking with SVD
+
+Exercises for week 37, September 6-10
+
+Exercise 1: Adding Ridge and Lasso Regression
+
+x = np.random.rand(100)
+y = 2.0+5*x*x+0.1*np.random.randn(100)
+
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 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.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+$$ 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.
+
+Exercise: Linear Regression for a two-dimensional function
+
+
+$$
+\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*}
+$$
+
+
+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()
+
+$$ MSE(\hat{y},\hat{\tilde{y}}) = \frac{1}{n}
+\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2,
+$$
+
+
+
+$$
+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},
+$$
+
+
+
+$$
+\bar{y} = \frac{1}{n} \sum_{i=0}^{n - 1} y_i.
+$$
+
+
+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
+
+
+def create_X(x, y, n ):
+ if len(x.shape) > 1:
+ x = np.ravel(x)
+ y = np.ravel(y)
+
+ N = len(x)
+ l = int((n+1)*(n+2)/2) # Number of elements in beta
+ X = np.ones((N,l))
+
+ for i in range(1,n+1):
+ q = int((i)*(i+1)/2)
+ for k in range(i+1):
+ X[:,q+k] = (x**(i-k))*(y**k)
+
+ return X
+
+
+# Making meshgrid of datapoints and compute Franke's function
+n = 5
+N = 1000
+x = np.sort(np.random.uniform(0, 1, N))
+y = np.sort(np.random.uniform(0, 1, N))
+z = FrankeFunction(x, y)
+X = create_X(x, y, n=n)
+
@@ -1420,7 +1443,72 @@ where \( \min(x_j) \) and \( \max(x_j) \) return the minimum and maximum value o
-Simple preprocessing examples, Franke function and regression
+Testing the Means Squared Error as function of Complexity
+One of
+the aims is to reproduce Figure 2.11 of Hastie et al.
+We will also use Ridge and Lasso regression.
+
+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)
+
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()
+
+
+More preprocessing examples, Franke function and regression
-Friday September 3
+Mathematical Interpretation of Ordinary Least Squares
@@ -1539,12 +1627,15 @@ Lasso and Ridge regression
-Friday September 12
+Friday September 3
Linking with SVD
+
+Exercises for week 37, September 6-10
+
+Exercise 1: Adding Ridge and Lasso Regression
+
+x = np.random.rand(100)
+y = 2.0+5*x*x+0.1*np.random.randn(100)
+
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 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.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
Exercise: Linear Regression for a two-dimensional function
+
+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()
+
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
+
+
+def create_X(x, y, n ):
+ if len(x.shape) > 1:
+ x = np.ravel(x)
+ y = np.ravel(y)
+
+ N = len(x)
+ l = int((n+1)*(n+2)/2) # Number of elements in beta
+ X = np.ones((N,l))
+
+ for i in range(1,n+1):
+ q = int((i)*(i+1)/2)
+ for k in range(i+1):
+ X[:,q+k] = (x**(i-k))*(y**k)
+
+ return X
+
+
+# Making meshgrid of datapoints and compute Franke's function
+n = 5
+N = 1000
+x = np.sort(np.random.uniform(0, 1, N))
+y = np.sort(np.random.uniform(0, 1, N))
+z = FrankeFunction(x, y)
+X = create_X(x, y, n=n)
+
@@ -1425,7 +1448,72 @@ where \( \min(x_j) \) and \( \max(x_j) \) return the minimum and maximum value o
-Simple preprocessing examples, Franke function and regression
+Testing the Means Squared Error as function of Complexity
+One of
+the aims is to reproduce Figure 2.11 of Hastie et al.
+We will also use Ridge and Lasso regression.
+
+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)
+
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()
+
+
+More preprocessing examples, Franke function and regression
-Friday September 3
+Mathematical Interpretation of Ordinary Least Squares
@@ -1544,12 +1632,15 @@ Lasso and Ridge regression
-Friday September 12
+Friday September 3
Linking with SVD
+
+Exercises for week 37, September 6-10
+
+Exercise 1: Adding Ridge and Lasso Regression
+
+x = np.random.rand(100)
+y = 2.0+5*x*x+0.1*np.random.randn(100)
+
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 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.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
Exercise: Linear Regression for a two-dimensional function
+
+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()
+
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
+
+
+def create_X(x, y, n ):
+ if len(x.shape) > 1:
+ x = np.ravel(x)
+ y = np.ravel(y)
+
+ N = len(x)
+ l = int((n+1)*(n+2)/2) # Number of elements in beta
+ X = np.ones((N,l))
+
+ for i in range(1,n+1):
+ q = int((i)*(i+1)/2)
+ for k in range(i+1):
+ X[:,q+k] = (x**(i-k))*(y**k)
+
+ return X
+
+
+# Making meshgrid of datapoints and compute Franke's function
+n = 5
+N = 1000
+x = np.sort(np.random.uniform(0, 1, N))
+y = np.sort(np.random.uniform(0, 1, N))
+z = FrankeFunction(x, y)
+X = create_X(x, y, n=n)
+