update

2023-08-28 10:35:01 +02:00
parent b4954d5316
commit 707c8e08f3
15 changed files with 1891 additions and 1772 deletions
@@ -1137,8 +1137,8 @@ doconce format html week35.do.txt --no_mako -->
 <h3>Reading recommendations:<a class="headerlink" href="#reading-recommendations" title="Permalink to this headline">¶</a></h3>
 <ol class="simple">
 <li><p>See lecture notes for week 35 at <a class="reference external" href="https://compphysics.github.io/MachineLearning/doc/web/course.html">https://compphysics.github.io/MachineLearning/doc/web/course.html</a></p></li>
-<li><p>Goodfellow, Bengio and Courville, Deep Learning, chapter 2 on linear algebra and sections 3.1-3.10 on elements of statistics</p></li>
-<li><p>Hastie, Tibshirani and Friedman, The elements of statistical learning, sections 3.1-3.4</p></li>
+<li><p>Goodfellow, Bengio and Courville, Deep Learning, chapter 2 on linear algebra and sections 3.1-3.10 on elements of statistics (background)</p></li>
+<li><p>Hastie, Tibshirani and Friedman, The elements of statistical learning, sections 3.1-3.4 (on relevance for the discussion of linear regression).</p></li>
 </ol>
 </div>
 </div>
@@ -1278,14 +1278,13 @@ will treat <span class="math notranslate nohighlight">\(y_i\)</span> as our exac
 \]</div>
 <p>We note also that since our design matrix is defined as <span class="math notranslate nohighlight">\(\boldsymbol{X}\in
 {\mathbb{R}}^{n\times p}\)</span>, the product <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X} \in
-{\mathbb{R}}^{p\times p}\)</span>.  In the above case we have that <span class="math notranslate nohighlight">\(p \ll n\)</span>,
-in our case <span class="math notranslate nohighlight">\(p=5\)</span> meaning that we end up with inverting a small
+{\mathbb{R}}^{p\times p}\)</span>.  In most cases we have that <span class="math notranslate nohighlight">\(p \ll n\)</span>. In our example case below we have <span class="math notranslate nohighlight">\(p=5\)</span> meaning. We end up with inverting a small
 <span class="math notranslate nohighlight">\(5\times 5\)</span> matrix. This is a rather common situation, in many cases we end up with low-dimensional
 matrices to invert. The methods discussed here and for many other
 supervised learning algorithms like classification with logistic
 regression or support vector machines, exhibit dimensionalities which
 allow for the usage of direct linear algebra methods such as <strong>LU</strong> decomposition or <strong>Singular Value Decomposition</strong> (SVD) for finding the inverse of the matrix
-<span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>.</p>
+<span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>.  This is discussed on Thursday this week.</p>
 <p><strong>Small question</strong>: Do you think the example we have at hand here (the nuclear binding energies) can lead to problems in inverting the matrix  <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>? What kind of problems can we expect?</p>
 </div>
 <div class="section" id="some-useful-matrix-and-vector-expressions">
@@ -1453,7 +1452,7 @@ C(\boldsymbol{\beta})=\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\bold
 \]</div>
 <p>We note that the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> does not depend on the unknown parameters defined by the vector <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>.
 We are now interested in minimizing the cost function with respect to the unknown parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>.</p>
-<p>The mean squared error is a scalar and if we use the results from the last example, we define a new vector</p>
+<p>The mean squared error is a scalar and if we use the results from example three above, we can define a new vector</p>
 <div class="math notranslate nohighlight">
 \[
 \boldsymbol{w}=\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta},
@@ -1553,8 +1552,8 @@ We assume our data can represented by a fourth-order polynomial. For the <span c
 \[
 \tilde{y}_i = \beta_0+\beta_1x_i+\beta_2x_i^2+\beta_3x_i^3+\beta_4x_i^4.
 \]</div>
