lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2021-09-13 14:00:16 +02:00
parent 1af325e7fb
commit 30d664a82b
9 changed files with 2392 additions and 942 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/BookChapters/chapter1.dlog
+27
View File
@@ -7,3 +7,30 @@ Translating doconce text in chapter1.do.txt to ipynb
Failed to remove ans_at_end environment
Failed to remove sol_at_end environment
output in chapter1.ipynb
Translating doconce text in chapter1.do.txt to ipynb
*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
Failed to remove ans_at_end environment
Failed to remove sol_at_end environment
output in chapter1.ipynb
Translating doconce text in chapter1.do.txt to ipynb
*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
Failed to remove ans_at_end environment
Failed to remove sol_at_end environment
output in chapter1.ipynb
Translating doconce text in chapter1.do.txt to ipynb
*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
Failed to remove ans_at_end environment
Failed to remove sol_at_end environment
output in chapter1.ipynb
doc/BookChapters/chapter1.do.txt
+159 -140
View File
@@ -7,38 +7,37 @@
Our emphasis throughout this series of lectures
is on understanding the mathematical aspects of
different algorithms used in the fields of data analysis and machine learning.
Our emphasis throughout this series of lectures is on understanding
the mathematical aspects of different algorithms used in the fields of
data analysis and machine learning.
However, where possible we will emphasize the
importance of using available software. We start thus with a hands-on
and top-down approach to machine learning. The aim is thus to start with
relevant data or data we have produced
and use these to introduce statistical data analysis
concepts and machine learning algorithms before we delve into the
algorithms themselves. The examples we will use in the beginning, start with simple
polynomials with random noise added. We will use the Python
software package "Scikit-Learn":"http://scikit-learn.org/stable/" and
introduce various machine learning algorithms to make fits of
the data and predictions. We move thereafter to more interesting
cases such as data from say experiments (below we will look at experimental nuclear binding energies as an example).
These are examples where we can easily set up the data and
then use machine learning algorithms included in for example
_Scikit-Learn_.
However, where possible we will emphasize the importance of using
available software. We start thus with a hands-on and top-down
approach to machine learning. The aim is thus to start with relevant
data or data we have produced and use these to introduce statistical
data analysis concepts and machine learning algorithms before we delve
into the algorithms themselves. The examples we will use in the
beginning, start with simple polynomials with random noise added. We
will use the Python software package
"Scikit-Learn":"http://scikit-learn.org/stable/" and introduce various
machine learning algorithms to make fits of the data and
predictions. We move thereafter to more interesting cases such as data
from say experiments (below we will look at experimental nuclear
binding energies as an example). These are examples where we can
easily set up the data and then use machine learning algorithms
included in for example _Scikit-Learn_.
These examples will serve us the purpose of getting
started. Furthermore, they allow us to catch more than two birds with
a stone. They will allow us to bring in some programming specific
topics and tools as well as showing the power of various Python
libraries for machine learning and statistical data analysis.
topics and tools as well as showing the power of various Python
libraries for machine learning and statistical data analysis.
Here, we will mainly focus on two
specific Python packages for Machine Learning, Scikit-Learn and
Tensorflow (see below for links etc). Moreover, the examples we
introduce will serve as inputs to many of our discussions later, as
well as allowing you to set up models and produce your own data and
get started with programming.
Here, we will mainly focus on two specific Python packages for Machine
Learning, Scikit-Learn and Tensorflow (see below for links etc).
Moreover, the examples we introduce will serve as inputs to many of
our discussions later, as well as allowing you to set up models and
produce your own data and get started with programming.
@@ -275,7 +274,17 @@ quality figures. Feel free to explore the extensive
"gallery":"https://matplotlib.org/gallery/index.html" of examples. In
this example we plot our original values of $x$ and $y$ as well as the
prediction _ypredict_ ($\tilde{y}$), which attempts at fitting our
data with a straight line.
data with a straight line. Note also that _Scikit-Learn_ requires a
matrix as input for the input values $x$ and $y$. In the above code we
have solved this by declaring $x$ and $y$ as arrays of dimension
$n\times 1$.
In the code here we have also made a new array for $x\in [0,1]$. Our
prediction is computed for these values, meaning that they were not
included in the data set used to *train* (or fit) the model.
This is a recurrring theme in machine learning and data analysis. We would like to train a model on a specific given data set.
Thereafter we wish to apply it to data which were not included in the training. Below we will encounter this again in the so-called *train-validate-test* spliting. We will typically split our data into different sets, oen for training, one for validation and finally, our data from the untouched test vault!
The Python code follows here.
!bc pycod
@@ -288,6 +297,7 @@ x = np.random.rand(100,1)
y = 2*x+np.random.randn(100,1)
linreg = LinearRegression()
linreg.fit(x,y)
# This is our new x-array to which we test our model
xnew = np.array([[0],[1]])
ypredict = linreg.predict(xnew)
@@ -394,8 +404,7 @@ plt.show()
Depending on the parameter in front of the normal distribution, we may
have a small or larger relative error. Try to play around with
different training data sets and study (graphically) the value of the
relative error. Note also that _Scikit-Learn_ requires a matrix as input for the input values $x$ and $y$. In the above code we have
solved this by declaring $x$ and $y$ as arrays of dimension $n\times 1$.
relative error.
As mentioned above, _Scikit-Learn_ has an impressive functionality.
We can for example extract the values of $\alpha$ and $\beta$ and
@@ -756,7 +765,7 @@ Energies = Masses['Ebinding']
print(Masses)
!ec
The next step, and we will define this mathematically later, is to set up the so-called _design matrix_. We will throughout call this matrix $\bm{X}$.
It has dimensionality $p\times n$, where $n$ is the number of data points and $p$ are the so-called predictors. In our case here they are given by the number of polynomials in $A$ we wish to include in the fit.
It has dimensionality $n\times p$, where $n$ is the number of data points and $p$ are the so-called predictors. In our case here they are given by the number of polynomials in $A$ we wish to include in the fit.
!bc pycod
# Now we set up the design matrix X
X = np.zeros((len(A),5))
@@ -767,9 +776,12 @@ X[:,3] = A**(-1.0/3.0)
X[:,4] = A**(-1.0)
!ec
Note well that we have made life simple here. We perform a fit in terms of the number of nucleons only. A more sophisticated fit can be done by including an explicit dependence on the number of protons and neutrons in the asymmetry and Coulomb terms.
Note well that we have made life simple here. We perform a fit in
terms of the number of nucleons only. A more sophisticated fit can be
done by including an explicit dependence on the number of protons and
neutrons in the asymmetry and Coulomb terms. We leave this as an exercise to you the reader.
With _scikitlearn_ we are now ready to use linear regression and fit our data.
With _Scikit-Learn_ we are now ready to use linear regression and fit our data.
!bc pycod
clf = skl.LinearRegression().fit(X, Energies)
fity = clf.predict(X)
@@ -799,7 +811,7 @@ plt.show()
!ec
As a teaser, let us now see how we can do this with decision trees using _scikit-learn_. Later we will switch to so-called _random forests_!
As a teaser, let us now see how we can do this with decision trees using _Scikit-Learn_. Later we will switch to so-called _random forests_!
!bc pycod
@@ -835,6 +847,16 @@ print(Masses)
print(np.mean( (Energies-y_1)**2))
!ec
With a deeper and deeper tree level, we can almost reproduce every
single data point by increasing the max depth of the tree.
We can actually decide to make a decision tree which fits every single point.
As we will
see later, this has the benefit that we can really train a model which
traverses every single data point. However, the price we pay is that
we will easily overfit. That is, if we apply our model to unseen data,
we will most likely fail miserably in our attempt at making
predictions. As an exercise, try to make the tree level larger by adjusting the maximum depth variable. When printing out the predicition, you will note that the binding energy of every nucleus is accurately reproduced.
The _seaborn_ package allows us to visualize data in an efficient way. Note that we use _scikit-learn_'s multi-layer perceptron (or feed forward neural network)
functionality.
@@ -1270,12 +1292,12 @@ matrices as upper case boldfaced letters.
!bt
\[
\frac{\partial (\bm{b}^T\bm{a})}{\partial\bm{a}}=\bm{b},
\frac{\partial\bm{b}^T\bm{a}}{\partial\bm{a}}=\bm{b},
\]
!et
!bt
\[
\frac{\partial (\bm{a}^T\bm{A}\bm{a})}{\partial\bm{a}}=(\bm{A}+\bm{A}^T)\bm{a},
\frac{\partial\bm{a}^T\bm{A}\bm{a}}{\partial\bm{a}}=(\bm{A}+\bm{A}^T)\bm{a},
\]
!et
!bt
@@ -1289,6 +1311,30 @@ matrices as upper case boldfaced letters.
\]
!et
We can then compute the second derivative of the cost function, which in our case is the second derivative
of the means squared error. This leads to
!bt
\[
\frac{\partial^2 C(\bm{\beta})}{\partial \bm{\beta}^T\partial \bm{\beta}} =\frac{2}{n}\bm{X}^T\bm{X}.
\]
!et
This quantity defines was what is called the Hessian matrix (the second derivative of a function we want to optimize).
The Hessian matrix plays an important role and is defined for the mean squared error as
!bt
\[
\bm{H}=\bm{X}^T\bm{X}.
\]
!et
The Hessian matrix for ordinary least squares is also proportional to
the covariance matrix. As we will see in the chapter on Ridge and Lasso regression, This means that we can use the Singular Value Decomposition of a matrix to find
the eigenvalues of the covariance matrix and the Hessian matrix in
terms of the singular values.
The residuals $\bm{\epsilon}$ are in turn given by
!bt
@@ -1546,9 +1592,8 @@ before, with the same initializations and declarations. We use also
_pandas_ again, rather extensively in order to organize our data.
The difference now is that we use _Scikit-Learn's_ regression tools
instead of our own matrix inversion implementation. Furthermore, we
sneak in _Ridge_ regression (to be discussed below) which includes a
hyperparameter $\lambda$, also to be explained below.
instead of our own matrix inversion implementation.
!bc pycod
@@ -1612,18 +1657,6 @@ print('Variance score: %.2f' % r2_score(Energies, ytilde))
print('Mean absolute error: %.2f' % mean_absolute_error(Energies, ytilde))
print(clf.coef_, clf.intercept_)
# The Ridge regression with a hyperparameter lambda = 0.1
_lambda = 0.1
clf_ridge = skl.Ridge(alpha=_lambda).fit(X, Energies)
yridge = clf_ridge.predict(X)
EoS['Eridge'] = yridge
# The mean squared error
print("Mean squared error: %.2f" % mean_squared_error(Energies, yridge))
# Explained variance score: 1 is perfect prediction
print('Variance score: %.2f' % r2_score(Energies, yridge))
# Mean absolute error
print('Mean absolute error: %.2f' % mean_absolute_error(Energies, yridge))
print(clf_ridge.coef_, clf_ridge.intercept_)
fig, ax = plt.subplots()
ax.set_xlabel(r'$\rho[\mathrm{fm}^{-3}]$')
@@ -1632,8 +1665,6 @@ ax.plot(EoS['Density'], EoS['Energy'], alpha=0.7, lw=2,
label='Theoretical data')
ax.plot(EoS['Density'], EoS['Eols'], alpha=0.7, lw=2, c='m',
label='OLS')
ax.plot(EoS['Density'], EoS['Eridge'], alpha=0.7, lw=2, c='g',
label='Ridge $\lambda = 0.1$')
ax.legend()
save_fig("EoSfitting")
plt.show()
@@ -1642,13 +1673,11 @@ plt.show()
The above simple polynomial in density $\rho$ gives an excellent fit
to the data.
We note also that there is a small deviation between the
standard OLS and the Ridge regression at higher densities. We discuss this in more detail
below.
===== Splitting our Data in Training and Test data =====
It is normal in essentially all Machine Learning studies to split the
data in a training set and a test set (sometimes also an additional
validation set). _Scikit-Learn_ has an own function for this. There
@@ -1660,6 +1689,81 @@ our discussion of the so-called _bias-variance_ tradeoff. Here we
limit ourselves to repeat the above equation of state fitting example
but now splitting the data into a training set and a test set.
Let us study some examples. The first code here takes a simple
one-dimensional second-order polynomial and we fit it to a
second-order polynomial. Depending on the strength of the added noise,
the various measures like the $R2$ score or the mean-squared error,
the fit becomes better or worse.
!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
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
x = np.random.rand(100)
y = 2.0+5*x*x+0.1*np.random.randn(100)
# The design matrix now as function of a given polynomial
X = np.zeros((len(x),3))
X[:,0] = 1.0
X[:,1] = x
X[:,2] = x**2
# 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)
# matrix inversion to find beta
beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train
print(beta)
# and then make the prediction
ytilde = X_train @ beta
print("Training R2")
print(R2(y_train,ytilde))
print("Training MSE")
print(MSE(y_train,ytilde))
ypredict = X_test @ beta
print("Test R2")
print(R2(y_test,ypredict))
print("Test MSE")
print(MSE(y_test,ypredict))
!ec
Alternatively, you could write your own test-train splitting function as shown here.
!bc pycod
# equivalently in numpy
def train_test_split_numpy(inputs, labels, train_size, test_size):
n_inputs = len(inputs)
inputs_shuffled = inputs.copy()
labels_shuffled = labels.copy()
np.random.shuffle(inputs_shuffled)
np.random.shuffle(labels_shuffled)
train_end = int(n_inputs*train_size)
X_train, X_test = inputs_shuffled[:train_end], inputs_shuffled[train_end:]
Y_train, Y_test = labels_shuffled[:train_end], labels_shuffled[train_end:]
return X_train, X_test, Y_train, Y_test
!ec
But since _scikit-learn_ has its own function for doing this and since
it interfaces easily with _tensorflow_ and other libraries, we
normally recommend using the latter functionality.
As another example, we apply the training and testing split to
to the above equation of state fitting example
but now splitting the data into a training set and a test set.
!bc pycod
import os
import numpy as np
@@ -1881,91 +1985,6 @@ plt.show()
!ec
===== Splitting our Data in Training and Test data =====
It is normal in essentially all Machine Learning studies to split the
data in a training set and a test set (sometimes also an additional
validation set). _Scikit-Learn_ has an own function for this. There
is no explicit recipe for how much data should be included as training
data and say test data. An accepted rule of thumb is to use
approximately $2/3$ to $4/5$ of the data as training data. We will
postpone a discussion of this splitting to the end of these notes and
our discussion of the so-called _bias-variance_ tradeoff. Here we
limit ourselves to repeat the above equation of state fitting example
but now splitting the data into a training set and a test set.
Let us study some examples. The first code here takes a simple
one-dimensional second-order polynomial and we fit it to a
second-order polynomial. Depending on the strength of the added noise,
the various measures like the $R2$ score or the mean-squared error,
the fit becomes better or worse.
!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
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
x = np.random.rand(100)
y = 2.0+5*x*x+0.1*np.random.randn(100)
# The design matrix now as function of a given polynomial
X = np.zeros((len(x),3))
X[:,0] = 1.0
X[:,1] = x
X[:,2] = x**2
# 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)
# matrix inversion to find beta
beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train
print(beta)
# and then make the prediction
ytilde = X_train @ beta
print("Training R2")
print(R2(y_train,ytilde))
print("Training MSE")
print(MSE(y_train,ytilde))
ypredict = X_test @ beta
print("Test R2")
print(R2(y_test,ypredict))
print("Test MSE")
print(MSE(y_test,ypredict))
!ec
Alternatively, you could write your own test-train splitting function as shown here.
!bc pycod
# equivalently in numpy
def train_test_split_numpy(inputs, labels, train_size, test_size):
n_inputs = len(inputs)
inputs_shuffled = inputs.copy()
labels_shuffled = labels.copy()
np.random.shuffle(inputs_shuffled)
np.random.shuffle(labels_shuffled)
train_end = int(n_inputs*train_size)
X_train, X_test = inputs_shuffled[:train_end], inputs_shuffled[train_end:]
Y_train, Y_test = labels_shuffled[:train_end], labels_shuffled[train_end:]
return X_train, X_test, Y_train, Y_test
!ec
But since _scikit-learn_ has its own function for doing this and since
it interfaces easily with _tensorflow_ and other libraries, we
normally recommend using the latter functionality.
===== Reducing the number of degrees of freedom, overarching view =====
doc/BookChapters/chapter2.do.txt
+56 -34
View File
@@ -7,11 +7,11 @@
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).
We have shown that in ordinary least squares the optimal parameters $\beta$ are given by
We have shown that in ordinary least squares (OLS) the optimal parameters $\beta$ are given by
!bt
\[
\hat{\bm{\beta}} = \left(\bm{X}^T\bm{X}\right)^{-1}\bm{X}^T\bm{y}.
\hat{\bm{\beta}}_{\mathrm{OLS}} = \left(\bm{X}^T\bm{X}\right)^{-1}\bm{X}^T\bm{y}.
\]
!et
@@ -39,7 +39,7 @@ We can rewrite
\]
!et
The matrix $\bm{A}$ has the important property that $\bm{A}^2=\bm{A}$. This is the definition of a projection matrix.
The matrix $\bm{A}$ has the important property that $\bm{A}^2=\bm{A}$. This is the definition of a "projection matrix":"https://en.wikipedia.org/wiki/Projection_matrix".
We can then interpret our optimal model $\tilde{\bm{y}}$ as being represented by an orthogonal projection of $\bm{y}$ onto a space defined by the column vectors of $\bm{X}$. In our case here the matrix $\bm{A}$ is a square matrix. If it is a general rectangular matrix we have an oblique projection matrix.
@@ -352,7 +352,7 @@ example
===== Code for SVD and Inversion of Matrices =====
How do we use the SVD to invert a matrix $\bm{X}^\bm{X}$ which is singular or near singular?
How do we use the SVD to invert a matrix $\bm{X}^T\bm{X}$ which is singular or near singular?
The simple answer is to use the linear algebra function for the pseudoinverse, that is
!bc pycod
#Ainv = np.linlag.pinv(A)
@@ -2057,7 +2057,7 @@ y = matrix( [4, 2, 3])
x = l1regls(X,y)
!ec
_More text will be added to this example._
===== Linking the regression analysis with a statistical interpretation =====
@@ -2279,7 +2279,9 @@ p(\bm{D}\vert\bm{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\lef
\]
!et
It is a conditional probability (see below) and reads as the likelihood of a domain of events $\bm{D}$ given a set of parameters $\bm{\beta}$.
It is a conditional probability (see below) and reads as the
likelihood of a domain of events $\bm{D}$ given a set of parameters
$\bm{\beta}$.
In statistics, maximum likelihood estimation (MLE) is a method of
@@ -2348,47 +2350,40 @@ Assume we have two domains of events $X=[x_0,x_1,\dots,x_{n-1}]$ and $Y=[y_0,y_1
We define also the likelihood for $X$ and $Y$ as $p(X)$ and $p(Y)$ respectively.
The likelihood of a specific event $x_i$ (or $y_i$) is then written as $p(X=x_i)$ or just $p(x_i)=p_i$.
!bblock Union of events is given by
The union of events is given by
!bt
\[
p(X \cup Y)= p(X)+p(Y)-p(X \cap Y).
\]
!et
!eblock
!bblock The product rule (aka joint probability) is given by
The product rule (aka joint probability) is given by
!bt
\[
p(X \cup Y)= p(X,Y)= p(X\vert Y)p(Y)=p(Y\vert X)p(X),
\]
!et
where we read $p(X\vert Y)$ as the likelihood of obtaining $X$ given $Y$.
!eblock
If we have independent events then $p(X,Y)=p(X)p(Y)$.
The marginal probability is defined in terms of only one of the set of variables $X,Y$. For a discrete probability we have
!bblock
!bt
\[
p(X)=\sum_{i=0}^{n-1}p(X,Y=y_i)=\sum_{i=0}^{n-1}p(X\vert Y=y_i)p(Y=y_i)=\sum_{i=0}^{n-1}p(X\vert y_i)p(y_i).
\]
!et
!eblock
The conditional probability, if $p(Y) > 0$, is
!bblock
!bt
\[
p(X\vert Y)= \frac{p(X,Y)}{p(Y)}=\frac{p(X,Y)}{\sum_{i=0}^{n-1}p(Y\vert X=x_i)p(x_i)}.
\]
!et
!eblock
@@ -2413,21 +2408,27 @@ evaluated for the observed data $Y$ and can be viewed as a function of
the parameter space represented by $X$. This function is not
necesseraly normalized and is normally called the likelihood function.
The function $p(X)$ on the right hand side is called the prior while the function on the left hand side is the called the posterior probability. The denominator on the right hand side serves as a normalization factor for the posterior distribution.
The function $p(X)$ on the right hand side is called the prior while
the function on the left hand side is the called the posterior
probability. The denominator on the right hand side serves as a
normalization factor for the posterior distribution.
Let us try to illustrate Bayes' theorem through an example.
Let us suppose that you are undergoing a series of mammography scans in
order to rule out possible breast cancer cases. We define the
Let us suppose that you are undergoing a series of mammography scans
in order to rule out possible breast cancer cases. We define the
sensitivity for a positive event by the variable $X$. It takes binary
values with $X=1$ representing a positive event and $X=0$ being a
negative event. We reserve $Y$ as a classification parameter for
either a negative or a positive breast cancer confirmation. (Short note on wordings: positive here means having breast cancer, although none of us would consider this being a positive thing).
either a negative or a positive breast cancer confirmation. (Short
note on wordings: positive here means having breast cancer, although
none of us would consider this being a positive thing).
We let $Y=1$ represent the the case of having breast cancer and $Y=0$ as not.
Let us assume that if you have breast cancer, the test will be positive with a probability of $0.8$, that is we have
Let us assume that if you have breast cancer, the test will be positive with a probability of $0.8$ (the numbers here are all made up),
that is we have
!bt
\[
@@ -2435,12 +2436,22 @@ p(X=1\vert Y=1) =0.8.
\]
!et
This obviously sounds scary since many would conclude that if the test is positive, there is a likelihood of $80\%$ for having cancer.
It is however not correct, as the following Bayesian analysis shows.
This obviously sounds scary since many would conclude that if the test
is positive, there is a likelihood of $80\%$ for having cancer. It is
however not correct, as the following Bayesian analysis shows. The correct question to pose is *what is the probability of having breast cancer in case of a positive test?*
We are thus interested in
!bt
\[
p(Y=1\vert X=1),
\]
!et
instead of $p(X=1\vert Y=1)$.
If we look at various national surveys on breast cancer, the general likelihood of developing breast cancer is a very small number.
Let us assume that the prior probability in the population as a whole is
If we look at various national surveys on breast cancer, the general
likelihood of developing breast cancer is a very small number. Let us
assume that the prior probability in the population as a whole is
!bt
\[
@@ -2448,19 +2459,24 @@ p(Y=1) =0.004.
\]
!et
We need also to account for the fact that the test may produce a false positive result (false alarm). Let us here assume that we have
We need also to account for the fact that the test may produce a false
positive result (false alarm). Let us here assume that we have
!bt
\[
p(X=1\vert Y=0) =0.1.
\]
!et
Using Bayes' theorem we can then find the posterior probability that the person has breast cancer in case of a positive test, that is we can compute
Using Bayes' theorem we can then find the posterior probability that
the person has breast cancer in case of a positive test, that is we
can compute
!bt
\[
p(Y=1\vert X=1)=\frac{p(X=1\vert Y=1)p(Y=1)}{p(X=1\vert Y=1)p(Y=1)+p(X=1\vert Y=0)p(Y=0)}=\frac{0.8\times 0.004}{0.8\times 0.004+0.1\times 0.996}=0.031.
\]
\begin{align}
p(Y=1\vert X=1)=\frac{p(X=1\vert Y=1)p(Y=1)}{p(X=1\vert Y=1)p(Y=1)+p(X=1\vert Y=0)p(Y=0)}&= \\
& \frac{0.8\times 0.004}{0.8\times 0.004+0.1\times 0.996}=0.031.
\end{align}
!et
That is, in case of a positive test, there is only a $3\%$ chance of having breast cancer!
@@ -2480,9 +2496,12 @@ We will play around with a study of the values for the optimal
parameters $\bm{\beta}$ using OLS, Ridge and Lasso regression. For
OLS, you will notice as function of the noise and polynomial degree,
that the parameters $\beta$ will fluctuate from order to order in the
polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS.
polynomial fit and that for larger and larger polynomial degrees of
freedom, the parameters will tend to increase in value for OLS.
For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
For Ridge and Lasso regression, the higher order parameters will
typically be reduced, providing thereby less fluctuations from one
order to another one.
!bc pycod
import numpy as np
@@ -2538,8 +2557,8 @@ for i in range(nlambdas):
# Compute the MSE and print it
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
MSELassoPredict[i] = MSE(y_test,ypredictLasso)
print(lmb,RegRidge.coef_)
print(lmb,RegLasso.coef_)
# print(lmb,RegRidge.coef_)
# print(lmb,RegLasso.coef_)
# Now plot the results
plt.figure()
plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
@@ -2552,10 +2571,13 @@ plt.show()
!ec
How can we understand this?
How can we understand this? _More text to be added_.
===== Linking Bayes' Theorem with Ridge and Lasso Regression =====
Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression.
For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\bm{D}$ (one-dimensional case)
doc/BookChapters/chapter3.do.txt
+157 -154
View File
@@ -427,160 +427,6 @@ plt.show()
!ec
===== Various steps in cross-validation =====
When the repetitive splitting of the data set is done randomly,
samples may accidently end up in a fast majority of the splits in
either training or test set. Such samples may have an unbalanced
influence on either model building or prediction evaluation. To avoid
this $k$-fold cross-validation structures the data splitting. The
samples are divided into $k$ more or less equally sized exhaustive and
mutually exclusive subsets. In turn (at each split) one of these
subsets plays the role of the test set while the union of the
remaining subsets constitutes the training set. Such a splitting
warrants a balanced representation of each sample in both training and
test set over the splits. Still the division into the $k$ subsets
involves a degree of randomness. This may be fully excluded when
choosing $k=n$. This particular case is referred to as leave-one-out
cross-validation (LOOCV).
* Define a range of interest for the penalty parameter.
* Divide the data set into training and test set comprising samples $\{1, \ldots, n\} \setminus i$ and $\{ i \}$, respectively.
* Fit the linear regression model by means of ridge estimation for each $\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\bm{\sigma}_{-i}^2(\lambda)$, as
!bt
\begin{align*}
\bm{\beta}_{-i}(\lambda) & = ( \bm{X}_{-i, \ast}^{T}
\bm{X}_{-i, \ast} + \lambda \bm{I}_{pp})^{-1}
\bm{X}_{-i, \ast}^{T} \bm{y}_{-i}
\end{align*}
!et
* Evaluate the prediction performance of these models on the test set by $\log\{L[y_i, \bm{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}$. Or, by the prediction error $|y_i - \bm{X}_{i, \ast} \bm{\beta}_{-i}(\lambda)|$, the relative error, the error squared or the R2 score function.
* Repeat the first three steps such that each sample plays the role of the test set once.
* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as
!bt
\begin{align*}
\frac{1}{n} \sum_{i = 1}^n \log\{L[y_i, \mathbf{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}.
\end{align*}
!et
For the various values of $k$
o shuffle the dataset randomly.
o Split the dataset into $k$ groups.
o For each unique group:
o Decide which group to use as set for test data
o Take the remaining groups as a training data set
o Fit a model on the training set and evaluate it on the test set
o Retain the evaluation score and discard the model
o Summarize the model using the sample of model evaluation scores
The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial.
!bc pycod
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import KFold
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import PolynomialFeatures
# A seed just to ensure that the random numbers are the same for every run.
# Useful for eventual debugging.
np.random.seed(3155)
# Generate the data.
nsamples = 100
x = np.random.randn(nsamples)
y = 3*x**2 + np.random.randn(nsamples)
## Cross-validation on Ridge regression using KFold only
# Decide degree on polynomial to fit
poly = PolynomialFeatures(degree = 6)
# Decide which values of lambda to use
nlambdas = 500
lambdas = np.logspace(-3, 5, nlambdas)
# Initialize a KFold instance
k = 5
kfold = KFold(n_splits = k)
# Perform the cross-validation to estimate MSE
scores_KFold = np.zeros((nlambdas, k))
i = 0
for lmb in lambdas:
ridge = Ridge(alpha = lmb)
j = 0
for train_inds, test_inds in kfold.split(x):
xtrain = x[train_inds]
ytrain = y[train_inds]
xtest = x[test_inds]
ytest = y[test_inds]
Xtrain = poly.fit_transform(xtrain[:, np.newaxis])
ridge.fit(Xtrain, ytrain[:, np.newaxis])
Xtest = poly.fit_transform(xtest[:, np.newaxis])
ypred = ridge.predict(Xtest)
scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred)
j += 1
i += 1
estimated_mse_KFold = np.mean(scores_KFold, axis = 1)
## Cross-validation using cross_val_score from sklearn along with KFold
# kfold is an instance initialized above as:
# kfold = KFold(n_splits = k)
estimated_mse_sklearn = np.zeros(nlambdas)
i = 0
for lmb in lambdas:
ridge = Ridge(alpha = lmb)
X = poly.fit_transform(x[:, np.newaxis])
estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold)
# cross_val_score return an array containing the estimated negative mse for every fold.
# we have to the the mean of every array in order to get an estimate of the mse of the model
estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
i += 1
## Plot and compare the slightly different ways to perform cross-validation
plt.figure()
plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold')
plt.xlabel('log10(lambda)')
plt.ylabel('mse')
plt.legend()
plt.show()
!ec
===== The bias-variance tradeoff =====
@@ -957,6 +803,163 @@ plt.show()
===== Cross-validation =====
When the repetitive splitting of the data set is done randomly,
samples may accidently end up in a fast majority of the splits in
either training or test set. Such samples may have an unbalanced
influence on either model building or prediction evaluation. To avoid
this $k$-fold cross-validation structures the data splitting. The
samples are divided into $k$ more or less equally sized exhaustive and
mutually exclusive subsets. In turn (at each split) one of these
subsets plays the role of the test set while the union of the
remaining subsets constitutes the training set. Such a splitting
warrants a balanced representation of each sample in both training and
test set over the splits. Still the division into the $k$ subsets
involves a degree of randomness. This may be fully excluded when
choosing $k=n$. This particular case is referred to as leave-one-out
cross-validation (LOOCV).
* Define a range of interest for the penalty parameter.
* Divide the data set into training and test set comprising samples $\{1, \ldots, n\} \setminus i$ and $\{ i \}$, respectively.
* Fit the linear regression model by means of ridge estimation for each $\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\bm{\sigma}_{-i}^2(\lambda)$, as
!bt
\begin{align*}
\bm{\beta}_{-i}(\lambda) & = ( \bm{X}_{-i, \ast}^{T}
\bm{X}_{-i, \ast} + \lambda \bm{I}_{pp})^{-1}
\bm{X}_{-i, \ast}^{T} \bm{y}_{-i}
\end{align*}
!et
* Evaluate the prediction performance of these models on the test set by $\log\{L[y_i, \bm{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}$. Or, by the prediction error $|y_i - \bm{X}_{i, \ast} \bm{\beta}_{-i}(\lambda)|$, the relative error, the error squared or the R2 score function.
* Repeat the first three steps such that each sample plays the role of the test set once.
* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as
!bt
\begin{align*}
\frac{1}{n} \sum_{i = 1}^n \log\{L[y_i, \mathbf{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}.
\end{align*}
!et
For the various values of $k$
o shuffle the dataset randomly.
o Split the dataset into $k$ groups.
o For each unique group:
o Decide which group to use as set for test data
o Take the remaining groups as a training data set
o Fit a model on the training set and evaluate it on the test set
o Retain the evaluation score and discard the model
o Summarize the model using the sample of model evaluation scores
The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial.
!bc pycod
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import KFold
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import PolynomialFeatures
# A seed just to ensure that the random numbers are the same for every run.
# Useful for eventual debugging.
np.random.seed(3155)
# Generate the data.
nsamples = 100
x = np.random.randn(nsamples)
y = 3*x**2 + np.random.randn(nsamples)
## Cross-validation on Ridge regression using KFold only
# Decide degree on polynomial to fit
poly = PolynomialFeatures(degree = 6)
# Decide which values of lambda to use
nlambdas = 500
lambdas = np.logspace(-3, 5, nlambdas)
# Initialize a KFold instance
k = 5
kfold = KFold(n_splits = k)
# Perform the cross-validation to estimate MSE
scores_KFold = np.zeros((nlambdas, k))
i = 0
for lmb in lambdas:
ridge = Ridge(alpha = lmb)
j = 0
for train_inds, test_inds in kfold.split(x):
xtrain = x[train_inds]
ytrain = y[train_inds]
xtest = x[test_inds]
ytest = y[test_inds]
Xtrain = poly.fit_transform(xtrain[:, np.newaxis])
ridge.fit(Xtrain, ytrain[:, np.newaxis])
Xtest = poly.fit_transform(xtest[:, np.newaxis])
ypred = ridge.predict(Xtest)
scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred)
j += 1
i += 1
estimated_mse_KFold = np.mean(scores_KFold, axis = 1)
## Cross-validation using cross_val_score from sklearn along with KFold
# kfold is an instance initialized above as:
# kfold = KFold(n_splits = k)
estimated_mse_sklearn = np.zeros(nlambdas)
i = 0
for lmb in lambdas:
ridge = Ridge(alpha = lmb)
X = poly.fit_transform(x[:, np.newaxis])
estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold)
# cross_val_score return an array containing the estimated negative mse for every fold.
# we have to the the mean of every array in order to get an estimate of the mse of the model
estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
i += 1
## Plot and compare the slightly different ways to perform cross-validation
plt.figure()
plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold')
plt.xlabel('log10(lambda)')
plt.ylabel('mse')
plt.legend()
plt.show()
!ec
More examples of the application of cross-validation follow here.
!bc pycod
# Common imports
import os
doc/BookChapters/chapter4.do.txt
+586 -126
View File
File diff suppressed because it is too large Load Diff
doc/LectureNotes/chapter1.ipynb
+200 -166
View File
@@ -13,38 +13,37 @@
"\n",
"\n",
"\n",
"Our emphasis throughout this series of lectures \n",
"is on understanding the mathematical aspects of\n",
"different algorithms used in the fields of data analysis and machine learning. \n",
"Our emphasis throughout this series of lectures is on understanding\n",
"the mathematical aspects of different algorithms used in the fields of\n",
"data analysis and machine learning.\n",
"\n",
"However, where possible we will emphasize the\n",
"importance of using available software. We start thus with a hands-on\n",
"and top-down approach to machine learning. The aim is thus to start with\n",
"relevant data or data we have produced \n",
"and use these to introduce statistical data analysis\n",
"concepts and machine learning algorithms before we delve into the\n",
"algorithms themselves. The examples we will use in the beginning, start with simple\n",
"polynomials with random noise added. We will use the Python\n",
"software package [Scikit-Learn](http://scikit-learn.org/stable/) and\n",
"introduce various machine learning algorithms to make fits of\n",
"the data and predictions. We move thereafter to more interesting\n",
"cases such as data from say experiments (below we will look at experimental nuclear binding energies as an example).\n",
"These are examples where we can easily set up the data and\n",
"then use machine learning algorithms included in for example\n",
"**Scikit-Learn**. \n",
"However, where possible we will emphasize the importance of using\n",
"available software. We start thus with a hands-on and top-down\n",
"approach to machine learning. The aim is thus to start with relevant\n",
"data or data we have produced and use these to introduce statistical\n",
"data analysis concepts and machine learning algorithms before we delve\n",
"into the algorithms themselves. The examples we will use in the\n",
"beginning, start with simple polynomials with random noise added. We\n",
"will use the Python software package\n",
"[Scikit-Learn](http://scikit-learn.org/stable/) and introduce various\n",
"machine learning algorithms to make fits of the data and\n",
"predictions. We move thereafter to more interesting cases such as data\n",
"from say experiments (below we will look at experimental nuclear\n",
"binding energies as an example). These are examples where we can\n",
"easily set up the data and then use machine learning algorithms\n",
"included in for example **Scikit-Learn**.\n",
"\n",
"These examples will serve us the purpose of getting\n",
"started. Furthermore, they allow us to catch more than two birds with\n",
"a stone. They will allow us to bring in some programming specific\n",
"topics and tools as well as showing the power of various Python \n",
"libraries for machine learning and statistical data analysis. \n",
"topics and tools as well as showing the power of various Python\n",
"libraries for machine learning and statistical data analysis.\n",
"\n",
"Here, we will mainly focus on two\n",
"specific Python packages for Machine Learning, Scikit-Learn and\n",
"Tensorflow (see below for links etc). Moreover, the examples we\n",
"introduce will serve as inputs to many of our discussions later, as\n",
"well as allowing you to set up models and produce your own data and\n",
"get started with programming.\n",
"Here, we will mainly focus on two specific Python packages for Machine\n",
"Learning, Scikit-Learn and Tensorflow (see below for links etc).\n",
"Moreover, the examples we introduce will serve as inputs to many of\n",
"our discussions later, as well as allowing you to set up models and\n",
"produce your own data and get started with programming.\n",
"\n",
"\n",
"\n",
@@ -286,7 +285,17 @@
"[gallery](https://matplotlib.org/gallery/index.html) of examples. In\n",
"this example we plot our original values of $x$ and $y$ as well as the\n",
"prediction **ypredict** ($\\tilde{y}$), which attempts at fitting our\n",
"data with a straight line.\n",
"data with a straight line. Note also that **Scikit-Learn** requires a\n",
"matrix as input for the input values $x$ and $y$. In the above code we\n",
"have solved this by declaring $x$ and $y$ as arrays of dimension\n",
"$n\\times 1$.\n",
"\n",
"In the code here we have also made a new array for $x\\in [0,1]$. Our\n",
"prediction is computed for these values, meaning that they were not\n",
"included in the data set used to *train* (or fit) the model.\n",
"This is a recurrring theme in machine learning and data analysis. We would like to train a model on a specific given data set.\n",
"Thereafter we wish to apply it to data which were not included in the training. Below we will encounter this again in the so-called *train-validate-test* spliting. We will typically split our data into different sets, oen for training, one for validation and finally, our data from the untouched test vault!\n",
"\n",
"\n",
"The Python code follows here."
]
@@ -311,6 +320,7 @@
"y = 2*x+np.random.randn(100,1)\n",
"linreg = LinearRegression()\n",
"linreg.fit(x,y)\n",
"# This is our new x-array to which we test our model\n",
"xnew = np.array([[0],[1]])\n",
"ypredict = linreg.predict(xnew)\n",
"\n",
@@ -460,8 +470,7 @@
"Depending on the parameter in front of the normal distribution, we may\n",
"have a small or larger relative error. Try to play around with\n",
"different training data sets and study (graphically) the value of the\n",
"relative error. Note also that **Scikit-Learn** requires a matrix as input for the input values $x$ and $y$. In the above code we have\n",
"solved this by declaring $x$ and $y$ as arrays of dimension $n\\times 1$.\n",
"relative error.\n",
"\n",
"As mentioned above, **Scikit-Learn** has an impressive functionality.\n",
"We can for example extract the values of $\\alpha$ and $\\beta$ and\n",
@@ -1039,7 +1048,7 @@
"metadata": {},
"source": [
"The next step, and we will define this mathematically later, is to set up the so-called **design matrix**. We will throughout call this matrix $\\boldsymbol{X}$.\n",
"It has dimensionality $p\\times n$, where $n$ is the number of data points and $p$ are the so-called predictors. In our case here they are given by the number of polynomials in $A$ we wish to include in the fit."
"It has dimensionality $n\\times p$, where $n$ is the number of data points and $p$ are the so-called predictors. In our case here they are given by the number of polynomials in $A$ we wish to include in the fit."
]
},
{
@@ -1064,9 +1073,12 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Note well that we have made life simple here. We perform a fit in terms of the number of nucleons only. A more sophisticated fit can be done by including an explicit dependence on the number of protons and neutrons in the asymmetry and Coulomb terms.\n",
"Note well that we have made life simple here. We perform a fit in\n",
"terms of the number of nucleons only. A more sophisticated fit can be\n",
"done by including an explicit dependence on the number of protons and\n",
"neutrons in the asymmetry and Coulomb terms. We leave this as an exercise to you the reader.\n",
"\n",
"With **scikitlearn** we are now ready to use linear regression and fit our data."
"With **Scikit-Learn** we are now ready to use linear regression and fit our data."
]
},
{
@@ -1124,7 +1136,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"As a teaser, let us now see how we can do this with decision trees using **scikit-learn**. Later we will switch to so-called **random forests**!"
"As a teaser, let us now see how we can do this with decision trees using **Scikit-Learn**. Later we will switch to so-called **random forests**!"
]
},
{
@@ -1172,6 +1184,17 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"With a deeper and deeper tree level, we can almost reproduce every\n",
"single data point by increasing the max depth of the tree.\n",
"We can actually decide to make a decision tree which fits every single point.\n",
"As we will\n",
"see later, this has the benefit that we can really train a model which\n",
"traverses every single data point. However, the price we pay is that\n",
"we will easily overfit. That is, if we apply our model to unseen data,\n",
"we will most likely fail miserably in our attempt at making\n",
"predictions. As an exercise, try to make the tree level larger by adjusting the maximum depth variable. When printing out the predicition, you will note that the binding energy of every nucleus is accurately reproduced.\n",
"\n",
"\n",
"The **seaborn** package allows us to visualize data in an efficient way. Note that we use **scikit-learn**'s multi-layer perceptron (or feed forward neural network) \n",
"functionality."
]
@@ -1977,6 +2000,47 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"We can then compute the second derivative of the cost function, which in our case is the second derivative\n",
"of the means squared error. This leads to"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial^2 C(\\boldsymbol{\\beta})}{\\partial \\boldsymbol{\\beta}^T\\partial \\boldsymbol{\\beta}} =\\frac{2}{n}\\boldsymbol{X}^T\\boldsymbol{X}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This quantity defines was what is called the Hessian matrix (the second derivative of a function we want to optimize).\n",
"\n",
"The Hessian matrix plays an important role and is defined for the mean squared error as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\boldsymbol{H}=\\boldsymbol{X}^T\\boldsymbol{X}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Hessian matrix for ordinary least squares is also proportional to\n",
"the covariance matrix. As we will see in the chapter on Ridge and Lasso regression, This means that we can use the Singular Value Decomposition of a matrix to find\n",
"the eigenvalues of the covariance matrix and the Hessian matrix in\n",
"terms of the singular values.\n",
"\n",
"\n",
"The residuals $\\boldsymbol{\\epsilon}$ are in turn given by"
]
},
@@ -2529,9 +2593,7 @@
"**pandas** again, rather extensively in order to organize our data.\n",
"\n",
"The difference now is that we use **Scikit-Learn's** regression tools\n",
"instead of our own matrix inversion implementation. Furthermore, we\n",
"sneak in **Ridge** regression (to be discussed below) which includes a\n",
"hyperparameter $\\lambda$, also to be explained below."
"instead of our own matrix inversion implementation."
]
},
{
@@ -2603,18 +2665,6 @@
"print('Mean absolute error: %.2f' % mean_absolute_error(Energies, ytilde))\n",
"print(clf.coef_, clf.intercept_)\n",
"\n",
"# The Ridge regression with a hyperparameter lambda = 0.1\n",
"_lambda = 0.1\n",
"clf_ridge = skl.Ridge(alpha=_lambda).fit(X, Energies)\n",
"yridge = clf_ridge.predict(X)\n",
"EoS['Eridge'] = yridge\n",
"# The mean squared error \n",
"print(\"Mean squared error: %.2f\" % mean_squared_error(Energies, yridge))\n",
"# Explained variance score: 1 is perfect prediction \n",
"print('Variance score: %.2f' % r2_score(Energies, yridge))\n",
"# Mean absolute error \n",
"print('Mean absolute error: %.2f' % mean_absolute_error(Energies, yridge))\n",
"print(clf_ridge.coef_, clf_ridge.intercept_)\n",
"\n",
"fig, ax = plt.subplots()\n",
"ax.set_xlabel(r'$\\rho[\\mathrm{fm}^{-3}]$')\n",
@@ -2623,8 +2673,6 @@
" label='Theoretical data')\n",
"ax.plot(EoS['Density'], EoS['Eols'], alpha=0.7, lw=2, c='m',\n",
" label='OLS')\n",
"ax.plot(EoS['Density'], EoS['Eridge'], alpha=0.7, lw=2, c='g',\n",
" label='Ridge $\\lambda = 0.1$')\n",
"ax.legend()\n",
"save_fig(\"EoSfitting\")\n",
"plt.show()"
@@ -2637,13 +2685,11 @@
"The above simple polynomial in density $\\rho$ gives an excellent fit\n",
"to the data. \n",
"\n",
"We note also that there is a small deviation between the\n",
"standard OLS and the Ridge regression at higher densities. We discuss this in more detail\n",
"below.\n",
"\n",
"\n",
"## Splitting our Data in Training and Test data\n",
"\n",
"\n",
"It is normal in essentially all Machine Learning studies to split the\n",
"data in a training set and a test set (sometimes also an additional\n",
"validation set). **Scikit-Learn** has an own function for this. There\n",
@@ -2653,6 +2699,107 @@
"postpone a discussion of this splitting to the end of these notes and\n",
"our discussion of the so-called **bias-variance** tradeoff. Here we\n",
"limit ourselves to repeat the above equation of state fitting example\n",
"but now splitting the data into a training set and a test set.\n",
"\n",
"Let us study some examples. The first code here takes a simple\n",
"one-dimensional second-order polynomial and we fit it to a\n",
"second-order polynomial. Depending on the strength of the added noise,\n",
"the various measures like the $R2$ score or the mean-squared error,\n",
"the fit becomes better or worse."
]
},
{
"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",
"\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",
"x = np.random.rand(100)\n",
"y = 2.0+5*x*x+0.1*np.random.randn(100)\n",
"\n",
"\n",
"# The design matrix now as function of a given polynomial\n",
"X = np.zeros((len(x),3))\n",
"X[:,0] = 1.0\n",
"X[:,1] = x\n",
"X[:,2] = x**2\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",
"# matrix inversion to find beta\n",
"beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train\n",
"print(beta)\n",
"# and then make the prediction\n",
"ytilde = X_train @ beta\n",
"print(\"Training R2\")\n",
"print(R2(y_train,ytilde))\n",
"print(\"Training MSE\")\n",
"print(MSE(y_train,ytilde))\n",
"ypredict = X_test @ beta\n",
"print(\"Test R2\")\n",
"print(R2(y_test,ypredict))\n",
"print(\"Test MSE\")\n",
"print(MSE(y_test,ypredict))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Alternatively, you could write your own test-train splitting function as shown here."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"# equivalently in numpy\n",
"def train_test_split_numpy(inputs, labels, train_size, test_size):\n",
" n_inputs = len(inputs)\n",
" inputs_shuffled = inputs.copy()\n",
" labels_shuffled = labels.copy()\n",
"\n",
" np.random.shuffle(inputs_shuffled)\n",
" np.random.shuffle(labels_shuffled)\n",
"\n",
" train_end = int(n_inputs*train_size)\n",
" X_train, X_test = inputs_shuffled[:train_end], inputs_shuffled[train_end:]\n",
" Y_train, Y_test = labels_shuffled[:train_end], labels_shuffled[train_end:]\n",
"\n",
" return X_train, X_test, Y_train, Y_test"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But since **scikit-learn** has its own function for doing this and since\n",
"it interfaces easily with **tensorflow** and other libraries, we\n",
"normally recommend using the latter functionality.\n",
"\n",
"\n",
"As another example, we apply the training and testing split to \n",
"to the above equation of state fitting example\n",
"but now splitting the data into a training set and a test set."
]
},
@@ -3050,119 +3197,6 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Splitting our Data in Training and Test data\n",
"\n",
"\n",
"It is normal in essentially all Machine Learning studies to split the\n",
"data in a training set and a test set (sometimes also an additional\n",
"validation set). **Scikit-Learn** has an own function for this. There\n",
"is no explicit recipe for how much data should be included as training\n",
"data and say test data. An accepted rule of thumb is to use\n",
"approximately $2/3$ to $4/5$ of the data as training data. We will\n",
"postpone a discussion of this splitting to the end of these notes and\n",
"our discussion of the so-called **bias-variance** tradeoff. Here we\n",
"limit ourselves to repeat the above equation of state fitting example\n",
"but now splitting the data into a training set and a test set.\n",
"\n",
"Let us study some examples. The first code here takes a simple\n",
"one-dimensional second-order polynomial and we fit it to a\n",
"second-order polynomial. Depending on the strength of the added noise,\n",
"the various measures like the $R2$ score or the mean-squared error,\n",
"the fit becomes better or worse."
]
},
{
"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",
"\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",
"x = np.random.rand(100)\n",
"y = 2.0+5*x*x+0.1*np.random.randn(100)\n",
"\n",
"\n",
"# The design matrix now as function of a given polynomial\n",
"X = np.zeros((len(x),3))\n",
"X[:,0] = 1.0\n",
"X[:,1] = x\n",
"X[:,2] = x**2\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",
"# matrix inversion to find beta\n",
"beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train\n",
"print(beta)\n",
"# and then make the prediction\n",
"ytilde = X_train @ beta\n",
"print(\"Training R2\")\n",
"print(R2(y_train,ytilde))\n",
"print(\"Training MSE\")\n",
"print(MSE(y_train,ytilde))\n",
"ypredict = X_test @ beta\n",
"print(\"Test R2\")\n",
"print(R2(y_test,ypredict))\n",
"print(\"Test MSE\")\n",
"print(MSE(y_test,ypredict))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Alternatively, you could write your own test-train splitting function as shown here."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"# equivalently in numpy\n",
"def train_test_split_numpy(inputs, labels, train_size, test_size):\n",
" n_inputs = len(inputs)\n",
" inputs_shuffled = inputs.copy()\n",
" labels_shuffled = labels.copy()\n",
"\n",
" np.random.shuffle(inputs_shuffled)\n",
" np.random.shuffle(labels_shuffled)\n",
"\n",
" train_end = int(n_inputs*train_size)\n",
" X_train, X_test = inputs_shuffled[:train_end], inputs_shuffled[train_end:]\n",
" Y_train, Y_test = labels_shuffled[:train_end], labels_shuffled[train_end:]\n",
"\n",
" return X_train, X_test, Y_train, Y_test"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But since **scikit-learn** has its own function for doing this and since\n",
"it interfaces easily with **tensorflow** and other libraries, we\n",
"normally recommend using the latter functionality.\n",
"\n",
"\n",
"\n",
"\n",
"## Reducing the number of degrees of freedom, overarching view\n",
"\n",
"Many Machine Learning problems involve thousands or even millions of\n",
doc/LectureNotes/chapter2.ipynb
+84 -26
View File
@@ -13,7 +13,7 @@
"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",
"We have shown that in ordinary least squares the optimal parameters $\\beta$ are given by"
"We have shown that in ordinary least squares (OLS) the optimal parameters $\\beta$ are given by"
]
},
{
@@ -21,7 +21,7 @@
"metadata": {},
"source": [
"$$\n",
"\\hat{\\boldsymbol{\\beta}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}.\n",
"\\hat{\\boldsymbol{\\beta}}_{\\mathrm{OLS}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}.\n",
"$$"
]
},
@@ -79,7 +79,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"The matrix $\\boldsymbol{A}$ has the important property that $\\boldsymbol{A}^2=\\boldsymbol{A}$. This is the definition of a projection matrix.\n",
"The matrix $\\boldsymbol{A}$ has the important property that $\\boldsymbol{A}^2=\\boldsymbol{A}$. This is the definition of a [projection matrix](https://en.wikipedia.org/wiki/Projection_matrix).\n",
"We can then interpret our optimal model $\\tilde{\\boldsymbol{y}}$ as being represented by an orthogonal projection of $\\boldsymbol{y}$ onto a space defined by the column vectors of $\\boldsymbol{X}$. In our case here the matrix $\\boldsymbol{A}$ is a square matrix. If it is a general rectangular matrix we have an oblique projection matrix.\n",
"\n",
"\n",
@@ -545,7 +545,7 @@
"\n",
"## Code for SVD and Inversion of Matrices\n",
"\n",
"How do we use the SVD to invert a matrix $\\boldsymbol{X}^\\boldsymbol{X}$ which is singular or near singular?\n",
"How do we use the SVD to invert a matrix $\\boldsymbol{X}^T\\boldsymbol{X}$ which is singular or near singular?\n",
"The simple answer is to use the linear algebra function for the pseudoinverse, that is"
]
},
@@ -3258,6 +3258,8 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"**More text will be added to this example.**\n",
"\n",
"## Linking the regression analysis with a statistical interpretation\n",
"\n",
"We will now couple the discussions of ordinary least squares, Ridge\n",
@@ -3624,7 +3626,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"It is a conditional probability (see below) and reads as the likelihood of a domain of events $\\boldsymbol{D}$ given a set of parameters $\\boldsymbol{\\beta}$.\n",
"It is a conditional probability (see below) and reads as the\n",
"likelihood of a domain of events $\\boldsymbol{D}$ given a set of parameters\n",
"$\\boldsymbol{\\beta}$.\n",
"\n",
"\n",
"In statistics, maximum likelihood estimation (MLE) is a method of\n",
@@ -3728,7 +3732,7 @@
"We define also the likelihood for $X$ and $Y$ as $p(X)$ and $p(Y)$ respectively.\n",
"The likelihood of a specific event $x_i$ (or $y_i$) is then written as $p(X=x_i)$ or just $p(x_i)=p_i$. \n",
"\n",
"**Union of events is given by.**"
"The union of events is given by"
]
},
{
@@ -3744,7 +3748,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"**The product rule (aka joint probability) is given by.**"
"The product rule (aka joint probability) is given by"
]
},
{
@@ -3762,8 +3766,6 @@
"source": [
"where we read $p(X\\vert Y)$ as the likelihood of obtaining $X$ given $Y$.\n",
"\n",
"\n",
"\n",
"If we have independent events then $p(X,Y)=p(X)p(Y)$.\n",
"\n",
"\n",
@@ -3840,21 +3842,27 @@
"the parameter space represented by $X$. This function is not\n",
"necesseraly normalized and is normally called the likelihood function.\n",
"\n",
"The function $p(X)$ on the right hand side is called the prior while the function on the left hand side is the called the posterior probability. The denominator on the right hand side serves as a normalization factor for the posterior distribution.\n",
"The function $p(X)$ on the right hand side is called the prior while\n",
"the function on the left hand side is the called the posterior\n",
"probability. The denominator on the right hand side serves as a\n",
"normalization factor for the posterior distribution.\n",
"\n",
"Let us try to illustrate Bayes' theorem through an example.\n",
"\n",
"\n",
"Let us suppose that you are undergoing a series of mammography scans in\n",
"order to rule out possible breast cancer cases. We define the\n",
"Let us suppose that you are undergoing a series of mammography scans\n",
"in order to rule out possible breast cancer cases. We define the\n",
"sensitivity for a positive event by the variable $X$. It takes binary\n",
"values with $X=1$ representing a positive event and $X=0$ being a\n",
"negative event. We reserve $Y$ as a classification parameter for\n",
"either a negative or a positive breast cancer confirmation. (Short note on wordings: positive here means having breast cancer, although none of us would consider this being a positive thing).\n",
"either a negative or a positive breast cancer confirmation. (Short\n",
"note on wordings: positive here means having breast cancer, although\n",
"none of us would consider this being a positive thing).\n",
"\n",
"We let $Y=1$ represent the the case of having breast cancer and $Y=0$ as not.\n",
"\n",
"Let us assume that if you have breast cancer, the test will be positive with a probability of $0.8$, that is we have"
"Let us assume that if you have breast cancer, the test will be positive with a probability of $0.8$ (the numbers here are all made up),\n",
"that is we have"
]
},
{
@@ -3870,12 +3878,32 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"This obviously sounds scary since many would conclude that if the test is positive, there is a likelihood of $80\\%$ for having cancer.\n",
"It is however not correct, as the following Bayesian analysis shows.\n",
"This obviously sounds scary since many would conclude that if the test\n",
"is positive, there is a likelihood of $80\\%$ for having cancer. It is\n",
"however not correct, as the following Bayesian analysis shows. The correct question to pose is *what is the probability of having breast cancer in case of a positive test?*\n",
"We are thus interested in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"p(Y=1\\vert X=1),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"instead of $p(X=1\\vert Y=1)$.\n",
"\n",
"\n",
"If we look at various national surveys on breast cancer, the general likelihood of developing breast cancer is a very small number.\n",
"Let us assume that the prior probability in the population as a whole is"
"\n",
"If we look at various national surveys on breast cancer, the general\n",
"likelihood of developing breast cancer is a very small number. Let us\n",
"assume that the prior probability in the population as a whole is"
]
},
{
@@ -3891,7 +3919,8 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"We need also to account for the fact that the test may produce a false positive result (false alarm). Let us here assume that we have"
"We need also to account for the fact that the test may produce a false\n",
"positive result (false alarm). Let us here assume that we have"
]
},
{
@@ -3907,15 +3936,38 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Using Bayes' theorem we can then find the posterior probability that the person has breast cancer in case of a positive test, that is we can compute"
"Using Bayes' theorem we can then find the posterior probability that\n",
"the person has breast cancer in case of a positive test, that is we\n",
"can compute"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto2\"></div>\n",
"\n",
"$$\n",
"p(Y=1\\vert X=1)=\\frac{p(X=1\\vert Y=1)p(Y=1)}{p(X=1\\vert Y=1)p(Y=1)+p(X=1\\vert Y=0)p(Y=0)}=\\frac{0.8\\times 0.004}{0.8\\times 0.004+0.1\\times 0.996}=0.031.\n",
"\\begin{equation}\n",
"p(Y=1\\vert X=1)=\\frac{p(X=1\\vert Y=1)p(Y=1)}{p(X=1\\vert Y=1)p(Y=1)+p(X=1\\vert Y=0)p(Y=0)}= \n",
"\\label{_auto2} \\tag{2}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto3\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
" \\frac{0.8\\times 0.004}{0.8\\times 0.004+0.1\\times 0.996}=0.031.\n",
"\\label{_auto3} \\tag{3}\n",
"\\end{equation}\n",
"$$"
]
},
@@ -3941,9 +3993,12 @@
"parameters $\\boldsymbol{\\beta}$ using OLS, Ridge and Lasso regression. For\n",
"OLS, you will notice as function of the noise and polynomial degree,\n",
"that the parameters $\\beta$ will fluctuate from order to order in the\n",
"polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS.\n",
"polynomial fit and that for larger and larger polynomial degrees of\n",
"freedom, the parameters will tend to increase in value for OLS.\n",
"\n",
"For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one."
"For Ridge and Lasso regression, the higher order parameters will\n",
"typically be reduced, providing thereby less fluctuations from one\n",
"order to another one."
]
},
{
@@ -4008,8 +4063,8 @@
" # Compute the MSE and print it\n",
" MSERidgePredict[i] = MSE(y_test,ypredictRidge)\n",
" MSELassoPredict[i] = MSE(y_test,ypredictLasso)\n",
" print(lmb,RegRidge.coef_)\n",
" print(lmb,RegLasso.coef_)\n",
"# print(lmb,RegRidge.coef_)\n",
"# print(lmb,RegLasso.coef_)\n",
"# Now plot the results\n",
"plt.figure()\n",
"plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')\n",
@@ -4024,10 +4079,13 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"How can we understand this? \n",
"How can we understand this? **More text to be added**.\n",
"\n",
"\n",
"\n",
"## Linking Bayes' Theorem with Ridge and Lasso Regression\n",
"\n",
"\n",
"Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression. \n",
"\n",
"For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\\boldsymbol{D}$ (one-dimensional case)"
doc/LectureNotes/chapter3.ipynb
+197 -190
View File
@@ -575,196 +575,6 @@
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Various steps in cross-validation\n",
"\n",
"When the repetitive splitting of the data set is done randomly,\n",
"samples may accidently end up in a fast majority of the splits in\n",
"either training or test set. Such samples may have an unbalanced\n",
"influence on either model building or prediction evaluation. To avoid\n",
"this $k$-fold cross-validation structures the data splitting. The\n",
"samples are divided into $k$ more or less equally sized exhaustive and\n",
"mutually exclusive subsets. In turn (at each split) one of these\n",
"subsets plays the role of the test set while the union of the\n",
"remaining subsets constitutes the training set. Such a splitting\n",
"warrants a balanced representation of each sample in both training and\n",
"test set over the splits. Still the division into the $k$ subsets\n",
"involves a degree of randomness. This may be fully excluded when\n",
"choosing $k=n$. This particular case is referred to as leave-one-out\n",
"cross-validation (LOOCV). \n",
"\n",
"\n",
"* Define a range of interest for the penalty parameter.\n",
"\n",
"* Divide the data set into training and test set comprising samples $\\{1, \\ldots, n\\} \\setminus i$ and $\\{ i \\}$, respectively.\n",
"\n",
"* Fit the linear regression model by means of ridge estimation for each $\\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\\boldsymbol{\\sigma}_{-i}^2(\\lambda)$, as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\boldsymbol{\\beta}_{-i}(\\lambda) & = ( \\boldsymbol{X}_{-i, \\ast}^{T}\n",
"\\boldsymbol{X}_{-i, \\ast} + \\lambda \\boldsymbol{I}_{pp})^{-1}\n",
"\\boldsymbol{X}_{-i, \\ast}^{T} \\boldsymbol{y}_{-i}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Evaluate the prediction performance of these models on the test set by $\\log\\{L[y_i, \\boldsymbol{X}_{i, \\ast}; \\boldsymbol{\\beta}_{-i}(\\lambda), \\boldsymbol{\\sigma}_{-i}^2(\\lambda)]\\}$. Or, by the prediction error $|y_i - \\boldsymbol{X}_{i, \\ast} \\boldsymbol{\\beta}_{-i}(\\lambda)|$, the relative error, the error squared or the R2 score function.\n",
"\n",
"* Repeat the first three steps such that each sample plays the role of the test set once.\n",
"\n",
"* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\frac{1}{n} \\sum_{i = 1}^n \\log\\{L[y_i, \\mathbf{X}_{i, \\ast}; \\boldsymbol{\\beta}_{-i}(\\lambda), \\boldsymbol{\\sigma}_{-i}^2(\\lambda)]\\}.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the various values of $k$\n",
"\n",
"1. shuffle the dataset randomly.\n",
"\n",
"2. Split the dataset into $k$ groups.\n",
"\n",
"3. For each unique group:\n",
"\n",
"a. Decide which group to use as set for test data\n",
"\n",
"b. Take the remaining groups as a training data set\n",
"\n",
"c. Fit a model on the training set and evaluate it on the test set\n",
"\n",
"d. Retain the evaluation score and discard the model\n",
"\n",
"\n",
"5. Summarize the model using the sample of model evaluation scores\n",
"\n",
"The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import KFold\n",
"from sklearn.linear_model import Ridge\n",
"from sklearn.model_selection import cross_val_score\n",
"from sklearn.preprocessing import PolynomialFeatures\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",
"# Generate the data.\n",
"nsamples = 100\n",
"x = np.random.randn(nsamples)\n",
"y = 3*x**2 + np.random.randn(nsamples)\n",
"\n",
"## Cross-validation on Ridge regression using KFold only\n",
"\n",
"# Decide degree on polynomial to fit\n",
"poly = PolynomialFeatures(degree = 6)\n",
"\n",
"# Decide which values of lambda to use\n",
"nlambdas = 500\n",
"lambdas = np.logspace(-3, 5, nlambdas)\n",
"\n",
"# Initialize a KFold instance\n",
"k = 5\n",
"kfold = KFold(n_splits = k)\n",
"\n",
"# Perform the cross-validation to estimate MSE\n",
"scores_KFold = np.zeros((nlambdas, k))\n",
"\n",
"i = 0\n",
"for lmb in lambdas:\n",
" ridge = Ridge(alpha = lmb)\n",
" j = 0\n",
" for train_inds, test_inds in kfold.split(x):\n",
" xtrain = x[train_inds]\n",
" ytrain = y[train_inds]\n",
"\n",
" xtest = x[test_inds]\n",
" ytest = y[test_inds]\n",
"\n",
" Xtrain = poly.fit_transform(xtrain[:, np.newaxis])\n",
" ridge.fit(Xtrain, ytrain[:, np.newaxis])\n",
"\n",
" Xtest = poly.fit_transform(xtest[:, np.newaxis])\n",
" ypred = ridge.predict(Xtest)\n",
"\n",
" scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred)\n",
"\n",
" j += 1\n",
" i += 1\n",
"\n",
"\n",
"estimated_mse_KFold = np.mean(scores_KFold, axis = 1)\n",
"\n",
"## Cross-validation using cross_val_score from sklearn along with KFold\n",
"\n",
"# kfold is an instance initialized above as:\n",
"# kfold = KFold(n_splits = k)\n",
"\n",
"estimated_mse_sklearn = np.zeros(nlambdas)\n",
"i = 0\n",
"for lmb in lambdas:\n",
" ridge = Ridge(alpha = lmb)\n",
"\n",
" X = poly.fit_transform(x[:, np.newaxis])\n",
" estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold)\n",
"\n",
" # cross_val_score return an array containing the estimated negative mse for every fold.\n",
" # we have to the the mean of every array in order to get an estimate of the mse of the model\n",
" estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)\n",
"\n",
" i += 1\n",
"\n",
"## Plot and compare the slightly different ways to perform cross-validation\n",
"\n",
"plt.figure()\n",
"\n",
"plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')\n",
"plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold')\n",
"\n",
"plt.xlabel('log10(lambda)')\n",
"plt.ylabel('mse')\n",
"\n",
"plt.legend()\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
@@ -1221,6 +1031,203 @@
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cross-validation\n",
"\n",
"When the repetitive splitting of the data set is done randomly,\n",
"samples may accidently end up in a fast majority of the splits in\n",
"either training or test set. Such samples may have an unbalanced\n",
"influence on either model building or prediction evaluation. To avoid\n",
"this $k$-fold cross-validation structures the data splitting. The\n",
"samples are divided into $k$ more or less equally sized exhaustive and\n",
"mutually exclusive subsets. In turn (at each split) one of these\n",
"subsets plays the role of the test set while the union of the\n",
"remaining subsets constitutes the training set. Such a splitting\n",
"warrants a balanced representation of each sample in both training and\n",
"test set over the splits. Still the division into the $k$ subsets\n",
"involves a degree of randomness. This may be fully excluded when\n",
"choosing $k=n$. This particular case is referred to as leave-one-out\n",
"cross-validation (LOOCV). \n",
"\n",
"\n",
"* Define a range of interest for the penalty parameter.\n",
"\n",
"* Divide the data set into training and test set comprising samples $\\{1, \\ldots, n\\} \\setminus i$ and $\\{ i \\}$, respectively.\n",
"\n",
"* Fit the linear regression model by means of ridge estimation for each $\\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\\boldsymbol{\\sigma}_{-i}^2(\\lambda)$, as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\boldsymbol{\\beta}_{-i}(\\lambda) & = ( \\boldsymbol{X}_{-i, \\ast}^{T}\n",
"\\boldsymbol{X}_{-i, \\ast} + \\lambda \\boldsymbol{I}_{pp})^{-1}\n",
"\\boldsymbol{X}_{-i, \\ast}^{T} \\boldsymbol{y}_{-i}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Evaluate the prediction performance of these models on the test set by $\\log\\{L[y_i, \\boldsymbol{X}_{i, \\ast}; \\boldsymbol{\\beta}_{-i}(\\lambda), \\boldsymbol{\\sigma}_{-i}^2(\\lambda)]\\}$. Or, by the prediction error $|y_i - \\boldsymbol{X}_{i, \\ast} \\boldsymbol{\\beta}_{-i}(\\lambda)|$, the relative error, the error squared or the R2 score function.\n",
"\n",
"* Repeat the first three steps such that each sample plays the role of the test set once.\n",
"\n",
"* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\frac{1}{n} \\sum_{i = 1}^n \\log\\{L[y_i, \\mathbf{X}_{i, \\ast}; \\boldsymbol{\\beta}_{-i}(\\lambda), \\boldsymbol{\\sigma}_{-i}^2(\\lambda)]\\}.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the various values of $k$\n",
"\n",
"1. shuffle the dataset randomly.\n",
"\n",
"2. Split the dataset into $k$ groups.\n",
"\n",
"3. For each unique group:\n",
"\n",
"a. Decide which group to use as set for test data\n",
"\n",
"b. Take the remaining groups as a training data set\n",
"\n",
"c. Fit a model on the training set and evaluate it on the test set\n",
"\n",
"d. Retain the evaluation score and discard the model\n",
"\n",
"\n",
"5. Summarize the model using the sample of model evaluation scores\n",
"\n",
"The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import KFold\n",
"from sklearn.linear_model import Ridge\n",
"from sklearn.model_selection import cross_val_score\n",
"from sklearn.preprocessing import PolynomialFeatures\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",
"# Generate the data.\n",
"nsamples = 100\n",
"x = np.random.randn(nsamples)\n",
"y = 3*x**2 + np.random.randn(nsamples)\n",
"\n",
"## Cross-validation on Ridge regression using KFold only\n",
"\n",
"# Decide degree on polynomial to fit\n",
"poly = PolynomialFeatures(degree = 6)\n",
"\n",
"# Decide which values of lambda to use\n",
"nlambdas = 500\n",
"lambdas = np.logspace(-3, 5, nlambdas)\n",
"\n",
"# Initialize a KFold instance\n",
"k = 5\n",
"kfold = KFold(n_splits = k)\n",
"\n",
"# Perform the cross-validation to estimate MSE\n",
"scores_KFold = np.zeros((nlambdas, k))\n",
"\n",
"i = 0\n",
"for lmb in lambdas:\n",
" ridge = Ridge(alpha = lmb)\n",
" j = 0\n",
" for train_inds, test_inds in kfold.split(x):\n",
" xtrain = x[train_inds]\n",
" ytrain = y[train_inds]\n",
"\n",
" xtest = x[test_inds]\n",
" ytest = y[test_inds]\n",
"\n",
" Xtrain = poly.fit_transform(xtrain[:, np.newaxis])\n",
" ridge.fit(Xtrain, ytrain[:, np.newaxis])\n",
"\n",
" Xtest = poly.fit_transform(xtest[:, np.newaxis])\n",
" ypred = ridge.predict(Xtest)\n",
"\n",
" scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred)\n",
"\n",
" j += 1\n",
" i += 1\n",
"\n",
"\n",
"estimated_mse_KFold = np.mean(scores_KFold, axis = 1)\n",
"\n",
"## Cross-validation using cross_val_score from sklearn along with KFold\n",
"\n",
"# kfold is an instance initialized above as:\n",
"# kfold = KFold(n_splits = k)\n",
"\n",
"estimated_mse_sklearn = np.zeros(nlambdas)\n",
"i = 0\n",
"for lmb in lambdas:\n",
" ridge = Ridge(alpha = lmb)\n",
"\n",
" X = poly.fit_transform(x[:, np.newaxis])\n",
" estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold)\n",
"\n",
" # cross_val_score return an array containing the estimated negative mse for every fold.\n",
" # we have to the the mean of every array in order to get an estimate of the mse of the model\n",
" estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)\n",
"\n",
" i += 1\n",
"\n",
"## Plot and compare the slightly different ways to perform cross-validation\n",
"\n",
"plt.figure()\n",
"\n",
"plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')\n",
"plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold')\n",
"\n",
"plt.xlabel('log10(lambda)')\n",
"plt.ylabel('mse')\n",
"\n",
"plt.legend()\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"More examples of the application of cross-validation follow here."
]
},
{
"cell_type": "code",
"execution_count": null,
doc/LectureNotes/chapter4.ipynb
+926 -106
View File
File diff suppressed because it is too large Load Diff