lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update book

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2023-09-14 15:27:57 +02:00
parent b5676b9e83
commit 59449d4e76
16 changed files with 3287 additions and 620 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/LectureNotes/_build/.doctrees/environment.pickle
BIN
View File
Binary file not shown.
doc/LectureNotes/_build/.doctrees/week37.doctree
BIN
View File
Binary file not shown.
doc/LectureNotes/_build/html/_images/week37_144_0.png
BIN
View File
Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 22 KiB

After

Width:  |  Height:  |  Size: 22 KiB

doc/LectureNotes/_build/html/_images/week37_162_3.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 29 KiB

doc/LectureNotes/_build/html/_images/week37_176_1.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 33 KiB

doc/LectureNotes/_build/html/_images/week37_208_1.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 33 KiB

doc/LectureNotes/_build/html/_sources/week37.ipynb
+814 -188
View File
File diff suppressed because it is too large Load Diff
doc/LectureNotes/_build/html/searchindex.js
+1 -1
View File
File diff suppressed because one or more lines are too long
doc/LectureNotes/_build/html/week37.html
+397 -27
View File
@@ -642,6 +642,11 @@ const thebe_selector_output = ".output, .cell_output"
The same example but now with cross-validation
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#notes-on-scaling-with-examples">
Notes on scaling with examples
</a>
</li>
</ul>
</nav>
@@ -933,6 +938,11 @@ const thebe_selector_output = ".output, .cell_output"
The same example but now with cross-validation
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#notes-on-scaling-with-examples">
Notes on scaling with examples
</a>
</li>
</ul>
</nav>
@@ -947,7 +957,7 @@ doconce format html week37.do.txt --no_mako -->
<!-- dom:TITLE: Week 37: Statitsitcal interpretations and Resampling Methods --><div class="tex2jax_ignore mathjax_ignore section" id="week-37-statitsitcal-interpretations-and-resampling-methods">
<h1>Week 37: Statitsitcal interpretations and Resampling Methods<a class="headerlink" href="#week-37-statitsitcal-interpretations-and-resampling-methods" title="Permalink to this headline"></a></h1>
<p><strong>Morten Hjorth-Jensen</strong>, Department of Physics, University of Oslo and Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University</p>
<p>Date: <strong>Sep 11, 2023</strong></p>
<p>Date: <strong>Sep 14, 2023</strong></p>
<p>Copyright 1999-2023, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license</p>
<!-- todo add link to videos and add link to Van Wieringens notes --><div class="section" id="plans-for-week-37">
<h2>Plans for week 37<a class="headerlink" href="#plans-for-week-37" title="Permalink to this headline"></a></h2>
@@ -956,17 +966,20 @@ doconce format html week37.do.txt --no_mako -->
<li><p>Lecture from last week on calculations of expectation values</p></li>
<li><p>Exercise for week 37</p></li>
<li><p>Work on project 1</p></li>
<li><p>See also additional note on scaling (jupyter-notebook) sent separately. This will be discussed during the first hour of each session.</p></li>
<li><p>See also additional note on scaling (jupyter-notebook) sent separately. This will be discussed during the first hour of each session. This note is added at the end of these slides.</p></li>
<li><p>For more discussions of Ridge regression and calculation of averages, <a class="reference external" href="https://arxiv.org/abs/1509.09169">Wessel van Wieringen’s</a> article is highly recommended.</p></li>
</ul>
<p><strong>Material for the lecture on Thursday September 7.</strong></p>
<ul class="simple">
<li><p>Statistical interpretation of Ridge and Lasso regression</p></li>
<li><p><a class="reference external" href="https://youtu.be/YOBBr_toYxc">Video of Lecture</a></p></li>
<li><p><a class="reference external" href="https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2023/NotesSep14.pdf">Whiteboard notes</a></p></li>
<li><p>Resampling techniques, Bootstrap and cross validation and bias-variance tradeoff</p></li>
<li><p>Reads and Videos:</p>
<li><p>Statistical interpretation of Ridge and Lasso regression</p></li>
<li><p>Readings and Videos:</p>
<ul>
<li><p>Hastie et al Chapter 7, here we recommend 7.1-7.5 and 7.10 (cross-validation) and 7.11 (bootstrap).</p></li>
<li><p><a class="reference external" href="https://www.youtube.com/watch?v=fSytzGwwBVw">Video on cross validation</a></p></li>
<li><p><a class="reference external" href="https://www.youtube.com/watch?v=Xz0x-8-cgaQ">Video on Bootstrapping</a></p></li>
<li><p><a class="reference external" href="https://www.youtube.com/watch?v=EuBBz3bI-aA">Video on bias-variance tradeoff</a></p></li>
</ul>
</li>
@@ -1073,13 +1086,13 @@ mean value <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbo
\begin{eqnarray*}
\mbox{Var}(\boldsymbol{\hat{\beta}}) &amp; = &amp; \mathbb{E} \{ [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})] [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})]^{T} \}
\\
&amp; = &amp; \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}]^{T} \}
&amp; = &amp; \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{y} - \boldsymbol{\beta}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{y} - \boldsymbol{\beta}]^{T} \}
\\
% &amp; = &amp; \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}]^{T} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
% &amp; = &amp; \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{y}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{y}]^{T} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
% \\
% &amp; = &amp; \mathbb{E} \{ (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} \, \mathbf{Y}^{T} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
% &amp; = &amp; \mathbb{E} \{ (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{y} \, \mathbf{y}^{T} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
% \\
&amp; = &amp; (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \mathbb{E} \{ \mathbf{Y} \, \mathbf{Y}^{T} \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
&amp; = &amp; (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \mathbb{E} \{ \mathbf{y} \, \mathbf{y}^{T} \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
\\
&amp; = &amp; (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \{ \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} + \sigma^2 \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
% \\
@@ -1091,7 +1104,7 @@ mean value <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbo
\, \, \, = \, \, \, \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1},
\end{eqnarray*}
\end{split}\]</div>
<p>where we have used that <span class="math notranslate nohighlight">\(\mathbb{E} (\mathbf{Y} \mathbf{Y}^{T}) =
<p>where we have used that <span class="math notranslate nohighlight">\(\mathbb{E} (\mathbf{y} \mathbf{y}^{T}) =
\mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} +
\sigma^2 \, \mathbf{I}_{nn}\)</span>. From <span class="math notranslate nohighlight">\(\mbox{Var}(\boldsymbol{\beta}) = \sigma^2
\, (\mathbf{X}^{T} \mathbf{X})^{-1}\)</span>, one obtains an estimate of the
@@ -1104,20 +1117,20 @@ when we employ Ridge regression, allowing us again to define a confidence interv
<p>It is rather straightforward to show that</p>
<div class="math notranslate nohighlight">
\[
\mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\boldsymbol{\beta}^{\mathrm{OLS}}.
\mathbb{E} \big[ \hat{\boldsymbol{\beta}}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\boldsymbol{\beta}.
\]</div>
<p>We see clearly that
<span class="math notranslate nohighlight">\(\mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big] \not= \boldsymbol{\beta}^{\mathrm{OLS}}\)</span> for any <span class="math notranslate nohighlight">\(\lambda &gt; 0\)</span>.</p>
<span class="math notranslate nohighlight">\(\mathbb{E} \big[ \hat{\boldsymbol{\beta}}^{\mathrm{Ridge}} \big] \not= \hat{\boldsymbol{\beta}}^{\mathrm{OLS}}\)</span> for any <span class="math notranslate nohighlight">\(\lambda &gt; 0\)</span>.</p>
<p>We can also compute the variance as</p>
<div class="math notranslate nohighlight">
\[
\mbox{Var}[\boldsymbol{\beta}^{\mathrm{Ridge}}]=\sigma^2[ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1} \mathbf{X}^{T} \mathbf{X} \{ [ \mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T},
\mbox{Var}[\hat{\boldsymbol{\beta}}^{\mathrm{Ridge}}]=\sigma^2[ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1} \mathbf{X}^{T} \mathbf{X} \{ [ \mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T},
\]</div>
<p>and it is easy to see that if the parameter <span class="math notranslate nohighlight">\(\lambda\)</span> goes to infinity then the variance of Ridge parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> goes to zero.</p>
<p>With this, we can compute the difference</p>
<div class="math notranslate nohighlight">
\[
\mbox{Var}[\boldsymbol{\beta}^{\mathrm{OLS}}]-\mbox{Var}(\boldsymbol{\beta}^{\mathrm{Ridge}})=\sigma^2 [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}[ 2\lambda\mathbf{I} + \lambda^2 (\mathbf{X}^{T} \mathbf{X})^{-1} ] \{ [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}.
\mbox{Var}[\hat{\boldsymbol{\beta}}^{\mathrm{OLS}}]-\mbox{Var}(\hat{\boldsymbol{\beta}}^{\mathrm{Ridge}})=\sigma^2 [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}[ 2\lambda\mathbf{I} + \lambda^2 (\mathbf{X}^{T} \mathbf{X})^{-1} ] \{ [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}.
\]</div>
<p>The difference is non-negative definite since each component of the
matrix product is non-negative definite.
@@ -1754,7 +1767,7 @@ theorem.</p>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Bootstrap Statistics :
original bias std. error
100.232 14.9426 100.231 0.148895
99.9713 14.943 99.972 0.149583
</pre></div>
</div>
</div>
@@ -2021,7 +2034,9 @@ Error: 0.021592704588021178
Bias^2: 0.010516485576646504
Var: 0.01107621901137467
0.021592704588021178 &gt;= 0.010516485576646504 + 0.01107621901137467 = 0.021592704588021174
Polynomial degree: 11
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Polynomial degree: 11
Error: 0.07160048164232538
Bias^2: 0.014436800088896381
Var: 0.05716368155342902
@@ -2038,7 +2053,7 @@ Var: 0.20867052175003364
0.2284246870217162 &gt;= 0.01975416527168255 + 0.20867052175003364 = 0.2284246870217162
</pre></div>
</div>
<img alt="_images/week37_162_2.png" src="_images/week37_162_2.png" />
<img alt="_images/week37_162_3.png" src="_images/week37_162_3.png" />
</div>
</div>
</div>
@@ -2458,12 +2473,12 @@ Mean squared error on test data: 238.16356503
Degree of polynomial: 20
Mean squared error on training data: 0.00140849
Mean squared error on test data: 1345.68592431
Degree of polynomial: 21
Mean squared error on training data: 0.00119699
Mean squared error on test data: 1836.21110005
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 22
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 21
Mean squared error on training data: 0.00119699
Mean squared error on test data: 1836.21110005
Degree of polynomial: 22
Mean squared error on training data: 0.00092904
Mean squared error on test data: 1182.64316482
Degree of polynomial: 23
@@ -2478,12 +2493,12 @@ Mean squared error on test data: 7697.35412147
Degree of polynomial: 26
Mean squared error on training data: 0.00075597
Mean squared error on test data: 1078.81597834
Degree of polynomial: 27
Mean squared error on training data: 0.00068088
Mean squared error on test data: 3189.20355156
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 28
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 27
Mean squared error on training data: 0.00068088
Mean squared error on test data: 3189.20355156
Degree of polynomial: 28
Mean squared error on training data: 0.00063364
Mean squared error on test data: 692.24085321
Degree of polynomial: 29
@@ -2491,9 +2506,9 @@ Mean squared error on training data: 0.00063862
Mean squared error on test data: 3073.63180447
</pre></div>
</div>
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/var/folders/td/3yk470mj5p931p9dtkk0y6jw0000gn/T/ipykernel_7768/626635268.py:73: RuntimeWarning: divide by zero encountered in log10
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/var/folders/td/3yk470mj5p931p9dtkk0y6jw0000gn/T/ipykernel_31560/626635268.py:73: RuntimeWarning: divide by zero encountered in log10
plt.plot(polynomial, np.log10(trainingerror), label=&#39;Training Error&#39;)
/var/folders/td/3yk470mj5p931p9dtkk0y6jw0000gn/T/ipykernel_7768/626635268.py:74: RuntimeWarning: divide by zero encountered in log10
/var/folders/td/3yk470mj5p931p9dtkk0y6jw0000gn/T/ipykernel_31560/626635268.py:74: RuntimeWarning: divide by zero encountered in log10
plt.plot(polynomial, np.log10(testerror), label=&#39;Test Error&#39;)
</pre></div>
</div>
@@ -2578,7 +2593,7 @@ Mean squared error on test data: 3073.63180447
</div>
</div>
<div class="cell_output docutils container">
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/var/folders/td/3yk470mj5p931p9dtkk0y6jw0000gn/T/ipykernel_7768/3817475779.py:63: RuntimeWarning: divide by zero encountered in log10
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/var/folders/td/3yk470mj5p931p9dtkk0y6jw0000gn/T/ipykernel_31560/3817475779.py:63: RuntimeWarning: divide by zero encountered in log10
plt.plot(polynomial, np.log10(estimated_mse_sklearn), label=&#39;Test Error&#39;)
</pre></div>
</div>
@@ -2586,6 +2601,361 @@ Mean squared error on test data: 3073.63180447
</div>
</div>
</div>
<div class="section" id="notes-on-scaling-with-examples">
<h2>Notes on scaling with examples<a class="headerlink" href="#notes-on-scaling-with-examples" title="Permalink to this headline"></a></h2>
<p>The programs here use both ordinrary least squares (OLS) and Ridge
regression with one value only for the hyperparameter <span class="math notranslate nohighlight">\(\lambda\)</span>. The
first example has no scaling and includes the intercept as well and we
are trying to fit a second-order polynomial. The second code takes out
the intercept and subtracts the mean values of each column of the
design matrix and the mean value of the outputs.</p>
<p>The third and final code uses <strong>Scikit-Learn</strong> as library in order to
calculate the optimal parameters for OLS and Ridge regression. Note
that it is highly recommended to not include the intercept in Ridge
and Lasso regression, in order to avoid penalizing the optimization by
the intercept. The second and third codes do thus not include the
intercept. In the second code we do the scaling ourselves while the
last code uses the standard scaler option included in <strong>Scikit-Learn</strong>, known as centering (where
we subtract the mean values).</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">sklearn.linear_model</span> <span class="kn">import</span> <span class="n">LinearRegression</span>
<span class="kn">from</span> <span class="nn">sklearn.preprocessing</span> <span class="kn">import</span> <span class="n">PolynomialFeatures</span>
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">train_test_split</span>
<span class="kn">from</span> <span class="nn">sklearn.preprocessing</span> <span class="kn">import</span> <span class="n">StandardScaler</span>
<span class="k">def</span> <span class="nf">MSE</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span><span class="n">y_model</span><span class="p">):</span>
<span class="n">n</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">size</span><span class="p">(</span><span class="n">y_model</span><span class="p">)</span>
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span><span class="o">-</span><span class="n">y_model</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">/</span><span class="n">n</span>
<span class="k">def</span> <span class="nf">OLS_fit_beta</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">):</span>
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">pinv</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="o">@</span> <span class="n">y</span>
<span class="k">def</span> <span class="nf">Ridge_fit_beta</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span><span class="n">L</span><span class="p">,</span><span class="n">d</span><span class="p">):</span>
<span class="n">I</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="n">d</span><span class="p">,</span><span class="n">d</span><span class="p">)</span>
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">pinv</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="o">+</span> <span class="n">L</span><span class="o">*</span><span class="n">I</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="o">@</span> <span class="n">y</span>
<span class="c1"># Same random numbers for each test.</span>
<span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">seed</span><span class="p">(</span><span class="mi">2018</span><span class="p">)</span>
<span class="n">n</span> <span class="o">=</span> <span class="mi">100</span>
<span class="n">d</span> <span class="o">=</span> <span class="mi">3</span>
<span class="c1"># hyperparameter lambda</span>
<span class="n">Lambda</span> <span class="o">=</span> <span class="mf">0.01</span>
<span class="c1"># Make data set, simple second-order polynomial</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="n">n</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="mf">0.5</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="mf">5.0</span><span class="o">*</span><span class="p">(</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)</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="n">n</span><span class="p">)</span>
<span class="c1"># The design matrix X includes the intercept and no scaling is made</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="n">d</span><span class="p">))</span>
<span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">d</span><span class="p">):</span>
<span class="n">X</span><span class="p">[:,</span> <span class="n">p</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span> <span class="o">**</span> <span class="p">(</span><span class="n">p</span><span class="p">)</span>
<span class="c1">#Split data, no scaling is used and we include the intercept</span>
<span class="n">X_train</span><span class="p">,</span> <span class="n">X_test</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span> <span class="n">y_test</span> <span class="o">=</span> <span class="n">train_test_split</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">test_size</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
<span class="c1">#Calculate beta, own code</span>
<span class="n">beta_OLS</span> <span class="o">=</span> <span class="n">OLS_fit_beta</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span> <span class="n">y_train</span><span class="p">)</span>
<span class="n">beta_Ridge</span> <span class="o">=</span> <span class="n">Ridge_fit_beta</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span><span class="n">Lambda</span><span class="p">,</span><span class="n">d</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">beta_OLS</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">beta_Ridge</span><span class="p">)</span>
<span class="c1">#predict value</span>
<span class="n">ytilde_test_OLS</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">@</span> <span class="n">beta_OLS</span>
<span class="n">ytilde_test_Ridge</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">@</span> <span class="n">beta_Ridge</span>
<span class="c1">#Calculate MSE</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot; &quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;test MSE of OLS:&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ytilde_test_OLS</span><span class="p">))</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot; &quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;test MSE of Ridge&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ytilde_test_Ridge</span><span class="p">))</span>
<span class="n">plt</span><span class="o">.</span><span class="n">scatter</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">label</span><span class="o">=</span><span class="s1">&#39;Data&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta_OLS</span><span class="p">,</span><span class="s1">&#39;*&#39;</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">&quot;OLS_Fit&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta_Ridge</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">&quot;Ridge_Fit&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">grid</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[1.79934087 0.47179152 5.01549939]
[1.79909592 0.47176716 5.01550546]
test MSE of OLS:
1.13943111290393
test MSE of Ridge
1.1395235273363686
</pre></div>
</div>
<img alt="_images/week37_176_1.png" src="_images/week37_176_1.png" />
</div>
</div>
<p>In this example we do not include the intercept and we scale the data by subtracting the mean values. This follows the discussion in the <a class="reference external" href="https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/chapter3.html#more-on-rescaling-data">lecture material</a>.
see also the weekly slides <a class="reference external" href="https://compphysics.github.io/MachineLearning/doc/pub/week36/html/._week36-bs029.html">for week 36</a>.
It is recommended whrn we use Ridge and Lasso regression to not include the intercept in the optimization process.</p>
<p>Before we discuss the code, we repeat some of the basic math from the slides of week 36.</p>
<p>Let us try to understand what this may imply mathematically when we
subtract the mean values, also known as <em>zero centering</em> or simply <em>centering</em>. For
simplicity, we will focus on ordinary regression, as done in the above example.</p>
<p>The cost/loss function for regression is</p>
<div class="math notranslate nohighlight">
\[
C(\beta_0, \beta_1, ... , \beta_{p-1}) = \frac{1}{n}\sum_{i=0}^{n} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij}\beta_j\right)^2,.
\]</div>
<p>Recall also that we use the squared value. This expression can lead to an
increased penalty for higher differences between predicted and
output/target values.</p>
<p>What we have done is to single out the <span class="math notranslate nohighlight">\(\beta_0\)</span> term in the
definition of the mean squared error (MSE). The design matrix <span class="math notranslate nohighlight">\(X\)</span>
does in this case not contain any intercept column. When we take the
derivative with respect to <span class="math notranslate nohighlight">\(\beta_0\)</span>, we want the derivative to obey</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial C}{\partial \beta_j} = 0,
\]</div>
<p>for all <span class="math notranslate nohighlight">\(j\)</span>. For <span class="math notranslate nohighlight">\(\beta_0\)</span> we have</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial C}{\partial \beta_0} = -\frac{2}{n}\sum_{i=0}^{n-1} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij} \beta_j\right).
\]</div>
<p>Multiplying away the constant <span class="math notranslate nohighlight">\(2/n\)</span>, we obtain</p>
<div class="math notranslate nohighlight">
\[
\sum_{i=0}^{n-1} \beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} \sum_{j=1}^{p-1} X_{ij} \beta_j.
\]</div>
<p>Let us specialize first to the case where we have only two parameters <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span>.
Our result for <span class="math notranslate nohighlight">\(\beta_0\)</span> simplifies then to</p>
<div class="math notranslate nohighlight">
\[
n\beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \beta_1.
\]</div>
<p>We obtain then</p>
<div class="math notranslate nohighlight">
\[
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \beta_1\frac{1}{n}\sum_{i=0}^{n-1} X_{i1}.
\]</div>
<p>If we define</p>
<div class="math notranslate nohighlight">
\[
\mu_{\boldsymbol{x}_1}=\frac{1}{n}\sum_{i=0}^{n-1} X_{i1},
\]</div>
<p>and the mean value of the outputs as</p>
<div class="math notranslate nohighlight">
\[
\mu_y=\frac{1}{n}\sum_{i=0}^{n-1}y_i,
\]</div>
<p>we have</p>
<div class="math notranslate nohighlight">
\[
\beta_0 = \mu_y - \beta_1\mu_{\boldsymbol{x}_1}.
\]</div>
<p>In the general case with more parameters than <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span>, we have</p>
<div class="math notranslate nohighlight">
\[
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \frac{1}{n}\sum_{i=0}^{n-1}\sum_{j=1}^{p-1} X_{ij}\beta_j.
\]</div>
<p>We can rewrite the latter equation as</p>
<div class="math notranslate nohighlight">
\[
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \sum_{j=1}^{p-1} \mu_{\boldsymbol{x}_j}\beta_j,
\]</div>
<p>where we have defined</p>
<div class="math notranslate nohighlight">
\[
\mu_{\boldsymbol{x}_j}=\frac{1}{n}\sum_{i=0}^{n-1} X_{ij},
\]</div>
<p>the mean value for all elements of the column vector <span class="math notranslate nohighlight">\(\boldsymbol{x}_j\)</span>.</p>
<p>Replacing <span class="math notranslate nohighlight">\(y_i\)</span> with <span class="math notranslate nohighlight">\(y_i - y_i - \overline{\boldsymbol{y}}\)</span> and centering also our design matrix results in a cost function (in vector-matrix disguise)</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{\beta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta}).
\]</div>
<p>If we minimize with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> we have then</p>
<div class="math notranslate nohighlight">
\[
\hat{\boldsymbol{\beta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}},
\]</div>
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\boldsymbol{y}}\)</span>
and <span class="math notranslate nohighlight">\(\tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj}\)</span>.</p>
<p>For Ridge regression we need to add <span class="math notranslate nohighlight">\(\lambda \boldsymbol{\beta}^T\boldsymbol{\beta}\)</span> to the cost function and get then</p>
<div class="math notranslate nohighlight">
\[
\hat{\boldsymbol{\beta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}.
\]</div>
<p>Now we try to implement this.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">seed</span><span class="p">(</span><span class="mi">2018</span><span class="p">)</span>
<span class="n">n</span> <span class="o">=</span> <span class="mi">100</span>
<span class="c1"># we do not include the intercept</span>
<span class="n">d</span> <span class="o">=</span> <span class="mi">2</span>
<span class="n">Lambda</span> <span class="o">=</span> <span class="mf">0.01</span>
<span class="c1"># Make data set.</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="n">n</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="mf">0.5</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="mf">5.0</span><span class="o">*</span><span class="p">(</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)</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="n">n</span><span class="p">)</span>
<span class="c1">#Design matrix X does not include the intercept. </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="n">d</span><span class="p">))</span>
<span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">d</span><span class="p">):</span>
<span class="n">X</span><span class="p">[:,</span> <span class="n">p</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span> <span class="o">**</span> <span class="p">(</span><span class="n">p</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
<span class="c1">#Split data in train and test</span>
<span class="n">X_train</span><span class="p">,</span> <span class="n">X_test</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span> <span class="n">y_test</span> <span class="o">=</span> <span class="n">train_test_split</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">test_size</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
<span class="c1"># Scale data by subtracting mean value,own implementation</span>
<span class="c1">#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable</span>
<span class="n">X_train_mean</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="c1">#Center by removing mean from each feature</span>
<span class="n">X_train_scaled</span> <span class="o">=</span> <span class="n">X_train</span> <span class="o">-</span> <span class="n">X_train_mean</span>
<span class="n">X_test_scaled</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">-</span> <span class="n">X_train_mean</span>
<span class="c1">#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note)</span>
<span class="n">y_scaler</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_train</span><span class="p">)</span>
<span class="n">y_train_scaled</span> <span class="o">=</span> <span class="n">y_train</span> <span class="o">-</span> <span class="n">y_scaler</span>
<span class="c1">#Calculate beta</span>
<span class="n">beta_OLS</span> <span class="o">=</span> <span class="n">OLS_fit_beta</span><span class="p">(</span><span class="n">X_train_scaled</span><span class="p">,</span> <span class="n">y_train_scaled</span><span class="p">)</span>
<span class="n">beta_Ridge</span> <span class="o">=</span> <span class="n">Ridge_fit_beta</span><span class="p">(</span><span class="n">X_train_scaled</span><span class="p">,</span> <span class="n">y_train_scaled</span><span class="p">,</span><span class="n">Lambda</span><span class="p">,</span><span class="n">d</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">beta_OLS</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">beta_Ridge</span><span class="p">)</span>
<span class="c1"># calculate intercepts and print them</span>
<span class="n">interceptOLS</span> <span class="o">=</span> <span class="n">y_scaler</span> <span class="o">-</span> <span class="n">X_train_mean</span> <span class="o">@</span> <span class="n">beta_OLS</span>
<span class="n">interceptRidge</span> <span class="o">=</span> <span class="n">y_scaler</span> <span class="o">-</span> <span class="n">X_train_mean</span> <span class="o">@</span> <span class="n">beta_Ridge</span>
<span class="nb">print</span><span class="p">(</span><span class="n">interceptOLS</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">interceptRidge</span><span class="p">)</span>
<span class="c1">#predict value with intercept</span>
<span class="n">ytilde_test_OLS</span> <span class="o">=</span> <span class="n">X_test_scaled</span> <span class="o">@</span> <span class="n">beta_OLS</span><span class="o">+</span><span class="n">y_scaler</span>
<span class="n">ytilde_test_Ridge</span> <span class="o">=</span> <span class="n">X_test_scaled</span> <span class="o">@</span> <span class="n">beta_Ridge</span><span class="o">+</span><span class="n">y_scaler</span>
<span class="c1">#Calculate MSE</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot; &quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;test MSE of OLS:&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ytilde_test_OLS</span><span class="p">))</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot; &quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;test MSE of Ridge&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ytilde_test_Ridge</span><span class="p">))</span>
<span class="n">plt</span><span class="o">.</span><span class="n">scatter</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">label</span><span class="o">=</span><span class="s1">&#39;Data&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta_OLS</span><span class="o">+</span><span class="n">interceptOLS</span><span class="p">,</span><span class="s1">&#39;*&#39;</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">&quot;OLS_Fit&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta_Ridge</span><span class="o">+</span><span class="n">interceptRidge</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">&quot;Ridge_Fit&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">grid</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[0.47179152 5.01549939]
[0.47176783 5.01542292]
1.7993408651198877
1.7995707762668065
test MSE of OLS:
1.1394311129039245
test MSE of Ridge
1.1395084586525954
</pre></div>
</div>
<img alt="_images/week37_208_1.png" src="_images/week37_208_1.png" />
</div>
</div>
<p>Finally, instead of using our own function we repeat the same example
using the <strong>standardscaler</strong> functionality of the library
<strong>Scikit-Learn</strong>. Here we limit ourselves to Ridge regression only.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">from</span> <span class="nn">sklearn</span> <span class="kn">import</span> <span class="n">linear_model</span>
<span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">seed</span><span class="p">(</span><span class="mi">2018</span><span class="p">)</span>
<span class="n">n</span> <span class="o">=</span> <span class="mi">10</span>
<span class="n">d</span> <span class="o">=</span> <span class="mi">2</span>
<span class="n">Lambda</span> <span class="o">=</span> <span class="mf">0.01</span>
<span class="c1"># Make data set.</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="n">n</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="mf">0.5</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="mf">5.0</span><span class="o">*</span><span class="p">(</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)</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="n">n</span><span class="p">)</span>
<span class="c1"># Design matrix X does not include the intercept. </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="n">n</span><span class="p">,</span> <span class="n">d</span><span class="p">))</span>
<span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">d</span><span class="p">):</span>
<span class="n">X</span><span class="p">[:,</span> <span class="n">p</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span> <span class="o">**</span> <span class="p">(</span><span class="n">p</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
<span class="c1">#Split data in train and test</span>
<span class="n">X_train</span><span class="p">,</span> <span class="n">X_test</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span> <span class="n">y_test</span> <span class="o">=</span> <span class="n">train_test_split</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">test_size</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
<span class="c1"># Scale data by subtracting mean value of the input using scikit-learn</span>
<span class="n">scaler</span> <span class="o">=</span> <span class="n">StandardScaler</span><span class="p">(</span><span class="n">with_std</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
<span class="n">scaler</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X_train</span><span class="p">)</span>
<span class="n">X_train_mean</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="n">X_train_scaled</span> <span class="o">=</span> <span class="n">scaler</span><span class="o">.</span><span class="n">transform</span><span class="p">(</span><span class="n">X_train</span><span class="p">)</span>
<span class="n">X_test_scaled</span> <span class="o">=</span> <span class="n">scaler</span><span class="o">.</span><span class="n">transform</span><span class="p">(</span><span class="n">X_test</span><span class="p">)</span>
<span class="c1"># We scale also the output, here by our own code</span>
<span class="n">y_scaler</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_train</span><span class="p">)</span>
<span class="n">y_train_scaled</span> <span class="o">=</span> <span class="n">y_train</span> <span class="o">-</span> <span class="n">y_scaler</span>
<span class="n">y_test_scaled</span> <span class="o">=</span> <span class="n">y_test</span><span class="o">-</span> <span class="n">y_scaler</span>
<span class="c1">#Calculate beta</span>
<span class="n">OLS</span> <span class="o">=</span> <span class="n">LinearRegression</span><span class="p">()</span>
<span class="n">betaOLS</span><span class="o">=</span><span class="n">OLS</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X_train_scaled</span><span class="p">,</span><span class="n">y_train_scaled</span><span class="p">)</span>
<span class="n">ypredictOLS</span> <span class="o">=</span> <span class="n">OLS</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test_scaled</span><span class="p">)</span>
<span class="n">linear_model</span><span class="o">.</span><span class="n">Ridge</span><span class="p">(</span><span class="n">Lambda</span><span class="p">)</span>
<span class="n">RegRidge</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X_train_scaled</span><span class="p">,</span><span class="n">y_train_scaled</span><span class="p">)</span>
<span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test_scaled</span><span class="p">)</span>
<span class="n">betaOLS</span> <span class="o">=</span> <span class="n">OLS</span><span class="o">.</span><span class="n">coef_</span>
<span class="n">betaRidge</span> <span class="o">=</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">coef_</span>
<span class="nb">print</span><span class="p">(</span><span class="n">betaOLS</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">betaRidge</span><span class="p">)</span>
<span class="n">interceptOLS</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_train</span><span class="p">)</span> <span class="o">-</span> <span class="n">X_train_mean</span> <span class="o">@</span> <span class="n">betaOLS</span>
<span class="n">interceptRidge</span> <span class="o">=</span> <span class="n">y_scaler</span> <span class="o">-</span> <span class="n">X_train_mean</span> <span class="o">@</span> <span class="n">betaRidge</span>
<span class="nb">print</span><span class="p">(</span><span class="n">interceptOLS</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">interceptRidge</span><span class="p">)</span>
<span class="c1">#predict value </span>
<span class="n">ytilde_test_Ridge</span> <span class="o">=</span> <span class="n">X_test_scaled</span> <span class="o">@</span> <span class="n">betaRidge</span><span class="o">+</span><span class="n">y_scaler</span>
<span class="n">ytilde_test_OLS</span> <span class="o">=</span> <span class="n">X_test_scaled</span> <span class="o">@</span> <span class="n">betaOLS</span><span class="o">+</span><span class="n">y_scaler</span>
<span class="c1">#Calculate MSE</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot; &quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;test MSE of OLS&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ytilde_test_OLS</span><span class="p">))</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot; &quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;test MSE of Ridge&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ytilde_test_Ridge</span><span class="p">))</span>
<span class="n">plt</span><span class="o">.</span><span class="n">scatter</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">label</span><span class="o">=</span><span class="s1">&#39;Data&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">X</span> <span class="o">@</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">coef_</span> <span class="o">+</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">intercept_</span> <span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">&quot;Ridge_Fit&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">grid</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</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 [11],</span> in <span class="ni">&lt;cell line: 34&gt;</span><span class="nt">()</span>
<span class="g g-Whitespace"> </span><span class="mi">32</span> <span class="n">ypredictOLS</span> <span class="o">=</span> <span class="n">OLS</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test_scaled</span><span class="p">)</span>
<span class="g g-Whitespace"> </span><span class="mi">33</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Ridge</span><span class="p">(</span><span class="n">Lambda</span><span class="p">)</span>
<span class="ne">---&gt; </span><span class="mi">34</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X_train_scaled</span><span class="p">,</span><span class="n">y_train_scaled</span><span class="p">)</span>
<span class="g g-Whitespace"> </span><span class="mi">35</span> <span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test_scaled</span><span class="p">)</span>
<span class="g g-Whitespace"> </span><span class="mi">36</span> <span class="n">betaOLS</span> <span class="o">=</span> <span class="n">OLS</span><span class="o">.</span><span class="n">coef_</span>
<span class="ne">NameError</span>: name &#39;RegRidge&#39; is not defined
</pre></div>
</div>
</div>
</div>
</div>
</div>
<script type="text/x-thebe-config">
doc/LectureNotes/_build/jupyter_execute/week37.ipynb
+907 -203
View File
File diff suppressed because one or more lines are too long
doc/LectureNotes/_build/jupyter_execute/week37.py
+354 -13
View File
@@ -8,7 +8,7 @@
# # Week 37: Statitsitcal interpretations and Resampling Methods
# **Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University
#
# Date: **Sep 11, 2023**
# Date: **Sep 14, 2023**
#
# Copyright 1999-2023, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license
#
@@ -24,23 +24,29 @@
#
# * Work on project 1
#
# * See also additional note on scaling (jupyter-notebook) sent separately. This will be discussed during the first hour of each session.
# * See also additional note on scaling (jupyter-notebook) sent separately. This will be discussed during the first hour of each session. This note is added at the end of these slides.
#
# * For more discussions of Ridge regression and calculation of averages, [Wessel van Wieringen's](https://arxiv.org/abs/1509.09169) article is highly recommended.
#
#
# **Material for the lecture on Thursday September 7.**
#
# * Statistical interpretation of Ridge and Lasso regression
# * [Video of Lecture](https://youtu.be/YOBBr_toYxc)
#
# * [Whiteboard notes](https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2023/NotesSep14.pdf)
#
# * Resampling techniques, Bootstrap and cross validation and bias-variance tradeoff
#
# * Reads and Videos:
# * Statistical interpretation of Ridge and Lasso regression
#
# * Readings and Videos:
#
# * Hastie et al Chapter 7, here we recommend 7.1-7.5 and 7.10 (cross-validation) and 7.11 (bootstrap).
#
# * [Video on cross validation](https://www.youtube.com/watch?v=fSytzGwwBVw)
#
# * [Video on Bootstrapping](https://www.youtube.com/watch?v=Xz0x-8-cgaQ)
#
# * [Video on bias-variance tradeoff](https://www.youtube.com/watch?v=EuBBz3bI-aA)
# ## Material from last week and relevant for the weekly exercises
@@ -153,13 +159,13 @@
# \begin{eqnarray*}
# \mbox{Var}(\boldsymbol{\hat{\beta}}) & = & \mathbb{E} \{ [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})] [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})]^{T} \}
# \\
# & = & \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}]^{T} \}
# & = & \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{y} - \boldsymbol{\beta}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{y} - \boldsymbol{\beta}]^{T} \}
# \\
# % & = & \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}]^{T} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
# % & = & \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{y}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{y}]^{T} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
# % \\
# % & = & \mathbb{E} \{ (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} \, \mathbf{Y}^{T} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
# % & = & \mathbb{E} \{ (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{y} \, \mathbf{y}^{T} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
# % \\
# & = & (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \mathbb{E} \{ \mathbf{Y} \, \mathbf{Y}^{T} \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
# & = & (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \mathbb{E} \{ \mathbf{y} \, \mathbf{y}^{T} \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
# \\
# & = & (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \{ \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} + \sigma^2 \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
# % \\
@@ -172,7 +178,7 @@
# \end{eqnarray*}
# $$
# where we have used that $\mathbb{E} (\mathbf{Y} \mathbf{Y}^{T}) =
# where we have used that $\mathbb{E} (\mathbf{y} \mathbf{y}^{T}) =
# \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} +
# \sigma^2 \, \mathbf{I}_{nn}$. From $\mbox{Var}(\boldsymbol{\beta}) = \sigma^2
# \, (\mathbf{X}^{T} \mathbf{X})^{-1}$, one obtains an estimate of the
@@ -187,16 +193,16 @@
# It is rather straightforward to show that
# $$
# \mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\boldsymbol{\beta}^{\mathrm{OLS}}.
# \mathbb{E} \big[ \hat{\boldsymbol{\beta}}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\boldsymbol{\beta}.
# $$
# We see clearly that
# $\mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big] \not= \boldsymbol{\beta}^{\mathrm{OLS}}$ for any $\lambda > 0$.
# $\mathbb{E} \big[ \hat{\boldsymbol{\beta}}^{\mathrm{Ridge}} \big] \not= \hat{\boldsymbol{\beta}}^{\mathrm{OLS}}$ for any $\lambda > 0$.
#
# We can also compute the variance as
# $$
# \mbox{Var}[\boldsymbol{\beta}^{\mathrm{Ridge}}]=\sigma^2[ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1} \mathbf{X}^{T} \mathbf{X} \{ [ \mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T},
# \mbox{Var}[\hat{\boldsymbol{\beta}}^{\mathrm{Ridge}}]=\sigma^2[ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1} \mathbf{X}^{T} \mathbf{X} \{ [ \mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T},
# $$
# and it is easy to see that if the parameter $\lambda$ goes to infinity then the variance of Ridge parameters $\boldsymbol{\beta}$ goes to zero.
@@ -204,7 +210,7 @@
# With this, we can compute the difference
# $$
# \mbox{Var}[\boldsymbol{\beta}^{\mathrm{OLS}}]-\mbox{Var}(\boldsymbol{\beta}^{\mathrm{Ridge}})=\sigma^2 [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}[ 2\lambda\mathbf{I} + \lambda^2 (\mathbf{X}^{T} \mathbf{X})^{-1} ] \{ [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}.
# \mbox{Var}[\hat{\boldsymbol{\beta}}^{\mathrm{OLS}}]-\mbox{Var}(\hat{\boldsymbol{\beta}}^{\mathrm{Ridge}})=\sigma^2 [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}[ 2\lambda\mathbf{I} + \lambda^2 (\mathbf{X}^{T} \mathbf{X})^{-1} ] \{ [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}.
# $$
# The difference is non-negative definite since each component of the
@@ -1513,3 +1519,338 @@ plt.ylabel('log10[MSE]')
plt.legend()
plt.show()
# ## Notes on scaling with examples
#
# The programs here use both ordinrary least squares (OLS) and Ridge
# regression with one value only for the hyperparameter $\lambda$. The
# first example has no scaling and includes the intercept as well and we
# are trying to fit a second-order polynomial. The second code takes out
# the intercept and subtracts the mean values of each column of the
# design matrix and the mean value of the outputs.
#
# The third and final code uses **Scikit-Learn** as library in order to
# calculate the optimal parameters for OLS and Ridge regression. Note
# that it is highly recommended to not include the intercept in Ridge
# and Lasso regression, in order to avoid penalizing the optimization by
# the intercept. The second and third codes do thus not include the
# intercept. In the second code we do the scaling ourselves while the
# last code uses the standard scaler option included in **Scikit-Learn**, known as centering (where
# we subtract the mean values).
# In[9]:
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
def MSE(y_data,y_model):
n = np.size(y_model)
return np.sum((y_data-y_model)**2)/n
def OLS_fit_beta(X, y):
return np.linalg.pinv(X.T @ X) @ X.T @ y
def Ridge_fit_beta(X, y,L,d):
I = np.eye(d,d)
return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y
# Same random numbers for each test.
np.random.seed(2018)
n = 100
d = 3
# hyperparameter lambda
Lambda = 0.01
# Make data set, simple second-order polynomial
x = np.linspace(-3, 3, n)
y = 2.0 + 0.5*x + 5.0*(x**2)+ np.random.randn(n)
# The design matrix X includes the intercept and no scaling is made
X = np.zeros((len(x), d))
for p in range(d):
X[:, p] = x ** (p)
#Split data, no scaling is used and we include the intercept
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
#Calculate beta, own code
beta_OLS = OLS_fit_beta(X_train, y_train)
beta_Ridge = Ridge_fit_beta(X_train, y_train,Lambda,d)
print(beta_OLS)
print(beta_Ridge)
#predict value
ytilde_test_OLS = X_test @ beta_OLS
ytilde_test_Ridge = X_test @ beta_Ridge
#Calculate MSE
print(" ")
print("test MSE of OLS:")
print(MSE(y_test,ytilde_test_OLS))
print(" ")
print("test MSE of Ridge")
print(MSE(y_test,ytilde_test_Ridge))
plt.scatter(x,y,label='Data')
plt.plot(x, X @ beta_OLS,'*', label="OLS_Fit")
plt.plot(x, X @ beta_Ridge, label="Ridge_Fit")
plt.grid()
plt.legend()
plt.show()
# In this example we do not include the intercept and we scale the data by subtracting the mean values. This follows the discussion in the [lecture material](https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/chapter3.html#more-on-rescaling-data).
# see also the weekly slides [for week 36](https://compphysics.github.io/MachineLearning/doc/pub/week36/html/._week36-bs029.html).
# It is recommended whrn we use Ridge and Lasso regression to not include the intercept in the optimization process.
#
# Before we discuss the code, we repeat some of the basic math from the slides of week 36.
#
# Let us try to understand what this may imply mathematically when we
# subtract the mean values, also known as *zero centering* or simply *centering*. For
# simplicity, we will focus on ordinary regression, as done in the above example.
#
# The cost/loss function for regression is
# $$
# C(\beta_0, \beta_1, ... , \beta_{p-1}) = \frac{1}{n}\sum_{i=0}^{n} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij}\beta_j\right)^2,.
# $$
# Recall also that we use the squared value. This expression can lead to an
# increased penalty for higher differences between predicted and
# output/target values.
#
# What we have done is to single out the $\beta_0$ term in the
# definition of the mean squared error (MSE). The design matrix $X$
# does in this case not contain any intercept column. When we take the
# derivative with respect to $\beta_0$, we want the derivative to obey
# $$
# \frac{\partial C}{\partial \beta_j} = 0,
# $$
# for all $j$. For $\beta_0$ we have
# $$
# \frac{\partial C}{\partial \beta_0} = -\frac{2}{n}\sum_{i=0}^{n-1} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij} \beta_j\right).
# $$
# Multiplying away the constant $2/n$, we obtain
# $$
# \sum_{i=0}^{n-1} \beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} \sum_{j=1}^{p-1} X_{ij} \beta_j.
# $$
# Let us specialize first to the case where we have only two parameters $\beta_0$ and $\beta_1$.
# Our result for $\beta_0$ simplifies then to
# $$
# n\beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \beta_1.
# $$
# We obtain then
# $$
# \beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \beta_1\frac{1}{n}\sum_{i=0}^{n-1} X_{i1}.
# $$
# If we define
# $$
# \mu_{\boldsymbol{x}_1}=\frac{1}{n}\sum_{i=0}^{n-1} X_{i1},
# $$
# and the mean value of the outputs as
# $$
# \mu_y=\frac{1}{n}\sum_{i=0}^{n-1}y_i,
# $$
# we have
# $$
# \beta_0 = \mu_y - \beta_1\mu_{\boldsymbol{x}_1}.
# $$
# In the general case with more parameters than $\beta_0$ and $\beta_1$, we have
# $$
# \beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \frac{1}{n}\sum_{i=0}^{n-1}\sum_{j=1}^{p-1} X_{ij}\beta_j.
# $$
# We can rewrite the latter equation as
# $$
# \beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \sum_{j=1}^{p-1} \mu_{\boldsymbol{x}_j}\beta_j,
# $$
# where we have defined
# $$
# \mu_{\boldsymbol{x}_j}=\frac{1}{n}\sum_{i=0}^{n-1} X_{ij},
# $$
# the mean value for all elements of the column vector $\boldsymbol{x}_j$.
#
# Replacing $y_i$ with $y_i - y_i - \overline{\boldsymbol{y}}$ and centering also our design matrix results in a cost function (in vector-matrix disguise)
# $$
# C(\boldsymbol{\beta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta}).
# $$
# If we minimize with respect to $\boldsymbol{\beta}$ we have then
# $$
# \hat{\boldsymbol{\beta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}},
# $$
# where $\boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\boldsymbol{y}}$
# and $\tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj}$.
#
# For Ridge regression we need to add $\lambda \boldsymbol{\beta}^T\boldsymbol{\beta}$ to the cost function and get then
# $$
# \hat{\boldsymbol{\beta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}.
# $$
# Now we try to implement this.
# In[10]:
np.random.seed(2018)
n = 100
# we do not include the intercept
d = 2
Lambda = 0.01
# Make data set.
x = np.linspace(-3, 3, n)
y = 2.0 + 0.5*x + 5.0*(x**2)+ np.random.randn(n)
#Design matrix X does not include the intercept.
X = np.zeros((len(x), d))
for p in range(d):
X[:, p] = x ** (p+1)
#Split data in train and test
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
# Scale data by subtracting mean value,own implementation
#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable
X_train_mean = np.mean(X_train,axis=0)
#Center by removing mean from each feature
X_train_scaled = X_train - X_train_mean
X_test_scaled = X_test - X_train_mean
#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note)
y_scaler = np.mean(y_train)
y_train_scaled = y_train - y_scaler
#Calculate beta
beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)
beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)
print(beta_OLS)
print(beta_Ridge)
# calculate intercepts and print them
interceptOLS = y_scaler - X_train_mean @ beta_OLS
interceptRidge = y_scaler - X_train_mean @ beta_Ridge
print(interceptOLS)
print(interceptRidge)
#predict value with intercept
ytilde_test_OLS = X_test_scaled @ beta_OLS+y_scaler
ytilde_test_Ridge = X_test_scaled @ beta_Ridge+y_scaler
#Calculate MSE
print(" ")
print("test MSE of OLS:")
print(MSE(y_test,ytilde_test_OLS))
print(" ")
print("test MSE of Ridge")
print(MSE(y_test,ytilde_test_Ridge))
plt.scatter(x,y,label='Data')
plt.plot(x, X @ beta_OLS+interceptOLS,'*', label="OLS_Fit")
plt.plot(x, X @ beta_Ridge+interceptRidge, label="Ridge_Fit")
plt.grid()
plt.legend()
plt.show()
# Finally, instead of using our own function we repeat the same example
# using the **standardscaler** functionality of the library
# **Scikit-Learn**. Here we limit ourselves to Ridge regression only.
# In[11]:
from sklearn import linear_model
np.random.seed(2018)
n = 10
d = 2
Lambda = 0.01
# Make data set.
x = np.linspace(-3, 3, n)
y = 2.0 + 0.5*x + 5.0*(x**2)+ np.random.randn(n)
# Design matrix X does not include the intercept.
X = np.zeros((n, d))
for p in range(d):
X[:, p] = x ** (p+1)
#Split data in train and test
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
# Scale data by subtracting mean value of the input using scikit-learn
scaler = StandardScaler(with_std=False)
scaler.fit(X_train)
X_train_mean = np.mean(X_train,axis=0)
X_train_scaled = scaler.transform(X_train)
X_test_scaled = scaler.transform(X_test)
# We scale also the output, here by our own code
y_scaler = np.mean(y_train)
y_train_scaled = y_train - y_scaler
y_test_scaled = y_test- y_scaler
#Calculate beta
OLS = LinearRegression()
betaOLS=OLS.fit(X_train_scaled,y_train_scaled)
ypredictOLS = OLS.predict(X_test_scaled)
linear_model.Ridge(Lambda)
RegRidge.fit(X_train_scaled,y_train_scaled)
ypredictRidge = RegRidge.predict(X_test_scaled)
betaOLS = OLS.coef_
betaRidge = RegRidge.coef_
print(betaOLS)
print(betaRidge)
interceptOLS = np.mean(y_train) - X_train_mean @ betaOLS
interceptRidge = y_scaler - X_train_mean @ betaRidge
print(interceptOLS)
print(interceptRidge)
#predict value
ytilde_test_Ridge = X_test_scaled @ betaRidge+y_scaler
ytilde_test_OLS = X_test_scaled @ betaOLS+y_scaler
#Calculate MSE
print(" ")
print("test MSE of OLS")
print(MSE(y_test,ytilde_test_OLS))
print(" ")
print("test MSE of Ridge")
print(MSE(y_test,ytilde_test_Ridge))
plt.scatter(x,y,label='Data')
plt.plot(x, X @ RegRidge.coef_ + RegRidge.intercept_ , label="Ridge_Fit")
plt.grid()
plt.legend()
plt.show()
doc/LectureNotes/_build/jupyter_execute/week37_144_0.png
BIN
View File
Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 22 KiB

After

Width:  |  Height:  |  Size: 22 KiB

doc/LectureNotes/_build/jupyter_execute/week37_162_3.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 29 KiB

doc/LectureNotes/_build/jupyter_execute/week37_176_1.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 33 KiB

doc/LectureNotes/_build/jupyter_execute/week37_208_1.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 33 KiB

doc/LectureNotes/week37.ipynb
+814 -188
View File
File diff suppressed because it is too large Load Diff