-<p>we have five predictors/features. The first is the intercept <span class="math notranslate nohighlight">\(\beta_0\)</span>. The other terms are <span class="math notranslate nohighlight">\(\beta_i\)</span> with <span class="math notranslate nohighlight">\(i=1,2,3,4\)</span>. Furthermore we have <span class="math notranslate nohighlight">\(n\)</span> entries for each predictor. It means that our design matrix is a
-<span class="math notranslate nohighlight">\(p\times n\)</span> matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
+<p>we have five predictors/features. The first is the intercept <span class="math notranslate nohighlight">\(\beta_0\)</span>. The other terms are <span class="math notranslate nohighlight">\(\beta_i\)</span> with <span class="math notranslate nohighlight">\(i=1,2,3,4\)</span>. Furthermore we have <span class="math notranslate nohighlight">\(n\)</span> entries for each predictor. It means that our design matrix is an
+<span class="math notranslate nohighlight">\(n\times p\)</span> matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
 </div>
 <div class="section" id="own-code-for-ordinary-least-squares">
 <h2>Own code for Ordinary Least Squares<a class="headerlink" href="#own-code-for-ordinary-least-squares" title="Permalink to this headline">¶</a></h2>
@@ -1562,25 +1561,24 @@ We assume our data can represented by a fourth-order polynomial. For the <span c
 <div class="cell docutils container">
 <div class="cell_input docutils container">
 <div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># matrix inversion to find beta</span>
+<span class="c1"># First we set up the data</span>
+<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
+<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">rand</span><span class="p">(</span><span class="mi">100</span><span class="p">)</span>
+<span class="n">y</span> <span class="o">=</span> <span class="mf">2.0</span><span class="o">+</span><span class="mi">5</span><span class="o">*</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="o">+</span><span class="mf">0.1</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">randn</span><span class="p">(</span><span class="mi">100</span><span class="p">)</span>
+<span class="c1"># and then the design matrix X including the intercept</span>
+<span class="c1">#  The design matrix now as function of a fourth-order polynomial</span>
+<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">),</span><span class="mi">5</span><span class="p">))</span>
+<span class="n">X</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>
+<span class="n">X</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span>
+<span class="n">X</span><span class="p">[:,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="mi">2</span>
+<span class="n">X</span><span class="p">[:,</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="mi">3</span>
+<span class="n">X</span><span class="p">[:,</span><span class="mi">4</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="mi">4</span>
 <span class="n">beta</span> <span class="o">=</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X</span><span class="p">)</span> <span class="o">@</span> <span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="p">)</span> <span class="o">@</span> <span class="n">y</span>
 <span class="c1"># and then make the prediction</span>
 <span class="n">ytilde</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta</span>
 </pre></div>
 </div>
 </div>
-<div class="cell_output docutils container">
-<div class="output traceback highlight-ipythontb notranslate"><div class="highlight"><pre><span></span><span class="gt">---------------------------------------------------------------------------</span>
-<span class="ne">NameError</span><span class="g g-Whitespace">                                 </span>Traceback (most recent call last)
-<span class="nn">Input In [1],</span> in <span class="ni">&lt;cell line: 2&gt;</span><span class="nt">()</span>
-<span class="g g-Whitespace">      </span><span class="mi">1</span> <span class="c1"># matrix inversion to find beta</span>
-<span class="ne">----&gt; </span><span class="mi">2</span> <span class="n">beta</span> <span class="o">=</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X</span><span class="p">)</span> <span class="o">@</span> <span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="p">)</span> <span class="o">@</span> <span class="n">y</span>
-<span class="g g-Whitespace">      </span><span class="mi">3</span> <span class="c1"># and then make the prediction</span>
-<span class="g g-Whitespace">      </span><span class="mi">4</span> <span class="n">ytilde</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta</span>
-
-<span class="ne">NameError</span>: name &#39;np&#39; is not defined
-</pre></div>
-</div>
-</div>
 </div>
 <p>Alternatively, you can use the least squares functionality in <strong>Numpy</strong> as</p>
 <div class="cell docutils container">
@@ -1590,6 +1588,17 @@ We assume our data can represented by a fourth-order polynomial. For the <span c
 </pre></div>
 </div>
 </div>
