updating typos
This commit is contained in:
Morten Hjorth-Jensen
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -2,7 +2,7 @@
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "cfebae03",
|
||||
"id": "1ff4a4e9",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -14,7 +14,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "353ebb07",
|
||||
"id": "0569f7a9",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -27,25 +27,34 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c02815dd",
|
||||
"id": "6f34e507",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"## Overarching aims of the exercises this week\n",
|
||||
"\n",
|
||||
"This exercise deals with various mean values and variances in linear regression method (here it may be useful to look up chapter 3, equation (3.8) of [Trevor Hastie, Robert Tibshirani, Jerome H. Friedman, The Elements of Statistical Learning, Springer](https://www.springer.com/gp/book/9780387848570)). The exercise is also a part of project 1 and can be reused in the theory part of the project.\n",
|
||||
"This exercise deals with various mean values and variances in linear\n",
|
||||
"regression method (here it may be useful to look up chapter 3,\n",
|
||||
"equation (3.8) of [Trevor Hastie, Robert Tibshirani, Jerome\n",
|
||||
"H. Friedman, The Elements of Statistical Learning,\n",
|
||||
"Springer](https://www.springer.com/gp/book/9780387848570)). The\n",
|
||||
"exercise is also a part of project 1 and can be reused in the theory\n",
|
||||
"part of the project.\n",
|
||||
"\n",
|
||||
"For more discussions on Ridge regression and calculation of expectation values, [Wessel van Wieringen's](https://arxiv.org/abs/1509.09169) article is highly recommended.\n",
|
||||
"For more discussions on Ridge regression and calculation of\n",
|
||||
"expectation values, [Wessel van\n",
|
||||
"Wieringen's](https://arxiv.org/abs/1509.09169) article is highly\n",
|
||||
"recommended.\n",
|
||||
"\n",
|
||||
"The assumption we have made is \n",
|
||||
"that there exists a continuous function $f(\\boldsymbol{x})$ and a normal distributed error $\\boldsymbol{\\varepsilon}\\sim N(0, \\sigma^2)$\n",
|
||||
"which describes our data"
|
||||
"The assumption we have made is that there exists a continuous function\n",
|
||||
"$f(\\boldsymbol{x})$ and a normal distributed error $\\boldsymbol{\\varepsilon}\\sim N(0,\n",
|
||||
"\\sigma^2)$ which describes our data"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "fb152043",
|
||||
"id": "3ff574a4",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -57,7 +66,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "9aba7a46",
|
||||
"id": "096a3165",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -68,7 +77,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "be52c097",
|
||||
"id": "cc9fb258",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -80,7 +89,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c54979e5",
|
||||
"id": "daf4eb5d",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -90,7 +99,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4b1eb811",
|
||||
"id": "3b491679",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -102,7 +111,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "036a5440",
|
||||
"id": "93a49739",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -114,7 +123,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "f7fb8094",
|
||||
"id": "a3fa80db",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -125,7 +134,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d7d92995",
|
||||
"id": "e6e2defd",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -137,7 +146,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "29eb9701",
|
||||
"id": "9aa19afa",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -150,7 +159,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "aa106511",
|
||||
"id": "5fa82a47",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -162,7 +171,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "58be1091",
|
||||
"id": "33099abb",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -172,7 +181,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4cb58e41",
|
||||
"id": "5d03f001",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -184,7 +193,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a04df28d",
|
||||
"id": "94bd151d",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -195,7 +204,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c095c96a",
|
||||
"id": "f0b84941",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -207,20 +216,20 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "5d425a79",
|
||||
"id": "60370a63",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"$$\n",
|
||||
"\\mathbb{E} \\big[ \\boldsymbol{\\beta}^{\\mathrm{Ridge}} \\big]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\\n",
|
||||
"\\mathbb{E} \\big[ \\hat{\\boldsymbol{\\beta}}^{\\mathrm{Ridge}} \\big]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\\n",
|
||||
"\\boldsymbol{\\beta}^{\\mathrm{OLS}}.\n",
|
||||
"$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "fb89cfab",
|
||||
"id": "af251217",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -233,19 +242,19 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "aea8801c",
|
||||
"id": "02e978da",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"$$\n",
|
||||
"\\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},\n",
|
||||
"\\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},\n",
|
||||
"$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "0e885693",
|
||||
"id": "3656a896",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -438,11 +438,20 @@ doconce format html exercisesweek37.do.txt -->
|
||||
<p>Date: <strong>Deadline is Sunday September 17 at midnight</strong></p>
|
||||
<div class="section" id="overarching-aims-of-the-exercises-this-week">
|
||||
<h2>Overarching aims of the exercises this week<a class="headerlink" href="#overarching-aims-of-the-exercises-this-week" title="Permalink to this headline">¶</a></h2>
|
||||
<p>This exercise deals with various mean values and variances in linear regression method (here it may be useful to look up chapter 3, equation (3.8) of <a class="reference external" href="https://www.springer.com/gp/book/9780387848570">Trevor Hastie, Robert Tibshirani, Jerome H. Friedman, The Elements of Statistical Learning, Springer</a>). The exercise is also a part of project 1 and can be reused in the theory part of the project.</p>
|
||||
<p>For more discussions on Ridge regression and calculation of expectation values, <a class="reference external" href="https://arxiv.org/abs/1509.09169">Wessel van Wieringen’s</a> article is highly recommended.</p>
|
||||
<p>The assumption we have made is
|
||||
that there exists a continuous function <span class="math notranslate nohighlight">\(f(\boldsymbol{x})\)</span> and a normal distributed error <span class="math notranslate nohighlight">\(\boldsymbol{\varepsilon}\sim N(0, \sigma^2)\)</span>
|
||||
which describes our data</p>
|
||||
<p>This exercise deals with various mean values and variances in linear
|
||||
regression method (here it may be useful to look up chapter 3,
|
||||
equation (3.8) of <a class="reference external" href="https://www.springer.com/gp/book/9780387848570">Trevor Hastie, Robert Tibshirani, Jerome
|
||||
H. Friedman, The Elements of Statistical Learning,
|
||||
Springer</a>). The
|
||||
exercise is also a part of project 1 and can be reused in the theory
|
||||
part of the project.</p>
|
||||
<p>For more discussions on Ridge regression and calculation of
|
||||
expectation values, <a class="reference external" href="https://arxiv.org/abs/1509.09169">Wessel van
|
||||
Wieringen’s</a> article is highly
|
||||
recommended.</p>
|
||||
<p>The assumption we have made is that there exists a continuous function
|
||||
<span class="math notranslate nohighlight">\(f(\boldsymbol{x})\)</span> and a normal distributed error <span class="math notranslate nohighlight">\(\boldsymbol{\varepsilon}\sim N(0,
|
||||
\sigma^2)\)</span> which describes our data</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\boldsymbol{y} = f(\boldsymbol{x})+\boldsymbol{\varepsilon}
|
||||
@@ -488,7 +497,7 @@ A given parameter <span class="math notranslate nohighlight">\(\beta_j\)</span>
|
||||
<p>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})\
|
||||
\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}^{\mathrm{OLS}}.
|
||||
\]</div>
|
||||
<p>We see clearly that
|
||||
@@ -496,7 +505,7 @@ A given parameter <span class="math notranslate nohighlight">\(\beta_j\)</span>
|
||||
<p>Show also that the variance is</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>
|
||||
</div>
|
||||
|
||||
File diff suppressed because one or more lines are too long
@@ -1646,7 +1646,7 @@ Since we are not using <strong>Scikit-Learn</strong> here we can define our own
|
||||
</div>
|
||||
</div>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.996535469511469
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.9952714635635723
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
@@ -1663,7 +1663,7 @@ Since we are not using <strong>Scikit-Learn</strong> here we can define our own
|
||||
</div>
|
||||
</div>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.0077750600806588575
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.009039197557276654
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
@@ -1678,23 +1678,31 @@ Since we are not using <strong>Scikit-Learn</strong> here we can define our own
|
||||
</div>
|
||||
</div>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[0.05751737 0.02644748 0.02533184 0.02193035 0.02521356 0.00712899
|
||||
0.00857688 0.0101339 0.01097227 0.01132947 0.01644452 0.01655455
|
||||
0.00771046 0.05177021 0.01374437 0.04208274 0.01667873 0.00228986
|
||||
0.00749579 0.00852965 0.01628823 0.0341899 0.01178253 0.00626339
|
||||
0.00514913 0.00791019 0.00363986 0.00471966 0.00390805 0.00910241
|
||||
0.00232255 0.02670606 0.03224523 0.01637147 0.00918914 0.02705154
|
||||
0.00321694 0.01904583 0.0177181 0.00538287 0.02291121 0.01028255
|
||||
0.04566103 0.00954257 0.00011551 0.0283589 0.00953728 0.03399521
|
||||
0.01205286 0.02251524 0.00353655 0.02344605 0.05288681 0.01950574
|
||||
0.00330907 0.01492855 0.01472766 0.03952093 0.01491299 0.00187154
|
||||
0.02897168 0.00072037 0.00895315 0.02775293 0.01547992 0.04117184
|
||||
0.0116726 0.03330335 0.01515966 0.0070381 0.01317074 0.00654702
|
||||
0.00329264 0.02502101 0.00226569 0.05029894 0.01929643 0.03028379
|
||||
0.053994 0.03413302 0.01854824 0.00393744 0.0600658 0.01855624
|
||||
0.0590702 0.01826743 0.01039549 0.02151219 0.00928016 0.03536062
|
||||
0.00503218 0.08526352 0.00506765 0.02609885 0.04771105 0.0010059
|
||||
0.00659545 0.00143188 0.01489692 0.08237141]
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[6.11430166e-02 7.59428035e-03 7.64295978e-03 5.34442654e-02
|
||||
2.53861506e-02 7.46543658e-02 2.79339933e-02 1.32946226e-02
|
||||
3.80320028e-03 2.40474117e-02 4.36028477e-02 3.24382528e-02
|
||||
1.26022420e-02 5.06848116e-03 1.96793569e-02 1.04826413e-04
|
||||
9.25580178e-03 1.80413812e-02 1.22519792e-02 1.01505087e-02
|
||||
2.67180442e-02 6.71316703e-02 7.68415876e-03 3.33174939e-02
|
||||
1.43570786e-02 8.43459933e-03 4.11046407e-03 9.10453697e-04
|
||||
4.33142087e-03 7.64507671e-02 2.19696497e-02 1.18460525e-02
|
||||
4.35665906e-02 3.62818572e-02 1.61018347e-03 3.78738205e-02
|
||||
4.57031517e-03 3.11279984e-02 1.07745809e-02 7.09848429e-03
|
||||
2.30006341e-03 4.62535124e-02 6.70493970e-02 6.26624267e-03
|
||||
4.23399553e-02 1.88690267e-02 7.46184816e-03 1.75934097e-02
|
||||
3.60820168e-02 1.09322328e-02 3.16749353e-02 2.30626515e-02
|
||||
5.76196895e-04 1.51000420e-02 1.95456434e-03 1.22138768e-02
|
||||
2.52945186e-03 2.12804129e-03 3.32695433e-03 1.90886410e-02
|
||||
1.84291049e-02 1.78104004e-02 6.34498455e-02 1.50032818e-02
|
||||
3.43525003e-02 1.52033394e-03 3.95921014e-02 7.67148008e-02
|
||||
2.40893740e-02 6.75638425e-02 1.22105472e-02 7.72706640e-03
|
||||
3.00721135e-02 8.77081499e-02 1.40684484e-02 4.57081718e-02
|
||||
3.51190686e-02 7.42092776e-02 4.03105744e-02 1.08994018e-02
|
||||
1.37071826e-02 3.69106251e-02 1.94966662e-03 1.71112140e-02
|
||||
3.12669602e-03 3.01610677e-02 1.72139847e-02 2.75703556e-02
|
||||
2.23333773e-02 2.76892410e-03 8.05886306e-02 1.15545946e-01
|
||||
4.88029792e-02 4.59648068e-02 1.98750794e-02 2.69178261e-03
|
||||
1.07925220e-02 1.13151043e-03 5.00620642e-03 3.69871367e-03]
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
@@ -1763,15 +1771,15 @@ but now splitting the data into a training set and a test set.</p>
|
||||
</div>
|
||||
</div>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[ 1.95976303 0.48974963 2.82202765 3.79224555 -2.14248744]
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[ 2.04038959 -0.81343612 7.94304103 -3.83770135 1.70073443]
|
||||
Training R2
|
||||
0.9947613223847728
|
||||
0.9954441212986554
|
||||
Training MSE
|
||||
0.010596520060222926
|
||||
0.008653750612762521
|
||||
Test R2
|
||||
0.989645675273948
|
||||
0.99255249100346
|
||||
Test MSE
|
||||
0.01889099091713886
|
||||
0.016970019460234076
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
|
||||
@@ -519,6 +519,16 @@ const thebe_selector_output = ".output, .cell_output"
|
||||
Material for lecture Thursday September 7
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#important-technicalities-more-on-rescaling-data">
|
||||
Important technicalities: More on Rescaling data
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#test-function-for-what-happens-with-ols-ridge-and-lasso">
|
||||
Test Function for what happens with OLS, Ridge and Lasso
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#linking-the-regression-analysis-with-a-statistical-interpretation">
|
||||
Linking the regression analysis with a statistical interpretation
|
||||
@@ -602,16 +612,6 @@ const thebe_selector_output = ".output, .cell_output"
|
||||
Bayes’ Theorem and Ridge and Lasso Regression
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#test-function-for-what-happens-with-ols-ridge-and-lasso">
|
||||
Test Function for what happens with OLS, Ridge and Lasso
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#invoking-bayes-theorem">
|
||||
Invoking Bayes’ theorem
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#ridge-and-bayes">
|
||||
Ridge and Bayes
|
||||
@@ -790,6 +790,16 @@ const thebe_selector_output = ".output, .cell_output"
|
||||
Material for lecture Thursday September 7
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#important-technicalities-more-on-rescaling-data">
|
||||
Important technicalities: More on Rescaling data
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#test-function-for-what-happens-with-ols-ridge-and-lasso">
|
||||
Test Function for what happens with OLS, Ridge and Lasso
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#linking-the-regression-analysis-with-a-statistical-interpretation">
|
||||
Linking the regression analysis with a statistical interpretation
|
||||
@@ -873,16 +883,6 @@ const thebe_selector_output = ".output, .cell_output"
|
||||
Bayes’ Theorem and Ridge and Lasso Regression
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#test-function-for-what-happens-with-ols-ridge-and-lasso">
|
||||
Test Function for what happens with OLS, Ridge and Lasso
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#invoking-bayes-theorem">
|
||||
Invoking Bayes’ theorem
|
||||
</a>
|
||||
</li>
|
||||
<li class="toc-h2 nav-item toc-entry">
|
||||
<a class="reference internal nav-link" href="#ridge-and-bayes">
|
||||
Ridge and Bayes
|
||||
@@ -920,16 +920,19 @@ doconce format html week36.do.txt --no_mako -->
|
||||
</li>
|
||||
<li><p>Material for the lecture on Thursday September 7</p>
|
||||
<ul>
|
||||
<li><p>Linear Regression and links with Statistics, Resampling methods</p></li>
|
||||
<li><p>Recommended Reading: Goodfellow et al chapter 3 on probability theory, see URL:””</p></li>
|
||||
<li><p>Technicalities related to scaling and other issues with data handling</p></li>
|
||||
<li><p>Linear Regression and links with Statistics</p></li>
|
||||
<li><p><a class="reference external" href="https://www.deeplearningbook.org/">Recommended Reading: Goodfellow et al chapter 3 on probability theory</a></p></li>
|
||||
<li><p>See also Murphy, sections 2.4 (Gaussian distributions) and 3.2 (Bayesian Statistics, basis)</p></li>
|
||||
<li><p><a class="reference external" href="https://youtu.be/Kc20CFK0z7Y">Video of lecture</a></p></li>
|
||||
<li><p><a class="reference external" href="https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2023/NotesSep7.pdf">Whiteboard notes</a></p></li>
|
||||
</ul>
|
||||
</li>
|
||||
</ul>
|
||||
</div>
|
||||
<div class="section" id="material-for-the-active-learning-sessions-tuesday-and-wednesday">
|
||||
<h2>Material for the active learning sessions Tuesday and Wednesday<a class="headerlink" href="#material-for-the-active-learning-sessions-tuesday-and-wednesday" title="Permalink to this headline">¶</a></h2>
|
||||
<p>The material here contains a summary from last Week and discussion of SVD, Ridge and Lasso regression with examples</p>
|
||||
<p>The material here contains a summary from last week and discussion of SVD, Ridge and Lasso regression with examples</p>
|
||||
</div>
|
||||
<div class="section" id="linear-regression-and-the-svd">
|
||||
<h2>Linear Regression and the SVD<a class="headerlink" href="#linear-regression-and-the-svd" title="Permalink to this headline">¶</a></h2>
|
||||
@@ -1700,6 +1703,531 @@ C(\boldsymbol{\beta})=(4-2\beta_0)^2+(2-\beta_1)^2+\lambda(\vert\beta_0\vert+\ve
|
||||
<div class="section" id="material-for-lecture-thursday-september-7">
|
||||
<h2>Material for lecture Thursday September 7<a class="headerlink" href="#material-for-lecture-thursday-september-7" title="Permalink to this headline">¶</a></h2>
|
||||
</div>
|
||||
<div class="section" id="important-technicalities-more-on-rescaling-data">
|
||||
<h2>Important technicalities: More on Rescaling data<a class="headerlink" href="#important-technicalities-more-on-rescaling-data" title="Permalink to this headline">¶</a></h2>
|
||||
<p>When you are comparing your own code with for example <strong>Scikit-Learn</strong>’s
|
||||
library, there are some technicalities to keep in mind. The examples
|
||||
here demonstrate some of these aspects with potential pitfalls.</p>
|
||||
<p>The discussion here focuses on the role of the intercept, how we can
|
||||
set up the design matrix, what scaling we should use and other topics
|
||||
which tend confuse us.</p>
|
||||
<p>The intercept can be interpreted as the expected value of our
|
||||
target/output variables when all other predictors are set to zero.
|
||||
Thus, if we cannot assume that the expected outputs/targets are zero
|
||||
when all predictors are zero (the columns in the design matrix), it
|
||||
may be a bad idea to implement a model which penalizes the intercept.
|
||||
Furthermore, in for example Ridge and Lasso regression, the default solutions
|
||||
from the library <strong>Scikit-Learn</strong> (when not shrinking <span class="math notranslate nohighlight">\(\beta_0\)</span>) for the unknown parameters
|
||||
<span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>, are derived under the assumption that both <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> and
|
||||
<span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> are zero centered, that is we subtract the mean values.</p>
|
||||
<p>If our predictors represent different scales, then it is important to
|
||||
standardize the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> by subtracting the mean of each
|
||||
column from the corresponding column and dividing the column with its
|
||||
standard deviation. Most machine learning libraries do this as a default. This means that if you compare your code with the results from a given library,
|
||||
the results may differ.</p>
|
||||
<p>The
|
||||
<a class="reference external" href="https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html">Standardscaler</a>
|
||||
function in <strong>Scikit-Learn</strong> does this for us. For the data sets we
|
||||
have been studying in our various examples, the data are in many cases
|
||||
already scaled and there is no need to scale them. You as a user of different machine learning algorithms, should always perform a
|
||||
survey of your data, with a critical assessment of them in case you need to scale the data.</p>
|
||||
<p>If you need to scale the data, not doing so will give an <em>unfair</em>
|
||||
penalization of the parameters since their magnitude depends on the
|
||||
scale of their corresponding predictor.</p>
|
||||
<p>Suppose as an example that you
|
||||
you have an input variable given by the heights of different persons.
|
||||
Human height might be measured in inches or meters or
|
||||
kilometers. If measured in kilometers, a standard linear regression
|
||||
model with this predictor would probably give a much bigger
|
||||
coefficient term, than if measured in millimeters.
|
||||
This can clearly lead to problems in evaluating the cost/loss functions.</p>
|
||||
<p>Keep in mind that when you transform your data set before training a model, the same transformation needs to be done
|
||||
on your eventual new data set before making a prediction. If we translate this into a Python code, it would could be implemented as</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="sd">"""</span>
|
||||
<span class="sd">#Model training, we compute the mean value of y and X</span>
|
||||
<span class="sd">y_train_mean = np.mean(y_train)</span>
|
||||
<span class="sd">X_train_mean = np.mean(X_train,axis=0)</span>
|
||||
<span class="sd">X_train = X_train - X_train_mean</span>
|
||||
<span class="sd">y_train = y_train - y_train_mean</span>
|
||||
|
||||
<span class="sd"># The we fit our model with the training data</span>
|
||||
<span class="sd">trained_model = some_model.fit(X_train,y_train)</span>
|
||||
|
||||
|
||||
<span class="sd">#Model prediction, we need also to transform our data set used for the prediction.</span>
|
||||
<span class="sd">X_test = X_test - X_train_mean #Use mean from training data</span>
|
||||
<span class="sd">y_pred = trained_model(X_test)</span>
|
||||
<span class="sd">y_pred = y_pred + y_train_mean</span>
|
||||
<span class="sd">"""</span>
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
</div>
|
||||
<p>Let us try to understand what this may imply mathematically when we
|
||||
subtract the mean values, also known as <em>zero 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>What does this mean? And why do we insist on all this? Let us look at some examples.</p>
|
||||
<p>This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (<em>code example thanks to Øyvind Sigmundson Schøyen</em>). Here our scaling of the data is done by subtracting the mean values only.
|
||||
Note also that we do not split the data into training and test.</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">numpy</span> <span class="k">as</span> <span class="nn">np</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">from</span> <span class="nn">sklearn.linear_model</span> <span class="kn">import</span> <span class="n">LinearRegression</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">2021</span><span class="p">)</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">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="n">true_beta</span> <span class="o">=</span> <span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mf">0.5</span><span class="p">,</span> <span class="mf">3.7</span><span class="p">]</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="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">11</span><span class="p">)</span>
|
||||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span>
|
||||
<span class="n">np</span><span class="o">.</span><span class="n">asarray</span><span class="p">([</span><span class="n">x</span> <span class="o">**</span> <span class="n">p</span> <span class="o">*</span> <span class="n">b</span> <span class="k">for</span> <span class="n">p</span><span class="p">,</span> <span class="n">b</span> <span class="ow">in</span> <span class="nb">enumerate</span><span class="p">(</span><span class="n">true_beta</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="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">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">))</span>
|
||||
|
||||
<span class="n">degree</span> <span class="o">=</span> <span class="mi">3</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">degree</span><span class="p">))</span>
|
||||
|
||||
<span class="c1"># Include the intercept in the design matrix</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">degree</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="n">p</span>
|
||||
|
||||
<span class="n">beta</span> <span class="o">=</span> <span class="n">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="c1"># Intercept is included in the design matrix</span>
|
||||
<span class="n">skl</span> <span class="o">=</span> <span class="n">LinearRegression</span><span class="p">(</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span><span class="o">.</span><span class="n">fit</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="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"True beta: </span><span class="si">{</span><span class="n">true_beta</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Fitted beta: </span><span class="si">{</span><span class="n">beta</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Sklearn fitted beta: </span><span class="si">{</span><span class="n">skl</span><span class="o">.</span><span class="n">coef_</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||||
<span class="n">ypredictOwn</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta</span>
|
||||
<span class="n">ypredictSKL</span> <span class="o">=</span> <span class="n">skl</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"MSE with intercept column"</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</span><span class="p">,</span><span class="n">ypredictOwn</span><span class="p">))</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"MSE with intercept column from SKL"</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</span><span class="p">,</span><span class="n">ypredictSKL</span><span class="p">))</span>
|
||||
|
||||
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</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="s2">"Data"</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</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">"Fit"</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">skl</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X</span><span class="p">),</span> <span class="n">label</span><span class="o">=</span><span class="s2">"Sklearn (fit_intercept=False)"</span><span class="p">)</span>
|
||||
|
||||
|
||||
<span class="c1"># Do not include the intercept in the design matrix</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">degree</span> <span class="o">-</span> <span class="mi">1</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">degree</span> <span class="o">-</span> <span class="mi">1</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"># Intercept is not included in the design matrix</span>
|
||||
<span class="n">skl</span> <span class="o">=</span> <span class="n">LinearRegression</span><span class="p">(</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span><span class="o">.</span><span class="n">fit</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="c1"># Use centered values for X and y when computing coefficients</span>
|
||||
<span class="n">y_offset</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">average</span><span class="p">(</span><span class="n">y</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_offset</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">average</span><span class="p">(</span><span class="n">X</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">beta</span> <span class="o">=</span> <span class="n">fit_beta</span><span class="p">(</span><span class="n">X</span> <span class="o">-</span> <span class="n">X_offset</span><span class="p">,</span> <span class="n">y</span> <span class="o">-</span> <span class="n">y_offset</span><span class="p">)</span>
|
||||
<span class="n">intercept</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_offset</span> <span class="o">-</span> <span class="n">X_offset</span> <span class="o">@</span> <span class="n">beta</span><span class="p">)</span>
|
||||
|
||||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Manual intercept: </span><span class="si">{</span><span class="n">intercept</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Fitted beta (wiothout intercept): </span><span class="si">{</span><span class="n">beta</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Sklearn intercept: </span><span class="si">{</span><span class="n">skl</span><span class="o">.</span><span class="n">intercept_</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Sklearn fitted beta (without intercept): </span><span class="si">{</span><span class="n">skl</span><span class="o">.</span><span class="n">coef_</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||||
<span class="n">ypredictOwn</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta</span>
|
||||
<span class="n">ypredictSKL</span> <span class="o">=</span> <span class="n">skl</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"MSE with Manual intercept"</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</span><span class="p">,</span><span class="n">ypredictOwn</span><span class="o">+</span><span class="n">intercept</span><span class="p">))</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"MSE with Sklearn intercept"</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</span><span class="p">,</span><span class="n">ypredictSKL</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</span> <span class="o">+</span> <span class="n">intercept</span><span class="p">,</span> <span class="s2">"--"</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">"Fit (manual intercept)"</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">skl</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X</span><span class="p">),</span> <span class="s2">"--"</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">"Sklearn (fit_intercept=True)"</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>
|
||||
<p>The intercept is the value of our output/target variable
|
||||
when all our features are zero and our function crosses the <span class="math notranslate nohighlight">\(y\)</span>-axis (for a one-dimensional case).</p>
|
||||
<p>Printing the MSE, we see first that both methods give the same MSE, as
|
||||
they should. However, when we move to for example Ridge regression,
|
||||
the way we treat the intercept may give a larger or smaller MSE,
|
||||
meaning that the MSE can be penalized by the value of the
|
||||
intercept. Not including the intercept in the fit, means that the
|
||||
regularization term does not include <span class="math notranslate nohighlight">\(\beta_0\)</span>. For different values
|
||||
of <span class="math notranslate nohighlight">\(\lambda\)</span>, this may lead to different MSE values.</p>
|
||||
<p>To remind the reader, the regularization term, with the intercept in Ridge regression, is given by</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\lambda \vert\vert \boldsymbol{\beta} \vert\vert_2^2 = \lambda \sum_{j=0}^{p-1}\beta_j^2,
|
||||
\]</div>
|
||||
<p>but when we take out the intercept, this equation becomes</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\lambda \vert\vert \boldsymbol{\beta} \vert\vert_2^2 = \lambda \sum_{j=1}^{p-1}\beta_j^2.
|
||||
\]</div>
|
||||
<p>For Lasso regression we have</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\lambda \vert\vert \boldsymbol{\beta} \vert\vert_1 = \lambda \sum_{j=1}^{p-1}\vert\beta_j\vert.
|
||||
\]</div>
|
||||
<p>It means that, when scaling the design matrix and the outputs/targets,
|
||||
by subtracting the mean values, we have an optimization problem which
|
||||
is not penalized by the intercept. The MSE value can then be smaller
|
||||
since it focuses only on the remaining quantities. If we however bring
|
||||
back the intercept, we will get a MSE which then contains the
|
||||
intercept.</p>
|
||||
<p>Armed with this wisdom, we attempt first to simply set the intercept equal to <strong>False</strong> in our implementation of Ridge regression for our well-known vanilla data set.</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">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||||
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</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">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</span> <span class="kn">import</span> <span class="n">linear_model</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="c1"># A seed just to ensure that the random numbers are the same for every run.</span>
|
||||
<span class="c1"># Useful for eventual debugging.</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">3155</span><span class="p">)</span>
|
||||
|
||||
<span class="n">n</span> <span class="o">=</span> <span class="mi">100</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="n">n</span><span class="p">)</span>
|
||||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</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="p">)</span> <span class="o">+</span> <span class="mf">1.5</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</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="mi">2</span><span class="p">)</span>
|
||||
|
||||
<span class="n">Maxpolydegree</span> <span class="o">=</span> <span class="mi">20</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">Maxpolydegree</span><span class="p">))</span>
|
||||
<span class="c1">#We include explicitely the intercept column</span>
|
||||
<span class="k">for</span> <span class="n">degree</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">Maxpolydegree</span><span class="p">):</span>
|
||||
<span class="n">X</span><span class="p">[:,</span><span class="n">degree</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="n">degree</span>
|
||||
<span class="c1"># We split the data in test and training data</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="n">p</span> <span class="o">=</span> <span class="n">Maxpolydegree</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">p</span><span class="p">,</span><span class="n">p</span><span class="p">)</span>
|
||||
<span class="c1"># Decide which values of lambda to use</span>
|
||||
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">6</span>
|
||||
<span class="n">MSEOwnRidgePredict</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">nlambdas</span><span class="p">)</span>
|
||||
<span class="n">MSERidgePredict</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">nlambdas</span><span class="p">)</span>
|
||||
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
|
||||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
|
||||
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
|
||||
<span class="n">OwnRidgeBeta</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">pinv</span><span class="p">(</span><span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">+</span><span class="n">lmb</span><span class="o">*</span><span class="n">I</span><span class="p">)</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y_train</span>
|
||||
<span class="c1"># Note: we include the intercept column and no scaling</span>
|
||||
<span class="n">RegRidge</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Ridge</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</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</span><span class="p">,</span><span class="n">y_train</span><span class="p">)</span>
|
||||
<span class="c1"># and then make the prediction</span>
|
||||
<span class="n">ytildeOwnRidge</span> <span class="o">=</span> <span class="n">X_train</span> <span class="o">@</span> <span class="n">OwnRidgeBeta</span>
|
||||
<span class="n">ypredictOwnRidge</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">@</span> <span class="n">OwnRidgeBeta</span>
|
||||
<span class="n">ytildeRidge</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_train</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</span><span class="p">)</span>
|
||||
<span class="n">MSEOwnRidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictOwnRidge</span><span class="p">)</span>
|
||||
<span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictRidge</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Beta values for own Ridge implementation"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">OwnRidgeBeta</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Beta values for Scikit-Learn Ridge implementation"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">RegRidge</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s2">"MSE values for own Ridge implementation"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">MSEOwnRidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s2">"MSE values for Scikit-Learn Ridge implementation"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
|
||||
|
||||
<span class="c1"># Now plot the results</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</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">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSEOwnRidgePredict</span><span class="p">,</span> <span class="s1">'r'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE own Ridge Test'</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">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSERidgePredict</span><span class="p">,</span> <span class="s1">'g'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Ridge Test'</span><span class="p">)</span>
|
||||
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'log10(lambda)'</span><span class="p">)</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'MSE'</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>
|
||||
<p>The results here agree when we force <strong>Scikit-Learn</strong>’s Ridge function to include the first column in our design matrix.
|
||||
We see that the results agree very well. Here we have thus explicitely included the intercept column in the design matrix.
|
||||
What happens if we do not include the intercept in our fit?
|
||||
Let us see how we can change this code by zero centering.</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">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||||
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</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">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</span> <span class="kn">import</span> <span class="n">linear_model</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="c1"># A seed just to ensure that the random numbers are the same for every run.</span>
|
||||
<span class="c1"># Useful for eventual debugging.</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">315</span><span class="p">)</span>
|
||||
|
||||
<span class="n">n</span> <span class="o">=</span> <span class="mi">100</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="n">n</span><span class="p">)</span>
|
||||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</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="p">)</span> <span class="o">+</span> <span class="mf">1.5</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</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="mi">2</span><span class="p">)</span>
|
||||
|
||||
<span class="n">Maxpolydegree</span> <span class="o">=</span> <span class="mi">20</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">Maxpolydegree</span><span class="o">-</span><span class="mi">1</span><span class="p">))</span>
|
||||
|
||||
<span class="k">for</span> <span class="n">degree</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="n">Maxpolydegree</span><span class="p">):</span> <span class="c1">#No intercept column</span>
|
||||
<span class="n">X</span><span class="p">[:,</span><span class="n">degree</span><span class="o">-</span><span class="mi">1</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">degree</span><span class="p">)</span>
|
||||
|
||||
<span class="c1"># We split the data in test and training data</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">#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)</span>
|
||||
<span class="c1">#Remove the intercept from the training data.</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">p</span> <span class="o">=</span> <span class="n">Maxpolydegree</span><span class="o">-</span><span class="mi">1</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">p</span><span class="p">,</span><span class="n">p</span><span class="p">)</span>
|
||||
<span class="c1"># Decide which values of lambda to use</span>
|
||||
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">6</span>
|
||||
<span class="n">MSEOwnRidgePredict</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">nlambdas</span><span class="p">)</span>
|
||||
<span class="n">MSERidgePredict</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">nlambdas</span><span class="p">)</span>
|
||||
|
||||
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
|
||||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
|
||||
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
|
||||
<span class="n">OwnRidgeBeta</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">pinv</span><span class="p">(</span><span class="n">X_train_scaled</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X_train_scaled</span><span class="o">+</span><span class="n">lmb</span><span class="o">*</span><span class="n">I</span><span class="p">)</span> <span class="o">@</span> <span class="n">X_train_scaled</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="p">(</span><span class="n">y_train_scaled</span><span class="p">)</span>
|
||||
<span class="n">intercept_</span> <span class="o">=</span> <span class="n">y_scaler</span> <span class="o">-</span> <span class="n">X_train_mean</span><span class="nd">@OwnRidgeBeta</span> <span class="c1">#The intercept can be shifted so the model can predict on uncentered data</span>
|
||||
<span class="c1">#Add intercept to prediction</span>
|
||||
<span class="n">ypredictOwnRidge</span> <span class="o">=</span> <span class="n">X_test_scaled</span> <span class="o">@</span> <span class="n">OwnRidgeBeta</span> <span class="o">+</span> <span class="n">y_scaler</span>
|
||||
<span class="n">RegRidge</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Ridge</span><span class="p">(</span><span class="n">lmb</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</span><span class="p">,</span><span class="n">y_train</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</span><span class="p">)</span>
|
||||
<span class="n">MSEOwnRidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictOwnRidge</span><span class="p">)</span>
|
||||
<span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictRidge</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Beta values for own Ridge implementation"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">OwnRidgeBeta</span><span class="p">)</span> <span class="c1">#Intercept is given by mean of target variable</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Beta values for Scikit-Learn Ridge implementation"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">RegRidge</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s1">'Intercept from own implementation:'</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">intercept_</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s1">'Intercept from Scikit-Learn Ridge implementation'</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">RegRidge</span><span class="o">.</span><span class="n">intercept_</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s2">"MSE values for own Ridge implementation"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">MSEOwnRidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s2">"MSE values for Scikit-Learn Ridge implementation"</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
|
||||
|
||||
|
||||
<span class="c1"># Now plot the results</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</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">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSEOwnRidgePredict</span><span class="p">,</span> <span class="s1">'b--'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE own Ridge Test'</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">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSERidgePredict</span><span class="p">,</span> <span class="s1">'g--'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE SL Ridge Test'</span><span class="p">)</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'log10(lambda)'</span><span class="p">)</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'MSE'</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>
|
||||
<p>We see here, when compared to the code which includes explicitely the
|
||||
intercept column, that our MSE value is actually smaller. This is
|
||||
because the regularization term does not include the intercept value
|
||||
<span class="math notranslate nohighlight">\(\beta_0\)</span> in the fitting. This applies to Lasso regularization as
|
||||
well. It means that our optimization is now done only with the
|
||||
centered matrix and/or vector that enter the fitting procedure.</p>
|
||||
</div>
|
||||
<div class="section" id="test-function-for-what-happens-with-ols-ridge-and-lasso">
|
||||
<h2>Test Function for what happens with OLS, Ridge and Lasso<a class="headerlink" href="#test-function-for-what-happens-with-ols-ridge-and-lasso" title="Permalink to this headline">¶</a></h2>
|
||||
<p>Hitherto we have discussed Ridge and Lasso regression in terms of a
|
||||
linear analysis. This may to many of you feel rather technical and
|
||||
perhaps not that intuitive. The question is whether we can develop a
|
||||
more intuitive way of understanding what Ridge and Lasso express.</p>
|
||||
<p>Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.</p>
|
||||
<p>We will play around with a study of the values for the optimal
|
||||
parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> using OLS, Ridge and Lasso regression. For
|
||||
OLS, you will notice as function of the noise and polynomial degree,
|
||||
that the parameters <span class="math notranslate nohighlight">\(\beta\)</span> will fluctuate from order to order in the
|
||||
polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS.</p>
|
||||
<p>For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.</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">numpy</span> <span class="k">as</span> <span class="nn">np</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">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</span> <span class="kn">import</span> <span class="n">linear_model</span>
|
||||
|
||||
<span class="k">def</span> <span class="nf">R2</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="k">return</span> <span class="mi">1</span> <span class="o">-</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">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">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_data</span><span class="p">))</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</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="c1"># Make data set.</span>
|
||||
<span class="n">n</span> <span class="o">=</span> <span class="mi">10000</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="n">n</span><span class="p">)</span>
|
||||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</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="p">)</span> <span class="o">+</span> <span class="mf">1.5</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</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="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="n">Maxpolydegree</span> <span class="o">=</span> <span class="mi">5</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">Maxpolydegree</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="k">for</span> <span class="n">polydegree</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="n">Maxpolydegree</span><span class="p">):</span>
|
||||
<span class="n">X</span><span class="p">[:,</span><span class="n">polydegree</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">polydegree</span><span class="p">)</span>
|
||||
|
||||
<span class="c1"># We split the data in test and training data</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"># matrix inversion to find beta</span>
|
||||
<span class="n">OLSbeta</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">pinv</span><span class="p">(</span><span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X_train</span><span class="p">)</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y_train</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">OLSbeta</span><span class="p">)</span>
|
||||
<span class="n">ypredictOLS</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">@</span> <span class="n">OLSbeta</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Test MSE OLS"</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">ypredictOLS</span><span class="p">))</span>
|
||||
<span class="c1"># Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn</span>
|
||||
<span class="c1"># Decide which values of lambda to use</span>
|
||||
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">4</span>
|
||||
<span class="n">MSERidgePredict</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">nlambdas</span><span class="p">)</span>
|
||||
<span class="n">MSELassoPredict</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">nlambdas</span><span class="p">)</span>
|
||||
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
|
||||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
|
||||
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
|
||||
<span class="c1"># Make the fit using Ridge and Lasso</span>
|
||||
<span class="n">RegRidge</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Ridge</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</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</span><span class="p">,</span><span class="n">y_train</span><span class="p">)</span>
|
||||
<span class="n">RegLasso</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Lasso</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
|
||||
<span class="n">RegLasso</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">y_train</span><span class="p">)</span>
|
||||
<span class="c1"># and then make the prediction</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</span><span class="p">)</span>
|
||||
<span class="n">ypredictLasso</span> <span class="o">=</span> <span class="n">RegLasso</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test</span><span class="p">)</span>
|
||||
<span class="c1"># Compute the MSE and print it</span>
|
||||
<span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictRidge</span><span class="p">)</span>
|
||||
<span class="n">MSELassoPredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictLasso</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">RegRidge</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">RegLasso</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||||
<span class="c1"># Now plot the results</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</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">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSERidgePredict</span><span class="p">,</span> <span class="s1">'b'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Ridge Test'</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">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSELassoPredict</span><span class="p">,</span> <span class="s1">'r'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Lasso Test'</span><span class="p">)</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'log10(lambda)'</span><span class="p">)</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'MSE'</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>
|
||||
<p>How can we understand this?</p>
|
||||
</div>
|
||||
<div class="section" id="linking-the-regression-analysis-with-a-statistical-interpretation">
|
||||
<h2>Linking the regression analysis with a statistical interpretation<a class="headerlink" href="#linking-the-regression-analysis-with-a-statistical-interpretation" title="Permalink to this headline">¶</a></h2>
|
||||
<p>We will now couple the discussions of ordinary least squares, Ridge
|
||||
@@ -2036,93 +2564,6 @@ p(Y=1\vert X=1)=\frac{p(X=1\vert Y=1)p(Y=1)}{p(X=1\vert Y=1)p(Y=1)+p(X=1\vert Y=
|
||||
</div>
|
||||
<div class="section" id="bayes-theorem-and-ridge-and-lasso-regression">
|
||||
<h2>Bayes’ Theorem and Ridge and Lasso Regression<a class="headerlink" href="#bayes-theorem-and-ridge-and-lasso-regression" title="Permalink to this headline">¶</a></h2>
|
||||
<p>Hitherto we have discussed Ridge and Lasso regression in terms of a
|
||||
linear analysis. This may to many of you feel rather technical and
|
||||
perhaps not that intuitive. The question is whether we can develop a
|
||||
more intuitive way of understanding what Ridge and Lasso express.</p>
|
||||
<p>Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.</p>
|
||||
</div>
|
||||
<div class="section" id="test-function-for-what-happens-with-ols-ridge-and-lasso">
|
||||
<h2>Test Function for what happens with OLS, Ridge and Lasso<a class="headerlink" href="#test-function-for-what-happens-with-ols-ridge-and-lasso" title="Permalink to this headline">¶</a></h2>
|
||||
<p>We will play around with a study of the values for the optimal
|
||||
parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> using OLS, Ridge and Lasso regression. For
|
||||
OLS, you will notice as function of the noise and polynomial degree,
|
||||
that the parameters <span class="math notranslate nohighlight">\(\beta\)</span> will fluctuate from order to order in the
|
||||
polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS.</p>
|
||||
<p>For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.</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">numpy</span> <span class="k">as</span> <span class="nn">np</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">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</span> <span class="kn">import</span> <span class="n">linear_model</span>
|
||||
|
||||
<span class="k">def</span> <span class="nf">R2</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="k">return</span> <span class="mi">1</span> <span class="o">-</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">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">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_data</span><span class="p">))</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</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="c1"># Make data set.</span>
|
||||
<span class="n">n</span> <span class="o">=</span> <span class="mi">10000</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="n">n</span><span class="p">)</span>
|
||||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</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="p">)</span> <span class="o">+</span> <span class="mf">1.5</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</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="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="n">Maxpolydegree</span> <span class="o">=</span> <span class="mi">5</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">Maxpolydegree</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="k">for</span> <span class="n">polydegree</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="n">Maxpolydegree</span><span class="p">):</span>
|
||||
<span class="k">for</span> <span class="n">degree</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">polydegree</span><span class="p">):</span>
|
||||
<span class="n">X</span><span class="p">[:,</span><span class="n">degree</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">degree</span><span class="p">)</span>
|
||||
|
||||
|
||||
<span class="c1"># We split the data in test and training data</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"># matrix inversion to find beta</span>
|
||||
<span class="n">OLSbeta</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">pinv</span><span class="p">(</span><span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X_train</span><span class="p">)</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y_train</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">OLSbeta</span><span class="p">)</span>
|
||||
<span class="n">ypredictOLS</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">@</span> <span class="n">OLSbeta</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Test MSE OLS"</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">ypredictOLS</span><span class="p">))</span>
|
||||
<span class="c1"># Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn</span>
|
||||
<span class="c1"># Decide which values of lambda to use</span>
|
||||
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">4</span>
|
||||
<span class="n">MSERidgePredict</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">nlambdas</span><span class="p">)</span>
|
||||
<span class="n">MSELassoPredict</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">nlambdas</span><span class="p">)</span>
|
||||
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
|
||||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
|
||||
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
|
||||
<span class="c1"># Make the fit using Ridge and Lasso</span>
|
||||
<span class="n">RegRidge</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Ridge</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</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</span><span class="p">,</span><span class="n">y_train</span><span class="p">)</span>
|
||||
<span class="n">RegLasso</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Lasso</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
|
||||
<span class="n">RegLasso</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">y_train</span><span class="p">)</span>
|
||||
<span class="c1"># and then make the prediction</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</span><span class="p">)</span>
|
||||
<span class="n">ypredictLasso</span> <span class="o">=</span> <span class="n">RegLasso</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test</span><span class="p">)</span>
|
||||
<span class="c1"># Compute the MSE and print it</span>
|
||||
<span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictRidge</span><span class="p">)</span>
|
||||
<span class="n">MSELassoPredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictLasso</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">RegRidge</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">RegLasso</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||||
<span class="c1"># Now plot the results</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</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">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSERidgePredict</span><span class="p">,</span> <span class="s1">'b'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Ridge Test'</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">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSELassoPredict</span><span class="p">,</span> <span class="s1">'r'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Lasso Test'</span><span class="p">)</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'log10(lambda)'</span><span class="p">)</span>
|
||||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'MSE'</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>
|
||||
<p>How can we understand this?</p>
|
||||
</div>
|
||||
<div class="section" id="invoking-bayes-theorem">
|
||||
<h2>Invoking Bayes’ theorem<a class="headerlink" href="#invoking-bayes-theorem" title="Permalink to this headline">¶</a></h2>
|
||||
<p>Using Bayes’ theorem we can gain a better intuition about Ridge and Lasso regression.</p>
|
||||
<p>For ordinary least squares we postulated that the maximum likelihood for the doamin of events <span class="math notranslate nohighlight">\(\boldsymbol{D}\)</span> (one-dimensional case)</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
|
||||
@@ -2,7 +2,7 @@
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "cfebae03",
|
||||
"id": "1ff4a4e9",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -14,7 +14,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "353ebb07",
|
||||
"id": "0569f7a9",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -27,25 +27,34 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c02815dd",
|
||||
"id": "6f34e507",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"## Overarching aims of the exercises this week\n",
|
||||
"\n",
|
||||
"This exercise deals with various mean values and variances in linear regression method (here it may be useful to look up chapter 3, equation (3.8) of [Trevor Hastie, Robert Tibshirani, Jerome H. Friedman, The Elements of Statistical Learning, Springer](https://www.springer.com/gp/book/9780387848570)). The exercise is also a part of project 1 and can be reused in the theory part of the project.\n",
|
||||
"This exercise deals with various mean values and variances in linear\n",
|
||||
"regression method (here it may be useful to look up chapter 3,\n",
|
||||
"equation (3.8) of [Trevor Hastie, Robert Tibshirani, Jerome\n",
|
||||
"H. Friedman, The Elements of Statistical Learning,\n",
|
||||
"Springer](https://www.springer.com/gp/book/9780387848570)). The\n",
|
||||
"exercise is also a part of project 1 and can be reused in the theory\n",
|
||||
"part of the project.\n",
|
||||
"\n",
|
||||
"For more discussions on Ridge regression and calculation of expectation values, [Wessel van Wieringen's](https://arxiv.org/abs/1509.09169) article is highly recommended.\n",
|
||||
"For more discussions on Ridge regression and calculation of\n",
|
||||
"expectation values, [Wessel van\n",
|
||||
"Wieringen's](https://arxiv.org/abs/1509.09169) article is highly\n",
|
||||
"recommended.\n",
|
||||
"\n",
|
||||
"The assumption we have made is \n",
|
||||
"that there exists a continuous function $f(\\boldsymbol{x})$ and a normal distributed error $\\boldsymbol{\\varepsilon}\\sim N(0, \\sigma^2)$\n",
|
||||
"which describes our data"
|
||||
"The assumption we have made is that there exists a continuous function\n",
|
||||
"$f(\\boldsymbol{x})$ and a normal distributed error $\\boldsymbol{\\varepsilon}\\sim N(0,\n",
|
||||
"\\sigma^2)$ which describes our data"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "fb152043",
|
||||
"id": "3ff574a4",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -57,7 +66,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "9aba7a46",
|
||||
"id": "096a3165",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -68,7 +77,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "be52c097",
|
||||
"id": "cc9fb258",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -80,7 +89,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c54979e5",
|
||||
"id": "daf4eb5d",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -90,7 +99,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4b1eb811",
|
||||
"id": "3b491679",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -102,7 +111,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "036a5440",
|
||||
"id": "93a49739",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -114,7 +123,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "f7fb8094",
|
||||
"id": "a3fa80db",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -125,7 +134,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d7d92995",
|
||||
"id": "e6e2defd",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -137,7 +146,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "29eb9701",
|
||||
"id": "9aa19afa",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -150,7 +159,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "aa106511",
|
||||
"id": "5fa82a47",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -162,7 +171,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "58be1091",
|
||||
"id": "33099abb",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -172,7 +181,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4cb58e41",
|
||||
"id": "5d03f001",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -184,7 +193,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a04df28d",
|
||||
"id": "94bd151d",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -195,7 +204,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c095c96a",
|
||||
"id": "f0b84941",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -207,20 +216,20 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "5d425a79",
|
||||
"id": "60370a63",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"$$\n",
|
||||
"\\mathbb{E} \\big[ \\boldsymbol{\\beta}^{\\mathrm{Ridge}} \\big]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\\n",
|
||||
"\\mathbb{E} \\big[ \\hat{\\boldsymbol{\\beta}}^{\\mathrm{Ridge}} \\big]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\\n",
|
||||
"\\boldsymbol{\\beta}^{\\mathrm{OLS}}.\n",
|
||||
"$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "fb89cfab",
|
||||
"id": "af251217",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -233,19 +242,19 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "aea8801c",
|
||||
"id": "02e978da",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"$$\n",
|
||||
"\\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},\n",
|
||||
"\\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},\n",
|
||||
"$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "0e885693",
|
||||
"id": "3656a896",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
|
||||
@@ -1519,7 +1519,7 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"0.996535469511469\n"
|
||||
"0.9952714635635723\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
@@ -1550,7 +1550,7 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"0.0077750600806588575\n"
|
||||
"0.009039197557276654\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
@@ -1585,23 +1585,31 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"[0.05751737 0.02644748 0.02533184 0.02193035 0.02521356 0.00712899\n",
|
||||
" 0.00857688 0.0101339 0.01097227 0.01132947 0.01644452 0.01655455\n",
|
||||
" 0.00771046 0.05177021 0.01374437 0.04208274 0.01667873 0.00228986\n",
|
||||
" 0.00749579 0.00852965 0.01628823 0.0341899 0.01178253 0.00626339\n",
|
||||
" 0.00514913 0.00791019 0.00363986 0.00471966 0.00390805 0.00910241\n",
|
||||
" 0.00232255 0.02670606 0.03224523 0.01637147 0.00918914 0.02705154\n",
|
||||
" 0.00321694 0.01904583 0.0177181 0.00538287 0.02291121 0.01028255\n",
|
||||
" 0.04566103 0.00954257 0.00011551 0.0283589 0.00953728 0.03399521\n",
|
||||
" 0.01205286 0.02251524 0.00353655 0.02344605 0.05288681 0.01950574\n",
|
||||
" 0.00330907 0.01492855 0.01472766 0.03952093 0.01491299 0.00187154\n",
|
||||
" 0.02897168 0.00072037 0.00895315 0.02775293 0.01547992 0.04117184\n",
|
||||
" 0.0116726 0.03330335 0.01515966 0.0070381 0.01317074 0.00654702\n",
|
||||
" 0.00329264 0.02502101 0.00226569 0.05029894 0.01929643 0.03028379\n",
|
||||
" 0.053994 0.03413302 0.01854824 0.00393744 0.0600658 0.01855624\n",
|
||||
" 0.0590702 0.01826743 0.01039549 0.02151219 0.00928016 0.03536062\n",
|
||||
" 0.00503218 0.08526352 0.00506765 0.02609885 0.04771105 0.0010059\n",
|
||||
" 0.00659545 0.00143188 0.01489692 0.08237141]\n"
|
||||
"[6.11430166e-02 7.59428035e-03 7.64295978e-03 5.34442654e-02\n",
|
||||
" 2.53861506e-02 7.46543658e-02 2.79339933e-02 1.32946226e-02\n",
|
||||
" 3.80320028e-03 2.40474117e-02 4.36028477e-02 3.24382528e-02\n",
|
||||
" 1.26022420e-02 5.06848116e-03 1.96793569e-02 1.04826413e-04\n",
|
||||
" 9.25580178e-03 1.80413812e-02 1.22519792e-02 1.01505087e-02\n",
|
||||
" 2.67180442e-02 6.71316703e-02 7.68415876e-03 3.33174939e-02\n",
|
||||
" 1.43570786e-02 8.43459933e-03 4.11046407e-03 9.10453697e-04\n",
|
||||
" 4.33142087e-03 7.64507671e-02 2.19696497e-02 1.18460525e-02\n",
|
||||
" 4.35665906e-02 3.62818572e-02 1.61018347e-03 3.78738205e-02\n",
|
||||
" 4.57031517e-03 3.11279984e-02 1.07745809e-02 7.09848429e-03\n",
|
||||
" 2.30006341e-03 4.62535124e-02 6.70493970e-02 6.26624267e-03\n",
|
||||
" 4.23399553e-02 1.88690267e-02 7.46184816e-03 1.75934097e-02\n",
|
||||
" 3.60820168e-02 1.09322328e-02 3.16749353e-02 2.30626515e-02\n",
|
||||
" 5.76196895e-04 1.51000420e-02 1.95456434e-03 1.22138768e-02\n",
|
||||
" 2.52945186e-03 2.12804129e-03 3.32695433e-03 1.90886410e-02\n",
|
||||
" 1.84291049e-02 1.78104004e-02 6.34498455e-02 1.50032818e-02\n",
|
||||
" 3.43525003e-02 1.52033394e-03 3.95921014e-02 7.67148008e-02\n",
|
||||
" 2.40893740e-02 6.75638425e-02 1.22105472e-02 7.72706640e-03\n",
|
||||
" 3.00721135e-02 8.77081499e-02 1.40684484e-02 4.57081718e-02\n",
|
||||
" 3.51190686e-02 7.42092776e-02 4.03105744e-02 1.08994018e-02\n",
|
||||
" 1.37071826e-02 3.69106251e-02 1.94966662e-03 1.71112140e-02\n",
|
||||
" 3.12669602e-03 3.01610677e-02 1.72139847e-02 2.75703556e-02\n",
|
||||
" 2.23333773e-02 2.76892410e-03 8.05886306e-02 1.15545946e-01\n",
|
||||
" 4.88029792e-02 4.59648068e-02 1.98750794e-02 2.69178261e-03\n",
|
||||
" 1.07925220e-02 1.13151043e-03 5.00620642e-03 3.69871367e-03]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
@@ -1655,15 +1663,15 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"[ 1.95976303 0.48974963 2.82202765 3.79224555 -2.14248744]\n",
|
||||
"[ 2.04038959 -0.81343612 7.94304103 -3.83770135 1.70073443]\n",
|
||||
"Training R2\n",
|
||||
"0.9947613223847728\n",
|
||||
"0.9954441212986554\n",
|
||||
"Training MSE\n",
|
||||
"0.010596520060222926\n",
|
||||
"0.008653750612762521\n",
|
||||
"Test R2\n",
|
||||
"0.989645675273948\n",
|
||||
"0.99255249100346\n",
|
||||
"Test MSE\n",
|
||||
"0.01889099091713886\n"
|
||||
"0.016970019460234076\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -22,15 +22,21 @@
|
||||
#
|
||||
# * Material for the lecture on Thursday September 7
|
||||
#
|
||||
# * Linear Regression and links with Statistics, Resampling methods
|
||||
# * Technicalities related to scaling and other issues with data handling
|
||||
#
|
||||
# * Recommended Reading: Goodfellow et al chapter 3 on probability theory, see URL:""
|
||||
# * Linear Regression and links with Statistics
|
||||
#
|
||||
# * [Recommended Reading: Goodfellow et al chapter 3 on probability theory](https://www.deeplearningbook.org/)
|
||||
#
|
||||
# * See also Murphy, sections 2.4 (Gaussian distributions) and 3.2 (Bayesian Statistics, basis)
|
||||
#
|
||||
# * [Video of lecture](https://youtu.be/Kc20CFK0z7Y)
|
||||
#
|
||||
# * [Whiteboard notes](https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2023/NotesSep7.pdf)
|
||||
|
||||
# ## Material for the active learning sessions Tuesday and Wednesday
|
||||
#
|
||||
# The material here contains a summary from last Week and discussion of SVD, Ridge and Lasso regression with examples
|
||||
# The material here contains a summary from last week and discussion of SVD, Ridge and Lasso regression with examples
|
||||
|
||||
# ## Linear Regression and the SVD
|
||||
#
|
||||
@@ -844,6 +850,567 @@ plt.show()
|
||||
|
||||
# ## Material for lecture Thursday September 7
|
||||
|
||||
# ## Important technicalities: More on Rescaling data
|
||||
#
|
||||
# When you are comparing your own code with for example **Scikit-Learn**'s
|
||||
# library, there are some technicalities to keep in mind. The examples
|
||||
# here demonstrate some of these aspects with potential pitfalls.
|
||||
#
|
||||
# The discussion here focuses on the role of the intercept, how we can
|
||||
# set up the design matrix, what scaling we should use and other topics
|
||||
# which tend confuse us.
|
||||
#
|
||||
# The intercept can be interpreted as the expected value of our
|
||||
# target/output variables when all other predictors are set to zero.
|
||||
# Thus, if we cannot assume that the expected outputs/targets are zero
|
||||
# when all predictors are zero (the columns in the design matrix), it
|
||||
# may be a bad idea to implement a model which penalizes the intercept.
|
||||
# Furthermore, in for example Ridge and Lasso regression, the default solutions
|
||||
# from the library **Scikit-Learn** (when not shrinking $\beta_0$) for the unknown parameters
|
||||
# $\boldsymbol{\beta}$, are derived under the assumption that both $\boldsymbol{y}$ and
|
||||
# $\boldsymbol{X}$ are zero centered, that is we subtract the mean values.
|
||||
#
|
||||
# If our predictors represent different scales, then it is important to
|
||||
# standardize the design matrix $\boldsymbol{X}$ by subtracting the mean of each
|
||||
# column from the corresponding column and dividing the column with its
|
||||
# standard deviation. Most machine learning libraries do this as a default. This means that if you compare your code with the results from a given library,
|
||||
# the results may differ.
|
||||
#
|
||||
# The
|
||||
# [Standardscaler](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html)
|
||||
# function in **Scikit-Learn** does this for us. For the data sets we
|
||||
# have been studying in our various examples, the data are in many cases
|
||||
# already scaled and there is no need to scale them. You as a user of different machine learning algorithms, should always perform a
|
||||
# survey of your data, with a critical assessment of them in case you need to scale the data.
|
||||
#
|
||||
# If you need to scale the data, not doing so will give an *unfair*
|
||||
# penalization of the parameters since their magnitude depends on the
|
||||
# scale of their corresponding predictor.
|
||||
#
|
||||
# Suppose as an example that you
|
||||
# you have an input variable given by the heights of different persons.
|
||||
# Human height might be measured in inches or meters or
|
||||
# kilometers. If measured in kilometers, a standard linear regression
|
||||
# model with this predictor would probably give a much bigger
|
||||
# coefficient term, than if measured in millimeters.
|
||||
# This can clearly lead to problems in evaluating the cost/loss functions.
|
||||
#
|
||||
# Keep in mind that when you transform your data set before training a model, the same transformation needs to be done
|
||||
# on your eventual new data set before making a prediction. If we translate this into a Python code, it would could be implemented as
|
||||
|
||||
# In[7]:
|
||||
|
||||
|
||||
"""
|
||||
#Model training, we compute the mean value of y and X
|
||||
y_train_mean = np.mean(y_train)
|
||||
X_train_mean = np.mean(X_train,axis=0)
|
||||
X_train = X_train - X_train_mean
|
||||
y_train = y_train - y_train_mean
|
||||
|
||||
# The we fit our model with the training data
|
||||
trained_model = some_model.fit(X_train,y_train)
|
||||
|
||||
|
||||
#Model prediction, we need also to transform our data set used for the prediction.
|
||||
X_test = X_test - X_train_mean #Use mean from training data
|
||||
y_pred = trained_model(X_test)
|
||||
y_pred = y_pred + y_train_mean
|
||||
"""
|
||||
|
||||
|
||||
# Let us try to understand what this may imply mathematically when we
|
||||
# subtract the mean values, also known as *zero 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}}.
|
||||
# $$
|
||||
|
||||
# What does this mean? And why do we insist on all this? Let us look at some examples.
|
||||
#
|
||||
# This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (*code example thanks to Øyvind Sigmundson Schøyen*). Here our scaling of the data is done by subtracting the mean values only.
|
||||
# Note also that we do not split the data into training and test.
|
||||
|
||||
# In[8]:
|
||||
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from sklearn.linear_model import LinearRegression
|
||||
|
||||
|
||||
np.random.seed(2021)
|
||||
|
||||
def MSE(y_data,y_model):
|
||||
n = np.size(y_model)
|
||||
return np.sum((y_data-y_model)**2)/n
|
||||
|
||||
|
||||
def fit_beta(X, y):
|
||||
return np.linalg.pinv(X.T @ X) @ X.T @ y
|
||||
|
||||
|
||||
true_beta = [2, 0.5, 3.7]
|
||||
|
||||
x = np.linspace(0, 1, 11)
|
||||
y = np.sum(
|
||||
np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
|
||||
) + 0.1 * np.random.normal(size=len(x))
|
||||
|
||||
degree = 3
|
||||
X = np.zeros((len(x), degree))
|
||||
|
||||
# Include the intercept in the design matrix
|
||||
for p in range(degree):
|
||||
X[:, p] = x ** p
|
||||
|
||||
beta = fit_beta(X, y)
|
||||
|
||||
# Intercept is included in the design matrix
|
||||
skl = LinearRegression(fit_intercept=False).fit(X, y)
|
||||
|
||||
print(f"True beta: {true_beta}")
|
||||
print(f"Fitted beta: {beta}")
|
||||
print(f"Sklearn fitted beta: {skl.coef_}")
|
||||
ypredictOwn = X @ beta
|
||||
ypredictSKL = skl.predict(X)
|
||||
print(f"MSE with intercept column")
|
||||
print(MSE(y,ypredictOwn))
|
||||
print(f"MSE with intercept column from SKL")
|
||||
print(MSE(y,ypredictSKL))
|
||||
|
||||
|
||||
plt.figure()
|
||||
plt.scatter(x, y, label="Data")
|
||||
plt.plot(x, X @ beta, label="Fit")
|
||||
plt.plot(x, skl.predict(X), label="Sklearn (fit_intercept=False)")
|
||||
|
||||
|
||||
# Do not include the intercept in the design matrix
|
||||
X = np.zeros((len(x), degree - 1))
|
||||
|
||||
for p in range(degree - 1):
|
||||
X[:, p] = x ** (p + 1)
|
||||
|
||||
# Intercept is not included in the design matrix
|
||||
skl = LinearRegression(fit_intercept=True).fit(X, y)
|
||||
|
||||
# Use centered values for X and y when computing coefficients
|
||||
y_offset = np.average(y, axis=0)
|
||||
X_offset = np.average(X, axis=0)
|
||||
|
||||
beta = fit_beta(X - X_offset, y - y_offset)
|
||||
intercept = np.mean(y_offset - X_offset @ beta)
|
||||
|
||||
print(f"Manual intercept: {intercept}")
|
||||
print(f"Fitted beta (wiothout intercept): {beta}")
|
||||
print(f"Sklearn intercept: {skl.intercept_}")
|
||||
print(f"Sklearn fitted beta (without intercept): {skl.coef_}")
|
||||
ypredictOwn = X @ beta
|
||||
ypredictSKL = skl.predict(X)
|
||||
print(f"MSE with Manual intercept")
|
||||
print(MSE(y,ypredictOwn+intercept))
|
||||
print(f"MSE with Sklearn intercept")
|
||||
print(MSE(y,ypredictSKL))
|
||||
|
||||
plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
|
||||
plt.plot(x, skl.predict(X), "--", label="Sklearn (fit_intercept=True)")
|
||||
plt.grid()
|
||||
plt.legend()
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
# The intercept is the value of our output/target variable
|
||||
# when all our features are zero and our function crosses the $y$-axis (for a one-dimensional case).
|
||||
#
|
||||
# Printing the MSE, we see first that both methods give the same MSE, as
|
||||
# they should. However, when we move to for example Ridge regression,
|
||||
# the way we treat the intercept may give a larger or smaller MSE,
|
||||
# meaning that the MSE can be penalized by the value of the
|
||||
# intercept. Not including the intercept in the fit, means that the
|
||||
# regularization term does not include $\beta_0$. For different values
|
||||
# of $\lambda$, this may lead to different MSE values.
|
||||
#
|
||||
# To remind the reader, the regularization term, with the intercept in Ridge regression, is given by
|
||||
|
||||
# $$
|
||||
# \lambda \vert\vert \boldsymbol{\beta} \vert\vert_2^2 = \lambda \sum_{j=0}^{p-1}\beta_j^2,
|
||||
# $$
|
||||
|
||||
# but when we take out the intercept, this equation becomes
|
||||
|
||||
# $$
|
||||
# \lambda \vert\vert \boldsymbol{\beta} \vert\vert_2^2 = \lambda \sum_{j=1}^{p-1}\beta_j^2.
|
||||
# $$
|
||||
|
||||
# For Lasso regression we have
|
||||
|
||||
# $$
|
||||
# \lambda \vert\vert \boldsymbol{\beta} \vert\vert_1 = \lambda \sum_{j=1}^{p-1}\vert\beta_j\vert.
|
||||
# $$
|
||||
|
||||
# It means that, when scaling the design matrix and the outputs/targets,
|
||||
# by subtracting the mean values, we have an optimization problem which
|
||||
# is not penalized by the intercept. The MSE value can then be smaller
|
||||
# since it focuses only on the remaining quantities. If we however bring
|
||||
# back the intercept, we will get a MSE which then contains the
|
||||
# intercept.
|
||||
#
|
||||
# Armed with this wisdom, we attempt first to simply set the intercept equal to **False** in our implementation of Ridge regression for our well-known vanilla data set.
|
||||
|
||||
# In[9]:
|
||||
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.model_selection import train_test_split
|
||||
from sklearn import linear_model
|
||||
|
||||
def MSE(y_data,y_model):
|
||||
n = np.size(y_model)
|
||||
return np.sum((y_data-y_model)**2)/n
|
||||
|
||||
|
||||
# A seed just to ensure that the random numbers are the same for every run.
|
||||
# Useful for eventual debugging.
|
||||
np.random.seed(3155)
|
||||
|
||||
n = 100
|
||||
x = np.random.rand(n)
|
||||
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
|
||||
|
||||
Maxpolydegree = 20
|
||||
X = np.zeros((n,Maxpolydegree))
|
||||
#We include explicitely the intercept column
|
||||
for degree in range(Maxpolydegree):
|
||||
X[:,degree] = x**degree
|
||||
# We split the data in test and training data
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
|
||||
p = Maxpolydegree
|
||||
I = np.eye(p,p)
|
||||
# Decide which values of lambda to use
|
||||
nlambdas = 6
|
||||
MSEOwnRidgePredict = np.zeros(nlambdas)
|
||||
MSERidgePredict = np.zeros(nlambdas)
|
||||
lambdas = np.logspace(-4, 2, nlambdas)
|
||||
for i in range(nlambdas):
|
||||
lmb = lambdas[i]
|
||||
OwnRidgeBeta = np.linalg.pinv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
|
||||
# Note: we include the intercept column and no scaling
|
||||
RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
# and then make the prediction
|
||||
ytildeOwnRidge = X_train @ OwnRidgeBeta
|
||||
ypredictOwnRidge = X_test @ OwnRidgeBeta
|
||||
ytildeRidge = RegRidge.predict(X_train)
|
||||
ypredictRidge = RegRidge.predict(X_test)
|
||||
MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
|
||||
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
|
||||
print("Beta values for own Ridge implementation")
|
||||
print(OwnRidgeBeta)
|
||||
print("Beta values for Scikit-Learn Ridge implementation")
|
||||
print(RegRidge.coef_)
|
||||
print("MSE values for own Ridge implementation")
|
||||
print(MSEOwnRidgePredict[i])
|
||||
print("MSE values for Scikit-Learn Ridge implementation")
|
||||
print(MSERidgePredict[i])
|
||||
|
||||
# Now plot the results
|
||||
plt.figure()
|
||||
plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'r', label = 'MSE own Ridge Test')
|
||||
plt.plot(np.log10(lambdas), MSERidgePredict, 'g', label = 'MSE Ridge Test')
|
||||
|
||||
plt.xlabel('log10(lambda)')
|
||||
plt.ylabel('MSE')
|
||||
plt.legend()
|
||||
plt.show()
|
||||
|
||||
|
||||
# The results here agree when we force **Scikit-Learn**'s Ridge function to include the first column in our design matrix.
|
||||
# We see that the results agree very well. Here we have thus explicitely included the intercept column in the design matrix.
|
||||
# What happens if we do not include the intercept in our fit?
|
||||
# Let us see how we can change this code by zero centering.
|
||||
|
||||
# In[10]:
|
||||
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.model_selection import train_test_split
|
||||
from sklearn import linear_model
|
||||
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
|
||||
# A seed just to ensure that the random numbers are the same for every run.
|
||||
# Useful for eventual debugging.
|
||||
np.random.seed(315)
|
||||
|
||||
n = 100
|
||||
x = np.random.rand(n)
|
||||
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
|
||||
|
||||
Maxpolydegree = 20
|
||||
X = np.zeros((n,Maxpolydegree-1))
|
||||
|
||||
for degree in range(1,Maxpolydegree): #No intercept column
|
||||
X[:,degree-1] = x**(degree)
|
||||
|
||||
# We split the data in test and training data
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
|
||||
#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)
|
||||
#Remove the intercept from the training data.
|
||||
y_scaler = np.mean(y_train)
|
||||
y_train_scaled = y_train - y_scaler
|
||||
|
||||
p = Maxpolydegree-1
|
||||
I = np.eye(p,p)
|
||||
# Decide which values of lambda to use
|
||||
nlambdas = 6
|
||||
MSEOwnRidgePredict = np.zeros(nlambdas)
|
||||
MSERidgePredict = np.zeros(nlambdas)
|
||||
|
||||
lambdas = np.logspace(-4, 2, nlambdas)
|
||||
for i in range(nlambdas):
|
||||
lmb = lambdas[i]
|
||||
OwnRidgeBeta = np.linalg.pinv(X_train_scaled.T @ X_train_scaled+lmb*I) @ X_train_scaled.T @ (y_train_scaled)
|
||||
intercept_ = y_scaler - X_train_mean@OwnRidgeBeta #The intercept can be shifted so the model can predict on uncentered data
|
||||
#Add intercept to prediction
|
||||
ypredictOwnRidge = X_test_scaled @ OwnRidgeBeta + y_scaler
|
||||
RegRidge = linear_model.Ridge(lmb)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
ypredictRidge = RegRidge.predict(X_test)
|
||||
MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
|
||||
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
|
||||
print("Beta values for own Ridge implementation")
|
||||
print(OwnRidgeBeta) #Intercept is given by mean of target variable
|
||||
print("Beta values for Scikit-Learn Ridge implementation")
|
||||
print(RegRidge.coef_)
|
||||
print('Intercept from own implementation:')
|
||||
print(intercept_)
|
||||
print('Intercept from Scikit-Learn Ridge implementation')
|
||||
print(RegRidge.intercept_)
|
||||
print("MSE values for own Ridge implementation")
|
||||
print(MSEOwnRidgePredict[i])
|
||||
print("MSE values for Scikit-Learn Ridge implementation")
|
||||
print(MSERidgePredict[i])
|
||||
|
||||
|
||||
# Now plot the results
|
||||
plt.figure()
|
||||
plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'b--', label = 'MSE own Ridge Test')
|
||||
plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
|
||||
plt.xlabel('log10(lambda)')
|
||||
plt.ylabel('MSE')
|
||||
plt.legend()
|
||||
plt.show()
|
||||
|
||||
|
||||
# We see here, when compared to the code which includes explicitely the
|
||||
# intercept column, that our MSE value is actually smaller. This is
|
||||
# because the regularization term does not include the intercept value
|
||||
# $\beta_0$ in the fitting. This applies to Lasso regularization as
|
||||
# well. It means that our optimization is now done only with the
|
||||
# centered matrix and/or vector that enter the fitting procedure.
|
||||
|
||||
# ## Test Function for what happens with OLS, Ridge and Lasso
|
||||
#
|
||||
# Hitherto we have discussed Ridge and Lasso regression in terms of a
|
||||
# linear analysis. This may to many of you feel rather technical and
|
||||
# perhaps not that intuitive. The question is whether we can develop a
|
||||
# more intuitive way of understanding what Ridge and Lasso express.
|
||||
#
|
||||
# Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
|
||||
#
|
||||
# We will play around with a study of the values for the optimal
|
||||
# parameters $\boldsymbol{\beta}$ using OLS, Ridge and Lasso regression. For
|
||||
# OLS, you will notice as function of the noise and polynomial degree,
|
||||
# that the parameters $\beta$ will fluctuate from order to order in the
|
||||
# polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS.
|
||||
#
|
||||
# For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
|
||||
|
||||
# In[11]:
|
||||
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.model_selection import train_test_split
|
||||
from sklearn import linear_model
|
||||
|
||||
def R2(y_data, y_model):
|
||||
return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
|
||||
def MSE(y_data,y_model):
|
||||
n = np.size(y_model)
|
||||
return np.sum((y_data-y_model)**2)/n
|
||||
|
||||
# Make data set.
|
||||
n = 10000
|
||||
x = np.random.rand(n)
|
||||
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.randn(n)
|
||||
|
||||
Maxpolydegree = 5
|
||||
X = np.zeros((len(x),Maxpolydegree))
|
||||
X[:,0] = 1.0
|
||||
|
||||
|
||||
for polydegree in range(1,Maxpolydegree):
|
||||
X[:,polydegree] = x**(polydegree)
|
||||
|
||||
# We split the data in test and training data
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
|
||||
# matrix inversion to find beta
|
||||
OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train
|
||||
print(OLSbeta)
|
||||
ypredictOLS = X_test @ OLSbeta
|
||||
print("Test MSE OLS")
|
||||
print(MSE(y_test,ypredictOLS))
|
||||
# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn
|
||||
# Decide which values of lambda to use
|
||||
nlambdas = 4
|
||||
MSERidgePredict = np.zeros(nlambdas)
|
||||
MSELassoPredict = np.zeros(nlambdas)
|
||||
lambdas = np.logspace(-3, 1, nlambdas)
|
||||
for i in range(nlambdas):
|
||||
lmb = lambdas[i]
|
||||
# Make the fit using Ridge and Lasso
|
||||
RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
|
||||
RegLasso.fit(X_train,y_train)
|
||||
# and then make the prediction
|
||||
ypredictRidge = RegRidge.predict(X_test)
|
||||
ypredictLasso = RegLasso.predict(X_test)
|
||||
# Compute the MSE and print it
|
||||
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
|
||||
MSELassoPredict[i] = MSE(y_test,ypredictLasso)
|
||||
print(lmb,RegRidge.coef_)
|
||||
print(lmb,RegLasso.coef_)
|
||||
# Now plot the results
|
||||
plt.figure()
|
||||
plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
|
||||
plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test')
|
||||
plt.xlabel('log10(lambda)')
|
||||
plt.ylabel('MSE')
|
||||
plt.legend()
|
||||
plt.show()
|
||||
|
||||
|
||||
# How can we understand this?
|
||||
|
||||
# ## Linking the regression analysis with a statistical interpretation
|
||||
#
|
||||
# We will now couple the discussions of ordinary least squares, Ridge
|
||||
@@ -1227,95 +1794,6 @@ plt.show()
|
||||
|
||||
# ## Bayes' Theorem and Ridge and Lasso Regression
|
||||
#
|
||||
# Hitherto we have discussed Ridge and Lasso regression in terms of a
|
||||
# linear analysis. This may to many of you feel rather technical and
|
||||
# perhaps not that intuitive. The question is whether we can develop a
|
||||
# more intuitive way of understanding what Ridge and Lasso express.
|
||||
#
|
||||
# Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
|
||||
|
||||
# ## Test Function for what happens with OLS, Ridge and Lasso
|
||||
#
|
||||
# We will play around with a study of the values for the optimal
|
||||
# parameters $\boldsymbol{\beta}$ using OLS, Ridge and Lasso regression. For
|
||||
# OLS, you will notice as function of the noise and polynomial degree,
|
||||
# that the parameters $\beta$ will fluctuate from order to order in the
|
||||
# polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS.
|
||||
#
|
||||
# For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
|
||||
|
||||
# In[7]:
|
||||
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.model_selection import train_test_split
|
||||
from sklearn import linear_model
|
||||
|
||||
def R2(y_data, y_model):
|
||||
return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
|
||||
def MSE(y_data,y_model):
|
||||
n = np.size(y_model)
|
||||
return np.sum((y_data-y_model)**2)/n
|
||||
|
||||
# Make data set.
|
||||
n = 10000
|
||||
x = np.random.rand(n)
|
||||
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.randn(n)
|
||||
|
||||
Maxpolydegree = 5
|
||||
X = np.zeros((len(x),Maxpolydegree))
|
||||
X[:,0] = 1.0
|
||||
|
||||
for polydegree in range(1, Maxpolydegree):
|
||||
for degree in range(polydegree):
|
||||
X[:,degree] = x**(degree)
|
||||
|
||||
|
||||
# We split the data in test and training data
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
|
||||
# matrix inversion to find beta
|
||||
OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train
|
||||
print(OLSbeta)
|
||||
ypredictOLS = X_test @ OLSbeta
|
||||
print("Test MSE OLS")
|
||||
print(MSE(y_test,ypredictOLS))
|
||||
# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn
|
||||
# Decide which values of lambda to use
|
||||
nlambdas = 4
|
||||
MSERidgePredict = np.zeros(nlambdas)
|
||||
MSELassoPredict = np.zeros(nlambdas)
|
||||
lambdas = np.logspace(-3, 1, nlambdas)
|
||||
for i in range(nlambdas):
|
||||
lmb = lambdas[i]
|
||||
# Make the fit using Ridge and Lasso
|
||||
RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
|
||||
RegLasso.fit(X_train,y_train)
|
||||
# and then make the prediction
|
||||
ypredictRidge = RegRidge.predict(X_test)
|
||||
ypredictLasso = RegLasso.predict(X_test)
|
||||
# Compute the MSE and print it
|
||||
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
|
||||
MSELassoPredict[i] = MSE(y_test,ypredictLasso)
|
||||
print(lmb,RegRidge.coef_)
|
||||
print(lmb,RegLasso.coef_)
|
||||
# Now plot the results
|
||||
plt.figure()
|
||||
plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
|
||||
plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test')
|
||||
plt.xlabel('log10(lambda)')
|
||||
plt.ylabel('MSE')
|
||||
plt.legend()
|
||||
plt.show()
|
||||
|
||||
|
||||
# How can we understand this?
|
||||
|
||||
# ## Invoking Bayes' theorem
|
||||
#
|
||||
# Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression.
|
||||
#
|
||||
# For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\boldsymbol{D}$ (one-dimensional case)
|
||||
|
||||
@@ -2,7 +2,7 @@
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "cfebae03",
|
||||
"id": "1ff4a4e9",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -14,7 +14,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "353ebb07",
|
||||
"id": "0569f7a9",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -27,25 +27,34 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c02815dd",
|
||||
"id": "6f34e507",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"## Overarching aims of the exercises this week\n",
|
||||
"\n",
|
||||
"This exercise deals with various mean values and variances in linear regression method (here it may be useful to look up chapter 3, equation (3.8) of [Trevor Hastie, Robert Tibshirani, Jerome H. Friedman, The Elements of Statistical Learning, Springer](https://www.springer.com/gp/book/9780387848570)). The exercise is also a part of project 1 and can be reused in the theory part of the project.\n",
|
||||
"This exercise deals with various mean values and variances in linear\n",
|
||||
"regression method (here it may be useful to look up chapter 3,\n",
|
||||
"equation (3.8) of [Trevor Hastie, Robert Tibshirani, Jerome\n",
|
||||
"H. Friedman, The Elements of Statistical Learning,\n",
|
||||
"Springer](https://www.springer.com/gp/book/9780387848570)). The\n",
|
||||
"exercise is also a part of project 1 and can be reused in the theory\n",
|
||||
"part of the project.\n",
|
||||
"\n",
|
||||
"For more discussions on Ridge regression and calculation of expectation values, [Wessel van Wieringen's](https://arxiv.org/abs/1509.09169) article is highly recommended.\n",
|
||||
"For more discussions on Ridge regression and calculation of\n",
|
||||
"expectation values, [Wessel van\n",
|
||||
"Wieringen's](https://arxiv.org/abs/1509.09169) article is highly\n",
|
||||
"recommended.\n",
|
||||
"\n",
|
||||
"The assumption we have made is \n",
|
||||
"that there exists a continuous function $f(\\boldsymbol{x})$ and a normal distributed error $\\boldsymbol{\\varepsilon}\\sim N(0, \\sigma^2)$\n",
|
||||
"which describes our data"
|
||||
"The assumption we have made is that there exists a continuous function\n",
|
||||
"$f(\\boldsymbol{x})$ and a normal distributed error $\\boldsymbol{\\varepsilon}\\sim N(0,\n",
|
||||
"\\sigma^2)$ which describes our data"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "fb152043",
|
||||
"id": "3ff574a4",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -57,7 +66,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "9aba7a46",
|
||||
"id": "096a3165",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -68,7 +77,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "be52c097",
|
||||
"id": "cc9fb258",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -80,7 +89,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c54979e5",
|
||||
"id": "daf4eb5d",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -90,7 +99,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4b1eb811",
|
||||
"id": "3b491679",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -102,7 +111,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "036a5440",
|
||||
"id": "93a49739",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -114,7 +123,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "f7fb8094",
|
||||
"id": "a3fa80db",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -125,7 +134,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d7d92995",
|
||||
"id": "e6e2defd",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -137,7 +146,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "29eb9701",
|
||||
"id": "9aa19afa",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -150,7 +159,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "aa106511",
|
||||
"id": "5fa82a47",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -162,7 +171,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "58be1091",
|
||||
"id": "33099abb",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -172,7 +181,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4cb58e41",
|
||||
"id": "5d03f001",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -184,7 +193,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a04df28d",
|
||||
"id": "94bd151d",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -195,7 +204,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c095c96a",
|
||||
"id": "f0b84941",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -207,20 +216,20 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "5d425a79",
|
||||
"id": "60370a63",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"$$\n",
|
||||
"\\mathbb{E} \\big[ \\boldsymbol{\\beta}^{\\mathrm{Ridge}} \\big]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\\n",
|
||||
"\\mathbb{E} \\big[ \\hat{\\boldsymbol{\\beta}}^{\\mathrm{Ridge}} \\big]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\\n",
|
||||
"\\boldsymbol{\\beta}^{\\mathrm{OLS}}.\n",
|
||||
"$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "fb89cfab",
|
||||
"id": "af251217",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
@@ -233,19 +242,19 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "aea8801c",
|
||||
"id": "02e978da",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"$$\n",
|
||||
"\\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},\n",
|
||||
"\\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},\n",
|
||||
"$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "0e885693",
|
||||
"id": "3656a896",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
|
||||
+1594
-768
File diff suppressed because it is too large
Load Diff
@@ -6,14 +6,24 @@ DATE: Deadline is Sunday September 17 at midnight
|
||||
===== Overarching aims of the exercises this week =====
|
||||
|
||||
|
||||
This exercise deals with various mean values and variances in linear regression method (here it may be useful to look up chapter 3, equation (3.8) of "Trevor Hastie, Robert Tibshirani, Jerome H. Friedman, The Elements of Statistical Learning, Springer":"https://www.springer.com/gp/book/9780387848570"). The exercise is also a part of project 1 and can be reused in the theory part of the project.
|
||||
This exercise deals with various mean values and variances in linear
|
||||
regression method (here it may be useful to look up chapter 3,
|
||||
equation (3.8) of "Trevor Hastie, Robert Tibshirani, Jerome
|
||||
H. Friedman, The Elements of Statistical Learning,
|
||||
Springer":"https://www.springer.com/gp/book/9780387848570"). The
|
||||
exercise is also a part of project 1 and can be reused in the theory
|
||||
part of the project.
|
||||
|
||||
For more discussions on Ridge regression and calculation of expectation values, "Wessel van Wieringen's":"https://arxiv.org/abs/1509.09169" article is highly recommended.
|
||||
For more discussions on Ridge regression and calculation of
|
||||
expectation values, "Wessel van
|
||||
Wieringen's":"https://arxiv.org/abs/1509.09169" article is highly
|
||||
recommended.
|
||||
|
||||
|
||||
The assumption we have made is
|
||||
that there exists a continuous function $f(\bm{x})$ and a normal distributed error $\bm{\varepsilon}\sim N(0, \sigma^2)$
|
||||
which describes our data
|
||||
The assumption we have made is that there exists a continuous function
|
||||
$f(\bm{x})$ and a normal distributed error $\bm{\varepsilon}\sim N(0,
|
||||
\sigma^2)$ which describes our data
|
||||
|
||||
!bt
|
||||
\[
|
||||
\bm{y} = f(\bm{x})+\bm{\varepsilon}
|
||||
@@ -70,7 +80,7 @@ A given parameter $\beta_j$ is given by the diagonal matrix element of the above
|
||||
Show that
|
||||
!bt
|
||||
\[
|
||||
\mathbb{E} \big[ \bm{\beta}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\
|
||||
\mathbb{E} \big[ \hat{\bm{\beta}}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\
|
||||
\bm{\beta}^{\mathrm{OLS}}.
|
||||
\]
|
||||
!et
|
||||
@@ -82,7 +92,7 @@ Show also that the variance is
|
||||
|
||||
!bt
|
||||
\[
|
||||
\mbox{Var}[\bm{\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{\bm{\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},
|
||||
\]
|
||||
!et
|
||||
and it is easy to see that if the parameter $\lambda$ goes to infinity then the variance of Ridge parameters $\bm{\beta}$ goes to zero.
|
||||
|
||||
Reference in New Issue
Block a user