diff --git a/doc/src/Regression/Regression.do.txt b/doc/src/Regression/Regression.do.txt
index d85574dd6..cad65b121 100644
--- a/doc/src/Regression/Regression.do.txt
+++ b/doc/src/Regression/Regression.do.txt
@@ -92,9 +92,9 @@ 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_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_1\\
+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,n-1}+\epsilon_{n-1}.\\
 \end{align}
@@ -125,6 +125,7 @@ The left-hand side of this equation forms know. Our error vector $\hat{\epsilon}
 !split
 ===== Optimizing our parameters  =====
 !bblock
+We have defined the matrix $\hat{X}$
 !bt
 \begin{align}
 y_0&=\beta_0x_{00}+\beta_1x_{01}+\beta_2x_{02}+\dots+\beta_{n-1}x_{0n-1}+\epsilon_0\\
@@ -133,7 +134,65 @@ y_2&=\beta_0x_{20}+\beta_1x_{21}+\beta_2x_{22}+\dots+\beta_{n-1}x_{2n-1}+\epsilo
 \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,n-1}+\epsilon_{n-1}.\\
+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}
 !et
+We well use this matrix to define the approximation $\hat{\tilde{y}}$ via the unknown quantity $\hat{\beta}$ as
+!bt
+\[
+\hat{\tilde{y}}= \hat{X}\hat{\beta},
+\]
+!et
+and in order to find the optimal parameters $\beta_i$ instead of solving the above linear algebra problem, we define a function which gives a measure of the spread between the values $y_i$ (which represent hopefully the exact values) and the parametrized values $\tilde{y}_i$, namely
+!bt
+\[
+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),
+\]
+!et
+or using the matrix $\hat{X}$ as
+!bt
+\[
+Q(\hat{\beta})=\left \hat{y}-\hat{X}\hat{\beta}\right)^T\left \hat{y}-\hat{X}\hat{\beta}\right).
+\]
+!et
+!eblock
+
+
+!split
+===== Interpretations and optimizing our parameters  =====
+!bblock
+The function 
+!bt
+\[
+Q(\hat{\beta})=\left \hat{y}-\hat{X}\hat{\beta}\right)^T\left \hat{y}-\hat{X}\hat{\beta}\right),
+\]
+!et
+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
+!bt
+\[
+y_{i}=\langle y_i \rangle = \beta_0x_{i,0}+\beta_1x_{i,1}+\beta_2x_{i,2}+\dots+\beta_{n-1}x_{i,n-1}+\epsilon_i,
+\]
+!et
+where $\langle y_i \rangle$ is the mean value. Keep in mind also that till now  we have treated $y_i$ as the exact value. Normally, the response (dependent or outcome) variable $y_i$ the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat $y_i$ as our exact value for the response variable.
+
+In order to find the parameters $\beta_i$ we will then minimize the spread of $Q(\hat{\beta})$ by requiring
+!bt
+\[
+\frac{\partial Q(\hat{\beta})}{\partial \beta_j} = \frac{ }{\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, 
+\]
+!et
+which results in
+!bt
+\[
+\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, 
+\]
+!et
+or in a matrix-vector form as
+!bt
+\[
+\frac{\partial Q(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{X}^T\left( \hat{y}-\hat{X}\hat{\beta}\right).  
+\]
+!et
+
+
 !eblock