+<div class="cell_output docutils container">
+<div class="output traceback highlight-ipythontb notranslate"><div class="highlight"><pre><span></span><span class="gt">---------------------------------------------------------------------------</span>
+<span class="ne">NameError</span><span class="g g-Whitespace">                                 </span>Traceback (most recent call last)
+<span class="nn">Input In [2],</span> in <span class="ni">&lt;cell line: 1&gt;</span><span class="nt">()</span>
+<span class="ne">----&gt; </span><span class="mi">1</span> <span class="n">fit</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">lstsq</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">Energies</span><span class="p">,</span> <span class="n">rcond</span> <span class="o">=</span><span class="kc">None</span><span class="p">)[</span><span class="mi">0</span><span class="p">]</span>
+<span class="g g-Whitespace">      </span><span class="mi">2</span> <span class="n">ytildenp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">fit</span><span class="p">,</span><span class="n">X</span><span class="o">.</span><span class="n">T</span><span class="p">)</span>
+
+<span class="ne">NameError</span>: name &#39;Energies&#39; is not defined
+</pre></div>
+</div>
+</div>
 </div>
 </div>
 <div class="section" id="adding-error-analysis-and-training-set-up">
@@ -28,9 +28,9 @@
 # 
 # 1. See lecture notes for week 35 at <https://compphysics.github.io/MachineLearning/doc/web/course.html>
 # 
-# 2. Goodfellow, Bengio and Courville, Deep Learning, chapter 2 on linear algebra and sections 3.1-3.10 on elements of statistics
+# 2. Goodfellow, Bengio and Courville, Deep Learning, chapter 2 on linear algebra and sections 3.1-3.10 on elements of statistics (background)
 # 
-# 3. Hastie, Tibshirani and Friedman, The elements of statistical learning, sections 3.1-3.4
+# 3. Hastie, Tibshirani and Friedman, The elements of statistical learning, sections 3.1-3.4 (on relevance for the discussion of linear regression).

 # ## Why Linear Regression (aka Ordinary Least Squares and family), repeat from last week
 # 
@@ -193,14 +193,13 @@

 # We note also that since our design matrix is defined as $\boldsymbol{X}\in
 # {\mathbb{R}}^{n\times p}$, the product $\boldsymbol{X}^T\boldsymbol{X} \in
-# {\mathbb{R}}^{p\times p}$.  In the above case we have that $p \ll n$,
-# in our case $p=5$ meaning that we end up with inverting a small
+# {\mathbb{R}}^{p\times p}$.  In most cases we have that $p \ll n$. In our example case below we have $p=5$ meaning. We end up with inverting a small
 # $5\times 5$ matrix. This is a rather common situation, in many cases we end up with low-dimensional
 # matrices to invert. The methods discussed here and for many other
 # supervised learning algorithms like classification with logistic
 # regression or support vector machines, exhibit dimensionalities which
 # allow for the usage of direct linear algebra methods such as **LU** decomposition or **Singular Value Decomposition** (SVD) for finding the inverse of the matrix
-# $\boldsymbol{X}^T\boldsymbol{X}$.
+# $\boldsymbol{X}^T\boldsymbol{X}$.  This is discussed on Thursday this week.
 # 
 # **Small question**: Do you think the example we have at hand here (the nuclear binding energies) can lead to problems in inverting the matrix  $\boldsymbol{X}^T\boldsymbol{X}$? What kind of problems can we expect?

@@ -390,7 +389,7 @@
 # We note that the design matrix $\boldsymbol{X}$ does not depend on the unknown parameters defined by the vector $\boldsymbol{\beta}$.
 # We are now interested in minimizing the cost function with respect to the unknown parameters $\boldsymbol{\beta}$.
 # 
-# The mean squared error is a scalar and if we use the results from the last example, we define a new vector
+# The mean squared error is a scalar and if we use the results from example three above, we can define a new vector

 # $$
 # \boldsymbol{w}=\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta},
@@ -502,8 +501,8 @@
 # \tilde{y}_i = \beta_0+\beta_1x_i+\beta_2x_i^2+\beta_3x_i^3+\beta_4x_i^4.
 # $$

-# we have five predictors/features. The first is the intercept $\beta_0$. The other terms are $\beta_i$ with $i=1,2,3,4$. Furthermore we have $n$ entries for each predictor. It means that our design matrix is a 
-# $p\times n$ matrix $\boldsymbol{X}$.
+# we have five predictors/features. The first is the intercept $\beta_0$. The other terms are $\beta_i$ with $i=1,2,3,4$. Furthermore we have $n$ entries for each predictor. It means that our design matrix is an 
+# $n\times p$ matrix $\boldsymbol{X}$.

 # ## Own code for Ordinary Least Squares
 # 
