From 65052fe79a53545f618a51fbf2fd7e0e56fb0d4b Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Mon, 6 Sep 2021 21:59:45 +0200 Subject: [PATCH] Delete week35_2020.do.txt --- doc/src/week35/week35_2020.do.txt | 1311 ----------------------------- 1 file changed, 1311 deletions(-) delete mode 100644 doc/src/week35/week35_2020.do.txt diff --git a/doc/src/week35/week35_2020.do.txt b/doc/src/week35/week35_2020.do.txt deleted file mode 100644 index cce129cf0..000000000 --- a/doc/src/week35/week35_2020.do.txt +++ /dev/null @@ -1,1311 +0,0 @@ -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. - -