From 6fe1e826984b114e5c7edb81aec03c7269dab893 Mon Sep 17 00:00:00 2001
From: Morten Hjorth-Jensen Before we proceed, we define also a function for making our plots. You can obviously avoid this and simply set up various matplotlib commands every time you need them. You may however find it convenient to collect all such commands in one function and simply call this function. Our next step is to read the data on experimental binding energies and
reorganize them as functions of the mass number \( A \), the number of
protons \( Z \) and neutrons \( N \) using pandas. Before we do this it is
@@ -1192,7 +1146,7 @@ Now it is time to dive more into the details of various methods. We will start w
Fitting a continuous function with linear parameterization in terms of the parameters \( \boldsymbol{\beta} \). Fitting a continuous function with linear parameterization in terms of the parameters \( \boldsymbol{\theta} \). 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 \( \boldsymbol{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 \( \boldsymbol{x} \) is called the independent variable, or the predictor variable or the explanatory variable, or simply just the inputs.
+The first variable \( y \) is called the the outcome or the response variable, or simply just the outputs.
The set of variables \( \boldsymbol{x} \) is called the independent variable, or the predictor variable or the explanatory variable, or simply just the inputs. We will throughout the course just use inputs and outputs as names. A regression model aims at finding a likelihood function \( p(\boldsymbol{y}\vert \boldsymbol{x}) \) or in the more traditional sense a function \( \boldsymbol{y}(\boldsymbol{x}) \), that is the conditional distribution for \( \boldsymbol{y} \) with a given \( \boldsymbol{x} \). The estimation of \( p(\boldsymbol{y}\vert \boldsymbol{x}) \) is made using a data set with The goal of the regression analysis is to extract/exploit relationship between \( \boldsymbol{y} \) and \( \boldsymbol{x} \) in order to infer specific dependencies, approximations to the likelihood functions, functional relationships and to make predictions, making fits and many other things.from pylab import plt, mpl
-plt.style.use('seaborn')
-mpl.rcParams['font.family'] = 'serif'
-
-def MakePlot(x,y, styles, labels, axlabels):
- plt.figure(figsize=(10,6))
- for i in range(len(x)):
- plt.plot(x[i], y[i], styles[i], label = labels[i])
- plt.xlabel(axlabels[0])
- plt.ylabel(axlabels[1])
- plt.legend(loc=0)
-
-Why Linear Regression (aka Ordinary Least Squares and family)
-
diff --git a/doc/pub/week34/html/._week34-bs046.html b/doc/pub/week34/html/._week34-bs046.html
index ddd27c066..e21b3fcf9 100644
--- a/doc/pub/week34/html/._week34-bs046.html
+++ b/doc/pub/week34/html/._week34-bs046.html
@@ -232,16 +232,7 @@ doconce format html week34.do.txt --html_style=bootstrap --pygments_html_style=d
('The $\\chi^2$ function', 2, None, 'the-chi-2-function'),
('The $\\chi^2$ function', 2, None, 'the-chi-2-function'),
('The $\\chi^2$ function', 2, None, 'the-chi-2-function'),
- ('The $\\chi^2$ function', 2, None, 'the-chi-2-function'),
- ('Fitting an Equation of State for Dense Nuclear Matter',
- 2,
- None,
- 'fitting-an-equation-of-state-for-dense-nuclear-matter'),
- ('The code', 2, None, 'the-code'),
- ('Splitting our Data in Training and Test data',
- 2,
- None,
- 'splitting-our-data-in-training-and-test-data')]}
+ ('The $\\chi^2$ function', 2, None, 'the-chi-2-function')]}
end of tocinfo -->
@@ -347,9 +338,6 @@ MathJax.Hub.Config({
Linear regression gives us a set of analytical equations for the parameters \( \beta_j \).
+Linear regression gives us a set of analytical equations for the parameters \( \theta_j \).
@@ -414,7 +402,7 @@ the linear regression model where \( \boldsymbol{\beta} = [\beta_0, \ldSince 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
$$ -y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \beta_j x_i^j+\epsilon_i, +y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \theta_j x_i^j+\epsilon_i, $$where \( \epsilon_i \) is the error in our approximation.
@@ -402,7 +390,7 @@ $$For every set of values \( y_i,x_i \) we have thus the corresponding set of equations
$$ \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\\ +y_0&=\theta_0+\theta_1x_0^1+\theta_2x_0^2+\dots+\theta_{n-1}x_0^{n-1}+\epsilon_0\\ +y_1&=\theta_0+\theta_1x_1^1+\theta_2x_1^2+\dots+\theta_{n-1}x_1^{n-1}+\epsilon_1\\ +y_2&=\theta_0+\theta_1x_2^1+\theta_2x_2^2+\dots+\theta_{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}.\\ +y_{n-1}&=\theta_0+\theta_1x_{n-1}^1+\theta_2x_{n-1}^2+\dots+\theta_{n-1}x_{n-1}^{n-1}+\epsilon_{n-1}.\\ \end{align*} $$ @@ -404,7 +392,7 @@ $$and
$$ -\boldsymbol{\beta} = [\beta_0,\beta_1, \beta_2,\dots, \beta_{n-1}]^T, +\boldsymbol{\theta} = [\theta_0,\theta_1, \theta_2,\dots, \theta_{n-1}]^T, $$and
@@ -394,7 +382,7 @@ $$we can rewrite our equations as
$$ -\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}. +\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\theta}+\boldsymbol{\epsilon}. $$The above design matrix is called a Vandermonde matrix.
@@ -427,7 +415,7 @@ $$and without loss of generality we rewrite again our equations as
$$ -\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}. +\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\theta}+\boldsymbol{\epsilon}. $$ -The left-hand side of this equation is kwown. Our error vector \( \boldsymbol{\epsilon} \) and the parameter vector \( \boldsymbol{\beta} \) are our unknow quantities. How can we obtain the optimal set of \( \beta_i \) values?
+The left-hand side of this equation is kwown. Our error vector \( \boldsymbol{\epsilon} \) and the parameter vector \( \boldsymbol{\theta} \) are our unknow quantities. How can we obtain the optimal set of \( \theta_i \) values?
@@ -412,7 +400,7 @@ $$We have defined the matrix \( \boldsymbol{X} \) via the equations
$$ \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\\ +y_0&=\theta_0x_{00}+\theta_1x_{01}+\theta_2x_{02}+\dots+\theta_{n-1}x_{0n-1}+\epsilon_0\\ +y_1&=\theta_0x_{10}+\theta_1x_{11}+\theta_2x_{12}+\dots+\theta_{n-1}x_{1n-1}+\epsilon_1\\ +y_2&=\theta_0x_{20}+\theta_1x_{21}+\theta_2x_{22}+\dots+\theta_{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\\ +y_{i}&=\theta_0x_{i0}+\theta_1x_{i1}+\theta_2x_{i2}+\dots+\theta_{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}.\\ +y_{n-1}&=\theta_0x_{n-1,0}+\theta_1x_{n-1,2}+\theta_2x_{n-1,2}+\dots+\theta_{n-1}x_{n-1,n-1}+\epsilon_{n-1}.\\ \end{align*} $$ @@ -411,7 +399,7 @@ our matrix as \( \boldsymbol{X}\in {\mathbb{R}}^{n\times p} \), with the predictWith \( \boldsymbol{\beta}\in {\mathbb{R}}^{p\times 1} \), it means that we will hereafter write our equations for the approximation as
+With \( \boldsymbol{\theta}\in {\mathbb{R}}^{p\times 1} \), it means that we will hereafter write our equations for the approximation as
$$ -\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\theta}, $$throughout these lectures.
@@ -489,7 +477,7 @@ $$With the above we use the design matrix to define the approximation \( \boldsymbol{\tilde{y}} \) via the unknown quantity \( \boldsymbol{\beta} \) as
+With the above we use the design matrix to define the approximation \( \boldsymbol{\tilde{y}} \) via the unknown quantity \( \boldsymbol{\theta} \) as
$$ -\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\theta}, $$ -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
+and in order to find the optimal parameters \( \theta_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
$$ -C(\boldsymbol{\beta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\}, +C(\boldsymbol{\theta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\}, $$or using the matrix \( \boldsymbol{X} \) and in a more compact matrix-vector notation as
$$ -C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}. +C(\boldsymbol{\theta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}. $$This function is one possible way to define the so-called cost function.
@@ -387,10 +375,10 @@ the function \( C \) as $$ -C(\boldsymbol{\beta})=\frac{1}{2n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2, +C(\boldsymbol{\theta})=\frac{1}{2n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2, $$ -since when taking the first derivative with respect to the unknown parameters \( \beta \), the factor of \( 2 \) cancels out.
+since when taking the first derivative with respect to the unknown parameters \( \theta \), the factor of \( 2 \) cancels out.
The function
$$ -C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}, +C(\boldsymbol{\theta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}, $$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
$$ -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, +y_{i}=\langle y_i \rangle = \theta_0x_{i,0}+\theta_1x_{i,1}+\theta_2x_{i,2}+\dots+\theta_{n-1}x_{i,n-1}+\epsilon_i, $$where \( \langle y_i \rangle \) is the mean value. Keep in mind also that @@ -388,25 +376,25 @@ 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(\boldsymbol{\beta}) \), that is we are going to solve the problem
+In order to find the parameters \( \theta_i \) we will then minimize the spread of \( C(\boldsymbol{\theta}) \), that is we are going to solve the problem
$$ -{\displaystyle \min_{\boldsymbol{\beta}\in -{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}. +{\displaystyle \min_{\boldsymbol{\theta}\in +{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}. $$In practical terms it means we will require
$$ -\frac{\partial C(\boldsymbol{\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, +\frac{\partial C(\boldsymbol{\theta})}{\partial \theta_j} = \frac{\partial }{\partial \theta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}\right)^2\right]=0, $$which results in
$$ -\frac{\partial C(\boldsymbol{\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, +\frac{\partial C(\boldsymbol{\theta})}{\partial \theta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_{ij}\left(y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}\right)\right]=0, $$or in a matrix-vector form as
$$ -\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right). +\frac{\partial C(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right). $$ @@ -437,7 +425,7 @@ $$We can rewrite
$$ -\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right), +\frac{\partial C(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right), $$as
$$ -\boldsymbol{X}^T\boldsymbol{y} = \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{X}^T\boldsymbol{y} = \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}, $$and if the matrix \( \boldsymbol{X}^T\boldsymbol{X} \) is invertible we have the solution
$$ -\boldsymbol{\beta} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}. +\boldsymbol{\theta} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}. $$We note also that since our design matrix is defined as \( \boldsymbol{X}\in @@ -427,8 +415,6 @@ allow for the usage of direct linear algebra methods such as LU decomposi
The residuals \( \boldsymbol{\epsilon} \) are in turn given by
$$ -\boldsymbol{\epsilon} = \boldsymbol{y}-\boldsymbol{\tilde{y}} = \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{\epsilon} = \boldsymbol{y}-\boldsymbol{\tilde{y}} = \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}, $$and with
$$ -\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)= 0, +\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)= 0, $$we have
$$ -\boldsymbol{X}^T\boldsymbol{\epsilon}=\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)= 0, +\boldsymbol{X}^T\boldsymbol{\epsilon}=\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)= 0, $$ -meaning that the solution for \( \boldsymbol{\beta} \) is the one which minimizes the residuals. Later we will link this with the maximum likelihood approach.
+meaning that the solution for \( \boldsymbol{\theta} \) is the one which minimizes the residuals. Later we will link this with the maximum likelihood approach.
@@ -410,9 +398,6 @@ $$It is rather straightforward to implement the matrix inversion and obtain the parameters \( \boldsymbol{\beta} \). After having defined the matrix \( \boldsymbol{X} \) we simply need to +
It is rather straightforward to implement the matrix inversion and obtain the parameters \( \boldsymbol{\theta} \). After having defined the matrix \( \boldsymbol{X} \) we simply need to write
@@ -475,10 +463,6 @@ plt.show()where the matrix \( \boldsymbol{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as matrix elements.
@@ -409,9 +397,6 @@ $$In order to find the parameters \( \beta_i \) we will then minimize the spread of \( \chi^2(\boldsymbol{\beta}) \) by requiring
+In order to find the parameters \( \theta_i \) we will then minimize the spread of \( \chi^2(\boldsymbol{\theta}) \) by requiring
$$ -\frac{\partial \chi^2(\boldsymbol{\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, +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_j} = \frac{\partial }{\partial \theta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(\frac{y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0, $$which results in
$$ -\frac{\partial \chi^2(\boldsymbol{\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, +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0, $$or in a matrix-vector form as
$$ -\frac{\partial \chi^2(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\beta}\right). +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\theta}\right). $$where we have defined the matrix \( \boldsymbol{A} =\boldsymbol{X}/\boldsymbol{\Sigma} \) with matrix elements \( a_{ij} = x_{ij}/\sigma_i \) and the vector \( \boldsymbol{b} \) with elements \( b_i = y_i/\sigma_i \).
@@ -405,9 +393,6 @@ $$We can rewrite
$$ -\frac{\partial \chi^2(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\beta}\right), +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\theta}\right), $$as
$$ -\boldsymbol{A}^T\boldsymbol{b} = \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\beta}, +\boldsymbol{A}^T\boldsymbol{b} = \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\theta}, $$and if the matrix \( \boldsymbol{A}^T\boldsymbol{A} \) is invertible we have the solution
$$ -\boldsymbol{\beta} =\left(\boldsymbol{A}^T\boldsymbol{A}\right)^{-1}\boldsymbol{A}^T\boldsymbol{b}. +\boldsymbol{\theta} =\left(\boldsymbol{A}^T\boldsymbol{A}\right)^{-1}\boldsymbol{A}^T\boldsymbol{b}. $$we have then the following expression for the parameters \( \beta_j \) (the matrix elements of \( \boldsymbol{H} \) are \( h_{ij} \))
+we have then the following expression for the parameters \( \theta_j \) (the matrix elements of \( \boldsymbol{H} \) are \( h_{ij} \))
$$ -\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} +\theta_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} $$ -We state without proof the expression for the uncertainty in the parameters \( \beta_j \) as (we leave this as an exercise)
+We state without proof the expression for the uncertainty in the parameters \( \theta_j \) as (we leave this as an exercise)
$$ -\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2, +\sigma^2(\theta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \theta_j}{\partial y_i}\right)^2, $$resulting in
$$ -\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}! +\sigma^2(\theta_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}! $$ @@ -406,9 +394,6 @@ $$The first step here is to approximate the function \( y \) with a first-order polynomial, that is we write
$$ -y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i. +y=y(x) \rightarrow y(x_i) \approx \theta_0+\theta_1 x_i. $$ -By computing the derivatives of \( \chi^2 \) with respect to \( \beta_0 \) and \( \beta_1 \) show that these are given by
+By computing the derivatives of \( \chi^2 \) with respect to \( \theta_0 \) and \( \theta_1 \) show that these are given by
$$ -\frac{\partial \chi^2(\boldsymbol{\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, +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_0} = -2\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(\frac{y_i-\theta_0-\theta_1x_{i}}{\sigma_i^2}\right)\right]=0, $$and
$$ -\frac{\partial \chi^2(\boldsymbol{\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. +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_1} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_i\left(\frac{y_i-\theta_0-\theta_1x_{i}}{\sigma_i^2}\right)\right]=0. $$ @@ -399,9 +387,6 @@ $$from pylab import plt, mpl
-plt.style.use('seaborn')
-mpl.rcParams['font.family'] = 'serif'
-
-def MakePlot(x,y, styles, labels, axlabels):
- plt.figure(figsize=(10,6))
- for i in range(len(x)):
- plt.plot(x[i], y[i], styles[i], label = labels[i])
- plt.xlabel(axlabels[0])
- plt.ylabel(axlabels[1])
- plt.legend(loc=0)
-
-Our next step is to read the data on experimental binding energies and
reorganize them as functions of the mass number \( A \), the number of
protons \( Z \) and neutrons \( N \) using pandas. Before we do this it is
@@ -2564,11 +2530,11 @@ Now it is time to dive more into the details of various methods. We will start w
Fitting a continuous function with linear parameterization in terms of the parameters \( \boldsymbol{\beta} \). Fitting a continuous function with linear parameterization in terms of the parameters \( \boldsymbol{\theta} \).
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 \( \boldsymbol{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 \( \boldsymbol{x} \) is called the independent variable, or the predictor variable or the explanatory variable, or simply just the inputs.
+The first variable \( y \) is called the the outcome or the response variable, or simply just the outputs.
The set of variables \( \boldsymbol{x} \) is called the independent variable, or the predictor variable or the explanatory variable, or simply just the inputs. We will throughout the course just use inputs and outputs as names. A regression model aims at finding a likelihood function \( p(\boldsymbol{y}\vert \boldsymbol{x}) \) or in the more traditional sense a function \( \boldsymbol{y}(\boldsymbol{x}) \), that is the conditional distribution for \( \boldsymbol{y} \) with a given \( \boldsymbol{x} \). The estimation of \( p(\boldsymbol{y}\vert \boldsymbol{x}) \) is made using a data set with
The goal of the regression analysis is to extract/exploit relationship between \( \boldsymbol{y} \) and \( \boldsymbol{x} \) in order to infer specific dependencies, approximations to the likelihood functions, functional relationships and to make predictions, making fits and many other things.Why Linear Regression (aka Ordinary Least Squares and family)
-
Linear regression gives us a set of analytical equations for the parameters \( \beta_j \).
+Linear regression gives us a set of analytical equations for the parameters \( \theta_j \).
@@ -2670,7 +2638,7 @@ so-called-
With the above we use the design matrix to define the approximation \( \boldsymbol{\tilde{y}} \) via the unknown quantity \( \boldsymbol{\beta} \) as
+With the above we use the design matrix to define the approximation \( \boldsymbol{\tilde{y}} \) via the unknown quantity \( \boldsymbol{\theta} \) as
$$
-\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\beta},
+\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\theta},
$$
-
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
+and in order to find the optimal parameters \( \theta_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
$$
-C(\boldsymbol{\beta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\},
+C(\boldsymbol{\theta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\},
$$
or using the matrix \( \boldsymbol{X} \) and in a more compact matrix-vector notation as
$$
-C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}.
+C(\boldsymbol{\theta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}.
$$
@@ -2979,11 +2947,11 @@ the function \( C \) as
$$
-C(\boldsymbol{\beta})=\frac{1}{2n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2,
+C(\boldsymbol{\theta})=\frac{1}{2n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2,
$$
-
since when taking the first derivative with respect to the unknown parameters \( \beta \), the factor of \( 2 \) cancels out.
+since when taking the first derivative with respect to the unknown parameters \( \theta \), the factor of \( 2 \) cancels out.
@@ -2996,7 +2964,7 @@ $$The function
$$
-C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\},
+C(\boldsymbol{\theta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\},
$$
@@ -3005,7 +2973,7 @@ When linking (see the discussion below) with the maximum likelihood approach bel
$$
-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,
+y_{i}=\langle y_i \rangle = \theta_0x_{i,0}+\theta_1x_{i,1}+\theta_2x_{i,2}+\dots+\theta_{n-1}x_{i,n-1}+\epsilon_i,
$$
@@ -3019,32 +2987,32 @@ 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(\boldsymbol{\beta}) \), that is we are going to solve the problem
+In order to find the parameters \( \theta_i \) we will then minimize the spread of \( C(\boldsymbol{\theta}) \), that is we are going to solve the problem
$$
-{\displaystyle \min_{\boldsymbol{\beta}\in
-{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}.
+{\displaystyle \min_{\boldsymbol{\theta}\in
+{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}.
$$
In practical terms it means we will require
$$
-\frac{\partial C(\boldsymbol{\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,
+\frac{\partial C(\boldsymbol{\theta})}{\partial \theta_j} = \frac{\partial }{\partial \theta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}\right)^2\right]=0,
$$
which results in
$$
-\frac{\partial C(\boldsymbol{\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,
+\frac{\partial C(\boldsymbol{\theta})}{\partial \theta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_{ij}\left(y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}\right)\right]=0,
$$
or in a matrix-vector form as
$$
-\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right).
+\frac{\partial C(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right).
$$
@@ -3058,21 +3026,21 @@ $$
We can rewrite
$$
-\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right),
+\frac{\partial C(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right),
$$
as
$$
-\boldsymbol{X}^T\boldsymbol{y} = \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta},
+\boldsymbol{X}^T\boldsymbol{y} = \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta},
$$
and if the matrix \( \boldsymbol{X}^T\boldsymbol{X} \) is invertible we have the solution
$$
-\boldsymbol{\beta} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
+\boldsymbol{\theta} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
$$
@@ -3105,25 +3073,25 @@ allow for the usage of direct linear algebra methods such as LU decomposi
The residuals \( \boldsymbol{\epsilon} \) are in turn given by
$$
-\boldsymbol{\epsilon} = \boldsymbol{y}-\boldsymbol{\tilde{y}} = \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta},
+\boldsymbol{\epsilon} = \boldsymbol{y}-\boldsymbol{\tilde{y}} = \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta},
$$
and with
$$
-\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)= 0,
+\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)= 0,
$$
we have
$$
-\boldsymbol{X}^T\boldsymbol{\epsilon}=\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)= 0,
+\boldsymbol{X}^T\boldsymbol{\epsilon}=\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)= 0,
$$
-
meaning that the solution for \( \boldsymbol{\beta} \) is the one which minimizes the residuals. Later we will link this with the maximum likelihood approach.
+meaning that the solution for \( \boldsymbol{\theta} \) is the one which minimizes the residuals. Later we will link this with the maximum likelihood approach.
Let us now return to our nuclear binding energies and simply code the above equations.
@@ -3132,7 +3100,7 @@ $$It is rather straightforward to implement the matrix inversion and obtain the parameters \( \boldsymbol{\beta} \). After having defined the matrix \( \boldsymbol{X} \) we simply need to +
It is rather straightforward to implement the matrix inversion and obtain the parameters \( \boldsymbol{\theta} \). After having defined the matrix \( \boldsymbol{X} \) we simply need to write
@@ -3353,7 +3321,7 @@ as
$$
-\chi^2(\boldsymbol{\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(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\frac{1}{\boldsymbol{\Sigma^2}}\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\},
+\chi^2(\boldsymbol{\theta})=\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(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\frac{1}{\boldsymbol{\Sigma^2}}\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\},
$$
@@ -3367,24 +3335,24 @@ $$
-
In order to find the parameters \( \beta_i \) we will then minimize the spread of \( \chi^2(\boldsymbol{\beta}) \) by requiring
+In order to find the parameters \( \theta_i \) we will then minimize the spread of \( \chi^2(\boldsymbol{\theta}) \) by requiring
$$
-\frac{\partial \chi^2(\boldsymbol{\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,
+\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_j} = \frac{\partial }{\partial \theta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(\frac{y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0,
$$
which results in
$$
-\frac{\partial \chi^2(\boldsymbol{\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,
+\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0,
$$
or in a matrix-vector form as
$$
-\frac{\partial \chi^2(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\beta}\right).
+\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\theta}\right).
$$
@@ -3401,21 +3369,21 @@ $$
We can rewrite
$$
-\frac{\partial \chi^2(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\beta}\right),
+\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\theta}\right),
$$
as
$$
-\boldsymbol{A}^T\boldsymbol{b} = \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\beta},
+\boldsymbol{A}^T\boldsymbol{b} = \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\theta},
$$
and if the matrix \( \boldsymbol{A}^T\boldsymbol{A} \) is invertible we have the solution
$$
-\boldsymbol{\beta} =\left(\boldsymbol{A}^T\boldsymbol{A}\right)^{-1}\boldsymbol{A}^T\boldsymbol{b}.
+\boldsymbol{\theta} =\left(\boldsymbol{A}^T\boldsymbol{A}\right)^{-1}\boldsymbol{A}^T\boldsymbol{b}.
$$
@@ -3434,24 +3402,24 @@ $$
$$
-
we have then the following expression for the parameters \( \beta_j \) (the matrix elements of \( \boldsymbol{H} \) are \( h_{ij} \))
+we have then the following expression for the parameters \( \theta_j \) (the matrix elements of \( \boldsymbol{H} \) are \( h_{ij} \))
$$
-\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}
+\theta_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}
$$
-
We state without proof the expression for the uncertainty in the parameters \( \beta_j \) as (we leave this as an exercise)
+We state without proof the expression for the uncertainty in the parameters \( \theta_j \) as (we leave this as an exercise)
$$
-\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2,
+\sigma^2(\theta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \theta_j}{\partial y_i}\right)^2,
$$
resulting in
$$
-\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}!
+\sigma^2(\theta_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}!
$$
@@ -3465,21 +3433,21 @@ $$
The first step here is to approximate the function \( y \) with a first-order polynomial, that is we write
$$
-y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i.
+y=y(x) \rightarrow y(x_i) \approx \theta_0+\theta_1 x_i.
$$
-
By computing the derivatives of \( \chi^2 \) with respect to \( \beta_0 \) and \( \beta_1 \) show that these are given by
+By computing the derivatives of \( \chi^2 \) with respect to \( \theta_0 \) and \( \theta_1 \) show that these are given by
$$
-\frac{\partial \chi^2(\boldsymbol{\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,
+\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_0} = -2\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(\frac{y_i-\theta_0-\theta_1x_{i}}{\sigma_i^2}\right)\right]=0,
$$
and
$$
-\frac{\partial \chi^2(\boldsymbol{\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.
+\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_1} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_i\left(\frac{y_i-\theta_0-\theta_1x_{i}}{\sigma_i^2}\right)\right]=0.
$$
@@ -3528,270 +3496,24 @@ $$
$$
-\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
+\theta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
$$
$$
-\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}.
+\theta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}.
$$
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 +unknown coefficients \( \theta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed next week.
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. 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. -
-# 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()
-
-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. -
-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. -
- - - -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))
-
-from pylab import plt, mpl
-plt.style.use('seaborn')
-mpl.rcParams['font.family'] = 'serif'
-
-def MakePlot(x,y, styles, labels, axlabels):
- plt.figure(figsize=(10,6))
- for i in range(len(x)):
- plt.plot(x[i], y[i], styles[i], label = labels[i])
- plt.xlabel(axlabels[0])
- plt.ylabel(axlabels[1])
- plt.legend(loc=0)
-
-Our next step is to read the data on experimental binding energies and
reorganize them as functions of the mass number \( A \), the number of
protons \( Z \) and neutrons \( N \) using pandas. Before we do this it is
@@ -2544,11 +2501,11 @@ Now it is time to dive more into the details of various methods. We will start w
Fitting a continuous function with linear parameterization in terms of the parameters \( \boldsymbol{\beta} \).
+Fitting a continuous function with linear parameterization in terms of the parameters \( \boldsymbol{\theta} \).
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 \( \boldsymbol{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 \( \boldsymbol{x} \) is called the independent variable, or the predictor variable or the explanatory variable, or simply just the inputs. +The first variable \( y \) is called the the outcome or the response variable, or simply just the outputs.
+The set of variables \( \boldsymbol{x} \) is called the independent variable, or the predictor variable or the explanatory variable, or simply just the inputs. We will throughout the course just use inputs and outputs as names.
+A regression model aims at finding a likelihood function \( p(\boldsymbol{y}\vert \boldsymbol{x}) \) or in the more traditional sense a function \( \boldsymbol{y}(\boldsymbol{x}) \), that is the conditional distribution for \( \boldsymbol{y} \) with a given \( \boldsymbol{x} \). The estimation of \( p(\boldsymbol{y}\vert \boldsymbol{x}) \) is made using a data set with
The goal of the regression analysis is to extract/exploit relationship between \( \boldsymbol{y} \) and \( \boldsymbol{x} \) in order to infer specific dependencies, approximations to the likelihood functions, functional relationships and to make predictions, making fits and many other things.
@@ -2600,11 +2559,11 @@ regression analysis is to explain \( \boldsymbol{y} \) in terms of 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 \( \boldsymbol{X} \) and \( \boldsymbol{y} \). This assumption gives rise to -the linear regression model where \( \boldsymbol{\beta} = [\beta_0, \ldots, -\beta_{p-1}]^{T} \) are the regression parameters. +the linear regression model where \( \boldsymbol{\theta} = [\theta_0, \ldots, +\theta_{p-1}]^{T} \) are the regression parameters. -Linear regression gives us a set of analytical equations for the parameters \( \beta_j \).
+Linear regression gives us a set of analytical equations for the parameters \( \theta_j \).
@@ -2644,7 +2603,7 @@ so-called Vandermonde matrix. @@ -2725,13 +2684,13 @@ of values \( y_i,x_i \) we can then generalize the equations to $$ \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\\ +y_0&=\theta_0x_{00}+\theta_1x_{01}+\theta_2x_{02}+\dots+\theta_{n-1}x_{0n-1}+\epsilon_0\\ +y_1&=\theta_0x_{10}+\theta_1x_{11}+\theta_2x_{12}+\dots+\theta_{n-1}x_{1n-1}+\epsilon_1\\ +y_2&=\theta_0x_{20}+\theta_1x_{21}+\theta_2x_{22}+\dots+\theta_{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\\ +y_{i}&=\theta_0x_{i0}+\theta_1x_{i1}+\theta_2x_{i2}+\dots+\theta_{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}.\\ +y_{n-1}&=\theta_0x_{n-1,0}+\theta_1x_{n-1,2}+\theta_2x_{n-1,2}+\dots+\theta_{n-1}x_{n-1,n-1}+\epsilon_{n-1}.\\ \end{align*} $$ @@ -2759,10 +2718,10 @@ $$and without loss of generality we rewrite again our equations as
$$ -\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}. +\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\theta}+\boldsymbol{\epsilon}. $$ -The left-hand side of this equation is kwown. Our error vector \( \boldsymbol{\epsilon} \) and the parameter vector \( \boldsymbol{\beta} \) are our unknow quantities. How can we obtain the optimal set of \( \beta_i \) values?
+The left-hand side of this equation is kwown. Our error vector \( \boldsymbol{\epsilon} \) and the parameter vector \( \boldsymbol{\theta} \) are our unknow quantities. How can we obtain the optimal set of \( \theta_i \) values?
@@ -2774,13 +2733,13 @@ $$We have defined the matrix \( \boldsymbol{X} \) via the equations
$$ \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\\ +y_0&=\theta_0x_{00}+\theta_1x_{01}+\theta_2x_{02}+\dots+\theta_{n-1}x_{0n-1}+\epsilon_0\\ +y_1&=\theta_0x_{10}+\theta_1x_{11}+\theta_2x_{12}+\dots+\theta_{n-1}x_{1n-1}+\epsilon_1\\ +y_2&=\theta_0x_{20}+\theta_1x_{21}+\theta_2x_{22}+\dots+\theta_{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\\ +y_{i}&=\theta_0x_{i0}+\theta_1x_{i1}+\theta_2x_{i2}+\dots+\theta_{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}.\\ +y_{n-1}&=\theta_0x_{n-1,0}+\theta_1x_{n-1,2}+\theta_2x_{n-1,2}+\dots+\theta_{n-1}x_{n-1,n-1}+\epsilon_{n-1}.\\ \end{align*} $$ @@ -2888,9 +2847,9 @@ display(DesignMatrix) -With \( \boldsymbol{\beta}\in {\mathbb{R}}^{p\times 1} \), it means that we will hereafter write our equations for the approximation as
+With \( \boldsymbol{\theta}\in {\mathbb{R}}^{p\times 1} \), it means that we will hereafter write our equations for the approximation as
$$ -\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\theta}, $$throughout these lectures.
@@ -2900,19 +2859,19 @@ $$-
With the above we use the design matrix to define the approximation \( \boldsymbol{\tilde{y}} \) via the unknown quantity \( \boldsymbol{\beta} \) as
+With the above we use the design matrix to define the approximation \( \boldsymbol{\tilde{y}} \) via the unknown quantity \( \boldsymbol{\theta} \) as
$$ -\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\theta}, $$ -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
+and in order to find the optimal parameters \( \theta_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
$$ -C(\boldsymbol{\beta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\}, +C(\boldsymbol{\theta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\}, $$or using the matrix \( \boldsymbol{X} \) and in a more compact matrix-vector notation as
$$ -C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}. +C(\boldsymbol{\theta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}. $$This function is one possible way to define the so-called cost function.
@@ -2922,10 +2881,10 @@ the function \( C \) as $$ -C(\boldsymbol{\beta})=\frac{1}{2n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2, +C(\boldsymbol{\theta})=\frac{1}{2n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2, $$ -since when taking the first derivative with respect to the unknown parameters \( \beta \), the factor of \( 2 \) cancels out.
+since when taking the first derivative with respect to the unknown parameters \( \theta \), the factor of \( 2 \) cancels out.
The function
$$ -C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}, +C(\boldsymbol{\theta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}, $$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
$$ -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, +y_{i}=\langle y_i \rangle = \theta_0x_{i,0}+\theta_1x_{i,1}+\theta_2x_{i,2}+\dots+\theta_{n-1}x_{i,n-1}+\epsilon_i, $$where \( \langle y_i \rangle \) is the mean value. Keep in mind also that @@ -2957,25 +2916,25 @@ 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(\boldsymbol{\beta}) \), that is we are going to solve the problem
+In order to find the parameters \( \theta_i \) we will then minimize the spread of \( C(\boldsymbol{\theta}) \), that is we are going to solve the problem
$$ -{\displaystyle \min_{\boldsymbol{\beta}\in -{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}. +{\displaystyle \min_{\boldsymbol{\theta}\in +{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}. $$In practical terms it means we will require
$$ -\frac{\partial C(\boldsymbol{\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, +\frac{\partial C(\boldsymbol{\theta})}{\partial \theta_j} = \frac{\partial }{\partial \theta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}\right)^2\right]=0, $$which results in
$$ -\frac{\partial C(\boldsymbol{\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, +\frac{\partial C(\boldsymbol{\theta})}{\partial \theta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_{ij}\left(y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}\right)\right]=0, $$or in a matrix-vector form as
$$ -\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right). +\frac{\partial C(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right). $$ @@ -2987,17 +2946,17 @@ $$
We can rewrite
$$ -\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right), +\frac{\partial C(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right), $$as
$$ -\boldsymbol{X}^T\boldsymbol{y} = \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{X}^T\boldsymbol{y} = \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}, $$and if the matrix \( \boldsymbol{X}^T\boldsymbol{X} \) is invertible we have the solution
$$ -\boldsymbol{\beta} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}. +\boldsymbol{\theta} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}. $$We note also that since our design matrix is defined as \( \boldsymbol{X}\in @@ -3028,20 +2987,20 @@ allow for the usage of direct linear algebra methods such as LU decomposi
The residuals \( \boldsymbol{\epsilon} \) are in turn given by
$$ -\boldsymbol{\epsilon} = \boldsymbol{y}-\boldsymbol{\tilde{y}} = \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{\epsilon} = \boldsymbol{y}-\boldsymbol{\tilde{y}} = \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}, $$and with
$$ -\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)= 0, +\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)= 0, $$we have
$$ -\boldsymbol{X}^T\boldsymbol{\epsilon}=\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)= 0, +\boldsymbol{X}^T\boldsymbol{\epsilon}=\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)= 0, $$ -meaning that the solution for \( \boldsymbol{\beta} \) is the one which minimizes the residuals. Later we will link this with the maximum likelihood approach.
+meaning that the solution for \( \boldsymbol{\theta} \) is the one which minimizes the residuals. Later we will link this with the maximum likelihood approach.
@@ -3050,7 +3009,7 @@ $$It is rather straightforward to implement the matrix inversion and obtain the parameters \( \boldsymbol{\beta} \). After having defined the matrix \( \boldsymbol{X} \) we simply need to +
It is rather straightforward to implement the matrix inversion and obtain the parameters \( \boldsymbol{\theta} \). After having defined the matrix \( \boldsymbol{X} \) we simply need to write
@@ -3270,7 +3229,7 @@ as $$ -\chi^2(\boldsymbol{\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(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\frac{1}{\boldsymbol{\Sigma^2}}\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\}, +\chi^2(\boldsymbol{\theta})=\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(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\frac{1}{\boldsymbol{\Sigma^2}}\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\}, $$where the matrix \( \boldsymbol{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as matrix elements.
@@ -3283,19 +3242,19 @@ $$-
In order to find the parameters \( \beta_i \) we will then minimize the spread of \( \chi^2(\boldsymbol{\beta}) \) by requiring
+In order to find the parameters \( \theta_i \) we will then minimize the spread of \( \chi^2(\boldsymbol{\theta}) \) by requiring
$$ -\frac{\partial \chi^2(\boldsymbol{\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, +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_j} = \frac{\partial }{\partial \theta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(\frac{y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0, $$which results in
$$ -\frac{\partial \chi^2(\boldsymbol{\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, +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0, $$or in a matrix-vector form as
$$ -\frac{\partial \chi^2(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\beta}\right). +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\theta}\right). $$where we have defined the matrix \( \boldsymbol{A} =\boldsymbol{X}/\boldsymbol{\Sigma} \) with matrix elements \( a_{ij} = x_{ij}/\sigma_i \) and the vector \( \boldsymbol{b} \) with elements \( b_i = y_i/\sigma_i \).
@@ -3310,17 +3269,17 @@ $$We can rewrite
$$ -\frac{\partial \chi^2(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\beta}\right), +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\theta}\right), $$as
$$ -\boldsymbol{A}^T\boldsymbol{b} = \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\beta}, +\boldsymbol{A}^T\boldsymbol{b} = \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\theta}, $$and if the matrix \( \boldsymbol{A}^T\boldsymbol{A} \) is invertible we have the solution
$$ -\boldsymbol{\beta} =\left(\boldsymbol{A}^T\boldsymbol{A}\right)^{-1}\boldsymbol{A}^T\boldsymbol{b}. +\boldsymbol{\theta} =\left(\boldsymbol{A}^T\boldsymbol{A}\right)^{-1}\boldsymbol{A}^T\boldsymbol{b}. $$ @@ -3336,19 +3295,19 @@ $$ \boldsymbol{H} = \left(\boldsymbol{A}^T\boldsymbol{A}\right)^{-1}, $$ -we have then the following expression for the parameters \( \beta_j \) (the matrix elements of \( \boldsymbol{H} \) are \( h_{ij} \))
+we have then the following expression for the parameters \( \theta_j \) (the matrix elements of \( \boldsymbol{H} \) are \( h_{ij} \))
$$ -\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} +\theta_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} $$ -We state without proof the expression for the uncertainty in the parameters \( \beta_j \) as (we leave this as an exercise)
+We state without proof the expression for the uncertainty in the parameters \( \theta_j \) as (we leave this as an exercise)
$$ -\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2, +\sigma^2(\theta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \theta_j}{\partial y_i}\right)^2, $$resulting in
$$ -\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}! +\sigma^2(\theta_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}! $$ @@ -3360,17 +3319,17 @@ $$
The first step here is to approximate the function \( y \) with a first-order polynomial, that is we write
$$ -y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i. +y=y(x) \rightarrow y(x_i) \approx \theta_0+\theta_1 x_i. $$ -By computing the derivatives of \( \chi^2 \) with respect to \( \beta_0 \) and \( \beta_1 \) show that these are given by
+By computing the derivatives of \( \chi^2 \) with respect to \( \theta_0 \) and \( \theta_1 \) show that these are given by
$$ -\frac{\partial \chi^2(\boldsymbol{\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, +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_0} = -2\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(\frac{y_i-\theta_0-\theta_1x_{i}}{\sigma_i^2}\right)\right]=0, $$and
$$ -\frac{\partial \chi^2(\boldsymbol{\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. +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_1} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_i\left(\frac{y_i-\theta_0-\theta_1x_{i}}{\sigma_i^2}\right)\right]=0. $$ @@ -3411,266 +3370,22 @@ $$we obtain
$$ -\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}, +\theta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}, $$ $$ -\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}. +\theta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}. $$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 +unknown coefficients \( \theta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed next week.
-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. 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. -
- -# 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()
-
-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. -
- -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. -
- - - -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))
-
-Before we proceed, we define also a function for making our plots. You can obviously avoid this and simply set up various matplotlib commands every time you need them. You may however find it convenient to collect all such commands in one function and simply call this function.
- - -from pylab import plt, mpl
-plt.style.use('seaborn')
-mpl.rcParams['font.family'] = 'serif'
-
-def MakePlot(x,y, styles, labels, axlabels):
- plt.figure(figsize=(10,6))
- for i in range(len(x)):
- plt.plot(x[i], y[i], styles[i], label = labels[i])
- plt.xlabel(axlabels[0])
- plt.ylabel(axlabels[1])
- plt.legend(loc=0)
-
-Our next step is to read the data on experimental binding energies and
reorganize them as functions of the mass number \( A \), the number of
protons \( Z \) and neutrons \( N \) using pandas. Before we do this it is
@@ -2621,11 +2578,11 @@ Now it is time to dive more into the details of various methods. We will start w
Fitting a continuous function with linear parameterization in terms of the parameters \( \boldsymbol{\beta} \).
+Fitting a continuous function with linear parameterization in terms of the parameters \( \boldsymbol{\theta} \).
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 \( \boldsymbol{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 \( \boldsymbol{x} \) is called the independent variable, or the predictor variable or the explanatory variable, or simply just the inputs. +The first variable \( y \) is called the the outcome or the response variable, or simply just the outputs.
+The set of variables \( \boldsymbol{x} \) is called the independent variable, or the predictor variable or the explanatory variable, or simply just the inputs. We will throughout the course just use inputs and outputs as names.
+A regression model aims at finding a likelihood function \( p(\boldsymbol{y}\vert \boldsymbol{x}) \) or in the more traditional sense a function \( \boldsymbol{y}(\boldsymbol{x}) \), that is the conditional distribution for \( \boldsymbol{y} \) with a given \( \boldsymbol{x} \). The estimation of \( p(\boldsymbol{y}\vert \boldsymbol{x}) \) is made using a data set with
The goal of the regression analysis is to extract/exploit relationship between \( \boldsymbol{y} \) and \( \boldsymbol{x} \) in order to infer specific dependencies, approximations to the likelihood functions, functional relationships and to make predictions, making fits and many other things.
@@ -2677,11 +2636,11 @@ regression analysis is to explain \( \boldsymbol{y} \) in terms of 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 \( \boldsymbol{X} \) and \( \boldsymbol{y} \). This assumption gives rise to -the linear regression model where \( \boldsymbol{\beta} = [\beta_0, \ldots, -\beta_{p-1}]^{T} \) are the regression parameters. +the linear regression model where \( \boldsymbol{\theta} = [\theta_0, \ldots, +\theta_{p-1}]^{T} \) are the regression parameters. -Linear regression gives us a set of analytical equations for the parameters \( \beta_j \).
+Linear regression gives us a set of analytical equations for the parameters \( \theta_j \).
@@ -2721,7 +2680,7 @@ so-called Vandermonde matrix. @@ -2802,13 +2761,13 @@ of values \( y_i,x_i \) we can then generalize the equations to $$ \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\\ +y_0&=\theta_0x_{00}+\theta_1x_{01}+\theta_2x_{02}+\dots+\theta_{n-1}x_{0n-1}+\epsilon_0\\ +y_1&=\theta_0x_{10}+\theta_1x_{11}+\theta_2x_{12}+\dots+\theta_{n-1}x_{1n-1}+\epsilon_1\\ +y_2&=\theta_0x_{20}+\theta_1x_{21}+\theta_2x_{22}+\dots+\theta_{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\\ +y_{i}&=\theta_0x_{i0}+\theta_1x_{i1}+\theta_2x_{i2}+\dots+\theta_{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}.\\ +y_{n-1}&=\theta_0x_{n-1,0}+\theta_1x_{n-1,2}+\theta_2x_{n-1,2}+\dots+\theta_{n-1}x_{n-1,n-1}+\epsilon_{n-1}.\\ \end{align*} $$ @@ -2836,10 +2795,10 @@ $$and without loss of generality we rewrite again our equations as
$$ -\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}. +\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\theta}+\boldsymbol{\epsilon}. $$ -The left-hand side of this equation is kwown. Our error vector \( \boldsymbol{\epsilon} \) and the parameter vector \( \boldsymbol{\beta} \) are our unknow quantities. How can we obtain the optimal set of \( \beta_i \) values?
+The left-hand side of this equation is kwown. Our error vector \( \boldsymbol{\epsilon} \) and the parameter vector \( \boldsymbol{\theta} \) are our unknow quantities. How can we obtain the optimal set of \( \theta_i \) values?
@@ -2851,13 +2810,13 @@ $$We have defined the matrix \( \boldsymbol{X} \) via the equations
$$ \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\\ +y_0&=\theta_0x_{00}+\theta_1x_{01}+\theta_2x_{02}+\dots+\theta_{n-1}x_{0n-1}+\epsilon_0\\ +y_1&=\theta_0x_{10}+\theta_1x_{11}+\theta_2x_{12}+\dots+\theta_{n-1}x_{1n-1}+\epsilon_1\\ +y_2&=\theta_0x_{20}+\theta_1x_{21}+\theta_2x_{22}+\dots+\theta_{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\\ +y_{i}&=\theta_0x_{i0}+\theta_1x_{i1}+\theta_2x_{i2}+\dots+\theta_{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}.\\ +y_{n-1}&=\theta_0x_{n-1,0}+\theta_1x_{n-1,2}+\theta_2x_{n-1,2}+\dots+\theta_{n-1}x_{n-1,n-1}+\epsilon_{n-1}.\\ \end{align*} $$ @@ -2965,9 +2924,9 @@ display(DesignMatrix) -With \( \boldsymbol{\beta}\in {\mathbb{R}}^{p\times 1} \), it means that we will hereafter write our equations for the approximation as
+With \( \boldsymbol{\theta}\in {\mathbb{R}}^{p\times 1} \), it means that we will hereafter write our equations for the approximation as
$$ -\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\theta}, $$throughout these lectures.
@@ -2977,19 +2936,19 @@ $$-
With the above we use the design matrix to define the approximation \( \boldsymbol{\tilde{y}} \) via the unknown quantity \( \boldsymbol{\beta} \) as
+With the above we use the design matrix to define the approximation \( \boldsymbol{\tilde{y}} \) via the unknown quantity \( \boldsymbol{\theta} \) as
$$ -\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{\tilde{y}}= \boldsymbol{X}\boldsymbol{\theta}, $$ -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
+and in order to find the optimal parameters \( \theta_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
$$ -C(\boldsymbol{\beta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\}, +C(\boldsymbol{\theta})=\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\}, $$or using the matrix \( \boldsymbol{X} \) and in a more compact matrix-vector notation as
$$ -C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}. +C(\boldsymbol{\theta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}. $$This function is one possible way to define the so-called cost function.
@@ -2999,10 +2958,10 @@ the function \( C \) as $$ -C(\boldsymbol{\beta})=\frac{1}{2n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2, +C(\boldsymbol{\theta})=\frac{1}{2n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2, $$ -since when taking the first derivative with respect to the unknown parameters \( \beta \), the factor of \( 2 \) cancels out.
+since when taking the first derivative with respect to the unknown parameters \( \theta \), the factor of \( 2 \) cancels out.
The function
$$ -C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}, +C(\boldsymbol{\theta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}, $$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
$$ -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, +y_{i}=\langle y_i \rangle = \theta_0x_{i,0}+\theta_1x_{i,1}+\theta_2x_{i,2}+\dots+\theta_{n-1}x_{i,n-1}+\epsilon_i, $$where \( \langle y_i \rangle \) is the mean value. Keep in mind also that @@ -3034,25 +2993,25 @@ 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(\boldsymbol{\beta}) \), that is we are going to solve the problem
+In order to find the parameters \( \theta_i \) we will then minimize the spread of \( C(\boldsymbol{\theta}) \), that is we are going to solve the problem
$$ -{\displaystyle \min_{\boldsymbol{\beta}\in -{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}. +{\displaystyle \min_{\boldsymbol{\theta}\in +{\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)\right\}. $$In practical terms it means we will require
$$ -\frac{\partial C(\boldsymbol{\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, +\frac{\partial C(\boldsymbol{\theta})}{\partial \theta_j} = \frac{\partial }{\partial \theta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}\right)^2\right]=0, $$which results in
$$ -\frac{\partial C(\boldsymbol{\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, +\frac{\partial C(\boldsymbol{\theta})}{\partial \theta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_{ij}\left(y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}\right)\right]=0, $$or in a matrix-vector form as
$$ -\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right). +\frac{\partial C(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right). $$ @@ -3064,17 +3023,17 @@ $$
We can rewrite
$$ -\frac{\partial C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right), +\frac{\partial C(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right), $$as
$$ -\boldsymbol{X}^T\boldsymbol{y} = \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{X}^T\boldsymbol{y} = \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}, $$and if the matrix \( \boldsymbol{X}^T\boldsymbol{X} \) is invertible we have the solution
$$ -\boldsymbol{\beta} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}. +\boldsymbol{\theta} =\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}. $$We note also that since our design matrix is defined as \( \boldsymbol{X}\in @@ -3105,20 +3064,20 @@ allow for the usage of direct linear algebra methods such as LU decomposi
The residuals \( \boldsymbol{\epsilon} \) are in turn given by
$$ -\boldsymbol{\epsilon} = \boldsymbol{y}-\boldsymbol{\tilde{y}} = \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}, +\boldsymbol{\epsilon} = \boldsymbol{y}-\boldsymbol{\tilde{y}} = \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}, $$and with
$$ -\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)= 0, +\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)= 0, $$we have
$$ -\boldsymbol{X}^T\boldsymbol{\epsilon}=\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)= 0, +\boldsymbol{X}^T\boldsymbol{\epsilon}=\boldsymbol{X}^T\left( \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\right)= 0, $$ -meaning that the solution for \( \boldsymbol{\beta} \) is the one which minimizes the residuals. Later we will link this with the maximum likelihood approach.
+meaning that the solution for \( \boldsymbol{\theta} \) is the one which minimizes the residuals. Later we will link this with the maximum likelihood approach.
@@ -3127,7 +3086,7 @@ $$It is rather straightforward to implement the matrix inversion and obtain the parameters \( \boldsymbol{\beta} \). After having defined the matrix \( \boldsymbol{X} \) we simply need to +
It is rather straightforward to implement the matrix inversion and obtain the parameters \( \boldsymbol{\theta} \). After having defined the matrix \( \boldsymbol{X} \) we simply need to write
@@ -3347,7 +3306,7 @@ as $$ -\chi^2(\boldsymbol{\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(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\frac{1}{\boldsymbol{\Sigma^2}}\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\}, +\chi^2(\boldsymbol{\theta})=\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(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)^T\frac{1}{\boldsymbol{\Sigma^2}}\left(\boldsymbol{y}-\boldsymbol{\tilde{y}}\right)\right\}, $$where the matrix \( \boldsymbol{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as matrix elements.
@@ -3360,19 +3319,19 @@ $$-
In order to find the parameters \( \beta_i \) we will then minimize the spread of \( \chi^2(\boldsymbol{\beta}) \) by requiring
+In order to find the parameters \( \theta_i \) we will then minimize the spread of \( \chi^2(\boldsymbol{\theta}) \) by requiring
$$ -\frac{\partial \chi^2(\boldsymbol{\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, +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_j} = \frac{\partial }{\partial \theta_j}\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(\frac{y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0, $$which results in
$$ -\frac{\partial \chi^2(\boldsymbol{\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, +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_j} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\theta_0x_{i,0}-\theta_1x_{i,1}-\theta_2x_{i,2}-\dots-\theta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0, $$or in a matrix-vector form as
$$ -\frac{\partial \chi^2(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\beta}\right). +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\theta}\right). $$where we have defined the matrix \( \boldsymbol{A} =\boldsymbol{X}/\boldsymbol{\Sigma} \) with matrix elements \( a_{ij} = x_{ij}/\sigma_i \) and the vector \( \boldsymbol{b} \) with elements \( b_i = y_i/\sigma_i \).
@@ -3387,17 +3346,17 @@ $$We can rewrite
$$ -\frac{\partial \chi^2(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\beta}\right), +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 0 = \boldsymbol{A}^T\left( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\theta}\right), $$as
$$ -\boldsymbol{A}^T\boldsymbol{b} = \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\beta}, +\boldsymbol{A}^T\boldsymbol{b} = \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{\theta}, $$and if the matrix \( \boldsymbol{A}^T\boldsymbol{A} \) is invertible we have the solution
$$ -\boldsymbol{\beta} =\left(\boldsymbol{A}^T\boldsymbol{A}\right)^{-1}\boldsymbol{A}^T\boldsymbol{b}. +\boldsymbol{\theta} =\left(\boldsymbol{A}^T\boldsymbol{A}\right)^{-1}\boldsymbol{A}^T\boldsymbol{b}. $$ @@ -3413,19 +3372,19 @@ $$ \boldsymbol{H} = \left(\boldsymbol{A}^T\boldsymbol{A}\right)^{-1}, $$ -we have then the following expression for the parameters \( \beta_j \) (the matrix elements of \( \boldsymbol{H} \) are \( h_{ij} \))
+we have then the following expression for the parameters \( \theta_j \) (the matrix elements of \( \boldsymbol{H} \) are \( h_{ij} \))
$$ -\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} +\theta_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} $$ -We state without proof the expression for the uncertainty in the parameters \( \beta_j \) as (we leave this as an exercise)
+We state without proof the expression for the uncertainty in the parameters \( \theta_j \) as (we leave this as an exercise)
$$ -\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2, +\sigma^2(\theta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \theta_j}{\partial y_i}\right)^2, $$resulting in
$$ -\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}! +\sigma^2(\theta_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}! $$ @@ -3437,17 +3396,17 @@ $$
The first step here is to approximate the function \( y \) with a first-order polynomial, that is we write
$$ -y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i. +y=y(x) \rightarrow y(x_i) \approx \theta_0+\theta_1 x_i. $$ -By computing the derivatives of \( \chi^2 \) with respect to \( \beta_0 \) and \( \beta_1 \) show that these are given by
+By computing the derivatives of \( \chi^2 \) with respect to \( \theta_0 \) and \( \theta_1 \) show that these are given by
$$ -\frac{\partial \chi^2(\boldsymbol{\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, +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_0} = -2\left[ \frac{1}{n}\sum_{i=0}^{n-1}\left(\frac{y_i-\theta_0-\theta_1x_{i}}{\sigma_i^2}\right)\right]=0, $$and
$$ -\frac{\partial \chi^2(\boldsymbol{\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. +\frac{\partial \chi^2(\boldsymbol{\theta})}{\partial \theta_1} = -\frac{2}{n}\left[ \sum_{i=0}^{n-1}x_i\left(\frac{y_i-\theta_0-\theta_1x_{i}}{\sigma_i^2}\right)\right]=0. $$ @@ -3488,266 +3447,22 @@ $$we obtain
$$ -\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}, +\theta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}, $$ $$ -\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}. +\theta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}. $$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 +unknown coefficients \( \theta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed next week.
-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. 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. -
- -# 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()
-
-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. -
- -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. -
- - - -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))
-
-)cCkHG3UGvqxEI{)D?ApB9nDl!8`&P6dveoW3!fqchz+6>VmpJ0q
zUu9Odz3}&VA%(NucV6{PzAQqqkQ#wR5gOx;y@HG4u2@^Ah&MiC087kvWow@gwQJm(
znw(6I`1Q4 BjW(T?E2W4pljis`l
z>&F>0T`}GmO{hedk~y{CX_Ha!MZ4H)TPQ_ZP>yUuGa=X>8s3CImlAe{?%-&<{YAIG
zK$T8%HA7!|ZGDq2dk-#pD%Dlo$*3P=px+fax@kn=u8RuS{N?mgZdTScM!9?8dV^m6
zaT3Eb9TN{8x4}>LT*>t+`h}(qL*Cp0&e_r$$f-^uQ+Nq
zoHavU-yeI4u4IrrgN0g6DT+>=E*?%;t=PBywbm2?+n8#04w*Vyf09m;$KUVH*mUQL
zjDI4Peh$SXjKw1Bb;@(?ByUUQ$?ViB+U<6+2agNpCBt3yBrBPo225&jT>0dyZ|6+j
zTiNg`I~%(6G~+8hWpdne)*AD^{hAvk=NTS4bRwh**`_XXWmIw7K@oG1l8bF$lT)b3
zo`PX6Hn5QRmxH~~ q)UXn(a3a|aB(w`A;G)Fu{us10{5N(&Zk`&4~N|h
z&X2BXYhtlMZt;_JG13WfOD4uD$H(KICDA*^Pd21&-ZQlD=h(hqpT~Z;xq-nkhpCIk!;T!d
z$LhV)tamSm{$(&rejRW{##5hOMxYG_p3sHhQNCa~Kd-kEk&T~yiHuX`oglBFPqwrvcN_+k08wff8cMpZd%ac+tfCEjzU
zlvvr|yG+*ZZg+xVF0uaf_aylwJG-u0ztf1A@!^e;YQ09KEwQ9$bUZQh=S7R|Dl%O}
z&bQp=
h!(!H)KDysw9
zvMrtsgsoGMQ*r1R<`XGno