diff --git a/doc/pub/Regression/html/Regression-bs.html b/doc/pub/Regression/html/Regression-bs.html
index 618248fe1..57bafe39d 100644
--- a/doc/pub/Regression/html/Regression-bs.html
+++ b/doc/pub/Regression/html/Regression-bs.html
@@ -11,9 +11,9 @@ Automatically generated HTML file from DocOnce source
 <title>Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis</title>
 
 <!-- Bootstrap style: bootstrap -->
-<link href="http://netdna.bootstrapcdn.com/bootstrap/3.1.1/css/bootstrap.min.css" rel="stylesheet">
+<link href="https://netdna.bootstrapcdn.com/bootstrap/3.1.1/css/bootstrap.min.css" rel="stylesheet">
 <!-- not necessary
-<link href="http://netdna.bootstrapcdn.com/font-awesome/4.0.3/css/font-awesome.css" rel="stylesheet">
+<link href="https://netdna.bootstrapcdn.com/font-awesome/4.0.3/css/font-awesome.css" rel="stylesheet">
 -->
 
 <style type="text/css">
@@ -46,7 +46,8 @@ Automatically generated HTML file from DocOnce source
                2,
                None,
                '___sec2'),
-              ('Rewriting the fitting procedure as a linear algebra problem, follows',
+              ('Rewriting the fitting procedure as a linear algebra problem, '
+               'follows',
                2,
                None,
                '___sec3'),
@@ -88,13 +89,13 @@ end of tocinfo -->
 <script type="text/x-mathjax-config">
 MathJax.Hub.Config({
   TeX: {
-     equationNumbers: {  autoNumber: "none"  },
+     equationNumbers: {  autoNumber: "AMS"  },
      extensions: ["AMSmath.js", "AMSsymbols.js", "autobold.js", "color.js"]
   }
 });
 </script>
-<script type="text/javascript"
- src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
+<script type="text/javascript" async
+ src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
 </script>
 
 
@@ -116,24 +117,24 @@ MathJax.Hub.Config({
       <li class="dropdown">
         <a href="#" class="dropdown-toggle" data-toggle="dropdown">Contents <b class="caret"></b></a>
         <ul class="dropdown-menu">
-     <!-- navigation toc: --> <li><a href="._Regression-bs001.html#___sec0" style="font-size: 80%;">Regression analysis, overarching aims</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs002.html#___sec1" style="font-size: 80%;">General linear models</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs003.html#___sec2" style="font-size: 80%;">Rewriting the fitting procedure as a linear algebra problem</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs004.html#___sec3" style="font-size: 80%;">Rewriting the fitting procedure as a linear algebra problem, follows</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs005.html#___sec4" style="font-size: 80%;">Generalizing the fitting procedure as a linear algebra problem</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs006.html#___sec5" style="font-size: 80%;">Generalizing the fitting procedure as a linear algebra problem</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs007.html#___sec6" style="font-size: 80%;">Optimizing our parameters</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs008.html#___sec7" style="font-size: 80%;">Optimizing our parameters, more details</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs009.html#___sec8" style="font-size: 80%;">Interpretations and optimizing our parameters</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs010.html#___sec9" style="font-size: 80%;">Interpretations and optimizing our parameters</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs011.html#___sec10" style="font-size: 80%;">Interpretations and optimizing our parameters</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs012.html#___sec11" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs013.html#___sec12" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs014.html#___sec13" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs015.html#___sec14" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs016.html#___sec15" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs017.html#___sec16" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
-     <!-- navigation toc: --> <li><a href="._Regression-bs018.html#___sec17" style="font-size: 80%;">The singular value decompostion</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec0" style="font-size: 80%;">Regression analysis, overarching aims</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec1" style="font-size: 80%;">General linear models</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec2" style="font-size: 80%;">Rewriting the fitting procedure as a linear algebra problem</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec3" style="font-size: 80%;">Rewriting the fitting procedure as a linear algebra problem, follows</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec4" style="font-size: 80%;">Generalizing the fitting procedure as a linear algebra problem</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec5" style="font-size: 80%;">Generalizing the fitting procedure as a linear algebra problem</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec6" style="font-size: 80%;">Optimizing our parameters</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec7" style="font-size: 80%;">Optimizing our parameters, more details</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec8" style="font-size: 80%;">Interpretations and optimizing our parameters</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec9" style="font-size: 80%;">Interpretations and optimizing our parameters</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec10" style="font-size: 80%;">Interpretations and optimizing our parameters</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec11" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec12" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec13" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec14" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec15" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec16" style="font-size: 80%;">The \( \chi^2 \) function</a></li>
+     <!-- navigation toc: --> <li><a href="#___sec17" style="font-size: 80%;">The singular value decompostion</a></li>
 
         </ul>
       </li>
@@ -146,7 +147,6 @@ MathJax.Hub.Config({
 
 <p>&nbsp;</p><p>&nbsp;</p><p>&nbsp;</p> <!-- add vertical space -->
 
-<a name="part0000"></a>
 <!-- ------------------- main content ---------------------- -->
 
 
@@ -168,39 +168,523 @@ MathJax.Hub.Config({
 <center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
 <br>
 <p>
-<center><h4>Nov 2, 2017</h4></center> <!-- date -->
+<center><h4>Nov 21, 2017</h4></center> <!-- date -->
 <br>
 <p>
-
-
-<p><a href="._Regression-bs001.html" class="btn btn-primary btn-lg">Read &raquo;</a></p>
-
-
+<!-- potential-jumbotron-button -->
 </div> <!-- end jumbotron -->
 
+<!-- !split -->
+
+<h2 id="___sec0" class="anchor">Regression analysis, overarching aims  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+
 <p>
-<!-- navigation buttons at the bottom of the page -->
-<ul class="pagination">
-  <li class="active"><a href="._Regression-bs000.html">1</a></li>
-  <li><a href="._Regression-bs001.html">2</a></li>
-  <li><a href="._Regression-bs002.html">3</a></li>
-  <li><a href="._Regression-bs003.html">4</a></li>
-  <li><a href="._Regression-bs004.html">5</a></li>
-  <li><a href="._Regression-bs005.html">6</a></li>
-  <li><a href="._Regression-bs006.html">7</a></li>
-  <li><a href="._Regression-bs007.html">8</a></li>
-  <li><a href="._Regression-bs008.html">9</a></li>
-  <li><a href="._Regression-bs009.html">10</a></li>
-  <li><a href="">...</a></li>
-  <li><a href="._Regression-bs018.html">19</a></li>
-  <li><a href="._Regression-bs001.html">&raquo;</a></li>
+Regression modeling deals with the description of  the sampling distribution of a given random variable \( y \) varies as function of another variable or a set of such variables \( \hat{x} =[x_0, x_1,\dots, x_p]^T \). 
+The first variable is called the <b>dependent</b>, the <b>outcome</b> or the <b>response</b> variable while the set of variables \( \hat{x} \) is called the independent variable, or the predictor variable or the explanatory variable.
+
+<p>
+A regression model aims at finding a likelihood function \( p(y\vert \hat{x}) \), that is the conditional distribution for \( y \) with a given \( \hat{x} \). The estimation of  \( p(y\vert \hat{x}) \) is made using a data set with 
+
+<ul>
+<li> \( n \) cases \( i = 0, 1, 2, \dots, n-1 \)</li> 
+<li> Response (dependent or outcome) variable \( y_i \) with \( i = 0, 1, 2, \dots, n-1 \)</li> 
+<li> \( p \) Explanatory (independent or predictor) variables \( \hat{x}_i=[x_{i0}, x_{i1}, \dots, x_{ip}] \) with \( i = 0, 1, 2, \dots, n-1 \)</li>   
 </ul>
+
+ The goal of the regression analysis is to extract/exploit relationship between \( y_i \) and \( \hat{x}_i \) in or to infer causal dependencies, approximations to the likelihood functions, functional relationships and to make predictions .
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec1" class="anchor">General linear models  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+Before we proceed let us study a case from linear algebra where we aim at fitting a set of data \( \hat{y}=[y_0,y_1,\dots,y_{n-1}] \). We could think of these data as a result of an experiment or a complicated numerical experiment. These data are functions of a series of variables \( \hat{x}=[x_0,x_1,\dots,x_{n-1}] \), that is \( y_i = y(x_i) \) with \( i=0,1,2,\dots,n-1 \). The variables \( x_i \) could represent physical quantities like time, temperature, position etc. We assume that \( y(x) \) is a smooth function.
+
+<p>
+Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of \( y \) which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree \( n-1 \) with \( n \) points, that is
+$$
+y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \beta_i x_i^j+\epsilon_i,
+$$
+
+where \( \epsilon_i \) is the error in our approximation.
+
+<p>
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec2" class="anchor">Rewriting the fitting procedure as a linear algebra problem  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+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\\
+\dots & \dots \\
+y_{n-1}&=\beta_0+\beta_1x_{n-1}^1+\beta_2x_{n-1}^2+\dots+\beta_1x_{n-1}^{n-1}+\epsilon_{n-1}.\\
+\end{align*}
+$$
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec3" class="anchor">Rewriting the fitting procedure as a linear algebra problem, follows  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+Defining the vectors
+$$
+\hat{y} = [y_0,y_1, y_2,\dots, y_{n-1}]^T,
+$$
+
+$$
+\hat{\beta} = [\beta_0,\beta_1, \beta_2,\dots, \beta_{n-1}]^T,
+$$
+
+$$
+\hat{\epsilon} = [\epsilon_0,\epsilon_1, \epsilon_2,\dots, \epsilon_{n-1}]^T,
+$$
+
+and the matrix
+$$
+\hat{X}=
+\begin{bmatrix} 
+1& x_{0}^1 &x_{0}^2& \dots & \dots &x_{0}^{n-1}\\
+1& x_{1}^1 &x_{1}^2& \dots & \dots &x_{1}^{n-1}\\
+1& x_{2}^1 &x_{2}^2& \dots & \dots &x_{2}^{n-1}\\                      
+\dots& \dots &\dots& \dots & \dots &\dots\\
+1& x_{n-1}^1 &x_{n-1}^2& \dots & \dots &x_{n-1}^{n-1}\\
+\end{bmatrix} 
+$$
+
+we can rewrite our equations as
+$$
+\hat{y} = \hat{X}\hat{\beta}+\hat{\epsilon}.
+$$
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec4" class="anchor">Generalizing the fitting procedure as a linear algebra problem  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+We are obviously not limited to the above polynomial. We could replace the various powers of \( x \) with elements of Fourier series, that is, instead of \( x_i^j \) we could have \( \cos{(j x_i)} \) or \( \sin{(j x_i)} \), or time series or other orthogonal functions.
+For every set of values \( y_i,x_i \) we can then generalize the equations to 
+$$
+\begin{align*}
+y_0&=\beta_0x_{00}+\beta_1x_{01}+\beta_2x_{02}+\dots+\beta_{n-1}x_{0n-1}+\epsilon_0\\
+y_1&=\beta_0x_{10}+\beta_1x_{11}+\beta_2x_{12}+\dots+\beta_{n-1}x_{1n-1}+\epsilon_1\\
+y_2&=\beta_0x_{20}+\beta_1x_{21}+\beta_2x_{22}+\dots+\beta_{n-1}x_{2n-1}+\epsilon_2\\
+\dots & \dots \\
+y_{i}&=\beta_0x_{i0}+\beta_1x_{i1}+\beta_2x_{i2}+\dots+\beta_{n-1}x_{in-1}+\epsilon_i\\
+\dots & \dots \\
+y_{n-1}&=\beta_0x_{n-1,0}+\beta_1x_{n-1,2}+\beta_2x_{n-1,2}+\dots+\beta_1x_{n-1,n-1}+\epsilon_{n-1}.\\
+\end{align*}
+$$
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec5" class="anchor">Generalizing the fitting procedure as a linear algebra problem  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+We redefine in turn the matrix \( \hat{X} \) as
+$$
+\hat{X}=
+\begin{bmatrix} 
+x_{00}& x_{01} &x_{02}& \dots & \dots &x_{0,n-1}\\
+x_{10}& x_{11} &x_{12}& \dots & \dots &x_{1,n-1}\\
+x_{20}& x_{21} &x_{22}& \dots & \dots &x_{2,n-1}\\                      
+\dots& \dots &\dots& \dots & \dots &\dots\\
+x_{n-1,0}& x_{n-1,1} &x_{n-1,2}& \dots & \dots &x_{n-1,n-1}\\
+\end{bmatrix} 
+$$
+
+and without loss of generality we rewrite again  our equations as
+$$
+\hat{y} = \hat{X}\hat{\beta}+\hat{\epsilon}.
+$$
+
+The left-hand side of this equation forms know. Our error vector \( \hat{\epsilon} \) and the parameter vector \( \hat{\beta} \) are our unknow quantities. How can we obtain the optimal set of \( \beta_i \) values?
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec6" class="anchor">Optimizing our parameters  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+We have defined the matrix \( \hat{X} \)
+$$
+\begin{align*}
+y_0&=\beta_0x_{00}+\beta_1x_{01}+\beta_2x_{02}+\dots+\beta_{n-1}x_{0n-1}+\epsilon_0\\
+y_1&=\beta_0x_{10}+\beta_1x_{11}+\beta_2x_{12}+\dots+\beta_{n-1}x_{1n-1}+\epsilon_1\\
+y_2&=\beta_0x_{20}+\beta_1x_{21}+\beta_2x_{22}+\dots+\beta_{n-1}x_{2n-1}+\epsilon_1\\
+\dots & \dots \\
+y_{i}&=\beta_0x_{i0}+\beta_1x_{i1}+\beta_2x_{i2}+\dots+\beta_{n-1}x_{in-1}+\epsilon_1\\
+\dots & \dots \\
+y_{n-1}&=\beta_0x_{n-1,0}+\beta_1x_{n-1,2}+\beta_2x_{n-1,2}+\dots+\beta_1x_{n-1,n-1}+\epsilon_{n-1}.\\
+\end{align*}
+$$
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec7" class="anchor">Optimizing our parameters, more details  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+We well use this matrix to define the approximation \( \hat{\tilde{y}} \) via the unknown quantity \( \hat{\beta} \) as
+$$
+\hat{\tilde{y}}= \hat{X}\hat{\beta},
+$$
+
+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 parametrized values \( \tilde{y}_i \), namely
+$$
+Q(\hat{\beta})=\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\left(\hat{y}-\hat{\tilde{y}}\right)^T\left(\hat{y}-\hat{\tilde{y}}\right),
+$$
+
+or using the matrix \( \hat{X} \) as
+$$
+Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right).
+$$
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec8" class="anchor">Interpretations and optimizing our parameters  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+The function 
+$$
+Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right),
+$$
+
+can be linked to the variance of the quantity \( y_i \) if we interpret the latter as the mean value of for example a numerical  experiment. When linking 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,
+$$
+
+where \( \langle y_i \rangle \) is the mean value. Keep in mind also that till now  we have treated \( y_i \) as the exact value. Normally, the response (dependent or outcome) variable \( y_i \) the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat \( y_i \) as our exact value for the response variable.
+
+<p>
+In order to find the parameters \( \beta_i \) we will then minimize the spread of \( Q(\hat{\beta}) \) by requiring
+$$
+\frac{\partial Q(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \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, 
+$$
+
+which results in
+$$
+\frac{\partial Q(\hat{\beta})}{\partial \beta_j} = -2\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, 
+$$
+
+or in a matrix-vector form as
+$$
+\frac{\partial Q(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{X}^T\left( \hat{y}-\hat{X}\hat{\beta}\right).  
+$$
+
+<p>
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec9" class="anchor">Interpretations and optimizing our parameters  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+We can rewrite
+$$
+\frac{\partial Q(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{X}^T\left( \hat{y}-\hat{X}\hat{\beta}\right),  
+$$
+
+as
+$$
+\hat{X}^T\hat{y} = \hat{X}^T\hat{X}\hat{\beta},  
+$$
+
+and if the matrix \( \hat{X}^T\hat{X} \) is invertible we have the solution
+$$
+\hat{\beta} =\left(\hat{X}^T\hat{X}\right)^{-1}\hat{X}^T\hat{y}.
+$$
+
+<p>
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec10" class="anchor">Interpretations and optimizing our parameters  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+The residuals \( \hat{\epsilon} \) are in turn given by
+$$
+\hat{\epsilon} = \hat{y}-\hat{\tilde{y}} = \hat{y}-\hat{X}\hat{\beta},
+$$
+
+and with 
+$$
+\hat{X}^T\left( \hat{y}-\hat{X}\hat{\beta}\right)= 0, 
+$$
+
+we have
+$$
+\hat{X}^T\hat{\epsilon}=\hat{X}^T\left( \hat{y}-\hat{X}\hat{\beta}\right)= 0, 
+$$
+
+meaning that the solution for \( \hat{\beta} \) is the one which minimizes the residuals.  Later we will link this with the maximum likelihood approach.
+
+<p>
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec11" class="anchor">The \( \chi^2 \) function  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+
+<p>
+Normally, the response (dependent or outcome) variable \( y_i \) the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat \( y_i \) as our exact value for the response variable.
+
+<p>
+Introducing the standard deviation \( \sigma_i \) for each measurement \( y_i \), we define now the \( \chi^2 \) function as  
+$$
+\chi^2(\hat{\beta})=\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\left(\hat{y}-\hat{\tilde{y}}\right)^T\frac{1}{\hat{\Sigma^2}}\left(\hat{y}-\hat{\tilde{y}}\right),
+$$
+
+where the matrix \( \hat{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as matrix elements.
+
+<p>
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec12" class="anchor">The \( \chi^2 \) function  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+
+<p>
+In order to find the parameters \( \beta_i \) we will then minimize the spread of \( \chi^2(\hat{\beta}) \) by requiring
+$$
+\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \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, 
+$$
+
+which results in
+$$
+\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = -2\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, 
+$$
+
+or in a matrix-vector form as
+$$
+\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right).  
+$$
+
+where we have defined the matrix \( \hat{A} =\hat{X}/\hat{\Sigma} \) with matrix elements \( a_{ij} = x_{ij}/\sigma_i \) and the vector \( \hat{b} \) with elements \( b_i = y_i/\sigma_i \).
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec13" class="anchor">The \( \chi^2 \) function  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+
+<p>
+We can rewrite
+$$
+\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right),  
+$$
+
+as
+$$
+\hat{A}^T\hat{b} = \hat{A}^T\hat{A}\hat{\beta},  
+$$
+
+and if the matrix \( \hat{A}^T\hat{A} \) is invertible we have the solution
+$$
+\hat{\beta} =\left(\hat{A}^T\hat{A}\right)^{-1}\hat{A}^T\hat{b}.
+$$
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec14" class="anchor">The \( \chi^2 \) function  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+
+<p>
+If we then introduce the matrix
+$$
+\hat{H} =  \hat{A}^T\hat{A},
+$$
+
+we have then the following expression for the parameters \( \beta_j \) (the matrix elements of \( \hat{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}
+$$
+
+We state without proof the expression for the uncertainty  in the parameters \( \beta_j \) as
+$$
+\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_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}!
+$$
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec15" class="anchor">The \( \chi^2 \) function  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+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.
+$$
+
+By computing the derivatives of \( \chi^2 \) with respect to \( \beta_0 \) and \( \beta_1 \) show that these are given by
+$$
+\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0,
+$$
+
+and
+$$
+\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}x_i\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0.
+$$
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec16" class="anchor">The \( \chi^2 \) function  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+
+<p>
+We define then
+$$
+\gamma =  \sum_{i=0}^{1}\frac{1}{\sigma_i^2},
+$$
+
+
+$$
+\gamma_x =  \sum_{i=0}^{1}\frac{x_{i}}{\sigma_i^2},
+$$
+
+$$
+\gamma_y = \sum_{i=0}^{1}\left(\frac{y_i}{\sigma_i^2}\right),
+$$
+
+$$
+\gamma_{xx} =  \sum_{i=0}^{1}\frac{x_ix_{i}}{\sigma_i^2},
+$$
+
+$$
+\gamma_{xy} = \sum_{i=0}^{1}\frac{y_ix_{i}}{\sigma_i^2},
+$$
+
+and show that
+$$
+\beta_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}.
+$$
+
+<p>
+The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \).  A better approach is to use the Singular Value Decomposition (SVD) method discussed below.
+</div>
+</div>
+
+
+<p>
+<!-- !split -->
+
+<h2 id="___sec17" class="anchor">The singular value decompostion  </h2>
+<div class="panel panel-default">
+<div class="panel-body">
+<p> <!-- subsequent paragraphs come in larger fonts, so start with a paragraph -->
+How can we use the singular value decomposition to find the parameters \( \beta_j \)? More details will come. We first note that a general \( m\times n \) matrix \( \hat{A} \) can be written in terms of a diagonal matrix \( \hat{\Sigma} \) of dimensionality \( n\times n \) and two orthognal matrices \( \hat{U} \) and \( \hat{V} \), where the first has dimensionality \( m \times n \) and the last dimensionality \( n\times n \). We have then
+$$
+\hat{A} = \hat{U}\hat{\Sigma}\hat{V}
+$$
+</div>
+</div>
+
+
+<p>
+
 <!-- ------------------- end of main content --------------- -->
 
 </div>  <!-- end container -->
 <!-- include javascript, jQuery *first* -->
-<script src="http://ajax.googleapis.com/ajax/libs/jquery/1.10.2/jquery.min.js"></script>
-<script src="http://netdna.bootstrapcdn.com/bootstrap/3.0.0/js/bootstrap.min.js"></script>
+<script src="https://ajax.googleapis.com/ajax/libs/jquery/1.10.2/jquery.min.js"></script>
+<script src="https://netdna.bootstrapcdn.com/bootstrap/3.0.0/js/bootstrap.min.js"></script>
 
 <!-- Bootstrap footer
 <footer>
diff --git a/doc/pub/Regression/html/Regression-reveal.html b/doc/pub/Regression/html/Regression-reveal.html
index a35af3441..5653f2d53 100644
--- a/doc/pub/Regression/html/Regression-reveal.html
+++ b/doc/pub/Regression/html/Regression-reveal.html
@@ -1,3 +1,4 @@
+\
 <!DOCTYPE html>
 
 <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
@@ -116,8 +117,8 @@ MathJax.Hub.Config({
   }
 });
 </script>
-<script type="text/javascript"
- src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
+<script type="text/javascript" async
+ src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
 </script>
 
 
@@ -147,7 +148,7 @@ MathJax.Hub.Config({
 <center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
 <br>
 <p>&nbsp;<br>
-<center><h4>Nov 2, 2017</h4></center> <!-- date -->
+<center><h4>Nov 21, 2017</h4></center> <!-- date -->
 <br>
 <p>
 
diff --git a/doc/pub/Regression/html/Regression-solarized.html b/doc/pub/Regression/html/Regression-solarized.html
index 02c5cf4f8..79d8f0338 100644
--- a/doc/pub/Regression/html/Regression-solarized.html
+++ b/doc/pub/Regression/html/Regression-solarized.html
@@ -15,7 +15,7 @@ Automatically generated HTML file from DocOnce source
 <script src="https://cdn.rawgit.com/hplgit/doconce/master/bundled/html_styles/style_solarized_box/js/highlight.pack.js"></script>
 <script>hljs.initHighlightingOnLoad();</script>
 
-<link href="http://thomasf.github.io/solarized-css/solarized-light.min.css" rel="stylesheet">
+<link href="https://thomasf.github.io/solarized-css/solarized-light.min.css" rel="stylesheet">
 <style type="text/css">
 h1 {color: #b58900;}  /* yellow */
 /* h1 {color: #cb4b16;}  orange */
@@ -66,7 +66,8 @@ div { text-align: justify; text-justify: inter-word; }
                2,
                None,
                '___sec2'),
-              ('Rewriting the fitting procedure as a linear algebra problem, follows',
+              ('Rewriting the fitting procedure as a linear algebra problem, '
+               'follows',
                2,
                None,
                '___sec3'),
@@ -113,8 +114,8 @@ MathJax.Hub.Config({
   }
 });
 </script>
-<script type="text/javascript"
- src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
+<script type="text/javascript" async
+ src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
 </script>
 
 
@@ -140,7 +141,7 @@ MathJax.Hub.Config({
 <center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
 <br>
 <p>
-<center><h4>Nov 2, 2017</h4></center> <!-- date -->
+<center><h4>Nov 21, 2017</h4></center> <!-- date -->
 <br>
 <p>
 <!-- !split --><br><br><br><br><br><br><br><br><br><br>
diff --git a/doc/pub/Regression/html/Regression.html b/doc/pub/Regression/html/Regression.html
index 8a93ef294..497beb44d 100644
--- a/doc/pub/Regression/html/Regression.html
+++ b/doc/pub/Regression/html/Regression.html
@@ -71,7 +71,8 @@ div { text-align: justify; text-justify: inter-word; }
                2,
                None,
                '___sec2'),
-              ('Rewriting the fitting procedure as a linear algebra problem, follows',
+              ('Rewriting the fitting procedure as a linear algebra problem, '
+               'follows',
                2,
                None,
                '___sec3'),
@@ -118,8 +119,8 @@ MathJax.Hub.Config({
   }
 });
 </script>
-<script type="text/javascript"
- src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
+<script type="text/javascript" async
+ src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
 </script>
 
 
@@ -145,7 +146,7 @@ MathJax.Hub.Config({
 <center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
 <br>
 <p>
-<center><h4>Nov 2, 2017</h4></center> <!-- date -->
+<center><h4>Nov 21, 2017</h4></center> <!-- date -->
 <br>
 <p>
 <!-- !split --><br><br><br><br><br><br><br><br><br><br>
diff --git a/doc/pub/Regression/ipynb/Regression.ipynb b/doc/pub/Regression/ipynb/Regression.ipynb
index a415d58f8..9f7f21521 100644
--- a/doc/pub/Regression/ipynb/Regression.ipynb
+++ b/doc/pub/Regression/ipynb/Regression.ipynb
@@ -4,23 +4,918 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "<!-- dom:TITLE: Data Analysis and Machine Learning: Linear and more Advanced Regression Analysis -->\n",
-    "# Data Analysis and Machine Learning: Linear and more Advanced Regression Analysis\n",
+    "<!-- dom:TITLE: Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis -->\n",
+    "# Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis\n",
     "<!-- dom:AUTHOR: Morten Hjorth-Jensen at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University -->\n",
-    "<!-- Author: --> **Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
+    "<!-- Author: -->  \n",
+    "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
     "\n",
-    "Date: **Sep 20, 2017**\n",
+    "Date: **Nov 21, 2017**\n",
     "\n",
     "Copyright 1999-2017, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
     "\n",
     "\n",
     "\n",
     "\n",
-    "## Introduction"
+    "## Regression analysis, overarching aims\n",
+    "\n",
+    "Regression modeling deals with the description of  the sampling distribution of a given random variable $y$ varies as function of another variable or a set of such variables $\\hat{x} =[x_0, x_1,\\dots, x_p]^T$. \n",
+    "The first variable is called the **dependent**, the **outcome** or the **response** variable while the set of variables $\\hat{x}$ is called the independent variable, or the predictor variable or the explanatory variable. \n",
+    "\n",
+    "A regression model aims at finding a likelihood function $p(y\\vert \\hat{x})$, that is the conditional distribution for $y$ with a given $\\hat{x}$. The estimation of  $p(y\\vert \\hat{x})$ is made using a data set with \n",
+    "* $n$ cases $i = 0, 1, 2, \\dots, n-1$ \n",
+    "\n",
+    "* Response (dependent or outcome) variable $y_i$ with $i = 0, 1, 2, \\dots, n-1$ \n",
+    "\n",
+    "* $p$ Explanatory (independent or predictor) variables $\\hat{x}_i=[x_{i0}, x_{i1}, \\dots, x_{ip}]$ with $i = 0, 1, 2, \\dots, n-1$   \n",
+    "\n",
+    " The goal of the regression analysis is to extract/exploit relationship between $y_i$ and $\\hat{x}_i$ in or to infer causal dependencies, approximations to the likelihood functions, functional relationships and to make predictions .\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "## General linear models\n",
+    "Before we proceed let us study a case from linear algebra where we aim at fitting a set of data $\\hat{y}=[y_0,y_1,\\dots,y_{n-1}]$. We could think of these data as a result of an experiment or a complicated numerical experiment. These data are functions of a series of variables $\\hat{x}=[x_0,x_1,\\dots,x_{n-1}]$, that is $y_i = y(x_i)$ with $i=0,1,2,\\dots,n-1$. The variables $x_i$ could represent physical quantities like time, temperature, position etc. We assume that $y(x)$ is a smooth function. \n",
+    "\n",
+    "Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of $y$ which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree $n-1$ with $n$ points, that is"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "y=y(x) \\rightarrow y(x_i)=\\tilde{y}_i+\\epsilon_i=\\sum_{j=0}^{n-1} \\beta_i x_i^j+\\epsilon_i,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "where $\\epsilon_i$ is the error in our approximation.\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "## Rewriting the fitting procedure as a linear algebra problem\n",
+    "For every set of values $y_i,x_i$ we have thus the corresponding set of equations"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\begin{align*}\n",
+    "y_0&=\\beta_0+\\beta_1x_0^1+\\beta_2x_0^2+\\dots+\\beta_{n-1}x_0^{n-1}+\\epsilon_0\\\\\n",
+    "y_1&=\\beta_0+\\beta_1x_1^1+\\beta_2x_1^2+\\dots+\\beta_{n-1}x_1^{n-1}+\\epsilon_1\\\\\n",
+    "y_2&=\\beta_0+\\beta_1x_2^1+\\beta_2x_2^2+\\dots+\\beta_{n-1}x_2^{n-1}+\\epsilon_2\\\\\n",
+    "\\dots & \\dots \\\\\n",
+    "y_{n-1}&=\\beta_0+\\beta_1x_{n-1}^1+\\beta_2x_{n-1}^2+\\dots+\\beta_1x_{n-1}^{n-1}+\\epsilon_{n-1}.\\\\\n",
+    "\\end{align*}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Rewriting the fitting procedure as a linear algebra problem, follows\n",
+    "Defining the vectors"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "2\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "3\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{\\epsilon} = [\\epsilon_0,\\epsilon_1, \\epsilon_2,\\dots, \\epsilon_{n-1}]^T,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and the matrix"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{X}=\n",
+    "\\begin{bmatrix} \n",
+    "1& x_{0}^1 &x_{0}^2& \\dots & \\dots &x_{0}^{n-1}\\\\\n",
+    "1& x_{1}^1 &x_{1}^2& \\dots & \\dots &x_{1}^{n-1}\\\\\n",
+    "1& x_{2}^1 &x_{2}^2& \\dots & \\dots &x_{2}^{n-1}\\\\                      \n",
+    "\\dots& \\dots &\\dots& \\dots & \\dots &\\dots\\\\\n",
+    "1& x_{n-1}^1 &x_{n-1}^2& \\dots & \\dots &x_{n-1}^{n-1}\\\\\n",
+    "\\end{bmatrix}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "we can rewrite our equations as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{y} = \\hat{X}\\hat{\\beta}+\\hat{\\epsilon}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generalizing the fitting procedure as a linear algebra problem\n",
+    "We are obviously not limited to the above polynomial. We could replace the various powers of $x$ with elements of Fourier series, that is, instead of $x_i^j$ we could have $\\cos{(j x_i)}$ or $\\sin{(j x_i)}$, or time series or other orthogonal functions.\n",
+    "For every set of values $y_i,x_i$ we can then generalize the equations to"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\begin{align*}\n",
+    "y_0&=\\beta_0x_{00}+\\beta_1x_{01}+\\beta_2x_{02}+\\dots+\\beta_{n-1}x_{0n-1}+\\epsilon_0\\\\\n",
+    "y_1&=\\beta_0x_{10}+\\beta_1x_{11}+\\beta_2x_{12}+\\dots+\\beta_{n-1}x_{1n-1}+\\epsilon_1\\\\\n",
+    "y_2&=\\beta_0x_{20}+\\beta_1x_{21}+\\beta_2x_{22}+\\dots+\\beta_{n-1}x_{2n-1}+\\epsilon_2\\\\\n",
+    "\\dots & \\dots \\\\\n",
+    "y_{i}&=\\beta_0x_{i0}+\\beta_1x_{i1}+\\beta_2x_{i2}+\\dots+\\beta_{n-1}x_{in-1}+\\epsilon_i\\\\\n",
+    "\\dots & \\dots \\\\\n",
+    "y_{n-1}&=\\beta_0x_{n-1,0}+\\beta_1x_{n-1,2}+\\beta_2x_{n-1,2}+\\dots+\\beta_1x_{n-1,n-1}+\\epsilon_{n-1}.\\\\\n",
+    "\\end{align*}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generalizing the fitting procedure as a linear algebra problem\n",
+    "We redefine in turn the matrix $\\hat{X}$ as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{X}=\n",
+    "\\begin{bmatrix} \n",
+    "x_{00}& x_{01} &x_{02}& \\dots & \\dots &x_{0,n-1}\\\\\n",
+    "x_{10}& x_{11} &x_{12}& \\dots & \\dots &x_{1,n-1}\\\\\n",
+    "x_{20}& x_{21} &x_{22}& \\dots & \\dots &x_{2,n-1}\\\\                      \n",
+    "\\dots& \\dots &\\dots& \\dots & \\dots &\\dots\\\\\n",
+    "x_{n-1,0}& x_{n-1,1} &x_{n-1,2}& \\dots & \\dots &x_{n-1,n-1}\\\\\n",
+    "\\end{bmatrix}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and without loss of generality we rewrite again  our equations as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{y} = \\hat{X}\\hat{\\beta}+\\hat{\\epsilon}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The left-hand side of this equation forms know. Our error vector $\\hat{\\epsilon}$ and the parameter vector $\\hat{\\beta}$ are our unknow quantities. How can we obtain the optimal set of $\\beta_i$ values?\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "## Optimizing our parameters\n",
+    "We have defined the matrix $\\hat{X}$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\begin{align*}\n",
+    "y_0&=\\beta_0x_{00}+\\beta_1x_{01}+\\beta_2x_{02}+\\dots+\\beta_{n-1}x_{0n-1}+\\epsilon_0\\\\\n",
+    "y_1&=\\beta_0x_{10}+\\beta_1x_{11}+\\beta_2x_{12}+\\dots+\\beta_{n-1}x_{1n-1}+\\epsilon_1\\\\\n",
+    "y_2&=\\beta_0x_{20}+\\beta_1x_{21}+\\beta_2x_{22}+\\dots+\\beta_{n-1}x_{2n-1}+\\epsilon_1\\\\\n",
+    "\\dots & \\dots \\\\\n",
+    "y_{i}&=\\beta_0x_{i0}+\\beta_1x_{i1}+\\beta_2x_{i2}+\\dots+\\beta_{n-1}x_{in-1}+\\epsilon_1\\\\\n",
+    "\\dots & \\dots \\\\\n",
+    "y_{n-1}&=\\beta_0x_{n-1,0}+\\beta_1x_{n-1,2}+\\beta_2x_{n-1,2}+\\dots+\\beta_1x_{n-1,n-1}+\\epsilon_{n-1}.\\\\\n",
+    "\\end{align*}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Optimizing our parameters, more details\n",
+    "We well use this matrix to define the approximation $\\hat{\\tilde{y}}$ via the unknown quantity $\\hat{\\beta}$ as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{\\tilde{y}}= \\hat{X}\\hat{\\beta},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "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 parametrized values $\\tilde{y}_i$, namely"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "Q(\\hat{\\beta})=\\sum_{i=0}^{n-1}\\left(y_i-\\tilde{y}_i\\right)^2=\\left(\\hat{y}-\\hat{\\tilde{y}}\\right)^T\\left(\\hat{y}-\\hat{\\tilde{y}}\\right),\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "or using the matrix $\\hat{X}$ as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "Q(\\hat{\\beta})=\\left(\\hat{y}-\\hat{X}\\hat{\\beta}\\right)^T\\left(\\hat{y}-\\hat{X}\\hat{\\beta}\\right).\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Interpretations and optimizing our parameters\n",
+    "The function"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "Q(\\hat{\\beta})=\\left(\\hat{y}-\\hat{X}\\hat{\\beta}\\right)^T\\left(\\hat{y}-\\hat{X}\\hat{\\beta}\\right),\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "can be linked to the variance of the quantity $y_i$ if we interpret the latter as the mean value of for example a numerical  experiment. When linking below with the maximum likelihood approach below, we will indeed interpret $y_i$ as a mean value"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "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,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "where $\\langle y_i \\rangle$ is the mean value. Keep in mind also that till now  we have treated $y_i$ as the exact value. Normally, the response (dependent or outcome) variable $y_i$ the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat $y_i$ as our exact value for the response variable.\n",
+    "\n",
+    "In order to find the parameters $\\beta_i$ we will then minimize the spread of $Q(\\hat{\\beta})$ by requiring"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial Q(\\hat{\\beta})}{\\partial \\beta_j} = \\frac{\\partial }{\\partial \\beta_j}\\left[ \\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,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "which results in"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial Q(\\hat{\\beta})}{\\partial \\beta_j} = -2\\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,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "or in a matrix-vector form as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial Q(\\hat{\\beta})}{\\partial \\hat{\\beta}} = 0 = \\hat{X}^T\\left( \\hat{y}-\\hat{X}\\hat{\\beta}\\right).\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Interpretations and optimizing our parameters\n",
+    "We can rewrite"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial Q(\\hat{\\beta})}{\\partial \\hat{\\beta}} = 0 = \\hat{X}^T\\left( \\hat{y}-\\hat{X}\\hat{\\beta}\\right),\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{X}^T\\hat{y} = \\hat{X}^T\\hat{X}\\hat{\\beta},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and if the matrix $\\hat{X}^T\\hat{X}$ is invertible we have the solution"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{\\beta} =\\left(\\hat{X}^T\\hat{X}\\right)^{-1}\\hat{X}^T\\hat{y}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Interpretations and optimizing our parameters\n",
+    "The residuals $\\hat{\\epsilon}$ are in turn given by"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{\\epsilon} = \\hat{y}-\\hat{\\tilde{y}} = \\hat{y}-\\hat{X}\\hat{\\beta},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and with"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{X}^T\\left( \\hat{y}-\\hat{X}\\hat{\\beta}\\right)= 0,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "we have"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{X}^T\\hat{\\epsilon}=\\hat{X}^T\\left( \\hat{y}-\\hat{X}\\hat{\\beta}\\right)= 0,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "meaning that the solution for $\\hat{\\beta}$ is the one which minimizes the residuals.  Later we will link this with the maximum likelihood approach.\n",
+    "\n",
+    "\n",
+    "\n",
+    "## The $\\chi^2$ function\n",
+    "\n",
+    "Normally, the response (dependent or outcome) variable $y_i$ the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat $y_i$ as our exact value for the response variable.\n",
+    "\n",
+    "Introducing the standard deviation $\\sigma_i$ for each measurement $y_i$, we define now the $\\chi^2$ function as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\chi^2(\\hat{\\beta})=\\sum_{i=0}^{n-1}\\frac{\\left(y_i-\\tilde{y}_i\\right)^2}{\\sigma_i^2}=\\left(\\hat{y}-\\hat{\\tilde{y}}\\right)^T\\frac{1}{\\hat{\\Sigma^2}}\\left(\\hat{y}-\\hat{\\tilde{y}}\\right),\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "where the matrix $\\hat{\\Sigma}$ is a diagonal matrix with $\\sigma_i$ as matrix elements.\n",
+    "\n",
+    "\n",
+    "\n",
+    "## The $\\chi^2$ function\n",
+    "\n",
+    "In order to find the parameters $\\beta_i$ we will then minimize the spread of $\\chi^2(\\hat{\\beta})$ by requiring"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\beta_j} = \\frac{\\partial }{\\partial \\beta_j}\\left[ \\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,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "which results in"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\beta_j} = -2\\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,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "or in a matrix-vector form as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\hat{\\beta}} = 0 = \\hat{A}^T\\left( \\hat{b}-\\hat{A}\\hat{\\beta}\\right).\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "where we have defined the matrix $\\hat{A} =\\hat{X}/\\hat{\\Sigma}$ with matrix elements $a_{ij} = x_{ij}/\\sigma_i$ and the vector $\\hat{b}$ with elements $b_i = y_i/\\sigma_i$.\n",
+    "\n",
+    "\n",
+    "\n",
+    "## The $\\chi^2$ function\n",
+    "\n",
+    "We can rewrite"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\hat{\\beta}} = 0 = \\hat{A}^T\\left( \\hat{b}-\\hat{A}\\hat{\\beta}\\right),\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{A}^T\\hat{b} = \\hat{A}^T\\hat{A}\\hat{\\beta},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and if the matrix $\\hat{A}^T\\hat{A}$ is invertible we have the solution"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{\\beta} =\\left(\\hat{A}^T\\hat{A}\\right)^{-1}\\hat{A}^T\\hat{b}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The $\\chi^2$ function\n",
+    "\n",
+    "If we then introduce the matrix"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{H} =  \\hat{A}^T\\hat{A},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "we have then the following expression for the parameters $\\beta_j$ (the matrix elements of $\\hat{H}$ are $h_{ij}$)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\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}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We state without proof the expression for the uncertainty  in the parameters $\\beta_j$ as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\sigma^2(\\beta_j) = \\sum_{i=0}^{n-1}\\sigma_i^2\\left( \\frac{\\partial \\beta_j}{\\partial y_i}\\right)^2,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "resulting in"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\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}!\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The $\\chi^2$ function\n",
+    "The first step here is to approximate the function $y$ with a first-order polynomial, that is we write"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "y=y(x) \\rightarrow y(x_i) \\approx \\beta_0+\\beta_1 x_i.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "By computing the derivatives of $\\chi^2$ with respect to $\\beta_0$ and $\\beta_1$ show that these are given by"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\beta_0} = -2\\left[ \\sum_{i=0}^{1}\\left(\\frac{y_i-\\beta_0-\\beta_1x_{i}}{\\sigma_i^2}\\right)\\right]=0,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\beta_0} = -2\\left[ \\sum_{i=0}^{1}x_i\\left(\\frac{y_i-\\beta_0-\\beta_1x_{i}}{\\sigma_i^2}\\right)\\right]=0.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The $\\chi^2$ function\n",
+    "\n",
+    "We define then"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\gamma =  \\sum_{i=0}^{1}\\frac{1}{\\sigma_i^2},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "4\n",
+    "0\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "4\n",
+    "1\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "4\n",
+    "2\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\gamma_{xy} = \\sum_{i=0}^{1}\\frac{y_ix_{i}}{\\sigma_i^2},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and show that"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "4\n",
+    "4\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\beta_1 = \\frac{\\gamma_{xy}\\gamma-\\gamma_x\\gamma_y}{\\gamma\\gamma_{xx}-\\gamma_x^2}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients $\\beta_i$.  A better approach is to use the Singular Value Decomposition (SVD) method discussed below.\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "## The singular value decompostion\n",
+    "How can we use the singular value decomposition to find the parameters $\\beta_j$? More details will come. We first note that a general $m\\times n$ matrix $\\hat{A}$ can be written in terms of a diagonal matrix $\\hat{\\Sigma}$ of dimensionality $n\\times n$ and two orthognal matrices $\\hat{U}$ and $\\hat{V}$, where the first has dimensionality $m \\times n$ and the last dimensionality $n\\times n$. We have then"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{A} = \\hat{U}\\hat{\\Sigma}\\hat{V}\n",
+    "$$"
    ]
   }
  ],
  "metadata": {},
  "nbformat": 4,
- "nbformat_minor": 1
+ "nbformat_minor": 2
 }
diff --git a/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz b/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz
index 18c3a279a..139731c2f 100644
Binary files a/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz and b/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz differ
diff --git a/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf b/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf
index 345684570..d560a1849 100644
Binary files a/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf and b/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf differ
diff --git a/doc/pub/Regression/pdf/Regression-beamer.pdf b/doc/pub/Regression/pdf/Regression-beamer.pdf
index de93f8c9d..c344916dd 100644
Binary files a/doc/pub/Regression/pdf/Regression-beamer.pdf and b/doc/pub/Regression/pdf/Regression-beamer.pdf differ
diff --git a/doc/pub/Regression/pdf/Regression-minted.pdf b/doc/pub/Regression/pdf/Regression-minted.pdf
index b34a4e212..d710633b3 100644
Binary files a/doc/pub/Regression/pdf/Regression-minted.pdf and b/doc/pub/Regression/pdf/Regression-minted.pdf differ
diff --git a/doc/src/Regression/.ipynb_checkpoints/Regression-checkpoint.ipynb b/doc/src/Regression/.ipynb_checkpoints/Regression-checkpoint.ipynb
new file mode 100644
index 000000000..210a74508
--- /dev/null
+++ b/doc/src/Regression/.ipynb_checkpoints/Regression-checkpoint.ipynb
@@ -0,0 +1,939 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<!-- dom:TITLE: Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis -->\n",
+    "# Data Analysis and Machine Learning: Linear Regression and more Advanced Regression Analysis\n",
+    "<!-- dom:AUTHOR: Morten Hjorth-Jensen at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University -->\n",
+    "<!-- Author: -->  \n",
+    "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
+    "\n",
+    "Date: **Nov 21, 2017**\n",
+    "\n",
+    "Copyright 1999-2017, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "## Regression analysis, overarching aims\n",
+    "\n",
+    "Regression modeling deals with the description of  the sampling distribution of a given random variable $y$ varies as function of another variable or a set of such variables $\\hat{x} =[x_0, x_1,\\dots, x_p]^T$. \n",
+    "The first variable is called the **dependent**, the **outcome** or the **response** variable while the set of variables $\\hat{x}$ is called the independent variable, or the predictor variable or the explanatory variable. \n",
+    "\n",
+    "A regression model aims at finding a likelihood function $p(y\\vert \\hat{x})$, that is the conditional distribution for $y$ with a given $\\hat{x}$. The estimation of  $p(y\\vert \\hat{x})$ is made using a data set with \n",
+    "* $n$ cases $i = 0, 1, 2, \\dots, n-1$ \n",
+    "\n",
+    "* Response (dependent or outcome) variable $y_i$ with $i = 0, 1, 2, \\dots, n-1$ \n",
+    "\n",
+    "* $p$ Explanatory (independent or predictor) variables $\\hat{x}_i=[x_{i0}, x_{i1}, \\dots, x_{ip}]$ with $i = 0, 1, 2, \\dots, n-1$   \n",
+    "\n",
+    " The goal of the regression analysis is to extract/exploit relationship between $y_i$ and $\\hat{x}_i$ in or to infer causal dependencies, approximations to the likelihood functions, functional relationships and to make predictions .\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "## General linear models\n",
+    "Before we proceed let us study a case from linear algebra where we aim at fitting a set of data $\\hat{y}=[y_0,y_1,\\dots,y_{n-1}]$. We could think of these data as a result of an experiment or a complicated numerical experiment. These data are functions of a series of variables $\\hat{x}=[x_0,x_1,\\dots,x_{n-1}]$, that is $y_i = y(x_i)$ with $i=0,1,2,\\dots,n-1$. The variables $x_i$ could represent physical quantities like time, temperature, position etc. We assume that $y(x)$ is a smooth function. \n",
+    "\n",
+    "Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of $y$ which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree $n-1$ with $n$ points, that is"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "y=y(x) \\rightarrow y(x_i)=\\tilde{y}_i+\\epsilon_i=\\sum_{j=0}^{n-1} \\beta_i x_i^j+\\epsilon_i,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "where $\\epsilon_i$ is the error in our approximation.\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "## Rewriting the fitting procedure as a linear algebra problem\n",
+    "For every set of values $y_i,x_i$ we have thus the corresponding set of equations"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\begin{align*}\n",
+    "y_0&=\\beta_0+\\beta_1x_0^1+\\beta_2x_0^2+\\dots+\\beta_{n-1}x_0^{n-1}+\\epsilon_0\\\\\n",
+    "y_1&=\\beta_0+\\beta_1x_1^1+\\beta_2x_1^2+\\dots+\\beta_{n-1}x_1^{n-1}+\\epsilon_1\\\\\n",
+    "y_2&=\\beta_0+\\beta_1x_2^1+\\beta_2x_2^2+\\dots+\\beta_{n-1}x_2^{n-1}+\\epsilon_2\\\\\n",
+    "\\dots & \\dots \\\\\n",
+    "y_{n-1}&=\\beta_0+\\beta_1x_{n-1}^1+\\beta_2x_{n-1}^2+\\dots+\\beta_1x_{n-1}^{n-1}+\\epsilon_{n-1}.\\\\\n",
+    "\\end{align*}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Rewriting the fitting procedure as a linear algebra problem, follows\n",
+    "Defining the vectors"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "2\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "3\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{\\epsilon} = [\\epsilon_0,\\epsilon_1, \\epsilon_2,\\dots, \\epsilon_{n-1}]^T,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and the matrix"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{X}=\n",
+    "\\begin{bmatrix} \n",
+    "1& x_{0}^1 &x_{0}^2& \\dots & \\dots &x_{0}^{n-1}\\\\\n",
+    "1& x_{1}^1 &x_{1}^2& \\dots & \\dots &x_{1}^{n-1}\\\\\n",
+    "1& x_{2}^1 &x_{2}^2& \\dots & \\dots &x_{2}^{n-1}\\\\                      \n",
+    "\\dots& \\dots &\\dots& \\dots & \\dots &\\dots\\\\\n",
+    "1& x_{n-1}^1 &x_{n-1}^2& \\dots & \\dots &x_{n-1}^{n-1}\\\\\n",
+    "\\end{bmatrix}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "we can rewrite our equations as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{y} = \\hat{X}\\hat{\\beta}+\\hat{\\epsilon}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generalizing the fitting procedure as a linear algebra problem\n",
+    "We are obviously not limited to the above polynomial. We could replace the various powers of $x$ with elements of Fourier series, that is, instead of $x_i^j$ we could have $\\cos{(j x_i)}$ or $\\sin{(j x_i)}$, or time series or other orthogonal functions.\n",
+    "For every set of values $y_i,x_i$ we can then generalize the equations to"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\begin{align*}\n",
+    "y_0&=\\beta_0x_{00}+\\beta_1x_{01}+\\beta_2x_{02}+\\dots+\\beta_{n-1}x_{0n-1}+\\epsilon_0\\\\\n",
+    "y_1&=\\beta_0x_{10}+\\beta_1x_{11}+\\beta_2x_{12}+\\dots+\\beta_{n-1}x_{1n-1}+\\epsilon_1\\\\\n",
+    "y_2&=\\beta_0x_{20}+\\beta_1x_{21}+\\beta_2x_{22}+\\dots+\\beta_{n-1}x_{2n-1}+\\epsilon_2\\\\\n",
+    "\\dots & \\dots \\\\\n",
+    "y_{i}&=\\beta_0x_{i0}+\\beta_1x_{i1}+\\beta_2x_{i2}+\\dots+\\beta_{n-1}x_{in-1}+\\epsilon_i\\\\\n",
+    "\\dots & \\dots \\\\\n",
+    "y_{n-1}&=\\beta_0x_{n-1,0}+\\beta_1x_{n-1,2}+\\beta_2x_{n-1,2}+\\dots+\\beta_1x_{n-1,n-1}+\\epsilon_{n-1}.\\\\\n",
+    "\\end{align*}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generalizing the fitting procedure as a linear algebra problem\n",
+    "We redefine in turn the matrix $\\hat{X}$ as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{X}=\n",
+    "\\begin{bmatrix} \n",
+    "x_{00}& x_{01} &x_{02}& \\dots & \\dots &x_{0,n-1}\\\\\n",
+    "x_{10}& x_{11} &x_{12}& \\dots & \\dots &x_{1,n-1}\\\\\n",
+    "x_{20}& x_{21} &x_{22}& \\dots & \\dots &x_{2,n-1}\\\\                      \n",
+    "\\dots& \\dots &\\dots& \\dots & \\dots &\\dots\\\\\n",
+    "x_{n-1,0}& x_{n-1,1} &x_{n-1,2}& \\dots & \\dots &x_{n-1,n-1}\\\\\n",
+    "\\end{bmatrix}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and without loss of generality we rewrite again  our equations as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{y} = \\hat{X}\\hat{\\beta}+\\hat{\\epsilon}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The left-hand side of this equation forms know. Our error vector $\\hat{\\epsilon}$ and the parameter vector $\\hat{\\beta}$ are our unknow quantities. How can we obtain the optimal set of $\\beta_i$ values?\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "## Optimizing our parameters\n",
+    "We have defined the matrix $\\hat{X}$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\begin{align*}\n",
+    "y_0&=\\beta_0x_{00}+\\beta_1x_{01}+\\beta_2x_{02}+\\dots+\\beta_{n-1}x_{0n-1}+\\epsilon_0\\\\\n",
+    "y_1&=\\beta_0x_{10}+\\beta_1x_{11}+\\beta_2x_{12}+\\dots+\\beta_{n-1}x_{1n-1}+\\epsilon_1\\\\\n",
+    "y_2&=\\beta_0x_{20}+\\beta_1x_{21}+\\beta_2x_{22}+\\dots+\\beta_{n-1}x_{2n-1}+\\epsilon_1\\\\\n",
+    "\\dots & \\dots \\\\\n",
+    "y_{i}&=\\beta_0x_{i0}+\\beta_1x_{i1}+\\beta_2x_{i2}+\\dots+\\beta_{n-1}x_{in-1}+\\epsilon_1\\\\\n",
+    "\\dots & \\dots \\\\\n",
+    "y_{n-1}&=\\beta_0x_{n-1,0}+\\beta_1x_{n-1,2}+\\beta_2x_{n-1,2}+\\dots+\\beta_1x_{n-1,n-1}+\\epsilon_{n-1}.\\\\\n",
+    "\\end{align*}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Optimizing our parameters, more details\n",
+    "We well use this matrix to define the approximation $\\hat{\\tilde{y}}$ via the unknown quantity $\\hat{\\beta}$ as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{\\tilde{y}}= \\hat{X}\\hat{\\beta},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "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 parametrized values $\\tilde{y}_i$, namely"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "Q(\\hat{\\beta})=\\sum_{i=0}^{n-1}\\left(y_i-\\tilde{y}_i\\right)^2=\\left(\\hat{y}-\\hat{\\tilde{y}}\\right)^T\\left(\\hat{y}-\\hat{\\tilde{y}}\\right),\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "or using the matrix $\\hat{X}$ as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "Q(\\hat{\\beta})=\\left(\\hat{y}-\\hat{X}\\hat{\\beta}\\right)^T\\left(\\hat{y}-\\hat{X}\\hat{\\beta}\\right).\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Interpretations and optimizing our parameters\n",
+    "The function"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "Q(\\hat{\\beta})=\\left(\\hat{y}-\\hat{X}\\hat{\\beta}\\right)^T\\left(\\hat{y}-\\hat{X}\\hat{\\beta}\\right),\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "can be linked to the variance of the quantity $y_i$ if we interpret the latter as the mean value of for example a numerical  experiment. When linking below with the maximum likelihood approach below, we will indeed interpret $y_i$ as a mean value"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "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,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "where $\\langle y_i \\rangle$ is the mean value. Keep in mind also that till now  we have treated $y_i$ as the exact value. Normally, the response (dependent or outcome) variable $y_i$ the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat $y_i$ as our exact value for the response variable.\n",
+    "\n",
+    "In order to find the parameters $\\beta_i$ we will then minimize the spread of $Q(\\hat{\\beta})$ by requiring"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial Q(\\hat{\\beta})}{\\partial \\beta_j} = \\frac{\\partial }{\\partial \\beta_j}\\left[ \\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,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "which results in"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial Q(\\hat{\\beta})}{\\partial \\beta_j} = -2\\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,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "or in a matrix-vector form as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial Q(\\hat{\\beta})}{\\partial \\hat{\\beta}} = 0 = \\hat{X}^T\\left( \\hat{y}-\\hat{X}\\hat{\\beta}\\right).\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Interpretations and optimizing our parameters\n",
+    "We can rewrite"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial Q(\\hat{\\beta})}{\\partial \\hat{\\beta}} = 0 = \\hat{X}^T\\left( \\hat{y}-\\hat{X}\\hat{\\beta}\\right),\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{X}^T\\hat{y} = \\hat{X}^T\\hat{X}\\hat{\\beta},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and if the matrix $\\hat{X}^T\\hat{X}$ is invertible we have the solution"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{\\beta} =\\left(\\hat{X}^T\\hat{X}\\right)^{-1}\\hat{X}^T\\hat{y}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Interpretations and optimizing our parameters\n",
+    "The residuals $\\hat{\\epsilon}$ are in turn given by"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{\\epsilon} = \\hat{y}-\\hat{\\tilde{y}} = \\hat{y}-\\hat{X}\\hat{\\beta},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and with"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{X}^T\\left( \\hat{y}-\\hat{X}\\hat{\\beta}\\right)= 0,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "we have"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{X}^T\\hat{\\epsilon}=\\hat{X}^T\\left( \\hat{y}-\\hat{X}\\hat{\\beta}\\right)= 0,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "meaning that the solution for $\\hat{\\beta}$ is the one which minimizes the residuals.  Later we will link this with the maximum likelihood approach.\n",
+    "\n",
+    "\n",
+    "\n",
+    "## The $\\chi^2$ function\n",
+    "\n",
+    "Normally, the response (dependent or outcome) variable $y_i$ the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat $y_i$ as our exact value for the response variable.\n",
+    "\n",
+    "Introducing the standard deviation $\\sigma_i$ for each measurement $y_i$, we define now the $\\chi^2$ function as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\chi^2(\\hat{\\beta})=\\sum_{i=0}^{n-1}\\frac{\\left(y_i-\\tilde{y}_i\\right)^2}{\\sigma_i^2}=\\left(\\hat{y}-\\hat{\\tilde{y}}\\right)^T\\frac{1}{\\hat{\\Sigma^2}}\\left(\\hat{y}-\\hat{\\tilde{y}}\\right),\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "where the matrix $\\hat{\\Sigma}$ is a diagonal matrix with $\\sigma_i$ as matrix elements.\n",
+    "\n",
+    "\n",
+    "\n",
+    "## The $\\chi^2$ function\n",
+    "\n",
+    "In order to find the parameters $\\beta_i$ we will then minimize the spread of $\\chi^2(\\hat{\\beta})$ by requiring"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\beta_j} = \\frac{\\partial }{\\partial \\beta_j}\\left[ \\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,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "which results in"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\beta_j} = -2\\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,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "or in a matrix-vector form as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\hat{\\beta}} = 0 = \\hat{A}^T\\left( \\hat{b}-\\hat{A}\\hat{\\beta}\\right).\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "where we have defined the matrix $\\hat{A} =\\hat{X}/\\hat{\\Sigma}$ with matrix elements $a_{ij} = x_{ij}/\\sigma_i$ and the vector $\\hat{b}$ with elements $b_i = y_i/\\sigma_i$.\n",
+    "\n",
+    "\n",
+    "\n",
+    "## The $\\chi^2$ function\n",
+    "\n",
+    "We can rewrite"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\hat{\\beta}} = 0 = \\hat{A}^T\\left( \\hat{b}-\\hat{A}\\hat{\\beta}\\right),\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{A}^T\\hat{b} = \\hat{A}^T\\hat{A}\\hat{\\beta},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and if the matrix $\\hat{A}^T\\hat{A}$ is invertible we have the solution"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{\\beta} =\\left(\\hat{A}^T\\hat{A}\\right)^{-1}\\hat{A}^T\\hat{b}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The $\\chi^2$ function\n",
+    "\n",
+    "If we then introduce the matrix"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{H} =  \\hat{A}^T\\hat{A},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "we have then the following expression for the parameters $\\beta_j$ (the matrix elements of $\\hat{H}$ are $h_{ij}$)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\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}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We state without proof the expression for the uncertainty  in the parameters $\\beta_j$ as"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\sigma^2(\\beta_j) = \\sum_{i=0}^{n-1}\\sigma_i^2\\left( \\frac{\\partial \\beta_j}{\\partial y_i}\\right)^2,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "resulting in"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\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}!\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The $\\chi^2$ function\n",
+    "The first step here is to approximate the function $y$ with a first-order polynomial, that is we write"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "y=y(x) \\rightarrow y(x_i) \\approx \\beta_0+\\beta_1 x_i.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "By computing the derivatives of $\\chi^2$ with respect to $\\beta_0$ and $\\beta_1$ show that these are given by"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\beta_0} = -2\\left[ \\sum_{i=0}^{1}\\left(\\frac{y_i-\\beta_0-\\beta_1x_{i}}{\\sigma_i^2}\\right)\\right]=0,\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\frac{\\partial \\chi^2(\\hat{\\beta})}{\\partial \\beta_0} = -2\\left[ \\sum_{i=0}^{1}x_i\\left(\\frac{y_i-\\beta_0-\\beta_1x_{i}}{\\sigma_i^2}\\right)\\right]=0.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The $\\chi^2$ function\n",
+    "\n",
+    "We define then"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\gamma =  \\sum_{i=0}^{1}\\frac{1}{\\sigma_i^2},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "4\n",
+    "0\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "4\n",
+    "1\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "4\n",
+    "2\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\gamma_{xy} = \\sum_{i=0}^{1}\\frac{y_ix_{i}}{\\sigma_i^2},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and show that"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "4\n",
+    "4\n",
+    " \n",
+    "<\n",
+    "<\n",
+    "<\n",
+    "!\n",
+    "!\n",
+    "M\n",
+    "A\n",
+    "T\n",
+    "H\n",
+    "_\n",
+    "B\n",
+    "L\n",
+    "O\n",
+    "C\n",
+    "K"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\beta_1 = \\frac{\\gamma_{xy}\\gamma-\\gamma_x\\gamma_y}{\\gamma\\gamma_{xx}-\\gamma_x^2}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients $\\beta_i$.  A better approach is to use the Singular Value Decomposition (SVD) method discussed below.\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "## The singular value decompostion\n",
+    "How can we use the singular value decomposition to find the parameters $\\beta_j$? More details will come. We first note that a general $m\\times n$ matrix $\\hat{A}$ can be written in terms of a diagonal matrix $\\hat{\\Sigma}$ of dimensionality $n\\times n$ and two orthognal matrices $\\hat{U}$ and $\\hat{V}$, where the first has dimensionality $m \\times n$ and the last dimensionality $n\\times n$. We have then"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\hat{A} = \\hat{U}\\hat{\\Sigma}\\hat{V}\n",
+    "$$"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/doc/src/Regression/make.sh b/doc/src/Regression/make.sh
index 0671a98d3..d73ef27fc 100755
--- a/doc/src/Regression/make.sh
+++ b/doc/src/Regression/make.sh
@@ -44,10 +44,10 @@ system doconce split_html $html.html --method=space10
 # Bootstrap style
 html=${name}-bs
 system doconce format html $name --html_style=bootstrap --pygments_html_style=default --html_admon=bootstrap_panel --html_output=$html $opt
-system doconce split_html $html.html --method=split --pagination --nav_button=bottom
+#system doconce split_html $html.html --method=split --pagination --nav_button=bottom
 
 # IPython notebook
-#system doconce format ipynb $name $opt
+system doconce format ipynb $name $opt
 
 # LaTeX Beamer slides
 beamertheme=red_plain