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

updating week39

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2021-09-28 11:50:46 +02:00
parent 240f6a86cc
commit 233d259bf9
16 changed files with 1013 additions and 459 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/pub/week38/ipynb/week38.ipynb
+181 -127
View File
File diff suppressed because one or more lines are too long
doc/pub/week39/html/week39-bs.html
+54 -52
View File
@@ -48,6 +48,7 @@ Automatically generated HTML file from DocOnce source
None,
'searching-for-optimal-regularization-parameters-lambda'),
('Grid Search', 2, None, 'grid-search'),
('Randomized Grid Search', 2, None, 'randomized-grid-search'),
('Optimization, the central part of any Machine Learning '
'algortithm',
2,
@@ -218,56 +219,57 @@ MathJax.Hub.Config({
<!-- navigation toc: --> <li><a href="._week39-bs002.html#thursday-september-30" style="font-size: 80%;">Thursday September 30</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs003.html#searching-for-optimal-regularization-parameters-lambda" style="font-size: 80%;">Searching for Optimal Regularization Parameters \( \lambda \)</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs004.html#grid-search" style="font-size: 80%;">Grid Search</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs005.html#optimization-the-central-part-of-any-machine-learning-algortithm" style="font-size: 80%;">Optimization, the central part of any Machine Learning algortithm</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs006.html#revisiting-our-logistic-regression-case" style="font-size: 80%;">Revisiting our Logistic Regression case</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs007.html#the-equations-to-solve" style="font-size: 80%;">The equations to solve</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs008.html#solving-using-newton-raphson-s-method" style="font-size: 80%;">Solving using Newton-Raphson's method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs009.html#brief-reminder-on-newton-raphson-s-method" style="font-size: 80%;">Brief reminder on Newton-Raphson's method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs010.html#the-equations" style="font-size: 80%;">The equations</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs011.html#simple-geometric-interpretation" style="font-size: 80%;">Simple geometric interpretation</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs012.html#extending-to-more-than-one-variable" style="font-size: 80%;">Extending to more than one variable</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs013.html#steepest-descent" style="font-size: 80%;">Steepest descent</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs014.html#more-on-steepest-descent" style="font-size: 80%;">More on Steepest descent</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs015.html#the-ideal" style="font-size: 80%;">The ideal</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs016.html#the-sensitiveness-of-the-gradient-descent" style="font-size: 80%;">The sensitiveness of the gradient descent</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs017.html#convex-functions" style="font-size: 80%;">Convex functions</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs018.html#convex-function" style="font-size: 80%;">Convex function</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs019.html#conditions-on-convex-functions" style="font-size: 80%;">Conditions on convex functions</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs020.html#more-on-convex-functions" style="font-size: 80%;">More on convex functions</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs021.html#some-simple-problems" style="font-size: 80%;">Some simple problems</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs022.html#standard-steepest-descent" style="font-size: 80%;">Standard steepest descent</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs023.html#gradient-method" style="font-size: 80%;">Gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs025.html#steepest-descent-method" style="font-size: 80%;">Steepest descent method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs025.html#steepest-descent-method" style="font-size: 80%;">Steepest descent method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs026.html#final-expressions" style="font-size: 80%;">Final expressions</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs027.html#steepest-descent-example" style="font-size: 80%;">Steepest descent example</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs035.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs035.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs035.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs035.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs032.html#conjugate-gradient-method-and-iterations" style="font-size: 80%;">Conjugate gradient method and iterations</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs035.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs035.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs035.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs036.html#revisiting-our-first-homework" style="font-size: 80%;">Revisiting our first homework</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs041.html#gradient-descent-example" style="font-size: 80%;">Gradient descent example</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs038.html#the-derivative-of-the-cost-loss-function" style="font-size: 80%;">The derivative of the cost/loss function</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs039.html#the-hessian-matrix" style="font-size: 80%;">The Hessian matrix</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs040.html#simple-program" style="font-size: 80%;">Simple program</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs041.html#gradient-descent-example" style="font-size: 80%;">Gradient Descent Example</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs042.html#and-a-corresponding-example-using-_scikit-learn_" style="font-size: 80%;">And a corresponding example using <b>scikit-learn</b></a></li>
<!-- navigation toc: --> <li><a href="._week39-bs043.html#gradient-descent-and-ridge" style="font-size: 80%;">Gradient descent and Ridge</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs044.html#program-example-for-gradient-descent-with-ridge-regression" style="font-size: 80%;">Program example for gradient descent with Ridge Regression</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs045.html#using-gradient-descent-methods-limitations" style="font-size: 80%;">Using gradient descent methods, limitations</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs046.html#friday-october-1" style="font-size: 80%;">Friday October 1</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs047.html#stochastic-gradient-descent" style="font-size: 80%;">Stochastic Gradient Descent</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs048.html#computation-of-gradients" style="font-size: 80%;">Computation of gradients</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs049.html#sgd-example" style="font-size: 80%;">SGD example</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs050.html#the-gradient-step" style="font-size: 80%;">The gradient step</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs051.html#simple-example-code" style="font-size: 80%;">Simple example code</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs052.html#when-do-we-stop" style="font-size: 80%;">When do we stop?</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs053.html#slightly-different-approach" style="font-size: 80%;">Slightly different approach</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs054.html#program-for-stochastic-gradient" style="font-size: 80%;">Program for stochastic gradient</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs005.html#randomized-grid-search" style="font-size: 80%;">Randomized Grid Search</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs006.html#optimization-the-central-part-of-any-machine-learning-algortithm" style="font-size: 80%;">Optimization, the central part of any Machine Learning algortithm</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs007.html#revisiting-our-logistic-regression-case" style="font-size: 80%;">Revisiting our Logistic Regression case</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs008.html#the-equations-to-solve" style="font-size: 80%;">The equations to solve</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs009.html#solving-using-newton-raphson-s-method" style="font-size: 80%;">Solving using Newton-Raphson's method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs010.html#brief-reminder-on-newton-raphson-s-method" style="font-size: 80%;">Brief reminder on Newton-Raphson's method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs011.html#the-equations" style="font-size: 80%;">The equations</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs012.html#simple-geometric-interpretation" style="font-size: 80%;">Simple geometric interpretation</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs013.html#extending-to-more-than-one-variable" style="font-size: 80%;">Extending to more than one variable</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs014.html#steepest-descent" style="font-size: 80%;">Steepest descent</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs015.html#more-on-steepest-descent" style="font-size: 80%;">More on Steepest descent</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs016.html#the-ideal" style="font-size: 80%;">The ideal</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs017.html#the-sensitiveness-of-the-gradient-descent" style="font-size: 80%;">The sensitiveness of the gradient descent</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs018.html#convex-functions" style="font-size: 80%;">Convex functions</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs019.html#convex-function" style="font-size: 80%;">Convex function</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs020.html#conditions-on-convex-functions" style="font-size: 80%;">Conditions on convex functions</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs021.html#more-on-convex-functions" style="font-size: 80%;">More on convex functions</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs022.html#some-simple-problems" style="font-size: 80%;">Some simple problems</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs023.html#standard-steepest-descent" style="font-size: 80%;">Standard steepest descent</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs024.html#gradient-method" style="font-size: 80%;">Gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs026.html#steepest-descent-method" style="font-size: 80%;">Steepest descent method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs026.html#steepest-descent-method" style="font-size: 80%;">Steepest descent method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs027.html#final-expressions" style="font-size: 80%;">Final expressions</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs028.html#steepest-descent-example" style="font-size: 80%;">Steepest descent example</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs036.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs036.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs036.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs036.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs033.html#conjugate-gradient-method-and-iterations" style="font-size: 80%;">Conjugate gradient method and iterations</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs036.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs036.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs036.html#conjugate-gradient-method" style="font-size: 80%;">Conjugate gradient method</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs037.html#revisiting-our-first-homework" style="font-size: 80%;">Revisiting our first homework</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs042.html#gradient-descent-example" style="font-size: 80%;">Gradient descent example</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs039.html#the-derivative-of-the-cost-loss-function" style="font-size: 80%;">The derivative of the cost/loss function</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs040.html#the-hessian-matrix" style="font-size: 80%;">The Hessian matrix</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs041.html#simple-program" style="font-size: 80%;">Simple program</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs042.html#gradient-descent-example" style="font-size: 80%;">Gradient Descent Example</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs043.html#and-a-corresponding-example-using-_scikit-learn_" style="font-size: 80%;">And a corresponding example using <b>scikit-learn</b></a></li>
<!-- navigation toc: --> <li><a href="._week39-bs044.html#gradient-descent-and-ridge" style="font-size: 80%;">Gradient descent and Ridge</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs045.html#program-example-for-gradient-descent-with-ridge-regression" style="font-size: 80%;">Program example for gradient descent with Ridge Regression</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs046.html#using-gradient-descent-methods-limitations" style="font-size: 80%;">Using gradient descent methods, limitations</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs047.html#friday-october-1" style="font-size: 80%;">Friday October 1</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs048.html#stochastic-gradient-descent" style="font-size: 80%;">Stochastic Gradient Descent</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs049.html#computation-of-gradients" style="font-size: 80%;">Computation of gradients</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs050.html#sgd-example" style="font-size: 80%;">SGD example</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs051.html#the-gradient-step" style="font-size: 80%;">The gradient step</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs052.html#simple-example-code" style="font-size: 80%;">Simple example code</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs053.html#when-do-we-stop" style="font-size: 80%;">When do we stop?</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs054.html#slightly-different-approach" style="font-size: 80%;">Slightly different approach</a></li>
<!-- navigation toc: --> <li><a href="._week39-bs055.html#program-for-stochastic-gradient" style="font-size: 80%;">Program for stochastic gradient</a></li>
</ul>
</li>
@@ -302,7 +304,7 @@ MathJax.Hub.Config({
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Sep 27, 2021</h4></center> <!-- date -->
<center><h4>Sep 28, 2021</h4></center> <!-- date -->
<br>
<p>
@@ -326,7 +328,7 @@ MathJax.Hub.Config({
<li><a href="._week39-bs008.html">9</a></li>
<li><a href="._week39-bs009.html">10</a></li>
<li><a href="">...</a></li>
<li><a href="._week39-bs054.html">55</a></li>
<li><a href="._week39-bs055.html">56</a></li>
<li><a href="._week39-bs001.html">&raquo;</a></li>
</ul>
<!-- ------------------- end of main content --------------- -->
doc/pub/week39/html/week39-reveal.html
+133 -54
View File
@@ -148,7 +148,7 @@ MathJax.Hub.Config({
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>&nbsp;<br>
<center><h4>Sep 27, 2021</h4></center> <!-- date -->
<center><h4>Sep 28, 2021</h4></center> <!-- date -->
<br>
<p>
@@ -162,7 +162,7 @@ MathJax.Hub.Config({
<h2 id="plan-for-week-39">Plan for week 39 </h2>
<ul>
<p><li> Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1</li>
<p><li> Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1 and examples on how to implement Logistic Regression</li>
<p><li> Friday: Stochastic Gradient descent with examples and automatic differentiation</li>
<p><li> Reading recommendations:</li>
</ul>
@@ -196,44 +196,54 @@ which results in optimal scores (for example the smallest MSE or an \( R2=1 \)).
<!-- 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">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> KFold
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn.linear_model</span> <span style="color: #8B008B; font-weight: bold">import</span> Ridge
<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> cross_val_score
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn.preprocessing</span> <span style="color: #8B008B; font-weight: bold">import</span> PolynomialFeatures
<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">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>
np.random.seed(<span style="color: #B452CD">3155</span>)
<span style="color: #228B22"># Generate the data.</span>
<span style="color: #228B22"># Useful for eventual debugging.</span>
np.random.seed(<span style="color: #B452CD">315</span>)
n = <span style="color: #B452CD">100</span>
x = np.linspace(-<span style="color: #B452CD">3</span>, <span style="color: #B452CD">3</span>, n).reshape(-<span style="color: #B452CD">1</span>, <span style="color: #B452CD">1</span>)
y = np.exp(-x**<span style="color: #B452CD">2</span>) + <span style="color: #B452CD">1.5</span> * np.exp(-(x-<span style="color: #B452CD">2</span>)**<span style="color: #B452CD">2</span>)+ np.random.normal(<span style="color: #B452CD">0</span>, <span style="color: #B452CD">0.1</span>, x.shape)
<span style="color: #228B22"># Decide degree on polynomial to fit</span>
poly = PolynomialFeatures(degree = <span style="color: #B452CD">10</span>)
x = np.random.rand(n)
y = np.exp(-x**<span style="color: #B452CD">2</span>) + <span style="color: #B452CD">1.5</span> * np.exp(-(x-<span style="color: #B452CD">2</span>)**<span style="color: #B452CD">2</span>)
Maxpolydegree = <span style="color: #B452CD">5</span>
X = np.zeros((n,Maxpolydegree-<span style="color: #B452CD">1</span>))
<span style="color: #8B008B; font-weight: bold">for</span> degree <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">1</span>,Maxpolydegree): <span style="color: #228B22">#No intercept column</span>
X[:,degree-<span style="color: #B452CD">1</span>] = x**(degree)
<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"># Decide which values of lambda to use</span>
nlambdas = <span style="color: #B452CD">500</span>
lambdas = np.logspace(-<span style="color: #B452CD">3</span>, <span style="color: #B452CD">5</span>, nlambdas)
<span style="color: #228B22"># Initialize a KFold instance</span>
k = <span style="color: #B452CD">5</span>
kfold = KFold(n_splits = k)
estimated_mse_sklearn = np.zeros(nlambdas)
i = <span style="color: #B452CD">0</span>
<span style="color: #8B008B; font-weight: bold">for</span> lmb <span style="color: #8B008B">in</span> lambdas:
ridge = Ridge(alpha = lmb)
estimated_mse_folds = cross_val_score(ridge, x, y, scoring=<span style="color: #CD5555">&#39;neg_mean_squared_error&#39;</span>, cv=kfold)
estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
i += <span style="color: #B452CD">1</span>
MSERidgePredict = np.zeros(nlambdas)
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</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]
RegRidge = linear_model.Ridge(lmb)
RegRidge.fit(X_train,y_train)
ypredictRidge = RegRidge.predict(X_test)
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
<span style="color: #228B22"># Now plot the results</span>
plt.figure()
plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = <span style="color: #CD5555">&#39;cross_val_score&#39;</span>)
plt.plot(np.log10(lambdas), MSERidgePredict, <span style="color: #CD5555">&#39;g--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE SL 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 from this plot that the optimal MSE occurs for a value of \( \lambda\in [10,100] \).
In order to nail down the best value of \( \lambda \), we could in turn narrow down the search area.
Here we have performed a rather data greedy calculation as function of the regularization parameter \( \lambda \). There is no resampling here. The latter can easily be added by employing the function <b>RidgeCV</b> instead of just calling the <b>Ridge</b> function. For <b>RidgeCV</b> we need to passe the array of \( \lambda \) values.
By inspecting the figure we can in turn determine which is the optimal regularization parameter.
This becomes however less functional in the long run.
</section>
@@ -249,40 +259,109 @@ example here.
<!-- 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">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">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> KFold
<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.linear_model</span> <span style="color: #8B008B; font-weight: bold">import</span> Ridge
<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> cross_val_score
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn.preprocessing</span> <span style="color: #8B008B; font-weight: bold">import</span> PolynomialFeatures
<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> GridSearchCV
<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>
np.random.seed(<span style="color: #B452CD">3155</span>)
<span style="color: #228B22"># Generate the data.</span>
<span style="color: #228B22"># Useful for eventual debugging.</span>
np.random.seed(<span style="color: #B452CD">2021</span>)
n = <span style="color: #B452CD">100</span>
x = np.linspace(-<span style="color: #B452CD">3</span>, <span style="color: #B452CD">3</span>, n).reshape(-<span style="color: #B452CD">1</span>, <span style="color: #B452CD">1</span>)
y = np.exp(-x**<span style="color: #B452CD">2</span>) + <span style="color: #B452CD">1.5</span> * np.exp(-(x-<span style="color: #B452CD">2</span>)**<span style="color: #B452CD">2</span>)+ np.random.normal(<span style="color: #B452CD">0</span>, <span style="color: #B452CD">0.1</span>, x.shape)
<span style="color: #228B22"># Decide degree on polynomial to fit</span>
poly = PolynomialFeatures(degree = <span style="color: #B452CD">10</span>)
x = np.random.rand(n)
y = np.exp(-x**<span style="color: #B452CD">2</span>) + <span style="color: #B452CD">1.5</span> * np.exp(-(x-<span style="color: #B452CD">2</span>)**<span style="color: #B452CD">2</span>)
Maxpolydegree = <span style="color: #B452CD">5</span>
X = np.zeros((n,Maxpolydegree-<span style="color: #B452CD">1</span>))
<span style="color: #8B008B; font-weight: bold">for</span> degree <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">1</span>,Maxpolydegree): <span style="color: #228B22">#No intercept column</span>
X[:,degree-<span style="color: #B452CD">1</span>] = x**(degree)
<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"># Decide which values of lambda to use</span>
nlambdas = <span style="color: #B452CD">500</span>
lambdas = np.logspace(-<span style="color: #B452CD">3</span>, <span style="color: #B452CD">5</span>, nlambdas)
<span style="color: #228B22"># Initialize a KFold instance</span>
k = <span style="color: #B452CD">5</span>
kfold = KFold(n_splits = k)
estimated_mse_sklearn = np.zeros(nlambdas)
i = <span style="color: #B452CD">0</span>
<span style="color: #8B008B; font-weight: bold">for</span> lmb <span style="color: #8B008B">in</span> lambdas:
ridge = Ridge(alpha = lmb)
estimated_mse_folds = cross_val_score(ridge, x, y, scoring=<span style="color: #CD5555">&#39;neg_mean_squared_error&#39;</span>, cv=kfold)
estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
i += <span style="color: #B452CD">1</span>
plt.figure()
plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = <span style="color: #CD5555">&#39;cross_val_score&#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()
nlambdas = <span style="color: #B452CD">10</span>
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>, nlambdas)
<span style="color: #228B22"># create and fit a ridge regression model, testing each alpha</span>
model = Ridge()
gridsearch = GridSearchCV(estimator=model, param_grid=<span style="color: #658b00">dict</span>(alpha=lambdas))
gridsearch.fit(X_train, y_train)
<span style="color: #658b00">print</span>(gridsearch)
ypredictRidge = gridsearch.predict(X_test)
<span style="color: #228B22"># summarize the results of the grid search</span>
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;Best estimated lambda-value: {</span>gridsearch.best_estimator_.alpha<span style="color: #CD5555">}&quot;</span>)
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;MSE score: {</span>MSE(y_test,ypredictRidge)<span style="color: #CD5555">}&quot;</span>)
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;R2 score: {</span>R2(y_test,ypredictRidge)<span style="color: #CD5555">}&quot;</span>)
</pre></div>
<p>
By default the grid search function includes cross validation with five folds. The <a href="https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV" target="_blank">Scikit-Learn documentation</a> contains more information on how to set the different parameters.
</section>
<section>
<h2 id="randomized-grid-search">Randomized Grid Search </h2>
<p>
An alternative to the above manual grid set up, is to use a random
search where the parameters are tuned from a random distribution
(uniform below) for a fixed number of iterations. A model is
constructed and evaluated for each combination of chosen parameters.
We repeat the previous example but now with a random search.
<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">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">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.linear_model</span> <span style="color: #8B008B; font-weight: bold">import</span> Ridge
<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> GridSearchCV
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">scipy.stats</span> <span style="color: #8B008B; font-weight: bold">import</span> uniform <span style="color: #8B008B; font-weight: bold">as</span> randuniform
<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> RandomizedSearchCV
<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">2021</span>)
n = <span style="color: #B452CD">100</span>
x = np.random.rand(n)
y = np.exp(-x**<span style="color: #B452CD">2</span>) + <span style="color: #B452CD">1.5</span> * np.exp(-(x-<span style="color: #B452CD">2</span>)**<span style="color: #B452CD">2</span>)
Maxpolydegree = <span style="color: #B452CD">5</span>
X = np.zeros((n,Maxpolydegree-<span style="color: #B452CD">1</span>))
<span style="color: #8B008B; font-weight: bold">for</span> degree <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">1</span>,Maxpolydegree): <span style="color: #228B22">#No intercept column</span>
X[:,degree-<span style="color: #B452CD">1</span>] = x**(degree)
<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>)
param_grid = {<span style="color: #CD5555">&#39;alpha&#39;</span>: randuniform()}
<span style="color: #228B22"># create and fit a ridge regression model, testing each alpha</span>
model = Ridge()
gridsearch = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=<span style="color: #B452CD">100</span>)
gridsearch.fit(X_train, y_train)
<span style="color: #658b00">print</span>(gridsearch)
ypredictRidge = gridsearch.predict(X_test)
<span style="color: #228B22"># summarize the results of the grid search</span>
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;Best estimated lambda-value: {</span>gridsearch.best_estimator_.alpha<span style="color: #CD5555">}&quot;</span>)
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;MSE score: {</span>MSE(y_test,ypredictRidge)<span style="color: #CD5555">}&quot;</span>)
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;R2 score: {</span>R2(y_test,ypredictRidge)<span style="color: #CD5555">}&quot;</span>)
</pre></div>
</section>
doc/pub/week39/html/week39-solarized.html
+134 -54
View File
@@ -68,6 +68,7 @@ div { text-align: justify; text-justify: inter-word; }
None,
'searching-for-optimal-regularization-parameters-lambda'),
('Grid Search', 2, None, 'grid-search'),
('Randomized Grid Search', 2, None, 'randomized-grid-search'),
('Optimization, the central part of any Machine Learning '
'algortithm',
2,
@@ -238,7 +239,7 @@ MathJax.Hub.Config({
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Sep 27, 2021</h4></center> <!-- date -->
<center><h4>Sep 28, 2021</h4></center> <!-- date -->
<br>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
@@ -246,7 +247,7 @@ MathJax.Hub.Config({
<h2 id="plan-for-week-39">Plan for week 39 </h2>
<ul>
<li> Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1</li>
<li> Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1 and examples on how to implement Logistic Regression</li>
<li> Friday: Stochastic Gradient descent with examples and automatic differentiation</li>
<li> Reading recommendations:</li>
</ul>
@@ -279,44 +280,54 @@ which results in optimal scores (for example the smallest MSE or an \( R2=1 \)).
<!-- 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">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> KFold
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn.linear_model</span> <span style="color: #8B008B; font-weight: bold">import</span> Ridge
<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> cross_val_score
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn.preprocessing</span> <span style="color: #8B008B; font-weight: bold">import</span> PolynomialFeatures
<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">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>
np.random.seed(<span style="color: #B452CD">3155</span>)
<span style="color: #228B22"># Generate the data.</span>
<span style="color: #228B22"># Useful for eventual debugging.</span>
np.random.seed(<span style="color: #B452CD">315</span>)
n = <span style="color: #B452CD">100</span>
x = np.linspace(-<span style="color: #B452CD">3</span>, <span style="color: #B452CD">3</span>, n).reshape(-<span style="color: #B452CD">1</span>, <span style="color: #B452CD">1</span>)
y = np.exp(-x**<span style="color: #B452CD">2</span>) + <span style="color: #B452CD">1.5</span> * np.exp(-(x-<span style="color: #B452CD">2</span>)**<span style="color: #B452CD">2</span>)+ np.random.normal(<span style="color: #B452CD">0</span>, <span style="color: #B452CD">0.1</span>, x.shape)
<span style="color: #228B22"># Decide degree on polynomial to fit</span>
poly = PolynomialFeatures(degree = <span style="color: #B452CD">10</span>)
x = np.random.rand(n)
y = np.exp(-x**<span style="color: #B452CD">2</span>) + <span style="color: #B452CD">1.5</span> * np.exp(-(x-<span style="color: #B452CD">2</span>)**<span style="color: #B452CD">2</span>)
Maxpolydegree = <span style="color: #B452CD">5</span>
X = np.zeros((n,Maxpolydegree-<span style="color: #B452CD">1</span>))
<span style="color: #8B008B; font-weight: bold">for</span> degree <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">1</span>,Maxpolydegree): <span style="color: #228B22">#No intercept column</span>
X[:,degree-<span style="color: #B452CD">1</span>] = x**(degree)
<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"># Decide which values of lambda to use</span>
nlambdas = <span style="color: #B452CD">500</span>
lambdas = np.logspace(-<span style="color: #B452CD">3</span>, <span style="color: #B452CD">5</span>, nlambdas)
<span style="color: #228B22"># Initialize a KFold instance</span>
k = <span style="color: #B452CD">5</span>
kfold = KFold(n_splits = k)
estimated_mse_sklearn = np.zeros(nlambdas)
i = <span style="color: #B452CD">0</span>
<span style="color: #8B008B; font-weight: bold">for</span> lmb <span style="color: #8B008B">in</span> lambdas:
ridge = Ridge(alpha = lmb)
estimated_mse_folds = cross_val_score(ridge, x, y, scoring=<span style="color: #CD5555">&#39;neg_mean_squared_error&#39;</span>, cv=kfold)
estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
i += <span style="color: #B452CD">1</span>
MSERidgePredict = np.zeros(nlambdas)
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</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]
RegRidge = linear_model.Ridge(lmb)
RegRidge.fit(X_train,y_train)
ypredictRidge = RegRidge.predict(X_test)
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
<span style="color: #228B22"># Now plot the results</span>
plt.figure()
plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = <span style="color: #CD5555">&#39;cross_val_score&#39;</span>)
plt.plot(np.log10(lambdas), MSERidgePredict, <span style="color: #CD5555">&#39;g--&#39;</span>, label = <span style="color: #CD5555">&#39;MSE SL 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 from this plot that the optimal MSE occurs for a value of \( \lambda\in [10,100] \).
In order to nail down the best value of \( \lambda \), we could in turn narrow down the search area.
Here we have performed a rather data greedy calculation as function of the regularization parameter \( \lambda \). There is no resampling here. The latter can easily be added by employing the function <b>RidgeCV</b> instead of just calling the <b>Ridge</b> function. For <b>RidgeCV</b> we need to passe the array of \( \lambda \) values.
By inspecting the figure we can in turn determine which is the optimal regularization parameter.
This becomes however less functional in the long run.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
@@ -332,40 +343,109 @@ example here.
<!-- 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">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">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> KFold
<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.linear_model</span> <span style="color: #8B008B; font-weight: bold">import</span> Ridge
<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> cross_val_score
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn.preprocessing</span> <span style="color: #8B008B; font-weight: bold">import</span> PolynomialFeatures
<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> GridSearchCV
<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>
np.random.seed(<span style="color: #B452CD">3155</span>)
<span style="color: #228B22"># Generate the data.</span>
<span style="color: #228B22"># Useful for eventual debugging.</span>
np.random.seed(<span style="color: #B452CD">2021</span>)
n = <span style="color: #B452CD">100</span>
x = np.linspace(-<span style="color: #B452CD">3</span>, <span style="color: #B452CD">3</span>, n).reshape(-<span style="color: #B452CD">1</span>, <span style="color: #B452CD">1</span>)
y = np.exp(-x**<span style="color: #B452CD">2</span>) + <span style="color: #B452CD">1.5</span> * np.exp(-(x-<span style="color: #B452CD">2</span>)**<span style="color: #B452CD">2</span>)+ np.random.normal(<span style="color: #B452CD">0</span>, <span style="color: #B452CD">0.1</span>, x.shape)
<span style="color: #228B22"># Decide degree on polynomial to fit</span>
poly = PolynomialFeatures(degree = <span style="color: #B452CD">10</span>)
x = np.random.rand(n)
y = np.exp(-x**<span style="color: #B452CD">2</span>) + <span style="color: #B452CD">1.5</span> * np.exp(-(x-<span style="color: #B452CD">2</span>)**<span style="color: #B452CD">2</span>)
Maxpolydegree = <span style="color: #B452CD">5</span>
X = np.zeros((n,Maxpolydegree-<span style="color: #B452CD">1</span>))
<span style="color: #8B008B; font-weight: bold">for</span> degree <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">1</span>,Maxpolydegree): <span style="color: #228B22">#No intercept column</span>
X[:,degree-<span style="color: #B452CD">1</span>] = x**(degree)
<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"># Decide which values of lambda to use</span>
nlambdas = <span style="color: #B452CD">500</span>
lambdas = np.logspace(-<span style="color: #B452CD">3</span>, <span style="color: #B452CD">5</span>, nlambdas)
<span style="color: #228B22"># Initialize a KFold instance</span>
k = <span style="color: #B452CD">5</span>
kfold = KFold(n_splits = k)
estimated_mse_sklearn = np.zeros(nlambdas)
i = <span style="color: #B452CD">0</span>
<span style="color: #8B008B; font-weight: bold">for</span> lmb <span style="color: #8B008B">in</span> lambdas:
ridge = Ridge(alpha = lmb)
estimated_mse_folds = cross_val_score(ridge, x, y, scoring=<span style="color: #CD5555">&#39;neg_mean_squared_error&#39;</span>, cv=kfold)
estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
i += <span style="color: #B452CD">1</span>
plt.figure()
plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = <span style="color: #CD5555">&#39;cross_val_score&#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()
nlambdas = <span style="color: #B452CD">10</span>
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>, nlambdas)
<span style="color: #228B22"># create and fit a ridge regression model, testing each alpha</span>
model = Ridge()
gridsearch = GridSearchCV(estimator=model, param_grid=<span style="color: #658b00">dict</span>(alpha=lambdas))
gridsearch.fit(X_train, y_train)
<span style="color: #658b00">print</span>(gridsearch)
ypredictRidge = gridsearch.predict(X_test)
<span style="color: #228B22"># summarize the results of the grid search</span>
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;Best estimated lambda-value: {</span>gridsearch.best_estimator_.alpha<span style="color: #CD5555">}&quot;</span>)
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;MSE score: {</span>MSE(y_test,ypredictRidge)<span style="color: #CD5555">}&quot;</span>)
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;R2 score: {</span>R2(y_test,ypredictRidge)<span style="color: #CD5555">}&quot;</span>)
</pre></div>
<p>
By default the grid search function includes cross validation with five folds. The <a href="https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV" target="_blank">Scikit-Learn documentation</a> contains more information on how to set the different parameters.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="randomized-grid-search">Randomized Grid Search </h2>
<p>
An alternative to the above manual grid set up, is to use a random
search where the parameters are tuned from a random distribution
(uniform below) for a fixed number of iterations. A model is
constructed and evaluated for each combination of chosen parameters.
We repeat the previous example but now with a random search.
<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">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">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.linear_model</span> <span style="color: #8B008B; font-weight: bold">import</span> Ridge
<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> GridSearchCV
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">scipy.stats</span> <span style="color: #8B008B; font-weight: bold">import</span> uniform <span style="color: #8B008B; font-weight: bold">as</span> randuniform
<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> RandomizedSearchCV
<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">2021</span>)
n = <span style="color: #B452CD">100</span>
x = np.random.rand(n)
y = np.exp(-x**<span style="color: #B452CD">2</span>) + <span style="color: #B452CD">1.5</span> * np.exp(-(x-<span style="color: #B452CD">2</span>)**<span style="color: #B452CD">2</span>)
Maxpolydegree = <span style="color: #B452CD">5</span>
X = np.zeros((n,Maxpolydegree-<span style="color: #B452CD">1</span>))
<span style="color: #8B008B; font-weight: bold">for</span> degree <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">1</span>,Maxpolydegree): <span style="color: #228B22">#No intercept column</span>
X[:,degree-<span style="color: #B452CD">1</span>] = x**(degree)
<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>)
param_grid = {<span style="color: #CD5555">&#39;alpha&#39;</span>: randuniform()}
<span style="color: #228B22"># create and fit a ridge regression model, testing each alpha</span>
model = Ridge()
gridsearch = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=<span style="color: #B452CD">100</span>)
gridsearch.fit(X_train, y_train)
<span style="color: #658b00">print</span>(gridsearch)
ypredictRidge = gridsearch.predict(X_test)
<span style="color: #228B22"># summarize the results of the grid search</span>
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;Best estimated lambda-value: {</span>gridsearch.best_estimator_.alpha<span style="color: #CD5555">}&quot;</span>)
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;MSE score: {</span>MSE(y_test,ypredictRidge)<span style="color: #CD5555">}&quot;</span>)
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f&quot;R2 score: {</span>R2(y_test,ypredictRidge)<span style="color: #CD5555">}&quot;</span>)
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
doc/pub/week39/html/week39.html
+134 -54
View File
@@ -73,6 +73,7 @@ div { text-align: justify; text-justify: inter-word; }
None,
'searching-for-optimal-regularization-parameters-lambda'),
('Grid Search', 2, None, 'grid-search'),
('Randomized Grid Search', 2, None, 'randomized-grid-search'),
('Optimization, the central part of any Machine Learning '
'algortithm',
2,
@@ -243,7 +244,7 @@ MathJax.Hub.Config({
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Sep 27, 2021</h4></center> <!-- date -->
<center><h4>Sep 28, 2021</h4></center> <!-- date -->
<br>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
@@ -251,7 +252,7 @@ MathJax.Hub.Config({
<h2 id="plan-for-week-39">Plan for week 39 </h2>
<ul>
<li> Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1</li>
<li> Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1 and examples on how to implement Logistic Regression</li>
<li> Friday: Stochastic Gradient descent with examples and automatic differentiation</li>
<li> Reading recommendations:</li>
</ul>
@@ -284,44 +285,54 @@ which results in optimal scores (for example the smallest MSE or an \( R2=1 \)).
<!-- 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">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> KFold
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">sklearn.linear_model</span> <span style="color: #008000; font-weight: bold">import</span> Ridge
<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> cross_val_score
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">sklearn.preprocessing</span> <span style="color: #008000; font-weight: bold">import</span> PolynomialFeatures
<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">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>
np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>seed(<span style="color: #666666">3155</span>)
<span style="color: #408080; font-style: italic"># Generate the data.</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">315</span>)
n <span style="color: #666666">=</span> <span style="color: #666666">100</span>
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linspace(<span style="color: #666666">-3</span>, <span style="color: #666666">3</span>, n)<span style="color: #666666">.</span>reshape(<span style="color: #666666">-1</span>, <span style="color: #666666">1</span>)
y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>x<span style="color: #666666">**2</span>) <span style="color: #666666">+</span> <span style="color: #666666">1.5</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>(x<span style="color: #666666">-2</span>)<span style="color: #666666">**2</span>)<span style="color: #666666">+</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>normal(<span style="color: #666666">0</span>, <span style="color: #666666">0.1</span>, x<span style="color: #666666">.</span>shape)
<span style="color: #408080; font-style: italic"># Decide degree on polynomial to fit</span>
poly <span style="color: #666666">=</span> PolynomialFeatures(degree <span style="color: #666666">=</span> <span style="color: #666666">10</span>)
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(n)
y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>x<span style="color: #666666">**2</span>) <span style="color: #666666">+</span> <span style="color: #666666">1.5</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>(x<span style="color: #666666">-2</span>)<span style="color: #666666">**2</span>)
Maxpolydegree <span style="color: #666666">=</span> <span style="color: #666666">5</span>
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros((n,Maxpolydegree<span style="color: #666666">-1</span>))
<span style="color: #008000; font-weight: bold">for</span> degree <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(<span style="color: #666666">1</span>,Maxpolydegree): <span style="color: #408080; font-style: italic">#No intercept column</span>
X[:,degree<span style="color: #666666">-1</span>] <span style="color: #666666">=</span> x<span style="color: #666666">**</span>(degree)
<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"># Decide which values of lambda to use</span>
nlambdas <span style="color: #666666">=</span> <span style="color: #666666">500</span>
lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-3</span>, <span style="color: #666666">5</span>, nlambdas)
<span style="color: #408080; font-style: italic"># Initialize a KFold instance</span>
k <span style="color: #666666">=</span> <span style="color: #666666">5</span>
kfold <span style="color: #666666">=</span> KFold(n_splits <span style="color: #666666">=</span> k)
estimated_mse_sklearn <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(nlambdas)
i <span style="color: #666666">=</span> <span style="color: #666666">0</span>
<span style="color: #008000; font-weight: bold">for</span> lmb <span style="color: #AA22FF; font-weight: bold">in</span> lambdas:
ridge <span style="color: #666666">=</span> Ridge(alpha <span style="color: #666666">=</span> lmb)
estimated_mse_folds <span style="color: #666666">=</span> cross_val_score(ridge, x, y, scoring<span style="color: #666666">=</span><span style="color: #BA2121">&#39;neg_mean_squared_error&#39;</span>, cv<span style="color: #666666">=</span>kfold)
estimated_mse_sklearn[i] <span style="color: #666666">=</span> np<span style="color: #666666">.</span>mean(<span style="color: #666666">-</span>estimated_mse_folds)
i <span style="color: #666666">+=</span> <span style="color: #666666">1</span>
MSERidgePredict <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">2</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]
RegRidge <span style="color: #666666">=</span> linear_model<span style="color: #666666">.</span>Ridge(lmb)
RegRidge<span style="color: #666666">.</span>fit(X_train,y_train)
ypredictRidge <span style="color: #666666">=</span> RegRidge<span style="color: #666666">.</span>predict(X_test)
MSERidgePredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictRidge)
<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), estimated_mse_sklearn, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;cross_val_score&#39;</span>)
plt<span style="color: #666666">.</span>plot(np<span style="color: #666666">.</span>log10(lambdas), MSERidgePredict, <span style="color: #BA2121">&#39;g--&#39;</span>, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;MSE SL 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 from this plot that the optimal MSE occurs for a value of \( \lambda\in [10,100] \).
In order to nail down the best value of \( \lambda \), we could in turn narrow down the search area.
Here we have performed a rather data greedy calculation as function of the regularization parameter \( \lambda \). There is no resampling here. The latter can easily be added by employing the function <b>RidgeCV</b> instead of just calling the <b>Ridge</b> function. For <b>RidgeCV</b> we need to passe the array of \( \lambda \) values.
By inspecting the figure we can in turn determine which is the optimal regularization parameter.
This becomes however less functional in the long run.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
@@ -337,40 +348,109 @@ example here.
<!-- 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">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">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> KFold
<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.linear_model</span> <span style="color: #008000; font-weight: bold">import</span> Ridge
<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> cross_val_score
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">sklearn.preprocessing</span> <span style="color: #008000; font-weight: bold">import</span> PolynomialFeatures
<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> GridSearchCV
<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>
np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>seed(<span style="color: #666666">3155</span>)
<span style="color: #408080; font-style: italic"># Generate the data.</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">2021</span>)
n <span style="color: #666666">=</span> <span style="color: #666666">100</span>
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linspace(<span style="color: #666666">-3</span>, <span style="color: #666666">3</span>, n)<span style="color: #666666">.</span>reshape(<span style="color: #666666">-1</span>, <span style="color: #666666">1</span>)
y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>x<span style="color: #666666">**2</span>) <span style="color: #666666">+</span> <span style="color: #666666">1.5</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>(x<span style="color: #666666">-2</span>)<span style="color: #666666">**2</span>)<span style="color: #666666">+</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>normal(<span style="color: #666666">0</span>, <span style="color: #666666">0.1</span>, x<span style="color: #666666">.</span>shape)
<span style="color: #408080; font-style: italic"># Decide degree on polynomial to fit</span>
poly <span style="color: #666666">=</span> PolynomialFeatures(degree <span style="color: #666666">=</span> <span style="color: #666666">10</span>)
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(n)
y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>x<span style="color: #666666">**2</span>) <span style="color: #666666">+</span> <span style="color: #666666">1.5</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>(x<span style="color: #666666">-2</span>)<span style="color: #666666">**2</span>)
Maxpolydegree <span style="color: #666666">=</span> <span style="color: #666666">5</span>
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros((n,Maxpolydegree<span style="color: #666666">-1</span>))
<span style="color: #008000; font-weight: bold">for</span> degree <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(<span style="color: #666666">1</span>,Maxpolydegree): <span style="color: #408080; font-style: italic">#No intercept column</span>
X[:,degree<span style="color: #666666">-1</span>] <span style="color: #666666">=</span> x<span style="color: #666666">**</span>(degree)
<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"># Decide which values of lambda to use</span>
nlambdas <span style="color: #666666">=</span> <span style="color: #666666">500</span>
lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-3</span>, <span style="color: #666666">5</span>, nlambdas)
<span style="color: #408080; font-style: italic"># Initialize a KFold instance</span>
k <span style="color: #666666">=</span> <span style="color: #666666">5</span>
kfold <span style="color: #666666">=</span> KFold(n_splits <span style="color: #666666">=</span> k)
estimated_mse_sklearn <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(nlambdas)
i <span style="color: #666666">=</span> <span style="color: #666666">0</span>
<span style="color: #008000; font-weight: bold">for</span> lmb <span style="color: #AA22FF; font-weight: bold">in</span> lambdas:
ridge <span style="color: #666666">=</span> Ridge(alpha <span style="color: #666666">=</span> lmb)
estimated_mse_folds <span style="color: #666666">=</span> cross_val_score(ridge, x, y, scoring<span style="color: #666666">=</span><span style="color: #BA2121">&#39;neg_mean_squared_error&#39;</span>, cv<span style="color: #666666">=</span>kfold)
estimated_mse_sklearn[i] <span style="color: #666666">=</span> np<span style="color: #666666">.</span>mean(<span style="color: #666666">-</span>estimated_mse_folds)
i <span style="color: #666666">+=</span> <span style="color: #666666">1</span>
plt<span style="color: #666666">.</span>figure()
plt<span style="color: #666666">.</span>plot(np<span style="color: #666666">.</span>log10(lambdas), estimated_mse_sklearn, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;cross_val_score&#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()
nlambdas <span style="color: #666666">=</span> <span style="color: #666666">10</span>
lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-4</span>, <span style="color: #666666">2</span>, nlambdas)
<span style="color: #408080; font-style: italic"># create and fit a ridge regression model, testing each alpha</span>
model <span style="color: #666666">=</span> Ridge()
gridsearch <span style="color: #666666">=</span> GridSearchCV(estimator<span style="color: #666666">=</span>model, param_grid<span style="color: #666666">=</span><span style="color: #008000">dict</span>(alpha<span style="color: #666666">=</span>lambdas))
gridsearch<span style="color: #666666">.</span>fit(X_train, y_train)
<span style="color: #008000">print</span>(gridsearch)
ypredictRidge <span style="color: #666666">=</span> gridsearch<span style="color: #666666">.</span>predict(X_test)
<span style="color: #408080; font-style: italic"># summarize the results of the grid search</span>
<span style="color: #008000">print</span>(<span style="color: #BA2121">f&quot;Best estimated lambda-value: </span><span style="color: #BB6688; font-weight: bold">{</span>gridsearch<span style="color: #666666">.</span>best_estimator_<span style="color: #666666">.</span>alpha<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">&quot;</span>)
<span style="color: #008000">print</span>(<span style="color: #BA2121">f&quot;MSE score: </span><span style="color: #BB6688; font-weight: bold">{</span>MSE(y_test,ypredictRidge)<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">&quot;</span>)
<span style="color: #008000">print</span>(<span style="color: #BA2121">f&quot;R2 score: </span><span style="color: #BB6688; font-weight: bold">{</span>R2(y_test,ypredictRidge)<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">&quot;</span>)
</pre></div>
<p>
By default the grid search function includes cross validation with five folds. The <a href="https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV" target="_blank">Scikit-Learn documentation</a> contains more information on how to set the different parameters.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="randomized-grid-search">Randomized Grid Search </h2>
<p>
An alternative to the above manual grid set up, is to use a random
search where the parameters are tuned from a random distribution
(uniform below) for a fixed number of iterations. A model is
constructed and evaluated for each combination of chosen parameters.
We repeat the previous example but now with a random search.
<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">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">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.linear_model</span> <span style="color: #008000; font-weight: bold">import</span> Ridge
<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> GridSearchCV
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">scipy.stats</span> <span style="color: #008000; font-weight: bold">import</span> uniform <span style="color: #008000; font-weight: bold">as</span> randuniform
<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> RandomizedSearchCV
<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">2021</span>)
n <span style="color: #666666">=</span> <span style="color: #666666">100</span>
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(n)
y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>x<span style="color: #666666">**2</span>) <span style="color: #666666">+</span> <span style="color: #666666">1.5</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>(x<span style="color: #666666">-2</span>)<span style="color: #666666">**2</span>)
Maxpolydegree <span style="color: #666666">=</span> <span style="color: #666666">5</span>
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros((n,Maxpolydegree<span style="color: #666666">-1</span>))
<span style="color: #008000; font-weight: bold">for</span> degree <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(<span style="color: #666666">1</span>,Maxpolydegree): <span style="color: #408080; font-style: italic">#No intercept column</span>
X[:,degree<span style="color: #666666">-1</span>] <span style="color: #666666">=</span> x<span style="color: #666666">**</span>(degree)
<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>)
param_grid <span style="color: #666666">=</span> {<span style="color: #BA2121">&#39;alpha&#39;</span>: randuniform()}
<span style="color: #408080; font-style: italic"># create and fit a ridge regression model, testing each alpha</span>
model <span style="color: #666666">=</span> Ridge()
gridsearch <span style="color: #666666">=</span> RandomizedSearchCV(estimator<span style="color: #666666">=</span>model, param_distributions<span style="color: #666666">=</span>param_grid, n_iter<span style="color: #666666">=100</span>)
gridsearch<span style="color: #666666">.</span>fit(X_train, y_train)
<span style="color: #008000">print</span>(gridsearch)
ypredictRidge <span style="color: #666666">=</span> gridsearch<span style="color: #666666">.</span>predict(X_test)
<span style="color: #408080; font-style: italic"># summarize the results of the grid search</span>
<span style="color: #008000">print</span>(<span style="color: #BA2121">f&quot;Best estimated lambda-value: </span><span style="color: #BB6688; font-weight: bold">{</span>gridsearch<span style="color: #666666">.</span>best_estimator_<span style="color: #666666">.</span>alpha<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">&quot;</span>)
<span style="color: #008000">print</span>(<span style="color: #BA2121">f&quot;MSE score: </span><span style="color: #BB6688; font-weight: bold">{</span>MSE(y_test,ypredictRidge)<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">&quot;</span>)
<span style="color: #008000">print</span>(<span style="color: #BA2121">f&quot;R2 score: </span><span style="color: #BB6688; font-weight: bold">{</span>R2(y_test,ypredictRidge)<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">&quot;</span>)
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
doc/pub/week39/ipynb/ipynb-week39-src.tar.gz
BIN
View File
Binary file not shown.
doc/pub/week39/ipynb/week39.ipynb
+141 -54
View File
@@ -10,7 +10,7 @@
"<!-- Author: --> \n",
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
"\n",
"Date: **Sep 27, 2021**\n",
"Date: **Sep 28, 2021**\n",
"\n",
"Copyright 1999-2021, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
"\n",
@@ -18,7 +18,7 @@
"\n",
"## Plan for week 39\n",
"\n",
"* Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1\n",
"* Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1 and examples on how to implement Logistic Regression\n",
"\n",
"* Friday: Stochastic Gradient descent with examples and automatic differentiation\n",
"\n",
@@ -55,36 +55,45 @@
"%matplotlib inline\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import KFold\n",
"from sklearn.linear_model import Ridge\n",
"from sklearn.model_selection import cross_val_score\n",
"from sklearn.preprocessing import PolynomialFeatures\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn import linear_model\n",
"\n",
"def MSE(y_data,y_model):\n",
" n = np.size(y_model)\n",
" return np.sum((y_data-y_model)**2)/n\n",
"# A seed just to ensure that the random numbers are the same for every run.\n",
"np.random.seed(3155)\n",
"# Generate the data.\n",
"# Useful for eventual debugging.\n",
"np.random.seed(315)\n",
"\n",
"n = 100\n",
"x = np.linspace(-3, 3, n).reshape(-1, 1)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)\n",
"# Decide degree on polynomial to fit\n",
"poly = PolynomialFeatures(degree = 10)\n",
"x = np.random.rand(n)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)\n",
"\n",
"Maxpolydegree = 5\n",
"X = np.zeros((n,Maxpolydegree-1))\n",
"\n",
"for degree in range(1,Maxpolydegree): #No intercept column\n",
" X[:,degree-1] = x**(degree)\n",
"\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",
"# Decide which values of lambda to use\n",
"nlambdas = 500\n",
"lambdas = np.logspace(-3, 5, nlambdas)\n",
"# Initialize a KFold instance\n",
"k = 5\n",
"kfold = KFold(n_splits = k)\n",
"estimated_mse_sklearn = np.zeros(nlambdas)\n",
"i = 0\n",
"for lmb in lambdas:\n",
" ridge = Ridge(alpha = lmb)\n",
" estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)\n",
" estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)\n",
" i += 1\n",
"MSERidgePredict = np.zeros(nlambdas)\n",
"lambdas = np.logspace(-4, 2, nlambdas)\n",
"for i in range(nlambdas):\n",
" lmb = lambdas[i]\n",
" RegRidge = linear_model.Ridge(lmb)\n",
" RegRidge.fit(X_train,y_train)\n",
" ypredictRidge = RegRidge.predict(X_test)\n",
" MSERidgePredict[i] = MSE(y_test,ypredictRidge)\n",
"\n",
"# Now plot the results\n",
"plt.figure()\n",
"plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')\n",
"plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')\n",
"plt.xlabel('log10(lambda)')\n",
"plt.ylabel('MSE')\n",
"plt.legend()\n",
@@ -95,8 +104,10 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"We see from this plot that the optimal MSE occurs for a value of $\\lambda\\in [10,100]$.\n",
"In order to nail down the best value of $\\lambda$, we could in turn narrow down the search area.\n",
"Here we have performed a rather data greedy calculation as function of the regularization parameter $\\lambda$. There is no resampling here. The latter can easily be added by employing the function **RidgeCV** instead of just calling the **Ridge** function. For **RidgeCV** we need to passe the array of $\\lambda$ values.\n",
"By inspecting the figure we can in turn determine which is the optimal regularization parameter.\n",
"This becomes however less functional in the long run. \n",
"\n",
"\n",
"## Grid Search\n",
"\n",
@@ -116,40 +127,116 @@
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import KFold\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.linear_model import Ridge\n",
"from sklearn.model_selection import cross_val_score\n",
"from sklearn.preprocessing import PolynomialFeatures\n",
"from sklearn.model_selection import GridSearchCV\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",
"\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",
"# A seed just to ensure that the random numbers are the same for every run.\n",
"np.random.seed(3155)\n",
"# Generate the data.\n",
"# Useful for eventual debugging.\n",
"np.random.seed(2021)\n",
"\n",
"n = 100\n",
"x = np.linspace(-3, 3, n).reshape(-1, 1)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)\n",
"# Decide degree on polynomial to fit\n",
"poly = PolynomialFeatures(degree = 10)\n",
"x = np.random.rand(n)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)\n",
"\n",
"Maxpolydegree = 5\n",
"X = np.zeros((n,Maxpolydegree-1))\n",
"\n",
"for degree in range(1,Maxpolydegree): #No intercept column\n",
" X[:,degree-1] = x**(degree)\n",
"\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",
"# Decide which values of lambda to use\n",
"nlambdas = 500\n",
"lambdas = np.logspace(-3, 5, nlambdas)\n",
"# Initialize a KFold instance\n",
"k = 5\n",
"kfold = KFold(n_splits = k)\n",
"estimated_mse_sklearn = np.zeros(nlambdas)\n",
"i = 0\n",
"for lmb in lambdas:\n",
" ridge = Ridge(alpha = lmb)\n",
" estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)\n",
" estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)\n",
" i += 1\n",
"plt.figure()\n",
"plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')\n",
"plt.xlabel('log10(lambda)')\n",
"plt.ylabel('MSE')\n",
"plt.legend()\n",
"plt.show()"
"nlambdas = 10\n",
"lambdas = np.logspace(-4, 2, nlambdas)\n",
"# create and fit a ridge regression model, testing each alpha\n",
"model = Ridge()\n",
"gridsearch = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))\n",
"gridsearch.fit(X_train, y_train)\n",
"print(gridsearch)\n",
"ypredictRidge = gridsearch.predict(X_test)\n",
"# summarize the results of the grid search\n",
"print(f\"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}\")\n",
"print(f\"MSE score: {MSE(y_test,ypredictRidge)}\")\n",
"print(f\"R2 score: {R2(y_test,ypredictRidge)}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By default the grid search function includes cross validation with five folds. The [Scikit-Learn documentation](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV) contains more information on how to set the different parameters.\n",
"\n",
"## Randomized Grid Search\n",
"\n",
"An alternative to the above manual grid set up, is to use a random\n",
"search where the parameters are tuned from a random distribution\n",
"(uniform below) for a fixed number of iterations. A model is\n",
"constructed and evaluated for each combination of chosen parameters.\n",
"We repeat the previous example but now with a random search."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.linear_model import Ridge\n",
"from sklearn.model_selection import GridSearchCV\n",
"from scipy.stats import uniform as randuniform\n",
"from sklearn.model_selection import RandomizedSearchCV\n",
"\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",
"\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",
"# A seed just to ensure that the random numbers are the same for every run.\n",
"# Useful for eventual debugging.\n",
"np.random.seed(2021)\n",
"\n",
"n = 100\n",
"x = np.random.rand(n)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)\n",
"\n",
"Maxpolydegree = 5\n",
"X = np.zeros((n,Maxpolydegree-1))\n",
"\n",
"for degree in range(1,Maxpolydegree): #No intercept column\n",
" X[:,degree-1] = x**(degree)\n",
"\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",
"param_grid = {'alpha': randuniform()}\n",
"# create and fit a ridge regression model, testing each alpha\n",
"model = Ridge()\n",
"gridsearch = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=100)\n",
"gridsearch.fit(X_train, y_train)\n",
"print(gridsearch)\n",
"ypredictRidge = gridsearch.predict(X_test)\n",
"# summarize the results of the grid search\n",
"print(f\"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}\")\n",
"print(f\"MSE score: {MSE(y_test,ypredictRidge)}\")\n",
"print(f\"R2 score: {R2(y_test,ypredictRidge)}\")"
]
},
{
doc/src/week39/test1.py → doc/src/week39/codes/test1.py
View File
doc/src/week39/test2.py → doc/src/week39/codes/test2.py
View File
doc/src/week39/test3.py → doc/src/week39/codes/test3.py
View File
doc/src/week39/test4.py → doc/src/week39/codes/test4.py
View File
doc/src/week39/codes/test5.py
+43
View File
@@ -0,0 +1,43 @@
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge
from sklearn.model_selection import GridSearchCV
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(2021)
n = 100
x = np.random.rand(n)
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
Maxpolydegree = 5
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)
# Decide which values of lambda to use
nlambdas = 10
lambdas = np.logspace(-4, 2, nlambdas)
# create and fit a ridge regression model, testing each alpha
model = Ridge()
gridsearch = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))
gridsearch.fit(X_train, y_train)
print(gridsearch)
ypredictRidge = gridsearch.predict(X_test)
# summarize the results of the grid search
#print(gridsearch.best_score_)
print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}")
print(f"MSE score: {MSE(y_test,ypredictRidge)}")
print(f"R2 score: {R2(y_test,ypredictRidge)}")
doc/src/week39/test6.py → doc/src/week39/codes/test6.py
View File
doc/src/week39/codes/test7.py
+44
View File
@@ -0,0 +1,44 @@
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(315)
n = 100
x = np.random.rand(n)
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
Maxpolydegree = 5
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)
# Decide which values of lambda to use
nlambdas = 500
MSERidgePredict = np.zeros(nlambdas)
lambdas = np.logspace(-4, 2, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
RegRidge = linear_model.Ridge(lmb)
RegRidge.fit(X_train,y_train)
ypredictRidge = RegRidge.predict(X_test)
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
# Now plot the results
plt.figure()
plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
plt.xlabel('log10(lambda)')
plt.ylabel('MSE')
plt.legend()
plt.show()
doc/src/week39/test5.py
-61
View File
@@ -1,61 +0,0 @@
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import GridSearchCV
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 = 5
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 = 10
MSEOwnRidgePredict = np.zeros(nlambdas)
MSERidgePredict = np.zeros(nlambdas)
lambdas = np.logspace(-4, 2, nlambdas)
# create and fit a ridge regression model, testing each alpha
model = Ridge()
grid = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))
grid.fit(X_train_scaled, y_train_scaled)
print(grid)
ypredictRidge = grid.predict(X_test_scaled)
# summarize the results of the grid search
print(grid.best_score_)
print(grid.best_estimator_.alpha)
print(MSE(y_test,ypredictRidge))
doc/src/week39/week39.do.txt
+149 -3
View File
@@ -5,7 +5,7 @@ DATE: today
!split
===== Plan for week 39 =====
* Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1
* Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1 and examples on how to implement Logistic Regression
* Friday: Stochastic Gradient descent with examples and automatic differentiation
@@ -31,10 +31,57 @@ in evaluating the MSE as function of different $\lambda$ values.
Based on these calculations, one tries then to determine the value of the hyperparameter $\lambda$
which results in optimal scores (for example the smallest MSE or an $R2=1$).
!bc pycod
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(315)
n = 100
x = np.random.rand(n)
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
Maxpolydegree = 5
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)
# Decide which values of lambda to use
nlambdas = 500
MSERidgePredict = np.zeros(nlambdas)
lambdas = np.logspace(-4, 2, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
RegRidge = linear_model.Ridge(lmb)
RegRidge.fit(X_train,y_train)
ypredictRidge = RegRidge.predict(X_test)
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
# Now plot the results
plt.figure()
plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
plt.xlabel('log10(lambda)')
plt.ylabel('MSE')
plt.legend()
plt.show()
!ec
We see from this plot that the optimal MSE occurs for a value of $\lambda\in [10,100]$.
In order to nail down the best value of $\lambda$, we could in turn narrow down the search area.
Here we have performed a rather data greedy calculation as function of the regularization parameter $\lambda$. There is no resampling here. The latter can easily be added by employing the function _RidgeCV_ instead of just calling the _Ridge_ function. For _RidgeCV_ we need to passe the array of $\lambda$ values.
By inspecting the figure we can in turn determine which is the optimal regularization parameter.
This becomes however less functional in the long run.
!split
===== Grid Search =====
@@ -45,12 +92,111 @@ included with the library _Scikit-Learn_, as demonstrated for the same
example here.
!bc pycod
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge
from sklearn.model_selection import GridSearchCV
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(2021)
n = 100
x = np.random.rand(n)
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
Maxpolydegree = 5
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)
# Decide which values of lambda to use
nlambdas = 10
lambdas = np.logspace(-4, 2, nlambdas)
# create and fit a ridge regression model, testing each alpha
model = Ridge()
gridsearch = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))
gridsearch.fit(X_train, y_train)
print(gridsearch)
ypredictRidge = gridsearch.predict(X_test)
# summarize the results of the grid search
print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}")
print(f"MSE score: {MSE(y_test,ypredictRidge)}")
print(f"R2 score: {R2(y_test,ypredictRidge)}")
!ec
By default the grid search function includes cross validation with five folds. The "Scikit-Learn documentation":"https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV" contains more information on how to set the different parameters.
!split
===== Randomized Grid Search =====
An alternative to the above manual grid set up, is to use a random
search where the parameters are tuned from a random distribution
(uniform below) for a fixed number of iterations. A model is
constructed and evaluated for each combination of chosen parameters.
We repeat the previous example but now with a random search.
!bc pycod
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge
from sklearn.model_selection import GridSearchCV
from scipy.stats import uniform as randuniform
from sklearn.model_selection import RandomizedSearchCV
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(2021)
n = 100
x = np.random.rand(n)
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
Maxpolydegree = 5
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)
param_grid = {'alpha': randuniform()}
# create and fit a ridge regression model, testing each alpha
model = Ridge()
gridsearch = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=100)
gridsearch.fit(X_train, y_train)
print(gridsearch)
ypredictRidge = gridsearch.predict(X_test)
# summarize the results of the grid search
print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}")
print(f"MSE score: {MSE(y_test,ypredictRidge)}")
print(f"R2 score: {R2(y_test,ypredictRidge)}")
!ec
!split