@@ -513,6 +512,18 @@


 # matrix inversion to find beta
+# First we set up the data
+import numpy as np
+x = np.random.rand(100)
+y = 2.0+5*x*x+0.1*np.random.randn(100)
+# and then the design matrix X including the intercept
+#  The design matrix now as function of a fourth-order polynomial
+X = np.zeros((len(x),5))
+X[:,0] = 1.0
+X[:,1] = x
+X[:,2] = x**2
+X[:,3] = x**3
+X[:,4] = x**4
 beta = (np.linalg.inv(X.T @ X) @ X.T ) @ y
 # and then make the prediction
 ytilde = X @ beta
@@ -378,7 +378,18 @@ MathJax.Hub.Config({
      <div class="input_area">
        <div class="highlight" style="background: #f8f8f8">
  <pre style="line-height: 125%;"><span style="color: #408080; font-style: italic"># matrix inversion to find beta</span>
+<span style="color: #408080; font-style: italic"># First we set up the data</span>
 <span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
+x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(<span style="color: #666666">100</span>)
+y <span style="color: #666666">=</span> <span style="color: #666666">2.0+5*</span>x<span style="color: #666666">*</span>x<span style="color: #666666">+0.1*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(<span style="color: #666666">100</span>)
+<span style="color: #408080; font-style: italic"># and then the design matrix X including the intercept</span>
+<span style="color: #408080; font-style: italic">#  The design matrix now as function of a fourth-order polynomial</span>
+X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros((<span style="color: #008000">len</span>(x),<span style="color: #666666">5</span>))
+X[:,<span style="color: #666666">0</span>] <span style="color: #666666">=</span> <span style="color: #666666">1.0</span>
+X[:,<span style="color: #666666">1</span>] <span style="color: #666666">=</span> x
+X[:,<span style="color: #666666">2</span>] <span style="color: #666666">=</span> x<span style="color: #666666">**2</span>
+X[:,<span style="color: #666666">3</span>] <span style="color: #666666">=</span> x<span style="color: #666666">**3</span>
+X[:,<span style="color: #666666">4</span>] <span style="color: #666666">=</span> x<span style="color: #666666">**4</span>
 beta <span style="color: #666666">=</span> (np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X) <span style="color: #666666">@</span> X<span style="color: #666666">.</span>T ) <span style="color: #666666">@</span> y
 <span style="color: #408080; font-style: italic"># and then make the prediction</span>
 ytilde <span style="color: #666666">=</span> X <span style="color: #666666">@</span> beta
@@ -828,7 +828,18 @@ $$
      <div class="input_area">
        <div class="highlight" style="background: #eeeedd">
  <pre style="font-size: 80%; line-height: 125%;"><span style="color: #228B22"># matrix inversion to find beta</span>
+<span style="color: #228B22"># First we set up the data</span>
 <span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
+x = np.random.rand(<span style="color: #B452CD">100</span>)
+y = <span style="color: #B452CD">2.0</span>+<span style="color: #B452CD">5</span>*x*x+<span style="color: #B452CD">0.1</span>*np.random.randn(<span style="color: #B452CD">100</span>)
+<span style="color: #228B22"># and then the design matrix X including the intercept</span>
+<span style="color: #228B22">#  The design matrix now as function of a fourth-order polynomial</span>
+X = np.zeros((<span style="color: #658b00">len</span>(x),<span style="color: #B452CD">5</span>))
+X[:,<span style="color: #B452CD">0</span>] = <span style="color: #B452CD">1.0</span>
+X[:,<span style="color: #B452CD">1</span>] = x
+X[:,<span style="color: #B452CD">2</span>] = x**<span style="color: #B452CD">2</span>
+X[:,<span style="color: #B452CD">3</span>] = x**<span style="color: #B452CD">3</span>
+X[:,<span style="color: #B452CD">4</span>] = x**<span style="color: #B452CD">4</span>
 beta = (np.linalg.inv(X.T @ X) @ X.T ) @ y
 <span style="color: #228B22"># and then make the prediction</span>
 ytilde = X @ beta
@@ -838,7 +838,18 @@ $$
      <div class="input_area">
        <div class="highlight" style="background: #eeeedd">
  <pre style="line-height: 125%;"><span style="color: #228B22"># matrix inversion to find beta</span>
+<span style="color: #228B22"># First we set up the data</span>
 <span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
+x = np.random.rand(<span style="color: #B452CD">100</span>)
+y = <span style="color: #B452CD">2.0</span>+<span style="color: #B452CD">5</span>*x*x+<span style="color: #B452CD">0.1</span>*np.random.randn(<span style="color: #B452CD">100</span>)
+<span style="color: #228B22"># and then the design matrix X including the intercept</span>
+<span style="color: #228B22">#  The design matrix now as function of a fourth-order polynomial</span>
+X = np.zeros((<span style="color: #658b00">len</span>(x),<span style="color: #B452CD">5</span>))
+X[:,<span style="color: #B452CD">0</span>] = <span style="color: #B452CD">1.0</span>
+X[:,<span style="color: #B452CD">1</span>] = x
+X[:,<span style="color: #B452CD">2</span>] = x**<span style="color: #B452CD">2</span>
+X[:,<span style="color: #B452CD">3</span>] = x**<span style="color: #B452CD">3</span>
+X[:,<span style="color: #B452CD">4</span>] = x**<span style="color: #B452CD">4</span>
 beta = (np.linalg.inv(X.T @ X) @ X.T ) @ y
 <span style="color: #228B22"># and then make the prediction</span>
 ytilde = X @ beta
@@ -915,7 +915,18 @@ $$
      <div class="input_area">
        <div class="highlight" style="background: #f8f8f8">
  <pre style="line-height: 125%;"><span style="color: #408080; font-style: italic"># matrix inversion to find beta</span>
+<span style="color: #408080; font-style: italic"># First we set up the data</span>
 <span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
+x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(<span style="color: #666666">100</span>)
+y <span style="color: #666666">=</span> <span style="color: #666666">2.0+5*</span>x<span style="color: #666666">*</span>x<span style="color: #666666">+0.1*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(<span style="color: #666666">100</span>)
+<span style="color: #408080; font-style: italic"># and then the design matrix X including the intercept</span>
+<span style="color: #408080; font-style: italic">#  The design matrix now as function of a fourth-order polynomial</span>
+X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros((<span style="color: #008000">len</span>(x),<span style="color: #666666">5</span>))
+X[:,<span style="color: #666666">0</span>] <span style="color: #666666">=</span> <span style="color: #666666">1.0</span>
+X[:,<span style="color: #666666">1</span>] <span style="color: #666666">=</span> x
+X[:,<span style="color: #666666">2</span>] <span style="color: #666666">=</span> x<span style="color: #666666">**2</span>
+X[:,<span style="color: #666666">3</span>] <span style="color: #666666">=</span> x<span style="color: #666666">**3</span>
+X[:,<span style="color: #666666">4</span>] <span style="color: #666666">=</span> x<span style="color: #666666">**4</span>
 beta <span style="color: #666666">=</span> (np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X) <span style="color: #666666">@</span> X<span style="color: #666666">.</span>T ) <span style="color: #666666">@</span> y
 <span style="color: #408080; font-style: italic"># and then make the prediction</span>
 ytilde <span style="color: #666666">=</span> X <span style="color: #666666">@</span> beta
@@ -583,7 +583,18 @@ $n\times p$ matrix $\bm{X}$.
 It is rather straightforward to implement the matrix inversion and obtain the parameters $\bm{\beta}$. After having defined the matrix $\bm{X}$ and the outputs $\bm{y}$ we have 
 !bc pycod 
 # matrix inversion to find beta
+# First we set up the data
 import numpy as np
+x = np.random.rand(100)
+y = 2.0+5*x*x+0.1*np.random.randn(100)
+# and then the design matrix X including the intercept
+#  The design matrix now as function of a fourth-order polynomial
+X = np.zeros((len(x),5))
+X[:,0] = 1.0
+X[:,1] = x
+X[:,2] = x**2
+X[:,3] = x**3
+X[:,4] = x**4
 beta = (np.linalg.inv(X.T @ X) @ X.T ) @ y
 # and then make the prediction
 ytilde = X @ beta