1312 lines
43 KiB
Plaintext
1312 lines
43 KiB
Plaintext
TITLE: Week 35: Linear Regression and Review of Statistical Analysis and Probability Theory
|
||
AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
|
||
DATE: today
|
||
|
||
|
||
!split
|
||
===== Plans for week 35, August 24-28 =====
|
||
|
||
* Thursday: Introduction to ordinary Least Squares and derivation of basic equation
|
||
* Friday: Linear regression and statistical analysis and probability theory
|
||
|
||
|
||
!split
|
||
===== Thursday August 27 =====
|
||
|
||
"Video of Lecture":"https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h20/forelesningsvideoer/LectureAug27.mp4?vrtx=view-as-webpage".
|
||
|
||
|
||
!split
|
||
===== Why Linear Regression (aka Ordinary Least Squares and family) =====
|
||
|
||
Fitting a continuous function with linear parameterization in terms of the parameters $\bm{\beta}$.
|
||
* Method of choice for fitting a continuous function!
|
||
* Gives an excellent introduction to central Machine Learning features with _understandable pedagogical_ links to other methods like _Neural Networks_, _Support Vector Machines_ etc
|
||
* Analytical expression for the fitting parameters $\bm{\beta}$
|
||
* Analytical expressions for statistical propertiers like mean values, variances, confidence intervals and more
|
||
* Analytical relation with probabilistic interpretations
|
||
* Easy to introduce basic concepts like bias-variance tradeoff, cross-validation, resampling and regularization techniques and many other ML topics
|
||
* Easy to code! And links well with classification problems and logistic regression and neural networks
|
||
* Allows for _easy_ hands-on understanding of gradient descent methods
|
||
* and many more features
|
||
|
||
For more discussions of Ridge and Lasso regression, "Wessel van Wieringen's":"https://arxiv.org/abs/1509.09169" article is highly recommended.
|
||
Similarly, "Mehta et al's article":"https://arxiv.org/abs/1803.08823" is also recommended.
|
||
|
||
|
||
!split
|
||
===== Regression analysis, overarching aims =====
|
||
!bblock
|
||
|
||
Regression modeling deals with the description of the sampling distribution of a given random variable $y$ and how it varies as function of another variable or a set of such variables $\bm{x} =[x_0, x_1,\dots, x_{n-1}]^T$.
|
||
The first variable is called the _dependent_, the _outcome_ or the _response_ variable while the set of variables $\bm{x}$ is called the independent variable, or the predictor variable or the explanatory variable.
|
||
|
||
A regression model aims at finding a likelihood function $p(\bm{y}\vert \bm{x})$, that is the conditional distribution for $\bm{y}$ with a given $\bm{x}$. The estimation of $p(\bm{y}\vert \bm{x})$ is made using a data set with
|
||
* $n$ cases $i = 0, 1, 2, \dots, n-1$
|
||
* Response (target, dependent or outcome) variable $y_i$ with $i = 0, 1, 2, \dots, n-1$
|
||
* $p$ so-called explanatory (independent or predictor) variables $\bm{x}_i=[x_{i0}, x_{i1}, \dots, x_{ip-1}]$ with $i = 0, 1, 2, \dots, n-1$ and explanatory variables running from $0$ to $p-1$. See below for more explicit examples.
|
||
The goal of the regression analysis is to extract/exploit relationship between $\bm{y}$ and $\bm{x}$ in or to infer causal dependencies, approximations to the likelihood functions, functional relationships and to make predictions, making fits and many other things.
|
||
!eblock
|
||
|
||
!split
|
||
===== Regression analysis, overarching aims II =====
|
||
!bblock
|
||
|
||
|
||
Consider an experiment in which $p$ characteristics of $n$ samples are
|
||
measured. The data from this experiment, for various explanatory variables $p$ are normally represented by a matrix
|
||
$\mathbf{X}$.
|
||
|
||
The matrix $\mathbf{X}$ is called the *design
|
||
matrix*. Additional information of the samples is available in the
|
||
form of $\bm{y}$ (also as above). The variable $\bm{y}$ is
|
||
generally referred to as the *response variable*. The aim of
|
||
regression analysis is to explain $\bm{y}$ in terms of
|
||
$\bm{X}$ through a functional relationship like $y_i =
|
||
f(\mathbf{X}_{i,\ast})$. When no prior knowledge on the form of
|
||
$f(\cdot)$ is available, it is common to assume a linear relationship
|
||
between $\bm{X}$ and $\bm{y}$. This assumption gives rise to
|
||
the *linear regression model* where $\bm{\beta} = [\beta_0, \ldots,
|
||
\beta_{p-1}]^{T}$ are the *regression parameters*.
|
||
|
||
Linear regression gives us a set of analytical equations for the parameters $\beta_j$.
|
||
|
||
!eblock
|
||
|
||
|
||
|
||
!split
|
||
===== Examples =====
|
||
!bblock
|
||
In order to understand the relation among the predictors $p$, the set of data $n$ and the target (outcome, output etc) $\bm{y}$,
|
||
consider the model we discussed for describing nuclear binding energies.
|
||
|
||
There we assumed that we could parametrize the data using a polynomial approximation based on the liquid drop model.
|
||
Assuming
|
||
!bt
|
||
\[
|
||
BE(A) = a_0+a_1A+a_2A^{2/3}+a_3A^{-1/3}+a_4A^{-1},
|
||
\]
|
||
!et
|
||
we have five predictors, that is the intercept, the $A$ dependent term, the $A^{2/3}$ term and the $A^{-1/3}$ and $A^{-1}$ terms.
|
||
This gives $p=0,1,2,3,4$. Furthermore we have $n$ entries for each predictor. It means that our design matrix is a
|
||
$p\times n$ matrix $\bm{X}$.
|
||
|
||
Here the predictors are based on a model we have made. A popular data set which is widely encountered in ML applications is the
|
||
so-called "credit card default data from Taiwan":"https://www.sciencedirect.com/science/article/pii/S0957417407006719?via%3Dihub". The data set contains data on $n=30000$ credit card holders with predictors like gender, marital status, age, profession, education, etc. In total there are $24$ such predictors or attributes leading to a design matrix of dimensionality $24 \times 30000$. This is however a classification problem and we will come back to it when we discuss Logistic Regression.
|
||
|
||
!eblock
|
||
|
||
|
||
|
||
|
||
|
||
!split
|
||
===== General linear models =====
|
||
!bblock
|
||
Before we proceed let us study a case from linear algebra where we aim at fitting a set of data $\bm{y}=[y_0,y_1,\dots,y_{n-1}]$. We could think of these data as a result of an experiment or a complicated numerical experiment. These data are functions of a series of variables $\bm{x}=[x_0,x_1,\dots,x_{n-1}]$, that is $y_i = y(x_i)$ with $i=0,1,2,\dots,n-1$. The variables $x_i$ could represent physical quantities like time, temperature, position etc. We assume that $y(x)$ is a smooth function.
|
||
|
||
Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of $y$ which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree $n-1$ with $n$ points, that is
|
||
!bt
|
||
\[
|
||
y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \beta_j x_i^j+\epsilon_i,
|
||
\]
|
||
!et
|
||
where $\epsilon_i$ is the error in our approximation.
|
||
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Rewriting the fitting procedure as a linear algebra problem =====
|
||
!bblock
|
||
For every set of values $y_i,x_i$ we have thus the corresponding set of equations
|
||
!bt
|
||
\begin{align*}
|
||
y_0&=\beta_0+\beta_1x_0^1+\beta_2x_0^2+\dots+\beta_{n-1}x_0^{n-1}+\epsilon_0\\
|
||
y_1&=\beta_0+\beta_1x_1^1+\beta_2x_1^2+\dots+\beta_{n-1}x_1^{n-1}+\epsilon_1\\
|
||
y_2&=\beta_0+\beta_1x_2^1+\beta_2x_2^2+\dots+\beta_{n-1}x_2^{n-1}+\epsilon_2\\
|
||
\dots & \dots \\
|
||
y_{n-1}&=\beta_0+\beta_1x_{n-1}^1+\beta_2x_{n-1}^2+\dots+\beta_{n-1}x_{n-1}^{n-1}+\epsilon_{n-1}.\\
|
||
\end{align*}
|
||
!et
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Rewriting the fitting procedure as a linear algebra problem, more details =====
|
||
!bblock
|
||
Defining the vectors
|
||
!bt
|
||
\[
|
||
\bm{y} = [y_0,y_1, y_2,\dots, y_{n-1}]^T,
|
||
\]
|
||
!et
|
||
and
|
||
!bt
|
||
\[
|
||
\bm{\beta} = [\beta_0,\beta_1, \beta_2,\dots, \beta_{n-1}]^T,
|
||
\]
|
||
!et
|
||
and
|
||
!bt
|
||
\[
|
||
\bm{\epsilon} = [\epsilon_0,\epsilon_1, \epsilon_2,\dots, \epsilon_{n-1}]^T,
|
||
\]
|
||
!et
|
||
and the design matrix
|
||
!bt
|
||
\[
|
||
\bm{X}=
|
||
\begin{bmatrix}
|
||
1& x_{0}^1 &x_{0}^2& \dots & \dots &x_{0}^{n-1}\\
|
||
1& x_{1}^1 &x_{1}^2& \dots & \dots &x_{1}^{n-1}\\
|
||
1& x_{2}^1 &x_{2}^2& \dots & \dots &x_{2}^{n-1}\\
|
||
\dots& \dots &\dots& \dots & \dots &\dots\\
|
||
1& x_{n-1}^1 &x_{n-1}^2& \dots & \dots &x_{n-1}^{n-1}\\
|
||
\end{bmatrix}
|
||
\]
|
||
!et
|
||
we can rewrite our equations as
|
||
!bt
|
||
\[
|
||
\bm{y} = \bm{X}\bm{\beta}+\bm{\epsilon}.
|
||
\]
|
||
!et
|
||
The above design matrix is called a "Vandermonde matrix":"https://en.wikipedia.org/wiki/Vandermonde_matrix".
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Generalizing the fitting procedure as a linear algebra problem =====
|
||
!bblock
|
||
|
||
We are obviously not limited to the above polynomial expansions. We
|
||
could replace the various powers of $x$ with elements of Fourier
|
||
series or instead of $x_i^j$ we could have $\cos{(j x_i)}$ or $\sin{(j
|
||
x_i)}$, or time series or other orthogonal functions. For every set
|
||
of values $y_i,x_i$ we can then generalize the equations to
|
||
|
||
!bt
|
||
\begin{align*}
|
||
y_0&=\beta_0x_{00}+\beta_1x_{01}+\beta_2x_{02}+\dots+\beta_{n-1}x_{0n-1}+\epsilon_0\\
|
||
y_1&=\beta_0x_{10}+\beta_1x_{11}+\beta_2x_{12}+\dots+\beta_{n-1}x_{1n-1}+\epsilon_1\\
|
||
y_2&=\beta_0x_{20}+\beta_1x_{21}+\beta_2x_{22}+\dots+\beta_{n-1}x_{2n-1}+\epsilon_2\\
|
||
\dots & \dots \\
|
||
y_{i}&=\beta_0x_{i0}+\beta_1x_{i1}+\beta_2x_{i2}+\dots+\beta_{n-1}x_{in-1}+\epsilon_i\\
|
||
\dots & \dots \\
|
||
y_{n-1}&=\beta_0x_{n-1,0}+\beta_1x_{n-1,2}+\beta_2x_{n-1,2}+\dots+\beta_{n-1}x_{n-1,n-1}+\epsilon_{n-1}.\\
|
||
\end{align*}
|
||
!et
|
||
|
||
_Note that we have $p=n$ here. The matrix is symmetric. This is generally not the case!_
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Generalizing the fitting procedure as a linear algebra problem =====
|
||
!bblock
|
||
We redefine in turn the matrix $\bm{X}$ as
|
||
!bt
|
||
\[
|
||
\bm{X}=
|
||
\begin{bmatrix}
|
||
x_{00}& x_{01} &x_{02}& \dots & \dots &x_{0,n-1}\\
|
||
x_{10}& x_{11} &x_{12}& \dots & \dots &x_{1,n-1}\\
|
||
x_{20}& x_{21} &x_{22}& \dots & \dots &x_{2,n-1}\\
|
||
\dots& \dots &\dots& \dots & \dots &\dots\\
|
||
x_{n-1,0}& x_{n-1,1} &x_{n-1,2}& \dots & \dots &x_{n-1,n-1}\\
|
||
\end{bmatrix}
|
||
\]
|
||
!et
|
||
and without loss of generality we rewrite again our equations as
|
||
!bt
|
||
\[
|
||
\bm{y} = \bm{X}\bm{\beta}+\bm{\epsilon}.
|
||
\]
|
||
!et
|
||
The left-hand side of this equation is kwown. Our error vector $\bm{\epsilon}$ and the parameter vector $\bm{\beta}$ are our unknow quantities. How can we obtain the optimal set of $\beta_i$ values?
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Optimizing our parameters =====
|
||
!bblock
|
||
We have defined the matrix $\bm{X}$ via the equations
|
||
!bt
|
||
\begin{align*}
|
||
y_0&=\beta_0x_{00}+\beta_1x_{01}+\beta_2x_{02}+\dots+\beta_{n-1}x_{0n-1}+\epsilon_0\\
|
||
y_1&=\beta_0x_{10}+\beta_1x_{11}+\beta_2x_{12}+\dots+\beta_{n-1}x_{1n-1}+\epsilon_1\\
|
||
y_2&=\beta_0x_{20}+\beta_1x_{21}+\beta_2x_{22}+\dots+\beta_{n-1}x_{2n-1}+\epsilon_1\\
|
||
\dots & \dots \\
|
||
y_{i}&=\beta_0x_{i0}+\beta_1x_{i1}+\beta_2x_{i2}+\dots+\beta_{n-1}x_{in-1}+\epsilon_1\\
|
||
\dots & \dots \\
|
||
y_{n-1}&=\beta_0x_{n-1,0}+\beta_1x_{n-1,2}+\beta_2x_{n-1,2}+\dots+\beta_{n-1}x_{n-1,n-1}+\epsilon_{n-1}.\\
|
||
\end{align*}
|
||
!et
|
||
|
||
As we noted above, we stayed with a system with the design matrix
|
||
$\bm{X}\in {\mathbb{R}}^{n\times n}$, that is we have $p=n$. For reasons to come later (algorithmic arguments) we will hereafter define
|
||
our matrix as $\bm{X}\in {\mathbb{R}}^{n\times p}$, with the predictors refering to the column numbers and the entries $n$ being the row elements.
|
||
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Our model for the nuclear binding energies =====
|
||
|
||
In our "introductory notes":"https://compphysics.github.io/MachineLearning/doc/pub/How2ReadData/html/How2ReadData.html" we looked at the so-called "liquid drop model":"https://en.wikipedia.org/wiki/Semi-empirical_mass_formula". Let us remind ourselves about what we did by looking at the code.
|
||
|
||
We restate the parts of the code we are most interested in.
|
||
!bc pycod
|
||
# Common imports
|
||
import numpy as np
|
||
import pandas as pd
|
||
import matplotlib.pyplot as plt
|
||
from IPython.display import display
|
||
import os
|
||
|
||
# Where to save the figures and data files
|
||
PROJECT_ROOT_DIR = "Results"
|
||
FIGURE_ID = "Results/FigureFiles"
|
||
DATA_ID = "DataFiles/"
|
||
|
||
if not os.path.exists(PROJECT_ROOT_DIR):
|
||
os.mkdir(PROJECT_ROOT_DIR)
|
||
|
||
if not os.path.exists(FIGURE_ID):
|
||
os.makedirs(FIGURE_ID)
|
||
|
||
if not os.path.exists(DATA_ID):
|
||
os.makedirs(DATA_ID)
|
||
|
||
def image_path(fig_id):
|
||
return os.path.join(FIGURE_ID, fig_id)
|
||
|
||
def data_path(dat_id):
|
||
return os.path.join(DATA_ID, dat_id)
|
||
|
||
def save_fig(fig_id):
|
||
plt.savefig(image_path(fig_id) + ".png", format='png')
|
||
|
||
infile = open(data_path("MassEval2016.dat"),'r')
|
||
|
||
|
||
# Read the experimental data with Pandas
|
||
Masses = pd.read_fwf(infile, usecols=(2,3,4,6,11),
|
||
names=('N', 'Z', 'A', 'Element', 'Ebinding'),
|
||
widths=(1,3,5,5,5,1,3,4,1,13,11,11,9,1,2,11,9,1,3,1,12,11,1),
|
||
header=39,
|
||
index_col=False)
|
||
|
||
# Extrapolated values are indicated by '#' in place of the decimal place, so
|
||
# the Ebinding column won't be numeric. Coerce to float and drop these entries.
|
||
Masses['Ebinding'] = pd.to_numeric(Masses['Ebinding'], errors='coerce')
|
||
Masses = Masses.dropna()
|
||
# Convert from keV to MeV.
|
||
Masses['Ebinding'] /= 1000
|
||
|
||
# Group the DataFrame by nucleon number, A.
|
||
Masses = Masses.groupby('A')
|
||
# Find the rows of the grouped DataFrame with the maximum binding energy.
|
||
Masses = Masses.apply(lambda t: t[t.Ebinding==t.Ebinding.max()])
|
||
A = Masses['A']
|
||
Z = Masses['Z']
|
||
N = Masses['N']
|
||
Element = Masses['Element']
|
||
Energies = Masses['Ebinding']
|
||
|
||
# Now we set up the design matrix X
|
||
X = np.zeros((len(A),5))
|
||
X[:,0] = 1
|
||
X[:,1] = A
|
||
X[:,2] = A**(2.0/3.0)
|
||
X[:,3] = A**(-1.0/3.0)
|
||
X[:,4] = A**(-1.0)
|
||
# Then nice printout using pandas
|
||
DesignMatrix = pd.DataFrame(X)
|
||
DesignMatrix.index = A
|
||
DesignMatrix.columns = ['1', 'A', 'A^(2/3)', 'A^(-1/3)', '1/A']
|
||
display(DesignMatrix)
|
||
!ec
|
||
|
||
With $\bm{\beta}\in {\mathbb{R}}^{p\times 1}$, it means that we will hereafter write our equations for the approximation as
|
||
!bt
|
||
\[
|
||
\bm{\tilde{y}}= \bm{X}\bm{\beta},
|
||
\]
|
||
!et
|
||
throughout these lectures.
|
||
|
||
|
||
!split
|
||
===== Optimizing our parameters, more details =====
|
||
!bblock
|
||
With the above we use the design matrix to define the approximation $\bm{\tilde{y}}$ via the unknown quantity $\bm{\beta}$ as
|
||
!bt
|
||
\[
|
||
\bm{\tilde{y}}= \bm{X}\bm{\beta},
|
||
\]
|
||
!et
|
||
and in order to find the optimal parameters $\beta_i$ instead of solving the above linear algebra problem, we define a function which gives a measure of the spread between the values $y_i$ (which represent hopefully the exact values) and the parameterized values $\tilde{y}_i$, namely
|
||
!bt
|
||
\[
|
||
C(\bm{\beta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\bm{y}-\bm{\tilde{y}}\right)^T\left(\bm{y}-\bm{\tilde{y}}\right)\right\},
|
||
\]
|
||
!et
|
||
or using the matrix $\bm{X}$ and in a more compact matrix-vector notation as
|
||
!bt
|
||
\[
|
||
C(\bm{\beta})=\frac{1}{n}\left\{\left(\bm{y}-\bm{X}\bm{\beta}\right)^T\left(\bm{y}-\bm{X}\bm{\beta}\right)\right\}.
|
||
\]
|
||
!et
|
||
This function is one possible way to define the so-called cost function.
|
||
|
||
|
||
|
||
It is also common to define
|
||
the function $C$ as
|
||
|
||
!bt
|
||
\[
|
||
C(\bm{\beta})=\frac{1}{2n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2,
|
||
\]
|
||
!et
|
||
since when taking the first derivative with respect to the unknown parameters $\beta$, the factor of $2$ cancels out.
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Interpretations and optimizing our parameters =====
|
||
!bblock
|
||
|
||
The function
|
||
!bt
|
||
\[
|
||
C(\bm{\beta})=\frac{1}{n}\left\{\left(\bm{y}-\bm{X}\bm{\beta}\right)^T\left(\bm{y}-\bm{X}\bm{\beta}\right)\right\},
|
||
\]
|
||
!et
|
||
can be linked to the variance of the quantity $y_i$ if we interpret the latter as the mean value.
|
||
When linking (see the discussion below) with the maximum likelihood approach below, we will indeed interpret $y_i$ as a mean value
|
||
!bt
|
||
\[
|
||
y_{i}=\langle y_i \rangle = \beta_0x_{i,0}+\beta_1x_{i,1}+\beta_2x_{i,2}+\dots+\beta_{n-1}x_{i,n-1}+\epsilon_i,
|
||
\]
|
||
!et
|
||
|
||
where $\langle y_i \rangle$ is the mean value. Keep in mind also that
|
||
till now we have treated $y_i$ as the exact value. Normally, the
|
||
response (dependent or outcome) variable $y_i$ the outcome of a
|
||
numerical experiment or another type of experiment and is thus only an
|
||
approximation to the true value. It is then always accompanied by an
|
||
error estimate, often limited to a statistical error estimate given by
|
||
the standard deviation discussed earlier. In the discussion here we
|
||
will treat $y_i$ as our exact value for the response variable.
|
||
|
||
In order to find the parameters $\beta_i$ we will then minimize the spread of $C(\bm{\beta})$, that is we are going to solve the problem
|
||
!bt
|
||
\[
|
||
{\displaystyle \min_{\bm{\beta}\in
|
||
{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\bm{y}-\bm{X}\bm{\beta}\right)^T\left(\bm{y}-\bm{X}\bm{\beta}\right)\right\}.
|
||
\]
|
||
!et
|
||
In practical terms it means we will require
|
||
!bt
|
||
\[
|
||
\frac{\partial C(\bm{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}\right)^2\right]=0,
|
||
\]
|
||
!et
|
||
which results in
|
||
!bt
|
||
\[
|
||
\frac{\partial C(\bm{\beta})}{\partial \beta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_{ij}\left(y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}\right)\right]=0,
|
||
\]
|
||
!et
|
||
or in a matrix-vector form as
|
||
!bt
|
||
\[
|
||
\frac{\partial C(\bm{\beta})}{\partial \bm{\beta}} = 0 = \bm{X}^T\left( \bm{y}-\bm{X}\bm{\beta}\right).
|
||
\]
|
||
!et
|
||
|
||
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Interpretations and optimizing our parameters =====
|
||
!bblock
|
||
We can rewrite
|
||
!bt
|
||
\[
|
||
\frac{\partial C(\bm{\beta})}{\partial \bm{\beta}} = 0 = \bm{X}^T\left( \bm{y}-\bm{X}\bm{\beta}\right),
|
||
\]
|
||
!et
|
||
as
|
||
!bt
|
||
\[
|
||
\bm{X}^T\bm{y} = \bm{X}^T\bm{X}\bm{\beta},
|
||
\]
|
||
!et
|
||
and if the matrix $\bm{X}^T\bm{X}$ is invertible we have the solution
|
||
!bt
|
||
\[
|
||
\bm{\beta} =\left(\bm{X}^T\bm{X}\right)^{-1}\bm{X}^T\bm{y}.
|
||
\]
|
||
!et
|
||
|
||
We note also that since our design matrix is defined as $\bm{X}\in
|
||
{\mathbb{R}}^{n\times p}$, the product $\bm{X}^T\bm{X} \in
|
||
{\mathbb{R}}^{p\times p}$. In the above case we have that $p \ll n$,
|
||
in our case $p=5$ meaning that we end up with inverting a small
|
||
$5\times 5$ matrix. This is a rather common situation, in many cases we end up with low-dimensional
|
||
matrices to invert. The methods discussed here and for many other
|
||
supervised learning algorithms like classification with logistic
|
||
regression or support vector machines, exhibit dimensionalities which
|
||
allow for the usage of direct linear algebra methods such as _LU_ decomposition or _Singular Value Decomposition_ (SVD) for finding the inverse of the matrix
|
||
$\bm{X}^T\bm{X}$.
|
||
!eblock
|
||
|
||
!bblock
|
||
_Small question_: Do you think the example we have at hand here (the nuclear binding energies) can lead to problems in inverting the matrix $\bm{X}^T\bm{X}$? What kind of problems can we expect?
|
||
!eblock
|
||
|
||
!split
|
||
===== Some useful matrix and vector expressions =====
|
||
|
||
The following matrix and vector relation will be useful here and for the rest of the course. Vectors are always written as boldfaced lower case letters and
|
||
matrices as upper case boldfaced letters.
|
||
|
||
!bt
|
||
\[
|
||
\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},
|
||
\]
|
||
!et
|
||
!bt
|
||
\[
|
||
\frac{\partial tr(\bm{B}\bm{A})}{\partial \bm{A}} = \bm{B}^T,
|
||
\]
|
||
!et
|
||
!bt
|
||
\[
|
||
\frac{\partial \log{\vert\bm{A}\vert}}{\partial \bm{A}} = (\bm{A}^{-1})^T.
|
||
\]
|
||
!et
|
||
!split
|
||
===== Interpretations and optimizing our parameters =====
|
||
!bblock
|
||
The residuals $\bm{\epsilon}$ are in turn given by
|
||
!bt
|
||
\[
|
||
\bm{\epsilon} = \bm{y}-\bm{\tilde{y}} = \bm{y}-\bm{X}\bm{\beta},
|
||
\]
|
||
!et
|
||
and with
|
||
!bt
|
||
\[
|
||
\bm{X}^T\left( \bm{y}-\bm{X}\bm{\beta}\right)= 0,
|
||
\]
|
||
!et
|
||
we have
|
||
!bt
|
||
\[
|
||
\bm{X}^T\bm{\epsilon}=\bm{X}^T\left( \bm{y}-\bm{X}\bm{\beta}\right)= 0,
|
||
\]
|
||
!et
|
||
meaning that the solution for $\bm{\beta}$ is the one which minimizes the residuals. Later we will link this with the maximum likelihood approach.
|
||
|
||
!eblock
|
||
|
||
|
||
Let us now return to our nuclear binding energies and simply code the above equations.
|
||
|
||
!split
|
||
===== Own code for Ordinary Least Squares =====
|
||
|
||
It is rather straightforward to implement the matrix inversion and obtain the parameters $\bm{\beta}$. After having defined the matrix $\bm{X}$ we simply need to
|
||
write
|
||
!bc pycod
|
||
# matrix inversion to find beta
|
||
beta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Energies)
|
||
# and then make the prediction
|
||
ytilde = X @ beta
|
||
!ec
|
||
Alternatively, you can use the least squares functionality in _Numpy_ as
|
||
!bc pycod
|
||
fit = np.linalg.lstsq(X, Energies, rcond =None)[0]
|
||
ytildenp = np.dot(fit,X.T)
|
||
!ec
|
||
|
||
And finally we plot our fit with and compare with data
|
||
!bc pycod
|
||
Masses['Eapprox'] = ytilde
|
||
# Generate a plot comparing the experimental with the fitted values values.
|
||
fig, ax = plt.subplots()
|
||
ax.set_xlabel(r'$A = N + Z$')
|
||
ax.set_ylabel(r'$E_\mathrm{bind}\,/\mathrm{MeV}$')
|
||
ax.plot(Masses['A'], Masses['Ebinding'], alpha=0.7, lw=2,
|
||
label='Ame2016')
|
||
ax.plot(Masses['A'], Masses['Eapprox'], alpha=0.7, lw=2, c='m',
|
||
label='Fit')
|
||
ax.legend()
|
||
save_fig("Masses2016OLS")
|
||
plt.show()
|
||
!ec
|
||
|
||
!split
|
||
===== Adding error analysis and training set up =====
|
||
|
||
We can easily test our fit by computing the $R2$ score that we discussed in connection with the functionality of _Scikit-Learn_ in the introductory slides.
|
||
Since we are not using _Scikit-Learn_ here we can define our own $R2$ function as
|
||
!bc pycod
|
||
def R2(y_data, y_model):
|
||
return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
|
||
!ec
|
||
and we would be using it as
|
||
!bc pycod
|
||
print(R2(Energies,ytilde))
|
||
!ec
|
||
|
||
We can easily add our _MSE_ score as
|
||
!bc pycod
|
||
def MSE(y_data,y_model):
|
||
n = np.size(y_model)
|
||
return np.sum((y_data-y_model)**2)/n
|
||
|
||
print(MSE(Energies,ytilde))
|
||
!ec
|
||
and finally the relative error as
|
||
!bc pycod
|
||
def RelativeError(y_data,y_model):
|
||
return abs((y_data-y_model)/y_data)
|
||
print(RelativeError(Energies, ytilde))
|
||
!ec
|
||
|
||
|
||
|
||
!split
|
||
===== The $\chi^2$ function =====
|
||
!bblock
|
||
|
||
Normally, the response (dependent or outcome) variable $y_i$ is the
|
||
outcome of a numerical experiment or another type of experiment and is
|
||
thus only an approximation to the true value. It is then always
|
||
accompanied by an error estimate, often limited to a statistical error
|
||
estimate given by the standard deviation discussed earlier. In the
|
||
discussion here we will treat $y_i$ as our exact value for the
|
||
response variable.
|
||
|
||
Introducing the standard deviation $\sigma_i$ for each measurement
|
||
$y_i$, we define now the $\chi^2$ function (omitting the $1/n$ term)
|
||
as
|
||
|
||
!bt
|
||
\[
|
||
\chi^2(\bm{\beta})=\frac{1}{n}\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\frac{1}{n}\left\{\left(\bm{y}-\bm{\tilde{y}}\right)^T\frac{1}{\bm{\Sigma^2}}\left(\bm{y}-\bm{\tilde{y}}\right)\right\},
|
||
\]
|
||
!et
|
||
where the matrix $\bm{\Sigma}$ is a diagonal matrix with $\sigma_i$ as matrix elements.
|
||
|
||
!eblock
|
||
|
||
!split
|
||
===== The $\chi^2$ function =====
|
||
!bblock
|
||
|
||
In order to find the parameters $\beta_i$ we will then minimize the spread of $\chi^2(\bm{\beta})$ by requiring
|
||
!bt
|
||
\[
|
||
\frac{\partial \chi^2(\bm{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0,
|
||
\]
|
||
!et
|
||
which results in
|
||
!bt
|
||
\[
|
||
\frac{\partial \chi^2(\bm{\beta})}{\partial \beta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0,
|
||
\]
|
||
!et
|
||
or in a matrix-vector form as
|
||
!bt
|
||
\[
|
||
\frac{\partial \chi^2(\bm{\beta})}{\partial \bm{\beta}} = 0 = \bm{A}^T\left( \bm{b}-\bm{A}\bm{\beta}\right).
|
||
\]
|
||
!et
|
||
where we have defined the matrix $\bm{A} =\bm{X}/\bm{\Sigma}$ with matrix elements $a_{ij} = x_{ij}/\sigma_i$ and the vector $\bm{b}$ with elements $b_i = y_i/\sigma_i$.
|
||
!eblock
|
||
|
||
!split
|
||
===== The $\chi^2$ function =====
|
||
!bblock
|
||
|
||
We can rewrite
|
||
!bt
|
||
\[
|
||
\frac{\partial \chi^2(\bm{\beta})}{\partial \bm{\beta}} = 0 = \bm{A}^T\left( \bm{b}-\bm{A}\bm{\beta}\right),
|
||
\]
|
||
!et
|
||
as
|
||
!bt
|
||
\[
|
||
\bm{A}^T\bm{b} = \bm{A}^T\bm{A}\bm{\beta},
|
||
\]
|
||
!et
|
||
and if the matrix $\bm{A}^T\bm{A}$ is invertible we have the solution
|
||
!bt
|
||
\[
|
||
\bm{\beta} =\left(\bm{A}^T\bm{A}\right)^{-1}\bm{A}^T\bm{b}.
|
||
\]
|
||
!et
|
||
!eblock
|
||
|
||
!split
|
||
===== The $\chi^2$ function =====
|
||
!bblock
|
||
|
||
If we then introduce the matrix
|
||
!bt
|
||
\[
|
||
\bm{H} = \left(\bm{A}^T\bm{A}\right)^{-1},
|
||
\]
|
||
!et
|
||
we have then the following expression for the parameters $\beta_j$ (the matrix elements of $\bm{H}$ are $h_{ij}$)
|
||
!bt
|
||
\[
|
||
\beta_j = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}\frac{y_i}{\sigma_i}\frac{x_{ik}}{\sigma_i} = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}b_ia_{ik}
|
||
\]
|
||
!et
|
||
We state without proof the expression for the uncertainty in the parameters $\beta_j$ as (we leave this as an exercise)
|
||
!bt
|
||
\[
|
||
\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2,
|
||
\]
|
||
!et
|
||
resulting in
|
||
!bt
|
||
\[
|
||
\sigma^2(\beta_j) = \left(\sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}a_{ik}\right)\left(\sum_{l=0}^{p-1}h_{jl}\sum_{m=0}^{n-1}a_{ml}\right) = h_{jj}!
|
||
\]
|
||
!et
|
||
!eblock
|
||
|
||
!split
|
||
===== The $\chi^2$ function =====
|
||
!bblock
|
||
The first step here is to approximate the function $y$ with a first-order polynomial, that is we write
|
||
!bt
|
||
\[
|
||
y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i.
|
||
\]
|
||
!et
|
||
By computing the derivatives of $\chi^2$ with respect to $\beta_0$ and $\beta_1$ show that these are given by
|
||
!bt
|
||
\[
|
||
\frac{\partial \chi^2(\bm{\beta})}{\partial \beta_0} = -2\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0,
|
||
\]
|
||
!et
|
||
and
|
||
!bt
|
||
\[
|
||
\frac{\partial \chi^2(\bm{\beta})}{\partial \beta_1} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_i\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0.
|
||
\]
|
||
!et
|
||
!eblock
|
||
|
||
!split
|
||
===== The $\chi^2$ function =====
|
||
!bblock
|
||
|
||
For a linear fit (a first-order polynomial) we don't need to invert a matrix!!
|
||
Defining
|
||
!bt
|
||
\[
|
||
\gamma = \sum_{i=0}^{n-1}\frac{1}{\sigma_i^2},
|
||
\]
|
||
!et
|
||
|
||
!bt
|
||
\[
|
||
\gamma_x = \sum_{i=0}^{n-1}\frac{x_{i}}{\sigma_i^2},
|
||
\]
|
||
!et
|
||
|
||
!bt
|
||
\[
|
||
\gamma_y = \sum_{i=0}^{n-1}\left(\frac{y_i}{\sigma_i^2}\right),
|
||
\]
|
||
!et
|
||
|
||
!bt
|
||
\[
|
||
\gamma_{xx} = \sum_{i=0}^{n-1}\frac{x_ix_{i}}{\sigma_i^2},
|
||
\]
|
||
!et
|
||
|
||
!bt
|
||
\[
|
||
\gamma_{xy} = \sum_{i=0}^{n-1}\frac{y_ix_{i}}{\sigma_i^2},
|
||
\]
|
||
!et
|
||
|
||
we obtain
|
||
|
||
!bt
|
||
\[
|
||
\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
|
||
\]
|
||
!et
|
||
|
||
!bt
|
||
\[
|
||
\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}.
|
||
\]
|
||
!et
|
||
|
||
This approach (different linear and non-linear regression) suffers
|
||
often from both being underdetermined and overdetermined in the
|
||
unknown coefficients $\beta_i$. A better approach is to use the
|
||
Singular Value Decomposition (SVD) method discussed below. Or using
|
||
Lasso and Ridge regression. See below.
|
||
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Fitting an Equation of State for Dense Nuclear Matter =====
|
||
|
||
Before we continue, let us introduce yet another example. We are going to fit the
|
||
nuclear equation of state using results from many-body calculations.
|
||
The equation of state we have made available here, as function of
|
||
density, has been derived using modern nucleon-nucleon potentials with
|
||
"the addition of three-body
|
||
forces":"https://www.sciencedirect.com/science/article/pii/S0370157399001106". This
|
||
time the file is presented as a standard _csv_ file.
|
||
|
||
The beginning of the Python code here is similar to what you have seen
|
||
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.
|
||
|
||
!split
|
||
===== The code =====
|
||
|
||
!bc pycod
|
||
# Common imports
|
||
import os
|
||
import numpy as np
|
||
import pandas as pd
|
||
import matplotlib.pyplot as plt
|
||
import matplotlib.pyplot as plt
|
||
import sklearn.linear_model as skl
|
||
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
|
||
|
||
# Where to save the figures and data files
|
||
PROJECT_ROOT_DIR = "Results"
|
||
FIGURE_ID = "Results/FigureFiles"
|
||
DATA_ID = "DataFiles/"
|
||
|
||
if not os.path.exists(PROJECT_ROOT_DIR):
|
||
os.mkdir(PROJECT_ROOT_DIR)
|
||
|
||
if not os.path.exists(FIGURE_ID):
|
||
os.makedirs(FIGURE_ID)
|
||
|
||
if not os.path.exists(DATA_ID):
|
||
os.makedirs(DATA_ID)
|
||
|
||
def image_path(fig_id):
|
||
return os.path.join(FIGURE_ID, fig_id)
|
||
|
||
def data_path(dat_id):
|
||
return os.path.join(DATA_ID, dat_id)
|
||
|
||
def save_fig(fig_id):
|
||
plt.savefig(image_path(fig_id) + ".png", format='png')
|
||
|
||
infile = open(data_path("EoS.csv"),'r')
|
||
|
||
# Read the EoS data as csv file and organize the data into two arrays with density and energies
|
||
EoS = pd.read_csv(infile, names=('Density', 'Energy'))
|
||
EoS['Energy'] = pd.to_numeric(EoS['Energy'], errors='coerce')
|
||
EoS = EoS.dropna()
|
||
Energies = EoS['Energy']
|
||
Density = EoS['Density']
|
||
# The design matrix now as function of various polytrops
|
||
X = np.zeros((len(Density),4))
|
||
X[:,3] = Density**(4.0/3.0)
|
||
X[:,2] = Density
|
||
X[:,1] = Density**(2.0/3.0)
|
||
X[:,0] = 1
|
||
|
||
# We use now Scikit-Learn's linear regressor and ridge regressor
|
||
# OLS part
|
||
clf = skl.LinearRegression().fit(X, Energies)
|
||
ytilde = clf.predict(X)
|
||
EoS['Eols'] = ytilde
|
||
# The mean squared error
|
||
print("Mean squared error: %.2f" % mean_squared_error(Energies, ytilde))
|
||
# Explained variance score: 1 is perfect prediction
|
||
print('Variance score: %.2f' % r2_score(Energies, ytilde))
|
||
# Mean absolute error
|
||
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}]$')
|
||
ax.set_ylabel(r'Energy per particle')
|
||
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()
|
||
!ec
|
||
|
||
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.
|
||
|
||
|
||
!split
|
||
===== 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.
|
||
|
||
!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
|
||
# Where to save the figures and data files
|
||
PROJECT_ROOT_DIR = "Results"
|
||
FIGURE_ID = "Results/FigureFiles"
|
||
DATA_ID = "DataFiles/"
|
||
|
||
if not os.path.exists(PROJECT_ROOT_DIR):
|
||
os.mkdir(PROJECT_ROOT_DIR)
|
||
|
||
if not os.path.exists(FIGURE_ID):
|
||
os.makedirs(FIGURE_ID)
|
||
|
||
if not os.path.exists(DATA_ID):
|
||
os.makedirs(DATA_ID)
|
||
|
||
def image_path(fig_id):
|
||
return os.path.join(FIGURE_ID, fig_id)
|
||
|
||
def data_path(dat_id):
|
||
return os.path.join(DATA_ID, dat_id)
|
||
|
||
def save_fig(fig_id):
|
||
plt.savefig(image_path(fig_id) + ".png", format='png')
|
||
|
||
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
|
||
|
||
infile = open(data_path("EoS.csv"),'r')
|
||
|
||
# Read the EoS data as csv file and organized into two arrays with density and energies
|
||
EoS = pd.read_csv(infile, names=('Density', 'Energy'))
|
||
EoS['Energy'] = pd.to_numeric(EoS['Energy'], errors='coerce')
|
||
EoS = EoS.dropna()
|
||
Energies = EoS['Energy']
|
||
Density = EoS['Density']
|
||
# The design matrix now as function of various polytrops
|
||
X = np.zeros((len(Density),5))
|
||
X[:,0] = 1
|
||
X[:,1] = Density**(2.0/3.0)
|
||
X[:,2] = Density
|
||
X[:,3] = Density**(4.0/3.0)
|
||
X[:,4] = Density**(5.0/3.0)
|
||
# We split the data in test and training data
|
||
X_train, X_test, y_train, y_test = train_test_split(X, Energies, test_size=0.2)
|
||
# matrix inversion to find beta
|
||
beta = np.linalg.inv(X_train.T.dot(X_train)).dot(X_train.T).dot(y_train)
|
||
# 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
|
||
|
||
|
||
!split
|
||
===== The Boston housing data example =====
|
||
|
||
The Boston housing
|
||
data set was originally a part of UCI Machine Learning Repository
|
||
and has been removed now. The data set is now included in _Scikit-Learn_'s
|
||
library. There are 506 samples and 13 feature (predictor) variables
|
||
in this data set. The objective is to predict the value of prices of
|
||
the house using the features (predictors) listed here.
|
||
|
||
The features/predictors are
|
||
o CRIM: Per capita crime rate by town
|
||
o ZN: Proportion of residential land zoned for lots over 25000 square feet
|
||
o INDUS: Proportion of non-retail business acres per town
|
||
o CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
|
||
o NOX: Nitric oxide concentration (parts per 10 million)
|
||
o RM: Average number of rooms per dwelling
|
||
o AGE: Proportion of owner-occupied units built prior to 1940
|
||
o DIS: Weighted distances to five Boston employment centers
|
||
o RAD: Index of accessibility to radial highways
|
||
o TAX: Full-value property tax rate per USD10000
|
||
o B: $1000(Bk - 0.63)^2$, where $Bk$ is the proportion of [people of African American descent] by town
|
||
o LSTAT: Percentage of lower status of the population
|
||
o MEDV: Median value of owner-occupied homes in USD 1000s
|
||
|
||
!split
|
||
===== Housing data, the code =====
|
||
We start by importing the libraries
|
||
!bc pycod
|
||
import numpy as np
|
||
import matplotlib.pyplot as plt
|
||
|
||
import pandas as pd
|
||
import seaborn as sns
|
||
!ec
|
||
and load the Boston Housing DataSet from _Scikit-Learn_
|
||
|
||
|
||
!bc pycod
|
||
from sklearn.datasets import load_boston
|
||
|
||
boston_dataset = load_boston()
|
||
|
||
# boston_dataset is a dictionary
|
||
# let's check what it contains
|
||
boston_dataset.keys()
|
||
!ec
|
||
Then we invoke Pandas
|
||
!bc pycod
|
||
boston = pd.DataFrame(boston_dataset.data, columns=boston_dataset.feature_names)
|
||
boston.head()
|
||
boston['MEDV'] = boston_dataset.target
|
||
!ec
|
||
and preprocess the data
|
||
!bc pycod
|
||
# check for missing values in all the columns
|
||
boston.isnull().sum()
|
||
!ec
|
||
We can then visualize the data
|
||
!bc pycod
|
||
# set the size of the figure
|
||
sns.set(rc={'figure.figsize':(11.7,8.27)})
|
||
|
||
# plot a histogram showing the distribution of the target values
|
||
sns.distplot(boston['MEDV'], bins=30)
|
||
plt.show()
|
||
!ec
|
||
|
||
It is now useful to look at the correlation matrix
|
||
!bc pycod
|
||
# compute the pair wise correlation for all columns
|
||
correlation_matrix = boston.corr().round(2)
|
||
# use the heatmap function from seaborn to plot the correlation matrix
|
||
# annot = True to print the values inside the square
|
||
sns.heatmap(data=correlation_matrix, annot=True)
|
||
!ec
|
||
From the above coorelation plot we can see that _MEDV_ is strongly correlated to _LSTAT_ and _RM_. We see also that _RAD_ and _TAX_ are stronly correlated, but we don't include this in our features together to avoid multi-colinearity
|
||
|
||
!bc pycod
|
||
plt.figure(figsize=(20, 5))
|
||
|
||
features = ['LSTAT', 'RM']
|
||
target = boston['MEDV']
|
||
|
||
for i, col in enumerate(features):
|
||
plt.subplot(1, len(features) , i+1)
|
||
x = boston[col]
|
||
y = target
|
||
plt.scatter(x, y, marker='o')
|
||
plt.title(col)
|
||
plt.xlabel(col)
|
||
plt.ylabel('MEDV')
|
||
!ec
|
||
Now we start training our model
|
||
!bc pycod
|
||
X = pd.DataFrame(np.c_[boston['LSTAT'], boston['RM']], columns = ['LSTAT','RM'])
|
||
Y = boston['MEDV']
|
||
!ec
|
||
We split the data into training and test sets
|
||
|
||
!bc pycod
|
||
from sklearn.model_selection import train_test_split
|
||
|
||
# splits the training and test data set in 80% : 20%
|
||
# assign random_state to any value.This ensures consistency.
|
||
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2, random_state=5)
|
||
print(X_train.shape)
|
||
print(X_test.shape)
|
||
print(Y_train.shape)
|
||
print(Y_test.shape)
|
||
!ec
|
||
Then we use the linear regression functionality from _Scikit-Learn_
|
||
!bc pycod
|
||
from sklearn.linear_model import LinearRegression
|
||
from sklearn.metrics import mean_squared_error, r2_score
|
||
|
||
lin_model = LinearRegression()
|
||
lin_model.fit(X_train, Y_train)
|
||
|
||
# model evaluation for training set
|
||
|
||
y_train_predict = lin_model.predict(X_train)
|
||
rmse = (np.sqrt(mean_squared_error(Y_train, y_train_predict)))
|
||
r2 = r2_score(Y_train, y_train_predict)
|
||
|
||
print("The model performance for training set")
|
||
print("--------------------------------------")
|
||
print('RMSE is {}'.format(rmse))
|
||
print('R2 score is {}'.format(r2))
|
||
print("\n")
|
||
|
||
# model evaluation for testing set
|
||
|
||
y_test_predict = lin_model.predict(X_test)
|
||
# root mean square error of the model
|
||
rmse = (np.sqrt(mean_squared_error(Y_test, y_test_predict)))
|
||
|
||
# r-squared score of the model
|
||
r2 = r2_score(Y_test, y_test_predict)
|
||
|
||
print("The model performance for testing set")
|
||
print("--------------------------------------")
|
||
print('RMSE is {}'.format(rmse))
|
||
print('R2 score is {}'.format(r2))
|
||
!ec
|
||
|
||
!bc pycod
|
||
# plotting the y_test vs y_pred
|
||
# ideally should have been a straight line
|
||
plt.scatter(Y_test, y_test_predict)
|
||
plt.show()
|
||
!ec
|
||
|
||
|
||
!split
|
||
===== Reducing the number of degrees of freedom, overarching view =====
|
||
!bblock
|
||
|
||
Many Machine Learning problems involve thousands or even millions of
|
||
features for each training instance. Not only does this make training
|
||
extremely slow, it can also make it much harder to find a good
|
||
solution, as we will see. This problem is often referred to as the
|
||
curse of dimensionality. Fortunately, in real-world problems, it is
|
||
often possible to reduce the number of features considerably, turning
|
||
an intractable problem into a tractable one.
|
||
|
||
Later we will discuss some of the most popular dimensionality reduction
|
||
techniques: the principal component analysis (PCA), Kernel PCA, and
|
||
Locally Linear Embedding (LLE).
|
||
|
||
|
||
Principal component analysis and its various variants deal with the
|
||
problem of fitting a low-dimensional "affine
|
||
subspace":"https://en.wikipedia.org/wiki/Affine_space" to a set of of
|
||
data points in a high-dimensional space. With its family of methods it
|
||
is one of the most used tools in data modeling, compression and
|
||
visualization.
|
||
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Preprocessing our data =====
|
||
!bblock
|
||
|
||
Before we proceed however, we will discuss how to preprocess our
|
||
data. Till now and in connection with our previous examples we have
|
||
not met so many cases where we are too sensitive to the scaling of our
|
||
data. Normally the data may need a rescaling and/or may be sensitive
|
||
to extreme values. Scaling the data renders our inputs much more
|
||
suitable for the algorithms we want to employ.
|
||
|
||
_Scikit-Learn_ has several functions which allow us to rescale the
|
||
data, normally resulting in much better results in terms of various
|
||
accuracy scores. The _StandardScaler_ function in _Scikit-Learn_
|
||
ensures that for each feature/predictor we study the mean value is
|
||
zero and the variance is one (every column in the design/feature
|
||
matrix). This scaling has the drawback that it does not ensure that
|
||
we have a particular maximum or minimum in our data set. Another
|
||
function included in _Scikit-Learn_ is the _MinMaxScaler_ which
|
||
ensures that all features are exactly between $0$ and $1$. The
|
||
|
||
!split
|
||
===== More preprocessing =====
|
||
|
||
|
||
The _Normalizer_ scales each data
|
||
point such that the feature vector has a euclidean length of one. In other words, it
|
||
projects a data point on the circle (or sphere in the case of higher dimensions) with a
|
||
radius of 1. This means every data point is scaled by a different number (by the
|
||
inverse of it’s length).
|
||
This normalization is often used when only the direction (or angle) of the data matters,
|
||
not the length of the feature vector.
|
||
|
||
The _RobustScaler_ works similarly to the StandardScaler in that it
|
||
ensures statistical properties for each feature that guarantee that
|
||
they are on the same scale. However, the RobustScaler uses the median
|
||
and quartiles, instead of mean and variance. This makes the
|
||
RobustScaler ignore data points that are very different from the rest
|
||
(like measurement errors). These odd data points are also called
|
||
outliers, and might often lead to trouble for other scaling
|
||
techniques.
|
||
|
||
!eblock
|
||
|
||
!split
|
||
===== Simple preprocessing examples, Franke function and regression =====
|
||
|
||
!bc pycod
|
||
# Common imports
|
||
import os
|
||
import numpy as np
|
||
import pandas as pd
|
||
import matplotlib.pyplot as plt
|
||
import sklearn.linear_model as skl
|
||
from sklearn.metrics import mean_squared_error
|
||
from sklearn.model_selection import train_test_split
|
||
from sklearn.preprocessing import MinMaxScaler, StandardScaler, Normalizer
|
||
|
||
# Where to save the figures and data files
|
||
PROJECT_ROOT_DIR = "Results"
|
||
FIGURE_ID = "Results/FigureFiles"
|
||
DATA_ID = "DataFiles/"
|
||
|
||
if not os.path.exists(PROJECT_ROOT_DIR):
|
||
os.mkdir(PROJECT_ROOT_DIR)
|
||
|
||
if not os.path.exists(FIGURE_ID):
|
||
os.makedirs(FIGURE_ID)
|
||
|
||
if not os.path.exists(DATA_ID):
|
||
os.makedirs(DATA_ID)
|
||
|
||
def image_path(fig_id):
|
||
return os.path.join(FIGURE_ID, fig_id)
|
||
|
||
def data_path(dat_id):
|
||
return os.path.join(DATA_ID, dat_id)
|
||
|
||
def save_fig(fig_id):
|
||
plt.savefig(image_path(fig_id) + ".png", format='png')
|
||
|
||
|
||
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)
|
||
# split in training and test data
|
||
X_train, X_test, y_train, y_test = train_test_split(X,z,test_size=0.2)
|
||
|
||
|
||
clf = skl.LinearRegression().fit(X_train, y_train)
|
||
|
||
# The mean squared error and R2 score
|
||
print("MSE before scaling: {:.2f}".format(mean_squared_error(clf.predict(X_test), y_test)))
|
||
print("R2 score before scaling {:.2f}".format(clf.score(X_test,y_test)))
|
||
|
||
scaler = StandardScaler()
|
||
scaler.fit(X_train)
|
||
X_train_scaled = scaler.transform(X_train)
|
||
X_test_scaled = scaler.transform(X_test)
|
||
|
||
print("Feature min values before scaling:\n {}".format(X_train.min(axis=0)))
|
||
print("Feature max values before scaling:\n {}".format(X_train.max(axis=0)))
|
||
|
||
print("Feature min values after scaling:\n {}".format(X_train_scaled.min(axis=0)))
|
||
print("Feature max values after scaling:\n {}".format(X_train_scaled.max(axis=0)))
|
||
|
||
clf = skl.LinearRegression().fit(X_train_scaled, y_train)
|
||
|
||
|
||
print("MSE after scaling: {:.2f}".format(mean_squared_error(clf.predict(X_test_scaled), y_test)))
|
||
print("R2 score for scaled data: {:.2f}".format(clf.score(X_test_scaled,y_test)))
|
||
|
||
!ec
|
||
|
||
|
||
!split
|
||
===== Friday August 28 =====
|
||
|
||
"Video of Lecture":"https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h20/forelesningsvideoer/LectureAug28.mp4?vrtx=view-as-webpage" and "handwritten notes":"https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/NotesAugust28.pdf"
|
||
|
||
More material will be added here, see handwritten notes also.
|
||
|
||
|