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

updating lectures

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2021-09-09 09:24:06 +02:00
parent dc49a4c888
commit f4d8b8ceca
7 changed files with 1189 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/pub/week36/html/week36-bs.html
+48 -41
View File
@@ -98,6 +98,11 @@ Automatically generated HTML file from DocOnce source
2,
None,
'simple-code-for-solving-the-above-problem'),
('With Lasso Regression', 2, None, 'with-lasso-regression'),
('Another Example, now with a polynomial fit',
2,
None,
'another-example-now-with-a-polynomial-fit'),
('Linking the regression analysis with a statistical '
'interpretation',
2,
@@ -281,46 +286,48 @@ MathJax.Hub.Config({
<!-- navigation toc: --> <li><a href="._week36-bs021.html#lasso-case" style="font-size: 80%;">Lasso case</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs022.html#the-first-case" style="font-size: 80%;">The first Case</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs023.html#simple-code-for-solving-the-above-problem" style="font-size: 80%;">Simple code for solving the above problem</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs024.html#linking-the-regression-analysis-with-a-statistical-interpretation" style="font-size: 80%;">Linking the regression analysis with a statistical interpretation</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs025.html#assumptions-made" style="font-size: 80%;">Assumptions made</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs026.html#expectation-value-and-variance" style="font-size: 80%;">Expectation value and variance</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs027.html#expectation-value-and-variance-for-boldsymbol-beta" style="font-size: 80%;">Expectation value and variance for \( \boldsymbol{\beta} \)</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs028.html#deriving-ols-from-a-probability-distribution" style="font-size: 80%;">Deriving OLS from a probability distribution</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs029.html#independent-and-identically-distrubuted-iid" style="font-size: 80%;">Independent and Identically Distrubuted (iid)</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs030.html#maximum-likelihood-estimation-mle" style="font-size: 80%;">Maximum Likelihood Estimation (MLE)</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs031.html#a-new-cost-function" style="font-size: 80%;">A new Cost Function</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs032.html#more-basic-statistics-and-bayes-theorem" style="font-size: 80%;">More basic Statistics and Bayes' theorem</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs033.html#marginal-probability" style="font-size: 80%;">Marginal Probability</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs034.html#conditional-probability" style="font-size: 80%;">Conditional Probability</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs035.html#bayes-theorem" style="font-size: 80%;">Bayes' Theorem</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs036.html#frliday-september-10" style="font-size: 80%;">Frliday September 10</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs040.html#why-resampling-methods" style="font-size: 80%;">Why resampling methods</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs042.html#resampling-methods" style="font-size: 80%;">Resampling methods</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs039.html#resampling-approaches-can-be-computationally-expensive" style="font-size: 80%;">Resampling approaches can be computationally expensive</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs040.html#why-resampling-methods" style="font-size: 80%;">Why resampling methods ?</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs041.html#statistical-analysis" style="font-size: 80%;">Statistical analysis</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs042.html#resampling-methods" style="font-size: 80%;">Resampling methods</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs043.html#resampling-methods-jackknife-and-bootstrap" style="font-size: 80%;">Resampling methods: Jackknife and Bootstrap</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs044.html#resampling-methods-jackknife" style="font-size: 80%;">Resampling methods: Jackknife</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs045.html#jackknife-code-example" style="font-size: 80%;">Jackknife code example</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs046.html#resampling-methods-bootstrap" style="font-size: 80%;">Resampling methods: Bootstrap</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs047.html#resampling-methods-bootstrap-background" style="font-size: 80%;">Resampling methods: Bootstrap background</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs048.html#resampling-methods-more-bootstrap-background" style="font-size: 80%;">Resampling methods: More Bootstrap background</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs049.html#resampling-methods-bootstrap-approach" style="font-size: 80%;">Resampling methods: Bootstrap approach</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs050.html#resampling-methods-bootstrap-steps" style="font-size: 80%;">Resampling methods: Bootstrap steps</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs051.html#code-example-for-the-bootstrap-method" style="font-size: 80%;">Code example for the Bootstrap method</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs052.html#various-steps-in-cross-validation" style="font-size: 80%;">Various steps in cross-validation</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs053.html#how-to-set-up-the-cross-validation-for-ridge-and-or-lasso" style="font-size: 80%;">How to set up the cross-validation for Ridge and/or Lasso</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs054.html#cross-validation-in-brief" style="font-size: 80%;">Cross-validation in brief</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs055.html#code-example-for-cross-validation-and-k-fold-cross-validation" style="font-size: 80%;">Code Example for Cross-validation and \( k \)-fold Cross-validation</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs056.html#the-bias-variance-tradeoff" style="font-size: 80%;">The bias-variance tradeoff</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs057.html#example-code-for-bias-variance-tradeoff" style="font-size: 80%;">Example code for Bias-Variance tradeoff</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs058.html#understanding-what-happens" style="font-size: 80%;">Understanding what happens</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs059.html#summing-up" style="font-size: 80%;">Summing up</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs060.html#another-example-from-scikit-learn-s-repository" style="font-size: 80%;">Another Example from Scikit-Learn's Repository</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs061.html#more-examples-on-bootstrap-and-cross-validation-and-errors" style="font-size: 80%;">More examples on bootstrap and cross-validation and errors</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs062.html#the-same-example-but-now-with-cross-validation" style="font-size: 80%;">The same example but now with cross-validation</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs063.html#cross-validation-with-ridge" style="font-size: 80%;">Cross-validation with Ridge</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs024.html#with-lasso-regression" style="font-size: 80%;">With Lasso Regression</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs025.html#another-example-now-with-a-polynomial-fit" style="font-size: 80%;">Another Example, now with a polynomial fit</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs026.html#linking-the-regression-analysis-with-a-statistical-interpretation" style="font-size: 80%;">Linking the regression analysis with a statistical interpretation</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs027.html#assumptions-made" style="font-size: 80%;">Assumptions made</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs028.html#expectation-value-and-variance" style="font-size: 80%;">Expectation value and variance</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs029.html#expectation-value-and-variance-for-boldsymbol-beta" style="font-size: 80%;">Expectation value and variance for \( \boldsymbol{\beta} \)</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs030.html#deriving-ols-from-a-probability-distribution" style="font-size: 80%;">Deriving OLS from a probability distribution</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs031.html#independent-and-identically-distrubuted-iid" style="font-size: 80%;">Independent and Identically Distrubuted (iid)</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs032.html#maximum-likelihood-estimation-mle" style="font-size: 80%;">Maximum Likelihood Estimation (MLE)</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs033.html#a-new-cost-function" style="font-size: 80%;">A new Cost Function</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs034.html#more-basic-statistics-and-bayes-theorem" style="font-size: 80%;">More basic Statistics and Bayes' theorem</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs035.html#marginal-probability" style="font-size: 80%;">Marginal Probability</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs036.html#conditional-probability" style="font-size: 80%;">Conditional Probability</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs037.html#bayes-theorem" style="font-size: 80%;">Bayes' Theorem</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs038.html#frliday-september-10" style="font-size: 80%;">Frliday September 10</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs042.html#why-resampling-methods" style="font-size: 80%;">Why resampling methods</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs044.html#resampling-methods" style="font-size: 80%;">Resampling methods</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs041.html#resampling-approaches-can-be-computationally-expensive" style="font-size: 80%;">Resampling approaches can be computationally expensive</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs042.html#why-resampling-methods" style="font-size: 80%;">Why resampling methods ?</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs043.html#statistical-analysis" style="font-size: 80%;">Statistical analysis</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs044.html#resampling-methods" style="font-size: 80%;">Resampling methods</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs045.html#resampling-methods-jackknife-and-bootstrap" style="font-size: 80%;">Resampling methods: Jackknife and Bootstrap</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs046.html#resampling-methods-jackknife" style="font-size: 80%;">Resampling methods: Jackknife</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs047.html#jackknife-code-example" style="font-size: 80%;">Jackknife code example</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs048.html#resampling-methods-bootstrap" style="font-size: 80%;">Resampling methods: Bootstrap</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs049.html#resampling-methods-bootstrap-background" style="font-size: 80%;">Resampling methods: Bootstrap background</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs050.html#resampling-methods-more-bootstrap-background" style="font-size: 80%;">Resampling methods: More Bootstrap background</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs051.html#resampling-methods-bootstrap-approach" style="font-size: 80%;">Resampling methods: Bootstrap approach</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs052.html#resampling-methods-bootstrap-steps" style="font-size: 80%;">Resampling methods: Bootstrap steps</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs053.html#code-example-for-the-bootstrap-method" style="font-size: 80%;">Code example for the Bootstrap method</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs054.html#various-steps-in-cross-validation" style="font-size: 80%;">Various steps in cross-validation</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs055.html#how-to-set-up-the-cross-validation-for-ridge-and-or-lasso" style="font-size: 80%;">How to set up the cross-validation for Ridge and/or Lasso</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs056.html#cross-validation-in-brief" style="font-size: 80%;">Cross-validation in brief</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs057.html#code-example-for-cross-validation-and-k-fold-cross-validation" style="font-size: 80%;">Code Example for Cross-validation and \( k \)-fold Cross-validation</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs058.html#the-bias-variance-tradeoff" style="font-size: 80%;">The bias-variance tradeoff</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs059.html#example-code-for-bias-variance-tradeoff" style="font-size: 80%;">Example code for Bias-Variance tradeoff</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs060.html#understanding-what-happens" style="font-size: 80%;">Understanding what happens</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs061.html#summing-up" style="font-size: 80%;">Summing up</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs062.html#another-example-from-scikit-learn-s-repository" style="font-size: 80%;">Another Example from Scikit-Learn's Repository</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs063.html#more-examples-on-bootstrap-and-cross-validation-and-errors" style="font-size: 80%;">More examples on bootstrap and cross-validation and errors</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs064.html#the-same-example-but-now-with-cross-validation" style="font-size: 80%;">The same example but now with cross-validation</a></li>
<!-- navigation toc: --> <li><a href="._week36-bs065.html#cross-validation-with-ridge" style="font-size: 80%;">Cross-validation with Ridge</a></li>
</ul>
</li>
@@ -379,7 +386,7 @@ MathJax.Hub.Config({
<li><a href="._week36-bs008.html">9</a></li>
<li><a href="._week36-bs009.html">10</a></li>
<li><a href="">...</a></li>
<li><a href="._week36-bs063.html">64</a></li>
<li><a href="._week36-bs065.html">66</a></li>
<li><a href="._week36-bs001.html">&raquo;</a></li>
</ul>
<!-- ------------------- end of main content --------------- -->
doc/pub/week36/html/week36-reveal.html
+223
View File
@@ -165,6 +165,15 @@ MathJax.Hub.Config({
<p><li> Thursday: Summary from last week on SVD, Statistics, probability theory and linear regression</li>
<p><li> Friday: Linear Regression and links with Statistics, Resampling methods and presentation of first project.</li>
</ul>
<p>
Recommended Reading:
<ol>
<p><li> Lectures on Regression</li>
<p><li> Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1</li>
<p><li> Hastie et al chapter 3</li>
</ol>
</section>
@@ -864,6 +873,220 @@ Here we set up the OLS, Ridge and Lasso functionality in order to study the abov
<p>
First we study and compare the OLS and Ridge results. The next code compares all three methods.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%;"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">os</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">pandas</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">pd</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">matplotlib.pyplot</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">plt</span>
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">R2</span>(y_data, y_model):
<span style="color: #8B008B; font-weight: bold">return</span> <span style="color: #B452CD">1</span> - np.sum((y_data - y_model) ** <span style="color: #B452CD">2</span>) / np.sum((y_data - np.mean(y_data)) ** <span style="color: #B452CD">2</span>)
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">MSE</span>(y_data,y_model):
n = np.size(y_model)
<span style="color: #8B008B; font-weight: bold">return</span> np.sum((y_data-y_model)**<span style="color: #B452CD">2</span>)/n
<span style="color: #228B22"># A seed just to ensure that the random numbers are the same for every run.</span>
<span style="color: #228B22"># Useful for eventual debugging.</span>
X = np.array( [ [ <span style="color: #B452CD">2</span>, <span style="color: #B452CD">0</span>], [<span style="color: #B452CD">0</span>, <span style="color: #B452CD">1</span>], [<span style="color: #B452CD">0</span>,<span style="color: #B452CD">0</span>]])
y = np.array( [<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>, <span style="color: #B452CD">3</span>])
<span style="color: #228B22"># matrix inversion to find beta</span>
OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
<span style="color: #658b00">print</span>(OLSbeta)
<span style="color: #228B22"># and then make the prediction</span>
ytildeOLS = X @ OLSbeta
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Training MSE for OLS&quot;</span>)
<span style="color: #658b00">print</span>(MSE(y,ytildeOLS))
ypredictOLS = X @ OLSbeta
<span style="color: #228B22"># Repeat now for Ridge regression and various values of the regularization parameter</span>
I = np.eye(<span style="color: #B452CD">2</span>,<span style="color: #B452CD">2</span>)
<span style="color: #228B22"># Decide which values of lambda to use</span>
nlambdas = <span style="color: #B452CD">100</span>
MSEPredict = np.zeros(nlambdas)
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">4</span>, nlambdas)
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
<span style="color: #228B22"># print(Ridgebeta)</span>
<span style="color: #228B22"># and then make the prediction</span>
ypredictRidge = X @ Ridgebeta
MSEPredict[i] = MSE(y,ypredictRidge)
<span style="color: #228B22"># print(MSEPredict[i])</span>
<span style="color: #228B22"># Now plot the results</span>
plt.figure()
plt.plot(np.log10(lambdas), MSEPredict, <span style="color: #CD5555">&#39;r--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE Ridge Test&#39;</span>)
plt.xlabel(<span style="color: #CD5555">&#39;log10(lambda)&#39;</span>)
plt.ylabel(<span style="color: #CD5555">&#39;MSE&#39;</span>)
plt.legend()
plt.show()
</pre></div>
<p>
We see here that we reach a plateau. What is actually happening?
</section>
<section>
<h2 id="with-lasso-regression">With Lasso Regression </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%;"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">os</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">pandas</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">pd</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">matplotlib.pyplot</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">plt</span>
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn</span> <span style="color: #8B008B; font-weight: bold">import</span> linear_model
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">R2</span>(y_data, y_model):
<span style="color: #8B008B; font-weight: bold">return</span> <span style="color: #B452CD">1</span> - np.sum((y_data - y_model) ** <span style="color: #B452CD">2</span>) / np.sum((y_data - np.mean(y_data)) ** <span style="color: #B452CD">2</span>)
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">MSE</span>(y_data,y_model):
n = np.size(y_model)
<span style="color: #8B008B; font-weight: bold">return</span> np.sum((y_data-y_model)**<span style="color: #B452CD">2</span>)/n
<span style="color: #228B22"># A seed just to ensure that the random numbers are the same for every run.</span>
<span style="color: #228B22"># Useful for eventual debugging.</span>
X = np.array( [ [ <span style="color: #B452CD">2</span>, <span style="color: #B452CD">0</span>], [<span style="color: #B452CD">0</span>, <span style="color: #B452CD">1</span>], [<span style="color: #B452CD">0</span>,<span style="color: #B452CD">0</span>]])
y = np.array( [<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>, <span style="color: #B452CD">3</span>])
<span style="color: #228B22"># matrix inversion to find beta</span>
OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
<span style="color: #658b00">print</span>(OLSbeta)
<span style="color: #228B22"># and then make the prediction</span>
ytildeOLS = X @ OLSbeta
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Training MSE for OLS&quot;</span>)
<span style="color: #658b00">print</span>(MSE(y,ytildeOLS))
ypredictOLS = X @ OLSbeta
<span style="color: #228B22"># Repeat now for Ridge regression and various values of the regularization parameter</span>
I = np.eye(<span style="color: #B452CD">2</span>,<span style="color: #B452CD">2</span>)
<span style="color: #228B22"># Decide which values of lambda to use</span>
nlambdas = <span style="color: #B452CD">100</span>
MSERidgePredict = np.zeros(nlambdas)
MSELassoPredict = np.zeros(nlambdas)
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">4</span>, nlambdas)
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
<span style="color: #228B22"># print(Ridgebeta)</span>
<span style="color: #228B22"># and then make the prediction</span>
ypredictRidge = X @ Ridgebeta
MSERidgePredict[i] = MSE(y,ypredictRidge)
<span style="color: #228B22"># print(MSEPredict[i])</span>
RegLasso = linear_model.Lasso(lmb)
RegLasso.fit(X,y)
ypredictLasso = RegLasso.predict(X)
MSELassoPredict[i] = MSE(y,ypredictLasso)
<span style="color: #228B22"># Now plot the results</span>
plt.figure()
plt.plot(np.log10(lambdas), MSERidgePredict, <span style="color: #CD5555">&#39;r--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE Ridge Train&#39;</span>)
plt.plot(np.log10(lambdas), MSELassoPredict, <span style="color: #CD5555">&#39;r--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE Lasso Train&#39;</span>)
plt.xlabel(<span style="color: #CD5555">&#39;log10(lambda)&#39;</span>)
plt.ylabel(<span style="color: #CD5555">&#39;MSE&#39;</span>)
plt.legend()
plt.show()
</pre></div>
</section>
<section>
<h2 id="another-example-now-with-a-polynomial-fit">Another Example, now with a polynomial fit </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%;"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">os</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">pandas</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">pd</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">matplotlib.pyplot</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">plt</span>
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn.model_selection</span> <span style="color: #8B008B; font-weight: bold">import</span> train_test_split
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn</span> <span style="color: #8B008B; font-weight: bold">import</span> linear_model
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">R2</span>(y_data, y_model):
<span style="color: #8B008B; font-weight: bold">return</span> <span style="color: #B452CD">1</span> - np.sum((y_data - y_model) ** <span style="color: #B452CD">2</span>) / np.sum((y_data - np.mean(y_data)) ** <span style="color: #B452CD">2</span>)
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">MSE</span>(y_data,y_model):
n = np.size(y_model)
<span style="color: #8B008B; font-weight: bold">return</span> np.sum((y_data-y_model)**<span style="color: #B452CD">2</span>)/n
<span style="color: #228B22"># A seed just to ensure that the random numbers are the same for every run.</span>
<span style="color: #228B22"># Useful for eventual debugging.</span>
np.random.seed(<span style="color: #B452CD">3155</span>)
x = np.random.rand(<span style="color: #B452CD">100</span>)
y = <span style="color: #B452CD">2.0</span>+<span style="color: #B452CD">5</span>*x*x+<span style="color: #B452CD">0.1</span>*np.random.randn(<span style="color: #B452CD">100</span>)
<span style="color: #228B22"># number of features p (here degree of polynomial</span>
p = <span style="color: #B452CD">3</span>
<span style="color: #228B22"># The design matrix now as function of a given polynomial</span>
X = np.zeros((<span style="color: #658b00">len</span>(x),p))
X[:,<span style="color: #B452CD">0</span>] = <span style="color: #B452CD">1.0</span>
X[:,<span style="color: #B452CD">1</span>] = x
X[:,<span style="color: #B452CD">2</span>] = x*x
<span style="color: #228B22"># We split the data in test and training data</span>
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=<span style="color: #B452CD">0.2</span>)
<span style="color: #228B22"># matrix inversion to find beta</span>
OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train
<span style="color: #658b00">print</span>(OLSbeta)
<span style="color: #228B22"># and then make the prediction</span>
ytildeOLS = X_train @ OLSbeta
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Training R2 for OLS&quot;</span>)
<span style="color: #658b00">print</span>(R2(y_train,ytildeOLS))
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Training MSE for OLS&quot;</span>)
<span style="color: #658b00">print</span>(MSE(y_train,ytildeOLS))
ypredictOLS = X_test @ OLSbeta
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Test R2 for OLS&quot;</span>)
<span style="color: #658b00">print</span>(R2(y_test,ypredictOLS))
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Test MSE OLS&quot;</span>)
<span style="color: #658b00">print</span>(MSE(y_test,ypredictOLS))
<span style="color: #228B22"># Repeat now for Lasso and Ridge regression and various values of the regularization parameter</span>
I = np.eye(p,p)
<span style="color: #228B22"># Decide which values of lambda to use</span>
nlambdas = <span style="color: #B452CD">100</span>
MSEPredict = np.zeros(nlambdas)
MSETrain = np.zeros(nlambdas)
MSELassoPredict = np.zeros(nlambdas)
MSELassoTrain = np.zeros(nlambdas)
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">4</span>, nlambdas)
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
<span style="color: #228B22"># include lasso using Scikit-Learn</span>
RegLasso = linear_model.Lasso(lmb)
RegLasso.fit(X_train,y_train)
<span style="color: #228B22"># and then make the prediction</span>
ytildeRidge = X_train @ Ridgebeta
ypredictRidge = X_test @ Ridgebeta
ytildeLasso = RegLasso.predict(X_train)
ypredictLasso = RegLasso.predict(X_test)
MSEPredict[i] = MSE(y_test,ypredictRidge)
MSETrain[i] = MSE(y_train,ytildeRidge)
MSELassoPredict[i] = MSE(y_test,ypredictLasso)
MSELassoTrain[i] = MSE(y_train,ytildeLasso)
<span style="color: #228B22"># Now plot the results</span>
plt.figure()
plt.plot(np.log10(lambdas), MSETrain, label = <span style="color: #CD5555">&#39;MSE Ridge train&#39;</span>)
plt.plot(np.log10(lambdas), MSEPredict, <span style="color: #CD5555">&#39;r--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE Ridge Test&#39;</span>)
plt.plot(np.log10(lambdas), MSELassoTrain, label = <span style="color: #CD5555">&#39;MSE Lasso train&#39;</span>)
plt.plot(np.log10(lambdas), MSELassoPredict, <span style="color: #CD5555">&#39;r--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE Lasso Test&#39;</span>)
plt.xlabel(<span style="color: #CD5555">&#39;log10(lambda)&#39;</span>)
plt.ylabel(<span style="color: #CD5555">&#39;MSE&#39;</span>)
plt.legend()
plt.show()
</pre></div>
</section>
doc/pub/week36/html/week36-solarized.html
+225
View File
@@ -118,6 +118,11 @@ div { text-align: justify; text-justify: inter-word; }
2,
None,
'simple-code-for-solving-the-above-problem'),
('With Lasso Regression', 2, None, 'with-lasso-regression'),
('Another Example, now with a polynomial fit',
2,
None,
'another-example-now-with-a-polynomial-fit'),
('Linking the regression analysis with a statistical '
'interpretation',
2,
@@ -294,6 +299,14 @@ MathJax.Hub.Config({
<li> Friday: Linear Regression and links with Statistics, Resampling methods and presentation of first project.</li>
</ul>
Recommended Reading:
<ol>
<li> Lectures on Regression</li>
<li> Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1</li>
<li> Hastie et al chapter 3</li>
</ol>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="thursday-september-9">Thursday September 9 </h2>
@@ -879,6 +892,218 @@ Here we set up the OLS, Ridge and Lasso functionality in order to study the abov
<p>
First we study and compare the OLS and Ridge results. The next code compares all three methods.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%;"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">os</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">pandas</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">pd</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">matplotlib.pyplot</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">plt</span>
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">R2</span>(y_data, y_model):
<span style="color: #8B008B; font-weight: bold">return</span> <span style="color: #B452CD">1</span> - np.sum((y_data - y_model) ** <span style="color: #B452CD">2</span>) / np.sum((y_data - np.mean(y_data)) ** <span style="color: #B452CD">2</span>)
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">MSE</span>(y_data,y_model):
n = np.size(y_model)
<span style="color: #8B008B; font-weight: bold">return</span> np.sum((y_data-y_model)**<span style="color: #B452CD">2</span>)/n
<span style="color: #228B22"># A seed just to ensure that the random numbers are the same for every run.</span>
<span style="color: #228B22"># Useful for eventual debugging.</span>
X = np.array( [ [ <span style="color: #B452CD">2</span>, <span style="color: #B452CD">0</span>], [<span style="color: #B452CD">0</span>, <span style="color: #B452CD">1</span>], [<span style="color: #B452CD">0</span>,<span style="color: #B452CD">0</span>]])
y = np.array( [<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>, <span style="color: #B452CD">3</span>])
<span style="color: #228B22"># matrix inversion to find beta</span>
OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
<span style="color: #658b00">print</span>(OLSbeta)
<span style="color: #228B22"># and then make the prediction</span>
ytildeOLS = X @ OLSbeta
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Training MSE for OLS&quot;</span>)
<span style="color: #658b00">print</span>(MSE(y,ytildeOLS))
ypredictOLS = X @ OLSbeta
<span style="color: #228B22"># Repeat now for Ridge regression and various values of the regularization parameter</span>
I = np.eye(<span style="color: #B452CD">2</span>,<span style="color: #B452CD">2</span>)
<span style="color: #228B22"># Decide which values of lambda to use</span>
nlambdas = <span style="color: #B452CD">100</span>
MSEPredict = np.zeros(nlambdas)
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">4</span>, nlambdas)
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
<span style="color: #228B22"># print(Ridgebeta)</span>
<span style="color: #228B22"># and then make the prediction</span>
ypredictRidge = X @ Ridgebeta
MSEPredict[i] = MSE(y,ypredictRidge)
<span style="color: #228B22"># print(MSEPredict[i])</span>
<span style="color: #228B22"># Now plot the results</span>
plt.figure()
plt.plot(np.log10(lambdas), MSEPredict, <span style="color: #CD5555">&#39;r--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE Ridge Test&#39;</span>)
plt.xlabel(<span style="color: #CD5555">&#39;log10(lambda)&#39;</span>)
plt.ylabel(<span style="color: #CD5555">&#39;MSE&#39;</span>)
plt.legend()
plt.show()
</pre></div>
<p>
We see here that we reach a plateau. What is actually happening?
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="with-lasso-regression">With Lasso Regression </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%;"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">os</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">pandas</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">pd</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">matplotlib.pyplot</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">plt</span>
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn</span> <span style="color: #8B008B; font-weight: bold">import</span> linear_model
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">R2</span>(y_data, y_model):
<span style="color: #8B008B; font-weight: bold">return</span> <span style="color: #B452CD">1</span> - np.sum((y_data - y_model) ** <span style="color: #B452CD">2</span>) / np.sum((y_data - np.mean(y_data)) ** <span style="color: #B452CD">2</span>)
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">MSE</span>(y_data,y_model):
n = np.size(y_model)
<span style="color: #8B008B; font-weight: bold">return</span> np.sum((y_data-y_model)**<span style="color: #B452CD">2</span>)/n
<span style="color: #228B22"># A seed just to ensure that the random numbers are the same for every run.</span>
<span style="color: #228B22"># Useful for eventual debugging.</span>
X = np.array( [ [ <span style="color: #B452CD">2</span>, <span style="color: #B452CD">0</span>], [<span style="color: #B452CD">0</span>, <span style="color: #B452CD">1</span>], [<span style="color: #B452CD">0</span>,<span style="color: #B452CD">0</span>]])
y = np.array( [<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>, <span style="color: #B452CD">3</span>])
<span style="color: #228B22"># matrix inversion to find beta</span>
OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
<span style="color: #658b00">print</span>(OLSbeta)
<span style="color: #228B22"># and then make the prediction</span>
ytildeOLS = X @ OLSbeta
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Training MSE for OLS&quot;</span>)
<span style="color: #658b00">print</span>(MSE(y,ytildeOLS))
ypredictOLS = X @ OLSbeta
<span style="color: #228B22"># Repeat now for Ridge regression and various values of the regularization parameter</span>
I = np.eye(<span style="color: #B452CD">2</span>,<span style="color: #B452CD">2</span>)
<span style="color: #228B22"># Decide which values of lambda to use</span>
nlambdas = <span style="color: #B452CD">100</span>
MSERidgePredict = np.zeros(nlambdas)
MSELassoPredict = np.zeros(nlambdas)
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">4</span>, nlambdas)
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
<span style="color: #228B22"># print(Ridgebeta)</span>
<span style="color: #228B22"># and then make the prediction</span>
ypredictRidge = X @ Ridgebeta
MSERidgePredict[i] = MSE(y,ypredictRidge)
<span style="color: #228B22"># print(MSEPredict[i])</span>
RegLasso = linear_model.Lasso(lmb)
RegLasso.fit(X,y)
ypredictLasso = RegLasso.predict(X)
MSELassoPredict[i] = MSE(y,ypredictLasso)
<span style="color: #228B22"># Now plot the results</span>
plt.figure()
plt.plot(np.log10(lambdas), MSERidgePredict, <span style="color: #CD5555">&#39;r--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE Ridge Train&#39;</span>)
plt.plot(np.log10(lambdas), MSELassoPredict, <span style="color: #CD5555">&#39;r--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE Lasso Train&#39;</span>)
plt.xlabel(<span style="color: #CD5555">&#39;log10(lambda)&#39;</span>)
plt.ylabel(<span style="color: #CD5555">&#39;MSE&#39;</span>)
plt.legend()
plt.show()
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="another-example-now-with-a-polynomial-fit">Another Example, now with a polynomial fit </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%;"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">os</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">pandas</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">pd</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">matplotlib.pyplot</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">plt</span>
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn.model_selection</span> <span style="color: #8B008B; font-weight: bold">import</span> train_test_split
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn</span> <span style="color: #8B008B; font-weight: bold">import</span> linear_model
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">R2</span>(y_data, y_model):
<span style="color: #8B008B; font-weight: bold">return</span> <span style="color: #B452CD">1</span> - np.sum((y_data - y_model) ** <span style="color: #B452CD">2</span>) / np.sum((y_data - np.mean(y_data)) ** <span style="color: #B452CD">2</span>)
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">MSE</span>(y_data,y_model):
n = np.size(y_model)
<span style="color: #8B008B; font-weight: bold">return</span> np.sum((y_data-y_model)**<span style="color: #B452CD">2</span>)/n
<span style="color: #228B22"># A seed just to ensure that the random numbers are the same for every run.</span>
<span style="color: #228B22"># Useful for eventual debugging.</span>
np.random.seed(<span style="color: #B452CD">3155</span>)
x = np.random.rand(<span style="color: #B452CD">100</span>)
y = <span style="color: #B452CD">2.0</span>+<span style="color: #B452CD">5</span>*x*x+<span style="color: #B452CD">0.1</span>*np.random.randn(<span style="color: #B452CD">100</span>)
<span style="color: #228B22"># number of features p (here degree of polynomial</span>
p = <span style="color: #B452CD">3</span>
<span style="color: #228B22"># The design matrix now as function of a given polynomial</span>
X = np.zeros((<span style="color: #658b00">len</span>(x),p))
X[:,<span style="color: #B452CD">0</span>] = <span style="color: #B452CD">1.0</span>
X[:,<span style="color: #B452CD">1</span>] = x
X[:,<span style="color: #B452CD">2</span>] = x*x
<span style="color: #228B22"># We split the data in test and training data</span>
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=<span style="color: #B452CD">0.2</span>)
<span style="color: #228B22"># matrix inversion to find beta</span>
OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train
<span style="color: #658b00">print</span>(OLSbeta)
<span style="color: #228B22"># and then make the prediction</span>
ytildeOLS = X_train @ OLSbeta
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Training R2 for OLS&quot;</span>)
<span style="color: #658b00">print</span>(R2(y_train,ytildeOLS))
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Training MSE for OLS&quot;</span>)
<span style="color: #658b00">print</span>(MSE(y_train,ytildeOLS))
ypredictOLS = X_test @ OLSbeta
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Test R2 for OLS&quot;</span>)
<span style="color: #658b00">print</span>(R2(y_test,ypredictOLS))
<span style="color: #658b00">print</span>(<span style="color: #CD5555">&quot;Test MSE OLS&quot;</span>)
<span style="color: #658b00">print</span>(MSE(y_test,ypredictOLS))
<span style="color: #228B22"># Repeat now for Lasso and Ridge regression and various values of the regularization parameter</span>
I = np.eye(p,p)
<span style="color: #228B22"># Decide which values of lambda to use</span>
nlambdas = <span style="color: #B452CD">100</span>
MSEPredict = np.zeros(nlambdas)
MSETrain = np.zeros(nlambdas)
MSELassoPredict = np.zeros(nlambdas)
MSELassoTrain = np.zeros(nlambdas)
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">4</span>, nlambdas)
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
<span style="color: #228B22"># include lasso using Scikit-Learn</span>
RegLasso = linear_model.Lasso(lmb)
RegLasso.fit(X_train,y_train)
<span style="color: #228B22"># and then make the prediction</span>
ytildeRidge = X_train @ Ridgebeta
ypredictRidge = X_test @ Ridgebeta
ytildeLasso = RegLasso.predict(X_train)
ypredictLasso = RegLasso.predict(X_test)
MSEPredict[i] = MSE(y_test,ypredictRidge)
MSETrain[i] = MSE(y_train,ytildeRidge)
MSELassoPredict[i] = MSE(y_test,ypredictLasso)
MSELassoTrain[i] = MSE(y_train,ytildeLasso)
<span style="color: #228B22"># Now plot the results</span>
plt.figure()
plt.plot(np.log10(lambdas), MSETrain, label = <span style="color: #CD5555">&#39;MSE Ridge train&#39;</span>)
plt.plot(np.log10(lambdas), MSEPredict, <span style="color: #CD5555">&#39;r--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE Ridge Test&#39;</span>)
plt.plot(np.log10(lambdas), MSELassoTrain, label = <span style="color: #CD5555">&#39;MSE Lasso train&#39;</span>)
plt.plot(np.log10(lambdas), MSELassoPredict, <span style="color: #CD5555">&#39;r--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE Lasso Test&#39;</span>)
plt.xlabel(<span style="color: #CD5555">&#39;log10(lambda)&#39;</span>)
plt.ylabel(<span style="color: #CD5555">&#39;MSE&#39;</span>)
plt.legend()
plt.show()
</pre></div>
<p>
<!-- !split -->
doc/pub/week36/html/week36.html
+225
View File
@@ -123,6 +123,11 @@ div { text-align: justify; text-justify: inter-word; }
2,
None,
'simple-code-for-solving-the-above-problem'),
('With Lasso Regression', 2, None, 'with-lasso-regression'),
('Another Example, now with a polynomial fit',
2,
None,
'another-example-now-with-a-polynomial-fit'),
('Linking the regression analysis with a statistical '
'interpretation',
2,
@@ -299,6 +304,14 @@ MathJax.Hub.Config({
<li> Friday: Linear Regression and links with Statistics, Resampling methods and presentation of first project.</li>
</ul>
Recommended Reading:
<ol>
<li> Lectures on Regression</li>
<li> Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1</li>
<li> Hastie et al chapter 3</li>
</ol>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="thursday-september-9">Thursday September 9 </h2>
@@ -884,6 +897,218 @@ Here we set up the OLS, Ridge and Lasso functionality in order to study the abov
<p>
First we study and compare the OLS and Ridge results. The next code compares all three methods.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%;"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">os</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">pandas</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">pd</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">matplotlib.pyplot</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">plt</span>
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">R2</span>(y_data, y_model):
<span style="color: #008000; font-weight: bold">return</span> <span style="color: #666666">1</span> <span style="color: #666666">-</span> np<span style="color: #666666">.</span>sum((y_data <span style="color: #666666">-</span> y_model) <span style="color: #666666">**</span> <span style="color: #666666">2</span>) <span style="color: #666666">/</span> np<span style="color: #666666">.</span>sum((y_data <span style="color: #666666">-</span> np<span style="color: #666666">.</span>mean(y_data)) <span style="color: #666666">**</span> <span style="color: #666666">2</span>)
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">MSE</span>(y_data,y_model):
n <span style="color: #666666">=</span> np<span style="color: #666666">.</span>size(y_model)
<span style="color: #008000; font-weight: bold">return</span> np<span style="color: #666666">.</span>sum((y_data<span style="color: #666666">-</span>y_model)<span style="color: #666666">**2</span>)<span style="color: #666666">/</span>n
<span style="color: #408080; font-style: italic"># A seed just to ensure that the random numbers are the same for every run.</span>
<span style="color: #408080; font-style: italic"># Useful for eventual debugging.</span>
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>array( [ [ <span style="color: #666666">2</span>, <span style="color: #666666">0</span>], [<span style="color: #666666">0</span>, <span style="color: #666666">1</span>], [<span style="color: #666666">0</span>,<span style="color: #666666">0</span>]])
y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>array( [<span style="color: #666666">4</span>, <span style="color: #666666">2</span>, <span style="color: #666666">3</span>])
<span style="color: #408080; font-style: italic"># matrix inversion to find beta</span>
OLSbeta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X) <span style="color: #666666">@</span> X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y
<span style="color: #008000">print</span>(OLSbeta)
<span style="color: #408080; font-style: italic"># and then make the prediction</span>
ytildeOLS <span style="color: #666666">=</span> X <span style="color: #666666">@</span> OLSbeta
<span style="color: #008000">print</span>(<span style="color: #BA2121">&quot;Training MSE for OLS&quot;</span>)
<span style="color: #008000">print</span>(MSE(y,ytildeOLS))
ypredictOLS <span style="color: #666666">=</span> X <span style="color: #666666">@</span> OLSbeta
<span style="color: #408080; font-style: italic"># Repeat now for Ridge regression and various values of the regularization parameter</span>
I <span style="color: #666666">=</span> np<span style="color: #666666">.</span>eye(<span style="color: #666666">2</span>,<span style="color: #666666">2</span>)
<span style="color: #408080; font-style: italic"># Decide which values of lambda to use</span>
nlambdas <span style="color: #666666">=</span> <span style="color: #666666">100</span>
MSEPredict <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(nlambdas)
lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-4</span>, <span style="color: #666666">4</span>, nlambdas)
<span style="color: #008000; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(nlambdas):
lmb <span style="color: #666666">=</span> lambdas[i]
Ridgebeta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y
<span style="color: #408080; font-style: italic"># print(Ridgebeta)</span>
<span style="color: #408080; font-style: italic"># and then make the prediction</span>
ypredictRidge <span style="color: #666666">=</span> X <span style="color: #666666">@</span> Ridgebeta
MSEPredict[i] <span style="color: #666666">=</span> MSE(y,ypredictRidge)
<span style="color: #408080; font-style: italic"># print(MSEPredict[i])</span>
<span style="color: #408080; font-style: italic"># Now plot the results</span>
plt<span style="color: #666666">.</span>figure()
plt<span style="color: #666666">.</span>plot(np<span style="color: #666666">.</span>log10(lambdas), MSEPredict, <span style="color: #BA2121">&#39;r--&#39;</span>, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;MSE Ridge Test&#39;</span>)
plt<span style="color: #666666">.</span>xlabel(<span style="color: #BA2121">&#39;log10(lambda)&#39;</span>)
plt<span style="color: #666666">.</span>ylabel(<span style="color: #BA2121">&#39;MSE&#39;</span>)
plt<span style="color: #666666">.</span>legend()
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
We see here that we reach a plateau. What is actually happening?
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="with-lasso-regression">With Lasso Regression </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%;"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">os</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">pandas</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">pd</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">matplotlib.pyplot</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">plt</span>
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">sklearn</span> <span style="color: #008000; font-weight: bold">import</span> linear_model
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">R2</span>(y_data, y_model):
<span style="color: #008000; font-weight: bold">return</span> <span style="color: #666666">1</span> <span style="color: #666666">-</span> np<span style="color: #666666">.</span>sum((y_data <span style="color: #666666">-</span> y_model) <span style="color: #666666">**</span> <span style="color: #666666">2</span>) <span style="color: #666666">/</span> np<span style="color: #666666">.</span>sum((y_data <span style="color: #666666">-</span> np<span style="color: #666666">.</span>mean(y_data)) <span style="color: #666666">**</span> <span style="color: #666666">2</span>)
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">MSE</span>(y_data,y_model):
n <span style="color: #666666">=</span> np<span style="color: #666666">.</span>size(y_model)
<span style="color: #008000; font-weight: bold">return</span> np<span style="color: #666666">.</span>sum((y_data<span style="color: #666666">-</span>y_model)<span style="color: #666666">**2</span>)<span style="color: #666666">/</span>n
<span style="color: #408080; font-style: italic"># A seed just to ensure that the random numbers are the same for every run.</span>
<span style="color: #408080; font-style: italic"># Useful for eventual debugging.</span>
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>array( [ [ <span style="color: #666666">2</span>, <span style="color: #666666">0</span>], [<span style="color: #666666">0</span>, <span style="color: #666666">1</span>], [<span style="color: #666666">0</span>,<span style="color: #666666">0</span>]])
y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>array( [<span style="color: #666666">4</span>, <span style="color: #666666">2</span>, <span style="color: #666666">3</span>])
<span style="color: #408080; font-style: italic"># matrix inversion to find beta</span>
OLSbeta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X) <span style="color: #666666">@</span> X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y
<span style="color: #008000">print</span>(OLSbeta)
<span style="color: #408080; font-style: italic"># and then make the prediction</span>
ytildeOLS <span style="color: #666666">=</span> X <span style="color: #666666">@</span> OLSbeta
<span style="color: #008000">print</span>(<span style="color: #BA2121">&quot;Training MSE for OLS&quot;</span>)
<span style="color: #008000">print</span>(MSE(y,ytildeOLS))
ypredictOLS <span style="color: #666666">=</span> X <span style="color: #666666">@</span> OLSbeta
<span style="color: #408080; font-style: italic"># Repeat now for Ridge regression and various values of the regularization parameter</span>
I <span style="color: #666666">=</span> np<span style="color: #666666">.</span>eye(<span style="color: #666666">2</span>,<span style="color: #666666">2</span>)
<span style="color: #408080; font-style: italic"># Decide which values of lambda to use</span>
nlambdas <span style="color: #666666">=</span> <span style="color: #666666">100</span>
MSERidgePredict <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(nlambdas)
MSELassoPredict <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(nlambdas)
lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-4</span>, <span style="color: #666666">4</span>, nlambdas)
<span style="color: #008000; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(nlambdas):
lmb <span style="color: #666666">=</span> lambdas[i]
Ridgebeta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y
<span style="color: #408080; font-style: italic"># print(Ridgebeta)</span>
<span style="color: #408080; font-style: italic"># and then make the prediction</span>
ypredictRidge <span style="color: #666666">=</span> X <span style="color: #666666">@</span> Ridgebeta
MSERidgePredict[i] <span style="color: #666666">=</span> MSE(y,ypredictRidge)
<span style="color: #408080; font-style: italic"># print(MSEPredict[i])</span>
RegLasso <span style="color: #666666">=</span> linear_model<span style="color: #666666">.</span>Lasso(lmb)
RegLasso<span style="color: #666666">.</span>fit(X,y)
ypredictLasso <span style="color: #666666">=</span> RegLasso<span style="color: #666666">.</span>predict(X)
MSELassoPredict[i] <span style="color: #666666">=</span> MSE(y,ypredictLasso)
<span style="color: #408080; font-style: italic"># Now plot the results</span>
plt<span style="color: #666666">.</span>figure()
plt<span style="color: #666666">.</span>plot(np<span style="color: #666666">.</span>log10(lambdas), MSERidgePredict, <span style="color: #BA2121">&#39;r--&#39;</span>, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;MSE Ridge Train&#39;</span>)
plt<span style="color: #666666">.</span>plot(np<span style="color: #666666">.</span>log10(lambdas), MSELassoPredict, <span style="color: #BA2121">&#39;r--&#39;</span>, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;MSE Lasso Train&#39;</span>)
plt<span style="color: #666666">.</span>xlabel(<span style="color: #BA2121">&#39;log10(lambda)&#39;</span>)
plt<span style="color: #666666">.</span>ylabel(<span style="color: #BA2121">&#39;MSE&#39;</span>)
plt<span style="color: #666666">.</span>legend()
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="another-example-now-with-a-polynomial-fit">Another Example, now with a polynomial fit </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%;"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">os</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">pandas</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">pd</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">matplotlib.pyplot</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">plt</span>
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">sklearn.model_selection</span> <span style="color: #008000; font-weight: bold">import</span> train_test_split
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">sklearn</span> <span style="color: #008000; font-weight: bold">import</span> linear_model
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">R2</span>(y_data, y_model):
<span style="color: #008000; font-weight: bold">return</span> <span style="color: #666666">1</span> <span style="color: #666666">-</span> np<span style="color: #666666">.</span>sum((y_data <span style="color: #666666">-</span> y_model) <span style="color: #666666">**</span> <span style="color: #666666">2</span>) <span style="color: #666666">/</span> np<span style="color: #666666">.</span>sum((y_data <span style="color: #666666">-</span> np<span style="color: #666666">.</span>mean(y_data)) <span style="color: #666666">**</span> <span style="color: #666666">2</span>)
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">MSE</span>(y_data,y_model):
n <span style="color: #666666">=</span> np<span style="color: #666666">.</span>size(y_model)
<span style="color: #008000; font-weight: bold">return</span> np<span style="color: #666666">.</span>sum((y_data<span style="color: #666666">-</span>y_model)<span style="color: #666666">**2</span>)<span style="color: #666666">/</span>n
<span style="color: #408080; font-style: italic"># A seed just to ensure that the random numbers are the same for every run.</span>
<span style="color: #408080; font-style: italic"># Useful for eventual debugging.</span>
np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>seed(<span style="color: #666666">3155</span>)
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(<span style="color: #666666">100</span>)
y <span style="color: #666666">=</span> <span style="color: #666666">2.0+5*</span>x<span style="color: #666666">*</span>x<span style="color: #666666">+0.1*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(<span style="color: #666666">100</span>)
<span style="color: #408080; font-style: italic"># number of features p (here degree of polynomial</span>
p <span style="color: #666666">=</span> <span style="color: #666666">3</span>
<span style="color: #408080; font-style: italic"># The design matrix now as function of a given polynomial</span>
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros((<span style="color: #008000">len</span>(x),p))
X[:,<span style="color: #666666">0</span>] <span style="color: #666666">=</span> <span style="color: #666666">1.0</span>
X[:,<span style="color: #666666">1</span>] <span style="color: #666666">=</span> x
X[:,<span style="color: #666666">2</span>] <span style="color: #666666">=</span> x<span style="color: #666666">*</span>x
<span style="color: #408080; font-style: italic"># We split the data in test and training data</span>
X_train, X_test, y_train, y_test <span style="color: #666666">=</span> train_test_split(X, y, test_size<span style="color: #666666">=0.2</span>)
<span style="color: #408080; font-style: italic"># matrix inversion to find beta</span>
OLSbeta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X_train) <span style="color: #666666">@</span> X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y_train
<span style="color: #008000">print</span>(OLSbeta)
<span style="color: #408080; font-style: italic"># and then make the prediction</span>
ytildeOLS <span style="color: #666666">=</span> X_train <span style="color: #666666">@</span> OLSbeta
<span style="color: #008000">print</span>(<span style="color: #BA2121">&quot;Training R2 for OLS&quot;</span>)
<span style="color: #008000">print</span>(R2(y_train,ytildeOLS))
<span style="color: #008000">print</span>(<span style="color: #BA2121">&quot;Training MSE for OLS&quot;</span>)
<span style="color: #008000">print</span>(MSE(y_train,ytildeOLS))
ypredictOLS <span style="color: #666666">=</span> X_test <span style="color: #666666">@</span> OLSbeta
<span style="color: #008000">print</span>(<span style="color: #BA2121">&quot;Test R2 for OLS&quot;</span>)
<span style="color: #008000">print</span>(R2(y_test,ypredictOLS))
<span style="color: #008000">print</span>(<span style="color: #BA2121">&quot;Test MSE OLS&quot;</span>)
<span style="color: #008000">print</span>(MSE(y_test,ypredictOLS))
<span style="color: #408080; font-style: italic"># Repeat now for Lasso and Ridge regression and various values of the regularization parameter</span>
I <span style="color: #666666">=</span> np<span style="color: #666666">.</span>eye(p,p)
<span style="color: #408080; font-style: italic"># Decide which values of lambda to use</span>
nlambdas <span style="color: #666666">=</span> <span style="color: #666666">100</span>
MSEPredict <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(nlambdas)
MSETrain <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(nlambdas)
MSELassoPredict <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(nlambdas)
MSELassoTrain <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(nlambdas)
lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-4</span>, <span style="color: #666666">4</span>, nlambdas)
<span style="color: #008000; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(nlambdas):
lmb <span style="color: #666666">=</span> lambdas[i]
Ridgebeta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X_train<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y_train
<span style="color: #408080; font-style: italic"># include lasso using Scikit-Learn</span>
RegLasso <span style="color: #666666">=</span> linear_model<span style="color: #666666">.</span>Lasso(lmb)
RegLasso<span style="color: #666666">.</span>fit(X_train,y_train)
<span style="color: #408080; font-style: italic"># and then make the prediction</span>
ytildeRidge <span style="color: #666666">=</span> X_train <span style="color: #666666">@</span> Ridgebeta
ypredictRidge <span style="color: #666666">=</span> X_test <span style="color: #666666">@</span> Ridgebeta
ytildeLasso <span style="color: #666666">=</span> RegLasso<span style="color: #666666">.</span>predict(X_train)
ypredictLasso <span style="color: #666666">=</span> RegLasso<span style="color: #666666">.</span>predict(X_test)
MSEPredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictRidge)
MSETrain[i] <span style="color: #666666">=</span> MSE(y_train,ytildeRidge)
MSELassoPredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictLasso)
MSELassoTrain[i] <span style="color: #666666">=</span> MSE(y_train,ytildeLasso)
<span style="color: #408080; font-style: italic"># Now plot the results</span>
plt<span style="color: #666666">.</span>figure()
plt<span style="color: #666666">.</span>plot(np<span style="color: #666666">.</span>log10(lambdas), MSETrain, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;MSE Ridge train&#39;</span>)
plt<span style="color: #666666">.</span>plot(np<span style="color: #666666">.</span>log10(lambdas), MSEPredict, <span style="color: #BA2121">&#39;r--&#39;</span>, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;MSE Ridge Test&#39;</span>)
plt<span style="color: #666666">.</span>plot(np<span style="color: #666666">.</span>log10(lambdas), MSELassoTrain, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;MSE Lasso train&#39;</span>)
plt<span style="color: #666666">.</span>plot(np<span style="color: #666666">.</span>log10(lambdas), MSELassoPredict, <span style="color: #BA2121">&#39;r--&#39;</span>, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;MSE Lasso Test&#39;</span>)
plt<span style="color: #666666">.</span>xlabel(<span style="color: #BA2121">&#39;log10(lambda)&#39;</span>)
plt<span style="color: #666666">.</span>ylabel(<span style="color: #BA2121">&#39;MSE&#39;</span>)
plt<span style="color: #666666">.</span>legend()
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
<!-- !split -->
doc/pub/week36/ipynb/ipynb-week36-src.tar.gz
BIN
View File
Binary file not shown.
doc/pub/week36/ipynb/week36.ipynb
+251 -3
View File
@@ -23,6 +23,13 @@
"\n",
"* Friday: Linear Regression and links with Statistics, Resampling methods and presentation of first project.\n",
"\n",
"Recommended Reading:\n",
"1. Lectures on Regression\n",
"\n",
"2. Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1\n",
"\n",
"3. Hastie et al chapter 3\n",
"\n",
"## Thursday September 9\n",
"\n",
"\n",
@@ -1092,9 +1099,252 @@
"\n",
"Here we set up the OLS, Ridge and Lasso functionality in order to study the above example. Note that here we have opted for a set of values of $\\lambda$, meaning that we need to perform a search in order to find the optimal values.\n",
"\n",
"First we study and compare the OLS and Ridge results. The next code compares all three methods.\n",
"First we study and compare the OLS and Ridge results. The next code compares all three methods."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import os\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def R2(y_data, y_model):\n",
" return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n",
"def MSE(y_data,y_model):\n",
" n = np.size(y_model)\n",
" return np.sum((y_data-y_model)**2)/n\n",
"\n",
"\n",
"# A seed just to ensure that the random numbers are the same for every run.\n",
"# Useful for eventual debugging.\n",
"\n",
"X = np.array( [ [ 2, 0], [0, 1], [0,0]])\n",
"y = np.array( [4, 2, 3])\n",
"\n",
"\n",
"# matrix inversion to find beta\n",
"OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y\n",
"print(OLSbeta)\n",
"# and then make the prediction\n",
"ytildeOLS = X @ OLSbeta\n",
"print(\"Training MSE for OLS\")\n",
"print(MSE(y,ytildeOLS))\n",
"ypredictOLS = X @ OLSbeta\n",
"\n",
"# Repeat now for Ridge regression and various values of the regularization parameter\n",
"I = np.eye(2,2)\n",
"# Decide which values of lambda to use\n",
"nlambdas = 100\n",
"MSEPredict = np.zeros(nlambdas)\n",
"lambdas = np.logspace(-4, 4, nlambdas)\n",
"for i in range(nlambdas):\n",
" lmb = lambdas[i]\n",
" Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y\n",
"# print(Ridgebeta)\n",
" # and then make the prediction\n",
" ypredictRidge = X @ Ridgebeta\n",
" MSEPredict[i] = MSE(y,ypredictRidge)\n",
"# print(MSEPredict[i])\n",
" # Now plot the results\n",
"plt.figure()\n",
"plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')\n",
"plt.xlabel('log10(lambda)')\n",
"plt.ylabel('MSE')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see here that we reach a plateau. What is actually happening?\n",
"\n",
"\n",
"## With Lasso Regression"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import os\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from sklearn import linear_model\n",
"\n",
"def R2(y_data, y_model):\n",
" return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n",
"def MSE(y_data,y_model):\n",
" n = np.size(y_model)\n",
" return np.sum((y_data-y_model)**2)/n\n",
"\n",
"\n",
"# A seed just to ensure that the random numbers are the same for every run.\n",
"# Useful for eventual debugging.\n",
"\n",
"X = np.array( [ [ 2, 0], [0, 1], [0,0]])\n",
"y = np.array( [4, 2, 3])\n",
"\n",
"\n",
"# matrix inversion to find beta\n",
"OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y\n",
"print(OLSbeta)\n",
"# and then make the prediction\n",
"ytildeOLS = X @ OLSbeta\n",
"print(\"Training MSE for OLS\")\n",
"print(MSE(y,ytildeOLS))\n",
"ypredictOLS = X @ OLSbeta\n",
"\n",
"# Repeat now for Ridge regression and various values of the regularization parameter\n",
"I = np.eye(2,2)\n",
"# Decide which values of lambda to use\n",
"nlambdas = 100\n",
"MSERidgePredict = np.zeros(nlambdas)\n",
"MSELassoPredict = np.zeros(nlambdas)\n",
"lambdas = np.logspace(-4, 4, nlambdas)\n",
"for i in range(nlambdas):\n",
" lmb = lambdas[i]\n",
" Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y\n",
"# print(Ridgebeta)\n",
" # and then make the prediction\n",
" ypredictRidge = X @ Ridgebeta\n",
" MSERidgePredict[i] = MSE(y,ypredictRidge)\n",
"# print(MSEPredict[i])\n",
" RegLasso = linear_model.Lasso(lmb)\n",
" RegLasso.fit(X,y)\n",
" ypredictLasso = RegLasso.predict(X)\n",
" MSELassoPredict[i] = MSE(y,ypredictLasso)\n",
"# Now plot the results\n",
"plt.figure()\n",
"plt.plot(np.log10(lambdas), MSERidgePredict, 'r--', label = 'MSE Ridge Train')\n",
"plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Train')\n",
"plt.xlabel('log10(lambda)')\n",
"plt.ylabel('MSE')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Another Example, now with a polynomial fit"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import os\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn import linear_model\n",
"\n",
"def R2(y_data, y_model):\n",
" return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n",
"def MSE(y_data,y_model):\n",
" n = np.size(y_model)\n",
" return np.sum((y_data-y_model)**2)/n\n",
"\n",
"\n",
"# A seed just to ensure that the random numbers are the same for every run.\n",
"# Useful for eventual debugging.\n",
"np.random.seed(3155)\n",
"\n",
"x = np.random.rand(100)\n",
"y = 2.0+5*x*x+0.1*np.random.randn(100)\n",
"\n",
"# number of features p (here degree of polynomial\n",
"p = 3\n",
"# The design matrix now as function of a given polynomial\n",
"X = np.zeros((len(x),p))\n",
"X[:,0] = 1.0\n",
"X[:,1] = x\n",
"X[:,2] = x*x\n",
"# We split the data in test and training data\n",
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"\n",
"# matrix inversion to find beta\n",
"OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train\n",
"print(OLSbeta)\n",
"# and then make the prediction\n",
"ytildeOLS = X_train @ OLSbeta\n",
"print(\"Training R2 for OLS\")\n",
"print(R2(y_train,ytildeOLS))\n",
"print(\"Training MSE for OLS\")\n",
"print(MSE(y_train,ytildeOLS))\n",
"ypredictOLS = X_test @ OLSbeta\n",
"print(\"Test R2 for OLS\")\n",
"print(R2(y_test,ypredictOLS))\n",
"print(\"Test MSE OLS\")\n",
"print(MSE(y_test,ypredictOLS))\n",
"\n",
"# Repeat now for Lasso and Ridge regression and various values of the regularization parameter\n",
"I = np.eye(p,p)\n",
"# Decide which values of lambda to use\n",
"nlambdas = 100\n",
"MSEPredict = np.zeros(nlambdas)\n",
"MSETrain = np.zeros(nlambdas)\n",
"MSELassoPredict = np.zeros(nlambdas)\n",
"MSELassoTrain = np.zeros(nlambdas)\n",
"lambdas = np.logspace(-4, 4, nlambdas)\n",
"for i in range(nlambdas):\n",
" lmb = lambdas[i]\n",
" Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train\n",
" # include lasso using Scikit-Learn\n",
" RegLasso = linear_model.Lasso(lmb)\n",
" RegLasso.fit(X_train,y_train)\n",
" # and then make the prediction\n",
" ytildeRidge = X_train @ Ridgebeta\n",
" ypredictRidge = X_test @ Ridgebeta\n",
" ytildeLasso = RegLasso.predict(X_train)\n",
" ypredictLasso = RegLasso.predict(X_test)\n",
" MSEPredict[i] = MSE(y_test,ypredictRidge)\n",
" MSETrain[i] = MSE(y_train,ytildeRidge)\n",
" MSELassoPredict[i] = MSE(y_test,ypredictLasso)\n",
" MSELassoTrain[i] = MSE(y_train,ytildeLasso)\n",
"\n",
"# Now plot the results\n",
"plt.figure()\n",
"plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train')\n",
"plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')\n",
"plt.plot(np.log10(lambdas), MSELassoTrain, label = 'MSE Lasso train')\n",
"plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Test')\n",
"\n",
"plt.xlabel('log10(lambda)')\n",
"plt.ylabel('MSE')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- !split -->\n",
"## Linking the regression analysis with a statistical interpretation\n",
"\n",
@@ -1961,8 +2211,6 @@
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"from numpy import *\n",
"from numpy.random import randint, randn\n",
"from time import time\n",
doc/src/week36/week36.do.txt
+217
View File
@@ -9,6 +9,11 @@ DATE: today
* Thursday: Summary from last week on SVD, Statistics, probability theory and linear regression
* Friday: Linear Regression and links with Statistics, Resampling methods and presentation of first project.
Recommended Reading:
o Lectures on Regression
o Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1
o Hastie et al chapter 3
!split
===== Thursday September 9 =====
@@ -594,6 +599,218 @@ Here we set up the OLS, Ridge and Lasso functionality in order to study the abov
First we study and compare the OLS and Ridge results. The next code compares all three methods.
!bc pycod
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
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
# A seed just to ensure that the random numbers are the same for every run.
# Useful for eventual debugging.
X = np.array( [ [ 2, 0], [0, 1], [0,0]])
y = np.array( [4, 2, 3])
# matrix inversion to find beta
OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
print(OLSbeta)
# and then make the prediction
ytildeOLS = X @ OLSbeta
print("Training MSE for OLS")
print(MSE(y,ytildeOLS))
ypredictOLS = X @ OLSbeta
# Repeat now for Ridge regression and various values of the regularization parameter
I = np.eye(2,2)
# Decide which values of lambda to use
nlambdas = 100
MSEPredict = np.zeros(nlambdas)
lambdas = np.logspace(-4, 4, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
# print(Ridgebeta)
# and then make the prediction
ypredictRidge = X @ Ridgebeta
MSEPredict[i] = MSE(y,ypredictRidge)
# print(MSEPredict[i])
# Now plot the results
plt.figure()
plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')
plt.xlabel('log10(lambda)')
plt.ylabel('MSE')
plt.legend()
plt.show()
!ec
We see here that we reach a plateau. What is actually happening?
!split
===== With Lasso Regression =====
!bc pycod
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
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
# A seed just to ensure that the random numbers are the same for every run.
# Useful for eventual debugging.
X = np.array( [ [ 2, 0], [0, 1], [0,0]])
y = np.array( [4, 2, 3])
# matrix inversion to find beta
OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
print(OLSbeta)
# and then make the prediction
ytildeOLS = X @ OLSbeta
print("Training MSE for OLS")
print(MSE(y,ytildeOLS))
ypredictOLS = X @ OLSbeta
# Repeat now for Ridge regression and various values of the regularization parameter
I = np.eye(2,2)
# Decide which values of lambda to use
nlambdas = 100
MSERidgePredict = np.zeros(nlambdas)
MSELassoPredict = np.zeros(nlambdas)
lambdas = np.logspace(-4, 4, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
# print(Ridgebeta)
# and then make the prediction
ypredictRidge = X @ Ridgebeta
MSERidgePredict[i] = MSE(y,ypredictRidge)
# print(MSEPredict[i])
RegLasso = linear_model.Lasso(lmb)
RegLasso.fit(X,y)
ypredictLasso = RegLasso.predict(X)
MSELassoPredict[i] = MSE(y,ypredictLasso)
# Now plot the results
plt.figure()
plt.plot(np.log10(lambdas), MSERidgePredict, 'r--', label = 'MSE Ridge Train')
plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Train')
plt.xlabel('log10(lambda)')
plt.ylabel('MSE')
plt.legend()
plt.show()
!ec
!split
===== Another Example, now with a polynomial fit =====
!bc pycod
import os
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 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
# A seed just to ensure that the random numbers are the same for every run.
# Useful for eventual debugging.
np.random.seed(3155)
x = np.random.rand(100)
y = 2.0+5*x*x+0.1*np.random.randn(100)
# number of features p (here degree of polynomial
p = 3
# The design matrix now as function of a given polynomial
X = np.zeros((len(x),p))
X[:,0] = 1.0
X[:,1] = x
X[:,2] = x*x
# 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.inv(X_train.T @ X_train) @ X_train.T @ y_train
print(OLSbeta)
# and then make the prediction
ytildeOLS = X_train @ OLSbeta
print("Training R2 for OLS")
print(R2(y_train,ytildeOLS))
print("Training MSE for OLS")
print(MSE(y_train,ytildeOLS))
ypredictOLS = X_test @ OLSbeta
print("Test R2 for OLS")
print(R2(y_test,ypredictOLS))
print("Test MSE OLS")
print(MSE(y_test,ypredictOLS))
# Repeat now for Lasso and Ridge regression and various values of the regularization parameter
I = np.eye(p,p)
# Decide which values of lambda to use
nlambdas = 100
MSEPredict = np.zeros(nlambdas)
MSETrain = np.zeros(nlambdas)
MSELassoPredict = np.zeros(nlambdas)
MSELassoTrain = np.zeros(nlambdas)
lambdas = np.logspace(-4, 4, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
# include lasso using Scikit-Learn
RegLasso = linear_model.Lasso(lmb)
RegLasso.fit(X_train,y_train)
# and then make the prediction
ytildeRidge = X_train @ Ridgebeta
ypredictRidge = X_test @ Ridgebeta
ytildeLasso = RegLasso.predict(X_train)
ypredictLasso = RegLasso.predict(X_test)
MSEPredict[i] = MSE(y_test,ypredictRidge)
MSETrain[i] = MSE(y_train,ytildeRidge)
MSELassoPredict[i] = MSE(y_test,ypredictLasso)
MSELassoTrain[i] = MSE(y_train,ytildeLasso)
# Now plot the results
plt.figure()
plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train')
plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')
plt.plot(np.log10(lambdas), MSELassoTrain, label = 'MSE Lasso train')
plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Test')
plt.xlabel('log10(lambda)')
plt.ylabel('MSE')
plt.legend()
plt.show()
!ec
!split
===== Linking the regression analysis with a statistical interpretation =====