lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
4ba5ffd3b53cd4cfc0e814fd331ae02783f41fc7
FYS-STK4155/doc/pub/week38/html/week38.html
T
mhjensen dd955bf2b0 udpdates
2020-09-17 11:16:08 +02:00

1706 lines
110 KiB
HTML
Raw Blame History

<!--
Automatically generated HTML file from DocOnce source
(https://github.com/hplgit/doconce/)
-->
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="DocOnce: https://github.com/hplgit/doconce/" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" />
<meta name="description" content="Data Analysis and Machine Learning: Logistic Regression">
<title>Data Analysis and Machine Learning: Logistic Regression</title>
<style type="text/css">
/* bloodish style */
body {
font-family: Helvetica, Verdana, Arial, Sans-serif;
color: #404040;
background: #ffffff;
}
h1 { font-size: 1.8em; color: #8A0808; }
h2 { font-size: 1.6em; color: #8A0808; }
h3 { font-size: 1.4em; color: #8A0808; }
h4 { color: #8A0808; }
a { color: #8A0808; text-decoration:none; }
tt { font-family: "Courier New", Courier; }
/* pre style removed because it will interfer with pygments */
p { text-indent: 0px; }
hr { border: 0; width: 80%; border-bottom: 1px solid #aaa}
p.caption { width: 80%; font-style: normal; text-align: left; }
hr.figure { border: 0; width: 80%; border-bottom: 1px solid #aaa}
div { text-align: justify; text-justify: inter-word; }
</style>
</head>
<!-- tocinfo
{'highest level': 2,
'sections': [('Plans for week 38', 2, None, '___sec0'),
('Thursday September 17', 2, None, '___sec1'),
('Ridge and LASSO Regression, reminder', 2, None, '___sec2'),
('Various steps in cross-validation', 2, None, '___sec3'),
('How to set up the cross-validation for Ridge and/or Lasso',
2,
None,
'___sec4'),
('Cross-validation in brief', 2, None, '___sec5'),
('Code Example for Cross-validation and $k$-fold '
'Cross-validation',
2,
None,
'___sec6'),
('Bias-Variance tradeoff with Bootstrap', 2, None, '___sec7'),
("Another Example from Scikit-Learn's Repository",
2,
None,
'___sec8'),
('Cross-validation with Ridge', 2, None, '___sec9'),
('The Ising model', 2, None, '___sec10'),
('Reformulating the problem to suit regression',
2,
None,
'___sec11'),
('Linear regression', 2, None, '___sec12'),
('Singular Value decomposition', 2, None, '___sec13'),
('The one-dimensional Ising model', 2, None, '___sec14'),
('Ridge regression', 2, None, '___sec15'),
('LASSO regression', 2, None, '___sec16'),
('Performance as function of the regularization parameter',
2,
None,
'___sec17'),
('Finding the optimal value of $\\lambda$', 2, None, '___sec18'),
('Friday September 18: Intro to Logistic Regression',
2,
None,
'___sec19'),
('Logistic Regression', 2, None, '___sec20'),
('Classification problems', 2, None, '___sec21'),
('Optimization and Deep learning', 2, None, '___sec22'),
('Basics', 2, None, '___sec23'),
('Linear classifier', 2, None, '___sec24'),
('Some selected properties', 2, None, '___sec25'),
('The logistic function', 2, None, '___sec26'),
('Examples of likelihood functions used in logistic regression '
'and nueral networks',
2,
None,
'___sec27'),
('Two parameters', 2, None, '___sec28'),
('Maximum likelihood', 2, None, '___sec29'),
('The cost function rewritten', 2, None, '___sec30'),
('Minimizing the cross entropy', 2, None, '___sec31'),
('A more compact expression', 2, None, '___sec32'),
('Extending to more predictors', 2, None, '___sec33'),
('Including more classes', 2, None, '___sec34'),
('More classes', 2, None, '___sec35'),
('A simple classification problem', 2, None, '___sec36'),
('Cancer Data again now with Decision Trees and other Methods',
2,
None,
'___sec37'),
('Other measures in classification studies: Cancer Data again',
2,
None,
'___sec38')]}
end of tocinfo -->
<body>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
TeX: {
equationNumbers: { autoNumber: "AMS" },
extensions: ["AMSmath.js", "AMSsymbols.js", "autobold.js", "color.js"]
}
});
</script>
<script type="text/javascript" async
src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<!-- ------------------- main content ---------------------- -->
<center><h1>Data Analysis and Machine Learning: Logistic Regression</h1></center> <!-- document title -->
<p>
<!-- author(s): Morten Hjorth-Jensen -->
<center>
<b>Morten Hjorth-Jensen</b> [1, 2]
</center>
<p>
<!-- institution(s) -->
<center>[1] <b>Department of Physics, University of Oslo</b></center>
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Sep 17, 2020</h4></center> <!-- date -->
<br>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec0">Plans for week 38 </h2>
<ul>
<li> Thursday: Summary of regression methods and discussion of project 1. We revisit also cross-validation and bootstrap as resampling techniques with examples</li>
<li> Friday: Logistic Regression</li>
</ul>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec1">Thursday September 17 </h2>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec2">Ridge and LASSO Regression, reminder </h2>
<p>
The expression for the standard Mean Squared Error (MSE) which we used to define our cost function and the equations for the ordinary least squares (OLS) method, that is
our optimization problem is
$$
{\displaystyle \min_{\boldsymbol{\beta}\in {\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}.
$$
or we can state it as
$$
{\displaystyle \min_{\boldsymbol{\beta}\in
{\mathbb{R}}^{p}}}\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\vert\vert_2^2,
$$
where we have used the definition of a norm-2 vector, that is
$$
\vert\vert \boldsymbol{x}\vert\vert_2 = \sqrt{\sum_i x_i^2}.
$$
<p>
By minimizing the above equation with respect to the parameters
\( \boldsymbol{\beta} \) we could then obtain an analytical expression for the
parameters \( \boldsymbol{\beta} \). We can add a regularization parameter \( \lambda \) by
defining a new cost function to be optimized, that is
$$
{\displaystyle \min_{\boldsymbol{\beta}\in
{\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\vert\vert_2^2+\lambda\vert\vert \boldsymbol{\beta}\vert\vert_2^2
$$
<p>
which leads to the Ridge regression minimization problem where we
require that \( \vert\vert \boldsymbol{\beta}\vert\vert_2^2\le t \), where \( t \) is
a finite number larger than zero. By defining
$$
C(\boldsymbol{X},\boldsymbol{\beta})=\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\vert\vert_2^2+\lambda\vert\vert \boldsymbol{\beta}\vert\vert_1,
$$
<p>
we have a new optimization equation
$$
{\displaystyle \min_{\boldsymbol{\beta}\in
{\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\vert\vert_2^2+\lambda\vert\vert \boldsymbol{\beta}\vert\vert_1
$$
which leads to Lasso regression. Lasso stands for least absolute shrinkage and selection operator.
<p>
Here we have defined the norm-1 as
$$
\vert\vert \boldsymbol{x}\vert\vert_1 = \sum_i \vert x_i\vert.
$$
<p>
<!-- !split -->
<h2 id="___sec3">Various steps in cross-validation </h2>
<p>
When the repetitive splitting of the data set is done randomly,
samples may accidently end up in a fast majority of the splits in
either training or test set. Such samples may have an unbalanced
influence on either model building or prediction evaluation. To avoid
this \( k \)-fold cross-validation structures the data splitting. The
samples are divided into \( k \) more or less equally sized exhaustive and
mutually exclusive subsets. In turn (at each split) one of these
subsets plays the role of the test set while the union of the
remaining subsets constitutes the training set. Such a splitting
warrants a balanced representation of each sample in both training and
test set over the splits. Still the division into the \( k \) subsets
involves a degree of randomness. This may be fully excluded when
choosing \( k=n \). This particular case is referred to as leave-one-out
cross-validation (LOOCV).
<p>
<!-- !split -->
<h2 id="___sec4">How to set up the cross-validation for Ridge and/or Lasso </h2>
<ul>
<li> Define a range of interest for the penalty parameter.</li>
<li> Divide the data set into training and test set comprising samples \( \{1, \ldots, n\} \setminus i \) and \( \{ i \} \), respectively.</li>
<li> Fit the linear regression model by means of ridge estimation for each \( \lambda \) in the grid using the training set, and the corresponding estimate of the error variance \( \boldsymbol{\sigma}_{-i}^2(\lambda) \), as</li>
</ul>
$$
\begin{align*}
\boldsymbol{\beta}_{-i}(\lambda) & = ( \boldsymbol{X}_{-i, \ast}^{T}
\boldsymbol{X}_{-i, \ast} + \lambda \boldsymbol{I}_{pp})^{-1}
\boldsymbol{X}_{-i, \ast}^{T} \boldsymbol{y}_{-i}
\end{align*}
$$
<ul>
<li> Evaluate the prediction performance of these models on the test set by \( \log\{L[y_i, \boldsymbol{X}_{i, \ast}; \boldsymbol{\beta}_{-i}(\lambda), \boldsymbol{\sigma}_{-i}^2(\lambda)]\} \). Or, by the prediction error \( |y_i - \boldsymbol{X}_{i, \ast} \boldsymbol{\beta}_{-i}(\lambda)| \), the relative error, the error squared or the R2 score function.</li>
<li> Repeat the first three steps such that each sample plays the role of the test set once.</li>
<li> Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as</li>
</ul>
$$
\begin{align*}
\frac{1}{n} \sum_{i = 1}^n \log\{L[y_i, \mathbf{X}_{i, \ast}; \boldsymbol{\beta}_{-i}(\lambda), \boldsymbol{\sigma}_{-i}^2(\lambda)]\}.
\end{align*}
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec5">Cross-validation in brief </h2>
<p>
For the various values of \( k \)
<ol>
<li> shuffle the dataset randomly.</li>
<li> Split the dataset into \( k \) groups.</li>
<li> For each unique group:
<ol type="a"></li>
<li> Decide which group to use as set for test data</li>
<li> Take the remaining groups as a training data set</li>
<li> Fit a model on the training set and evaluate it on the test set</li>
<li> Retain the evaluation score and discard the model</li>
</ol>
<li> Summarize the model using the sample of model evaluation scores</li>
</ol>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec6">Code Example for Cross-validation and \( k \)-fold Cross-validation </h2>
<p>
The code here uses Ridge regression with cross-validation (CV) resampling and \( k \)-fold CV in order to fit a specific polynomial.
<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">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: #408080; font-style: italic"># A seed just to ensure that the random numbers are the same for every run.</span>
<span style="color: #408080; font-style: italic"># Useful for eventual debugging.</span>
np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>seed(<span style="color: #666666">3155</span>)
<span style="color: #408080; font-style: italic"># Generate the data.</span>
nsamples <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>randn(nsamples)
y <span style="color: #666666">=</span> <span style="color: #666666">3*</span>x<span style="color: #666666">**2</span> <span style="color: #666666">+</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(nsamples)
<span style="color: #408080; font-style: italic">## Cross-validation on Ridge regression using KFold only</span>
<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">6</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)
<span style="color: #408080; font-style: italic"># Perform the cross-validation to estimate MSE</span>
scores_KFold <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros((nlambdas, k))
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)
j <span style="color: #666666">=</span> <span style="color: #666666">0</span>
<span style="color: #008000; font-weight: bold">for</span> train_inds, test_inds <span style="color: #AA22FF; font-weight: bold">in</span> kfold<span style="color: #666666">.</span>split(x):
xtrain <span style="color: #666666">=</span> x[train_inds]
ytrain <span style="color: #666666">=</span> y[train_inds]
xtest <span style="color: #666666">=</span> x[test_inds]
ytest <span style="color: #666666">=</span> y[test_inds]
Xtrain <span style="color: #666666">=</span> poly<span style="color: #666666">.</span>fit_transform(xtrain[:, np<span style="color: #666666">.</span>newaxis])
ridge<span style="color: #666666">.</span>fit(Xtrain, ytrain[:, np<span style="color: #666666">.</span>newaxis])
Xtest <span style="color: #666666">=</span> poly<span style="color: #666666">.</span>fit_transform(xtest[:, np<span style="color: #666666">.</span>newaxis])
ypred <span style="color: #666666">=</span> ridge<span style="color: #666666">.</span>predict(Xtest)
scores_KFold[i,j] <span style="color: #666666">=</span> np<span style="color: #666666">.</span>sum((ypred <span style="color: #666666">-</span> ytest[:, np<span style="color: #666666">.</span>newaxis])<span style="color: #666666">**2</span>)<span style="color: #666666">/</span>np<span style="color: #666666">.</span>size(ypred)
j <span style="color: #666666">+=</span> <span style="color: #666666">1</span>
i <span style="color: #666666">+=</span> <span style="color: #666666">1</span>
estimated_mse_KFold <span style="color: #666666">=</span> np<span style="color: #666666">.</span>mean(scores_KFold, axis <span style="color: #666666">=</span> <span style="color: #666666">1</span>)
<span style="color: #408080; font-style: italic">## Cross-validation using cross_val_score from sklearn along with KFold</span>
<span style="color: #408080; font-style: italic"># kfold is an instance initialized above as:</span>
<span style="color: #408080; font-style: italic"># kfold = KFold(n_splits = k)</span>
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)
X <span style="color: #666666">=</span> poly<span style="color: #666666">.</span>fit_transform(x[:, np<span style="color: #666666">.</span>newaxis])
estimated_mse_folds <span style="color: #666666">=</span> cross_val_score(ridge, X, y[:, np<span style="color: #666666">.</span>newaxis], scoring<span style="color: #666666">=</span><span style="color: #BA2121">&#39;neg_mean_squared_error&#39;</span>, cv<span style="color: #666666">=</span>kfold)
<span style="color: #408080; font-style: italic"># cross_val_score return an array containing the estimated negative mse for every fold.</span>
<span style="color: #408080; font-style: italic"># we have to the the mean of every array in order to get an estimate of the mse of the model</span>
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>
<span style="color: #408080; font-style: italic">## Plot and compare the slightly different ways to perform cross-validation</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), estimated_mse_KFold, <span style="color: #BA2121">&#39;r--&#39;</span>, label <span style="color: #666666">=</span> <span style="color: #BA2121">&#39;KFold&#39;</span>)
plt<span style="color: #666666">.</span>xlabel(<span style="color: #BA2121">&#39;log10(lambda)&#39;</span>)
plt<span style="color: #666666">.</span>ylabel(<span style="color: #BA2121">&#39;mse&#39;</span>)
plt<span style="color: #666666">.</span>legend()
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec7">Bias-Variance tradeoff with Bootstrap </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">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">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.linear_model</span> <span style="color: #008000; font-weight: bold">import</span> LinearRegression, Ridge, Lasso
<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.pipeline</span> <span style="color: #008000; font-weight: bold">import</span> make_pipeline
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">sklearn.utils</span> <span style="color: #008000; font-weight: bold">import</span> resample
np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>seed(<span style="color: #666666">2018</span>)
n <span style="color: #666666">=</span> <span style="color: #666666">40</span>
n_boostraps <span style="color: #666666">=</span> <span style="color: #666666">100</span>
maxdegree <span style="color: #666666">=</span> <span style="color: #666666">14</span>
<span style="color: #408080; font-style: italic"># Make data set.</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)
error <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(maxdegree)
bias <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(maxdegree)
variance <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(maxdegree)
polydegree <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(maxdegree)
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: #008000; font-weight: bold">for</span> degree <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(maxdegree):
model <span style="color: #666666">=</span> make_pipeline(PolynomialFeatures(degree<span style="color: #666666">=</span>degree), LinearRegression(fit_intercept<span style="color: #666666">=</span><span style="color: #008000">False</span>))
y_pred <span style="color: #666666">=</span> np<span style="color: #666666">.</span>empty((y_test<span style="color: #666666">.</span>shape[<span style="color: #666666">0</span>], n_boostraps))
<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>(n_boostraps):
x_, y_ <span style="color: #666666">=</span> resample(x_train, y_train)
y_pred[:, i] <span style="color: #666666">=</span> model<span style="color: #666666">.</span>fit(x_, y_)<span style="color: #666666">.</span>predict(x_test)<span style="color: #666666">.</span>ravel()
polydegree[degree] <span style="color: #666666">=</span> degree
error[degree] <span style="color: #666666">=</span> np<span style="color: #666666">.</span>mean( np<span style="color: #666666">.</span>mean((y_test <span style="color: #666666">-</span> y_pred)<span style="color: #666666">**2</span>, axis<span style="color: #666666">=1</span>, keepdims<span style="color: #666666">=</span><span style="color: #008000">True</span>) )
bias[degree] <span style="color: #666666">=</span> np<span style="color: #666666">.</span>mean( (y_test <span style="color: #666666">-</span> np<span style="color: #666666">.</span>mean(y_pred, axis<span style="color: #666666">=1</span>, keepdims<span style="color: #666666">=</span><span style="color: #008000">True</span>))<span style="color: #666666">**2</span> )
variance[degree] <span style="color: #666666">=</span> np<span style="color: #666666">.</span>mean( np<span style="color: #666666">.</span>var(y_pred, axis<span style="color: #666666">=1</span>, keepdims<span style="color: #666666">=</span><span style="color: #008000">True</span>) )
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&#39;Polynomial degree:&#39;</span>, degree)
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&#39;Error:&#39;</span>, error[degree])
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&#39;Bias^2:&#39;</span>, bias[degree])
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&#39;Var:&#39;</span>, variance[degree])
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&#39;{} &gt;= {} + {} = {}&#39;</span><span style="color: #666666">.</span>format(error[degree], bias[degree], variance[degree], bias[degree]<span style="color: #666666">+</span>variance[degree]))
plt<span style="color: #666666">.</span>plot(polydegree, error, label<span style="color: #666666">=</span><span style="color: #BA2121">&#39;Error&#39;</span>)
plt<span style="color: #666666">.</span>plot(polydegree, bias, label<span style="color: #666666">=</span><span style="color: #BA2121">&#39;bias&#39;</span>)
plt<span style="color: #666666">.</span>plot(polydegree, variance, label<span style="color: #666666">=</span><span style="color: #BA2121">&#39;Variance&#39;</span>)
plt<span style="color: #666666">.</span>legend()
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec8">Another Example from Scikit-Learn's Repository </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;</span>
<span style="color: #BA2121; font-style: italic">============================</span>
<span style="color: #BA2121; font-style: italic">Underfitting vs. Overfitting</span>
<span style="color: #BA2121; font-style: italic">============================</span>
<span style="color: #BA2121; font-style: italic">This example demonstrates the problems of underfitting and overfitting and</span>
<span style="color: #BA2121; font-style: italic">how we can use linear regression with polynomial features to approximate</span>
<span style="color: #BA2121; font-style: italic">nonlinear functions. The plot shows the function that we want to approximate,</span>
<span style="color: #BA2121; font-style: italic">which is a part of the cosine function. In addition, the samples from the</span>
<span style="color: #BA2121; font-style: italic">real function and the approximations of different models are displayed. The</span>
<span style="color: #BA2121; font-style: italic">models have polynomial features of different degrees. We can see that a</span>
<span style="color: #BA2121; font-style: italic">linear function (polynomial with degree 1) is not sufficient to fit the</span>
<span style="color: #BA2121; font-style: italic">training samples. This is called **underfitting**. A polynomial of degree 4</span>
<span style="color: #BA2121; font-style: italic">approximates the true function almost perfectly. However, for higher degrees</span>
<span style="color: #BA2121; font-style: italic">the model will **overfit** the training data, i.e. it learns the noise of the</span>
<span style="color: #BA2121; font-style: italic">training data.</span>
<span style="color: #BA2121; font-style: italic">We evaluate quantitatively **overfitting** / **underfitting** by using</span>
<span style="color: #BA2121; font-style: italic">cross-validation. We calculate the mean squared error (MSE) on the validation</span>
<span style="color: #BA2121; font-style: italic">set, the higher, the less likely the model generalizes correctly from the</span>
<span style="color: #BA2121; font-style: italic">training data.</span>
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;</span>
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #19177C">__doc__</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.pipeline</span> <span style="color: #008000; font-weight: bold">import</span> Pipeline
<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.linear_model</span> <span style="color: #008000; font-weight: bold">import</span> LinearRegression
<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">def</span> <span style="color: #0000FF">true_fun</span>(X):
<span style="color: #008000; font-weight: bold">return</span> np<span style="color: #666666">.</span>cos(<span style="color: #666666">1.5</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>pi <span style="color: #666666">*</span> X)
np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>seed(<span style="color: #666666">0</span>)
n_samples <span style="color: #666666">=</span> <span style="color: #666666">30</span>
degrees <span style="color: #666666">=</span> [<span style="color: #666666">1</span>, <span style="color: #666666">4</span>, <span style="color: #666666">15</span>]
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>sort(np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(n_samples))
y <span style="color: #666666">=</span> true_fun(X) <span style="color: #666666">+</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(n_samples) <span style="color: #666666">*</span> <span style="color: #666666">0.1</span>
plt<span style="color: #666666">.</span>figure(figsize<span style="color: #666666">=</span>(<span style="color: #666666">14</span>, <span style="color: #666666">5</span>))
<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>(<span style="color: #008000">len</span>(degrees)):
ax <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>subplot(<span style="color: #666666">1</span>, <span style="color: #008000">len</span>(degrees), i <span style="color: #666666">+</span> <span style="color: #666666">1</span>)
plt<span style="color: #666666">.</span>setp(ax, xticks<span style="color: #666666">=</span>(), yticks<span style="color: #666666">=</span>())
polynomial_features <span style="color: #666666">=</span> PolynomialFeatures(degree<span style="color: #666666">=</span>degrees[i],
include_bias<span style="color: #666666">=</span><span style="color: #008000">False</span>)
linear_regression <span style="color: #666666">=</span> LinearRegression()
pipeline <span style="color: #666666">=</span> Pipeline([(<span style="color: #BA2121">&quot;polynomial_features&quot;</span>, polynomial_features),
(<span style="color: #BA2121">&quot;linear_regression&quot;</span>, linear_regression)])
pipeline<span style="color: #666666">.</span>fit(X[:, np<span style="color: #666666">.</span>newaxis], y)
<span style="color: #408080; font-style: italic"># Evaluate the models using crossvalidation</span>
scores <span style="color: #666666">=</span> cross_val_score(pipeline, X[:, np<span style="color: #666666">.</span>newaxis], y,
scoring<span style="color: #666666">=</span><span style="color: #BA2121">&quot;neg_mean_squared_error&quot;</span>, cv<span style="color: #666666">=10</span>)
X_test <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linspace(<span style="color: #666666">0</span>, <span style="color: #666666">1</span>, <span style="color: #666666">100</span>)
plt<span style="color: #666666">.</span>plot(X_test, pipeline<span style="color: #666666">.</span>predict(X_test[:, np<span style="color: #666666">.</span>newaxis]), label<span style="color: #666666">=</span><span style="color: #BA2121">&quot;Model&quot;</span>)
plt<span style="color: #666666">.</span>plot(X_test, true_fun(X_test), label<span style="color: #666666">=</span><span style="color: #BA2121">&quot;True function&quot;</span>)
plt<span style="color: #666666">.</span>scatter(X, y, edgecolor<span style="color: #666666">=</span><span style="color: #BA2121">&#39;b&#39;</span>, s<span style="color: #666666">=20</span>, label<span style="color: #666666">=</span><span style="color: #BA2121">&quot;Samples&quot;</span>)
plt<span style="color: #666666">.</span>xlabel(<span style="color: #BA2121">&quot;x&quot;</span>)
plt<span style="color: #666666">.</span>ylabel(<span style="color: #BA2121">&quot;y&quot;</span>)
plt<span style="color: #666666">.</span>xlim((<span style="color: #666666">0</span>, <span style="color: #666666">1</span>))
plt<span style="color: #666666">.</span>ylim((<span style="color: #666666">-2</span>, <span style="color: #666666">2</span>))
plt<span style="color: #666666">.</span>legend(loc<span style="color: #666666">=</span><span style="color: #BA2121">&quot;best&quot;</span>)
plt<span style="color: #666666">.</span>title(<span style="color: #BA2121">&quot;Degree {}</span><span style="color: #BB6622; font-weight: bold">\n</span><span style="color: #BA2121">MSE = {:.2e}(+/- {:.2e})&quot;</span><span style="color: #666666">.</span>format(
degrees[i], <span style="color: #666666">-</span>scores<span style="color: #666666">.</span>mean(), scores<span style="color: #666666">.</span>std()))
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec9">Cross-validation with Ridge </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">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.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: #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>
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>)
<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()
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec10">The Ising model </h2>
<p>
The one-dimensional Ising model with nearest neighbor interaction, no
external field and a constant coupling constant \( J \) is given by
$$
\begin{align}
H = -J \sum_{k}^L s_k s_{k + 1},
\label{_auto1}
\end{align}
$$
<p>
where \( s_i \in \{-1, 1\} \) and \( s_{N + 1} = s_1 \). The number of spins
in the system is determined by \( L \). For the one-dimensional system
there is no phase transition.
<p>
We will look at a system of \( L = 40 \) spins with a coupling constant of
\( J = 1 \). To get enough training data we will generate 10000 states
with their respective energies.
<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">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">mpl_toolkits.axes_grid1</span> <span style="color: #008000; font-weight: bold">import</span> make_axes_locatable
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">seaborn</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">sns</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">scipy.linalg</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">scl</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">import</span> <span style="color: #0000FF; font-weight: bold">tqdm</span>
sns<span style="color: #666666">.</span>set(color_codes<span style="color: #666666">=</span><span style="color: #008000">True</span>)
cmap_args<span style="color: #666666">=</span><span style="color: #008000">dict</span>(vmin<span style="color: #666666">=-1.</span>, vmax<span style="color: #666666">=1.</span>, cmap<span style="color: #666666">=</span><span style="color: #BA2121">&#39;seismic&#39;</span>)
L <span style="color: #666666">=</span> <span style="color: #666666">40</span>
n <span style="color: #666666">=</span> <span style="color: #008000">int</span>(<span style="color: #666666">1e4</span>)
spins <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>choice([<span style="color: #666666">-1</span>, <span style="color: #666666">1</span>], size<span style="color: #666666">=</span>(n, L))
J <span style="color: #666666">=</span> <span style="color: #666666">1.0</span>
energies <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(n)
<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>(n):
energies[i] <span style="color: #666666">=</span> <span style="color: #666666">-</span> J <span style="color: #666666">*</span> np<span style="color: #666666">.</span>dot(spins[i], np<span style="color: #666666">.</span>roll(spins[i], <span style="color: #666666">1</span>))
</pre></div>
<p>
Here we use ordinary least squares
regression to predict the energy for the nearest neighbor
one-dimensional Ising model on a ring, i.e., the endpoints wrap
around. We will use linear regression to fit a value for
the coupling constant to achieve this.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec11">Reformulating the problem to suit regression </h2>
<p>
A more general form for the one-dimensional Ising model is
$$
\begin{align}
H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
\label{_auto2}
\end{align}
$$
<p>
Here we allow for interactions beyond the nearest neighbors and a state dependent
coupling constant. This latter expression can be formulated as
a matrix-product
$$
\begin{align}
\boldsymbol{H} = \boldsymbol{X} J,
\label{_auto3}
\end{align}
$$
<p>
where \( X_{jk} = s_j s_k \) and \( J \) is a matrix which consists of the
elements \( -J_{jk} \). This form of writing the energy fits perfectly
with the form utilized in linear regression, that is
$$
\begin{align}
\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon},
\label{_auto4}
\end{align}
$$
<p>
We split the data in training and test data as discussed in the previous example
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros((n, L <span style="color: #666666">**</span> <span style="color: #666666">2</span>))
<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>(n):
X[i] <span style="color: #666666">=</span> np<span style="color: #666666">.</span>outer(spins[i], spins[i])<span style="color: #666666">.</span>ravel()
y <span style="color: #666666">=</span> energies
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>)
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec12">Linear regression </h2>
<p>
In the ordinary least squares method we choose the cost function
$$
\begin{align}
C(\boldsymbol{X}, \boldsymbol{\beta})= \frac{1}{n}\left\{(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})\right\}.
\label{_auto5}
\end{align}
$$
<p>
We then find the extremal point of \( C \) by taking the derivative with respect to \( \boldsymbol{\beta} \) as discussed above.
This yields the expression for \( \boldsymbol{\beta} \) to be
$$
\boldsymbol{\beta} = \frac{\boldsymbol{X}^T \boldsymbol{y}}{\boldsymbol{X}^T \boldsymbol{X}},
$$
<p>
which immediately imposes some requirements on \( \boldsymbol{X} \) as there must exist
an inverse of \( \boldsymbol{X}^T \boldsymbol{X} \). If the expression we are modeling contains an
intercept, i.e., a constant term, we must make sure that the
first column of \( \boldsymbol{X} \) consists of \( 1 \). We do this here
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>X_train_own <span style="color: #666666">=</span> np<span style="color: #666666">.</span>concatenate(
(np<span style="color: #666666">.</span>ones(<span style="color: #008000">len</span>(X_train))[:, np<span style="color: #666666">.</span>newaxis], X_train),
axis<span style="color: #666666">=1</span>
)
X_test_own <span style="color: #666666">=</span> np<span style="color: #666666">.</span>concatenate(
(np<span style="color: #666666">.</span>ones(<span style="color: #008000">len</span>(X_test))[:, np<span style="color: #666666">.</span>newaxis], X_test),
axis<span style="color: #666666">=1</span>
)
</pre></div>
<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">def</span> <span style="color: #0000FF">ols_inv</span>(x: np<span style="color: #666666">.</span>ndarray, y: np<span style="color: #666666">.</span>ndarray) <span style="color: #666666">-&gt;</span> np<span style="color: #666666">.</span>ndarray:
<span style="color: #008000; font-weight: bold">return</span> scl<span style="color: #666666">.</span>inv(x<span style="color: #666666">.</span>T @ x) @ (x<span style="color: #666666">.</span>T @ y)
beta <span style="color: #666666">=</span> ols_inv(X_train_own, y_train)
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec13">Singular Value decomposition </h2>
<p>
Doing the inversion directly turns out to be a bad idea since the matrix
\( \boldsymbol{X}^T\boldsymbol{X} \) is singular. An alternative approach is to use the <b>singular
value decomposition</b>. Using the definition of the Moore-Penrose
pseudoinverse we can write the equation for \( \boldsymbol{\beta} \) as
$$
\boldsymbol{\beta} = \boldsymbol{X}^{+}\boldsymbol{y},
$$
<p>
where the pseudoinverse of \( \boldsymbol{X} \) is given by
$$
\boldsymbol{X}^{+} = \frac{\boldsymbol{X}^T}{\boldsymbol{X}^T\boldsymbol{X}}.
$$
<p>
Using singular value decomposition we can decompose the matrix \( \boldsymbol{X} = \boldsymbol{U}\boldsymbol{\Sigma} \boldsymbol{V}^T \),
where \( \boldsymbol{U} \) and \( \boldsymbol{V} \) are orthogonal(unitary) matrices and \( \boldsymbol{\Sigma} \) contains the singular values (more details below).
where \( X^{+} = V\Sigma^{+} U^T \). This reduces the equation for
\( \omega \) to
$$
\begin{align}
\boldsymbol{\beta} = \boldsymbol{V}\boldsymbol{\Sigma}^{+} \boldsymbol{U}^T \boldsymbol{y}.
\label{_auto6}
\end{align}
$$
<p>
Note that solving this equation by actually doing the pseudoinverse
(which is what we will do) is not a good idea as this operation scales
as \( \mathcal{O}(n^3) \), where \( n \) is the number of elements in a
general matrix. Instead, doing \( QR \)-factorization and solving the
linear system as an equation would reduce this down to
\( \mathcal{O}(n^2) \) operations.
<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">def</span> <span style="color: #0000FF">ols_svd</span>(x: np<span style="color: #666666">.</span>ndarray, y: np<span style="color: #666666">.</span>ndarray) <span style="color: #666666">-&gt;</span> np<span style="color: #666666">.</span>ndarray:
u, s, v <span style="color: #666666">=</span> scl<span style="color: #666666">.</span>svd(x)
<span style="color: #008000; font-weight: bold">return</span> v<span style="color: #666666">.</span>T @ scl<span style="color: #666666">.</span>pinv(scl<span style="color: #666666">.</span>diagsvd(s, u<span style="color: #666666">.</span>shape[<span style="color: #666666">0</span>], v<span style="color: #666666">.</span>shape[<span style="color: #666666">0</span>])) @ u<span style="color: #666666">.</span>T @ y
</pre></div>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>beta <span style="color: #666666">=</span> ols_svd(X_train_own,y_train)
</pre></div>
<p>
When extracting the \( J \)-matrix we need to make sure that we remove the intercept, as is done here
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>J <span style="color: #666666">=</span> beta[<span style="color: #666666">1</span>:]<span style="color: #666666">.</span>reshape(L, L)
</pre></div>
<p>
A way of looking at the coefficients in \( J \) is to plot the matrices as images.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>fig <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>figure(figsize<span style="color: #666666">=</span>(<span style="color: #666666">20</span>, <span style="color: #666666">14</span>))
im <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>imshow(J, <span style="color: #666666">**</span>cmap_args)
plt<span style="color: #666666">.</span>title(<span style="color: #BA2121">&quot;OLS&quot;</span>, fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>xticks(fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>yticks(fontsize<span style="color: #666666">=18</span>)
cb <span style="color: #666666">=</span> fig<span style="color: #666666">.</span>colorbar(im)
cb<span style="color: #666666">.</span>ax<span style="color: #666666">.</span>set_yticklabels(cb<span style="color: #666666">.</span>ax<span style="color: #666666">.</span>get_yticklabels(), fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
It is interesting to note that OLS
considers both \( J_{j, j + 1} = -0.5 \) and \( J_{j, j - 1} = -0.5 \) as
valid matrix elements for \( J \).
In our discussion below on hyperparameters and Ridge and Lasso regression we will see that
this problem can be removed, partly and only with Lasso regression.
<p>
In this case our matrix inversion was actually possible. The obvious question now is what is the mathematics behind the SVD?
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec14">The one-dimensional Ising model </h2>
<p>
Let us bring back the Ising model again, but now with an additional
focus on Ridge and Lasso regression as well. We repeat some of the
basic parts of the Ising model and the setup of the training and test
data. The one-dimensional Ising model with nearest neighbor
interaction, no external field and a constant coupling constant \( J \) is
given by
$$
\begin{align}
H = -J \sum_{k}^L s_k s_{k + 1},
\label{_auto7}
\end{align}
$$
where \( s_i \in \{-1, 1\} \) and \( s_{N + 1} = s_1 \). The number of spins in the system is determined by \( L \). For the one-dimensional system there is no phase transition.
<p>
We will look at a system of \( L = 40 \) spins with a coupling constant of \( J = 1 \). To get enough training data we will generate 10000 states with their respective energies.
<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">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">mpl_toolkits.axes_grid1</span> <span style="color: #008000; font-weight: bold">import</span> make_axes_locatable
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">seaborn</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">sns</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">scipy.linalg</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">scl</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">import</span> <span style="color: #0000FF; font-weight: bold">sklearn.linear_model</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">skl</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">tqdm</span>
sns<span style="color: #666666">.</span>set(color_codes<span style="color: #666666">=</span><span style="color: #008000">True</span>)
cmap_args<span style="color: #666666">=</span><span style="color: #008000">dict</span>(vmin<span style="color: #666666">=-1.</span>, vmax<span style="color: #666666">=1.</span>, cmap<span style="color: #666666">=</span><span style="color: #BA2121">&#39;seismic&#39;</span>)
L <span style="color: #666666">=</span> <span style="color: #666666">40</span>
n <span style="color: #666666">=</span> <span style="color: #008000">int</span>(<span style="color: #666666">1e4</span>)
spins <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>choice([<span style="color: #666666">-1</span>, <span style="color: #666666">1</span>], size<span style="color: #666666">=</span>(n, L))
J <span style="color: #666666">=</span> <span style="color: #666666">1.0</span>
energies <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros(n)
<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>(n):
energies[i] <span style="color: #666666">=</span> <span style="color: #666666">-</span> J <span style="color: #666666">*</span> np<span style="color: #666666">.</span>dot(spins[i], np<span style="color: #666666">.</span>roll(spins[i], <span style="color: #666666">1</span>))
</pre></div>
<p>
A more general form for the one-dimensional Ising model is
$$
\begin{align}
H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
\label{_auto8}
\end{align}
$$
<p>
Here we allow for interactions beyond the nearest neighbors and a more
adaptive coupling matrix. This latter expression can be formulated as
a matrix-product on the form
$$
\begin{align}
H = X J,
\label{_auto9}
\end{align}
$$
<p>
where \( X_{jk} = s_j s_k \) and \( J \) is the matrix consisting of the
elements \( -J_{jk} \). This form of writing the energy fits perfectly
with the form utilized in linear regression, viz.
$$
\begin{align}
\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}.
\label{_auto10}
\end{align}
$$
We organize the data as we did above
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros((n, L <span style="color: #666666">**</span> <span style="color: #666666">2</span>))
<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>(n):
X[i] <span style="color: #666666">=</span> np<span style="color: #666666">.</span>outer(spins[i], spins[i])<span style="color: #666666">.</span>ravel()
y <span style="color: #666666">=</span> energies
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.96</span>)
X_train_own <span style="color: #666666">=</span> np<span style="color: #666666">.</span>concatenate(
(np<span style="color: #666666">.</span>ones(<span style="color: #008000">len</span>(X_train))[:, np<span style="color: #666666">.</span>newaxis], X_train),
axis<span style="color: #666666">=1</span>
)
X_test_own <span style="color: #666666">=</span> np<span style="color: #666666">.</span>concatenate(
(np<span style="color: #666666">.</span>ones(<span style="color: #008000">len</span>(X_test))[:, np<span style="color: #666666">.</span>newaxis], X_test),
axis<span style="color: #666666">=1</span>
)
</pre></div>
<p>
We will do all fitting with <b>Scikit-Learn</b>,
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>clf <span style="color: #666666">=</span> skl<span style="color: #666666">.</span>LinearRegression()<span style="color: #666666">.</span>fit(X_train, y_train)
</pre></div>
<p>
When extracting the \( J \)-matrix we make sure to remove the intercept
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>J_sk <span style="color: #666666">=</span> clf<span style="color: #666666">.</span>coef_<span style="color: #666666">.</span>reshape(L, L)
</pre></div>
<p>
And then we plot the results
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>fig <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>figure(figsize<span style="color: #666666">=</span>(<span style="color: #666666">20</span>, <span style="color: #666666">14</span>))
im <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>imshow(J_sk, <span style="color: #666666">**</span>cmap_args)
plt<span style="color: #666666">.</span>title(<span style="color: #BA2121">&quot;LinearRegression from Scikit-learn&quot;</span>, fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>xticks(fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>yticks(fontsize<span style="color: #666666">=18</span>)
cb <span style="color: #666666">=</span> fig<span style="color: #666666">.</span>colorbar(im)
cb<span style="color: #666666">.</span>ax<span style="color: #666666">.</span>set_yticklabels(cb<span style="color: #666666">.</span>ax<span style="color: #666666">.</span>get_yticklabels(), fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
The results perfectly with our previous discussion where we used our own code.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec15">Ridge regression </h2>
<p>
Having explored the ordinary least squares we move on to ridge
regression. In ridge regression we include a <b>regularizer</b>. This
involves a new cost function which leads to a new estimate for the
weights \( \boldsymbol{\beta} \). This results in a penalized regression problem. The
cost function is given by
$$
\begin{align}
C(\boldsymbol{X}, \boldsymbol{\beta}; \lambda) = (\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y}) + \lambda \boldsymbol{\beta}^T\boldsymbol{\beta}.
\label{_auto11}
\end{align}
$$
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>_lambda <span style="color: #666666">=</span> <span style="color: #666666">0.1</span>
clf_ridge <span style="color: #666666">=</span> skl<span style="color: #666666">.</span>Ridge(alpha<span style="color: #666666">=</span>_lambda)<span style="color: #666666">.</span>fit(X_train, y_train)
J_ridge_sk <span style="color: #666666">=</span> clf_ridge<span style="color: #666666">.</span>coef_<span style="color: #666666">.</span>reshape(L, L)
fig <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>figure(figsize<span style="color: #666666">=</span>(<span style="color: #666666">20</span>, <span style="color: #666666">14</span>))
im <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>imshow(J_ridge_sk, <span style="color: #666666">**</span>cmap_args)
plt<span style="color: #666666">.</span>title(<span style="color: #BA2121">&quot;Ridge from Scikit-learn&quot;</span>, fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>xticks(fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>yticks(fontsize<span style="color: #666666">=18</span>)
cb <span style="color: #666666">=</span> fig<span style="color: #666666">.</span>colorbar(im)
cb<span style="color: #666666">.</span>ax<span style="color: #666666">.</span>set_yticklabels(cb<span style="color: #666666">.</span>ax<span style="color: #666666">.</span>get_yticklabels(), fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec16">LASSO regression </h2>
<p>
In the <b>Least Absolute Shrinkage and Selection Operator</b> (LASSO)-method we get a third cost function.
$$
\begin{align}
C(\boldsymbol{X}, \boldsymbol{\beta}; \lambda) = (\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y}) + \lambda \sqrt{\boldsymbol{\beta}^T\boldsymbol{\beta}}.
\label{_auto12}
\end{align}
$$
<p>
Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from <b>Scikit-Learn</b>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>clf_lasso <span style="color: #666666">=</span> skl<span style="color: #666666">.</span>Lasso(alpha<span style="color: #666666">=</span>_lambda)<span style="color: #666666">.</span>fit(X_train, y_train)
J_lasso_sk <span style="color: #666666">=</span> clf_lasso<span style="color: #666666">.</span>coef_<span style="color: #666666">.</span>reshape(L, L)
fig <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>figure(figsize<span style="color: #666666">=</span>(<span style="color: #666666">20</span>, <span style="color: #666666">14</span>))
im <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>imshow(J_lasso_sk, <span style="color: #666666">**</span>cmap_args)
plt<span style="color: #666666">.</span>title(<span style="color: #BA2121">&quot;Lasso from Scikit-learn&quot;</span>, fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>xticks(fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>yticks(fontsize<span style="color: #666666">=18</span>)
cb <span style="color: #666666">=</span> fig<span style="color: #666666">.</span>colorbar(im)
cb<span style="color: #666666">.</span>ax<span style="color: #666666">.</span>set_yticklabels(cb<span style="color: #666666">.</span>ax<span style="color: #666666">.</span>get_yticklabels(), fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
It is quite striking how LASSO breaks the symmetry of the coupling
constant as opposed to ridge and OLS. We get a sparse solution with
\( J_{j, j + 1} = -1 \).
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec17">Performance as function of the regularization parameter </h2>
<p>
We see how the different models perform for a different set of values for \( \lambda \).
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-4</span>, <span style="color: #666666">5</span>, <span style="color: #666666">10</span>)
train_errors <span style="color: #666666">=</span> {
<span style="color: #BA2121">&quot;ols_sk&quot;</span>: np<span style="color: #666666">.</span>zeros(lambdas<span style="color: #666666">.</span>size),
<span style="color: #BA2121">&quot;ridge_sk&quot;</span>: np<span style="color: #666666">.</span>zeros(lambdas<span style="color: #666666">.</span>size),
<span style="color: #BA2121">&quot;lasso_sk&quot;</span>: np<span style="color: #666666">.</span>zeros(lambdas<span style="color: #666666">.</span>size)
}
test_errors <span style="color: #666666">=</span> {
<span style="color: #BA2121">&quot;ols_sk&quot;</span>: np<span style="color: #666666">.</span>zeros(lambdas<span style="color: #666666">.</span>size),
<span style="color: #BA2121">&quot;ridge_sk&quot;</span>: np<span style="color: #666666">.</span>zeros(lambdas<span style="color: #666666">.</span>size),
<span style="color: #BA2121">&quot;lasso_sk&quot;</span>: np<span style="color: #666666">.</span>zeros(lambdas<span style="color: #666666">.</span>size)
}
plot_counter <span style="color: #666666">=</span> <span style="color: #666666">1</span>
fig <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>figure(figsize<span style="color: #666666">=</span>(<span style="color: #666666">32</span>, <span style="color: #666666">54</span>))
<span style="color: #008000; font-weight: bold">for</span> i, _lambda <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">enumerate</span>(tqdm<span style="color: #666666">.</span>tqdm(lambdas)):
<span style="color: #008000; font-weight: bold">for</span> key, method <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">zip</span>(
[<span style="color: #BA2121">&quot;ols_sk&quot;</span>, <span style="color: #BA2121">&quot;ridge_sk&quot;</span>, <span style="color: #BA2121">&quot;lasso_sk&quot;</span>],
[skl<span style="color: #666666">.</span>LinearRegression(), skl<span style="color: #666666">.</span>Ridge(alpha<span style="color: #666666">=</span>_lambda), skl<span style="color: #666666">.</span>Lasso(alpha<span style="color: #666666">=</span>_lambda)]
):
method <span style="color: #666666">=</span> method<span style="color: #666666">.</span>fit(X_train, y_train)
train_errors[key][i] <span style="color: #666666">=</span> method<span style="color: #666666">.</span>score(X_train, y_train)
test_errors[key][i] <span style="color: #666666">=</span> method<span style="color: #666666">.</span>score(X_test, y_test)
omega <span style="color: #666666">=</span> method<span style="color: #666666">.</span>coef_<span style="color: #666666">.</span>reshape(L, L)
plt<span style="color: #666666">.</span>subplot(<span style="color: #666666">10</span>, <span style="color: #666666">5</span>, plot_counter)
plt<span style="color: #666666">.</span>imshow(omega, <span style="color: #666666">**</span>cmap_args)
plt<span style="color: #666666">.</span>title(<span style="color: #BA2121">r&quot;</span><span style="color: #BB6688; font-weight: bold">%s</span><span style="color: #BA2121">, $\lambda = </span><span style="color: #BB6688; font-weight: bold">%.4f</span><span style="color: #BA2121">$&quot;</span> <span style="color: #666666">%</span> (key, _lambda))
plot_counter <span style="color: #666666">+=</span> <span style="color: #666666">1</span>
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
We see that LASSO reaches a good solution for low
values of \( \lambda \), but will "wither" when we increase \( \lambda \) too
much. Ridge is more stable over a larger range of values for
\( \lambda \), but eventually also fades away.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec18">Finding the optimal value of \( \lambda \) </h2>
<p>
To determine which value of \( \lambda \) is best we plot the accuracy of
the models when predicting the training and the testing set. We expect
the accuracy of the training set to be quite good, but if the accuracy
of the testing set is much lower this tells us that we might be
subject to an overfit model. The ideal scenario is an accuracy on the
testing set that is close to the accuracy of the training set.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>fig <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>figure(figsize<span style="color: #666666">=</span>(<span style="color: #666666">20</span>, <span style="color: #666666">14</span>))
colors <span style="color: #666666">=</span> {
<span style="color: #BA2121">&quot;ols_sk&quot;</span>: <span style="color: #BA2121">&quot;r&quot;</span>,
<span style="color: #BA2121">&quot;ridge_sk&quot;</span>: <span style="color: #BA2121">&quot;y&quot;</span>,
<span style="color: #BA2121">&quot;lasso_sk&quot;</span>: <span style="color: #BA2121">&quot;c&quot;</span>
}
<span style="color: #008000; font-weight: bold">for</span> key <span style="color: #AA22FF; font-weight: bold">in</span> train_errors:
plt<span style="color: #666666">.</span>semilogx(
lambdas,
train_errors[key],
colors[key],
label<span style="color: #666666">=</span><span style="color: #BA2121">&quot;Train {0}&quot;</span><span style="color: #666666">.</span>format(key),
linewidth<span style="color: #666666">=4.0</span>
)
<span style="color: #008000; font-weight: bold">for</span> key <span style="color: #AA22FF; font-weight: bold">in</span> test_errors:
plt<span style="color: #666666">.</span>semilogx(
lambdas,
test_errors[key],
colors[key] <span style="color: #666666">+</span> <span style="color: #BA2121">&quot;--&quot;</span>,
label<span style="color: #666666">=</span><span style="color: #BA2121">&quot;Test {0}&quot;</span><span style="color: #666666">.</span>format(key),
linewidth<span style="color: #666666">=4.0</span>
)
plt<span style="color: #666666">.</span>legend(loc<span style="color: #666666">=</span><span style="color: #BA2121">&quot;best&quot;</span>, fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>xlabel(<span style="color: #BA2121">r&quot;$\lambda$&quot;</span>, fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>ylabel(<span style="color: #BA2121">r&quot;$R^2$&quot;</span>, fontsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>tick_params(labelsize<span style="color: #666666">=18</span>)
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
From the above figure we can see that LASSO with \( \lambda = 10^{-2} \)
achieves a very good accuracy on the test set. This by far surpasses the
other models for all values of \( \lambda \).
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec19">Friday September 18: Intro to Logistic Regression </h2>
<p>
<!-- !split -->
<h2 id="___sec20">Logistic Regression </h2>
<p>
In linear regression our main interest was centered on learning the
coefficients of a functional fit (say a polynomial) in order to be
able to predict the response of a continuous variable on some unseen
data. The fit to the continuous variable \( y_i \) is based on some
independent variables \( \hat{x}_i \). Linear regression resulted in
analytical expressions for standard ordinary Least Squares or Ridge
regression (in terms of matrices to invert) for several quantities,
ranging from the variance and thereby the confidence intervals of the
parameters \( \hat{\beta} \) to the mean squared error. If we can invert
the product of the design matrices, linear regression gives then a
simple recipe for fitting our data.
<p>
<!-- !split -->
<h2 id="___sec21">Classification problems </h2>
<p>
Classification problems, however, are concerned with outcomes taking
the form of discrete variables (i.e. categories). We may for example,
on the basis of DNA sequencing for a number of patients, like to find
out which mutations are important for a certain disease; or based on
scans of various patients' brains, figure out if there is a tumor or
not; or given a specific physical system, we'd like to identify its
state, say whether it is an ordered or disordered system (typical
situation in solid state physics); or classify the status of a
patient, whether she/he has a stroke or not and many other similar
situations.
<p>
The most common situation we encounter when we apply logistic
regression is that of two possible outcomes, normally denoted as a
binary outcome, true or false, positive or negative, success or
failure etc.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec22">Optimization and Deep learning </h2>
<p>
Logistic regression will also serve as our stepping stone towards
neural network algorithms and supervised deep learning. For logistic
learning, the minimization of the cost function leads to a non-linear
equation in the parameters \( \hat{\beta} \). The optimization of the
problem calls therefore for minimization algorithms. This forms the
bottle neck of all machine learning algorithms, namely how to find
reliable minima of a multi-variable function. This leads us to the
family of gradient descent methods. The latter are the working horses
of basically all modern machine learning algorithms.
<p>
We note also that many of the topics discussed here on logistic
regression are also commonly used in modern supervised Deep Learning
models, as we will see later.
<p>
<!-- !split -->
<h2 id="___sec23">Basics </h2>
<p>
We consider the case where the dependent variables, also called the
responses or the outcomes, \( y_i \) are discrete and only take values
from \( k=0,\dots,K-1 \) (i.e. \( K \) classes).
<p>
The goal is to predict the
output classes from the design matrix \( \hat{X}\in\mathbb{R}^{n\times p} \)
made of \( n \) samples, each of which carries \( p \) features or predictors. The
primary goal is to identify the classes to which new unseen samples
belong.
<p>
Let us specialize to the case of two classes only, with outputs
\( y_i=0 \) and \( y_i=1 \). Our outcomes could represent the status of a
credit card user that could default or not on her/his credit card
debt. That is
$$
y_i = \begin{bmatrix} 0 & \mathrm{no}\\ 1 & \mathrm{yes} \end{bmatrix}.
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec24">Linear classifier </h2>
<p>
Before moving to the logistic model, let us try to use our linear
regression model to classify these two outcomes. We could for example
fit a linear model to the default case if \( y_i > 0.5 \) and the no
default case \( y_i \leq 0.5 \).
<p>
We would then have our
weighted linear combination, namely
$$
\begin{equation}
\hat{y} = \hat{X}^T\hat{\beta} + \hat{\epsilon},
\label{_auto13}
\end{equation}
$$
where \( \hat{y} \) is a vector representing the possible outcomes, \( \hat{X} \) is our
\( n\times p \) design matrix and \( \hat{\beta} \) represents our estimators/predictors.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec25">Some selected properties </h2>
<p>
The main problem with our function is that it takes values on the
entire real axis. In the case of logistic regression, however, the
labels \( y_i \) are discrete variables. A typical example is the credit
card data discussed below here, where we can set the state of
defaulting the debt to \( y_i=1 \) and not to \( y_i=0 \) for one the persons
in the data set (see the full example below).
<p>
One simple way to get a discrete output is to have sign
functions that map the output of a linear regressor to values \( \{0,1\} \),
\( f(s_i)=sign(s_i)=1 \) if \( s_i\ge 0 \) and 0 if otherwise.
We will encounter this model in our first demonstration of neural networks. Historically it is called the &quot;perceptron" model in the machine learning
literature. This model is extremely simple. However, in many cases it is more
favorable to use a ``soft" classifier that outputs
the probability of a given category. This leads us to the logistic function.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec26">The logistic function </h2>
<p>
The perceptron is an example of a ``hard classification&quot; model. We
will encounter this model when we discuss neural networks as
well. Each datapoint is deterministically assigned to a category (i.e
\( y_i=0 \) or \( y_i=1 \)). In many cases, it is favorable to have a &quot;soft&quot;
classifier that outputs the probability of a given category rather
than a single value. For example, given \( x_i \), the classifier
outputs the probability of being in a category \( k \). Logistic regression
is the most common example of a so-called soft classifier. In logistic
regression, the probability that a data point \( x_i \)
belongs to a category \( y_i=\{0,1\} \) is given by the so-called logit function (or Sigmoid) which is meant to represent the likelihood for a given event,
$$
p(t) = \frac{1}{1+\mathrm \exp{-t}}=\frac{\exp{t}}{1+\mathrm \exp{t}}.
$$
Note that \( 1-p(t)= p(-t) \).
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec27">Examples of likelihood functions used in logistic regression and nueral networks </h2>
<p>
The following code plots the logistic function, the step function and other functions we will encounter from here and on.
<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: #BA2121; font-style: italic">&quot;&quot;&quot;The sigmoid function (or the logistic curve) is a</span>
<span style="color: #BA2121; font-style: italic">function that takes any real number, z, and outputs a number (0,1).</span>
<span style="color: #BA2121; font-style: italic">It is useful in neural networks for assigning weights on a relative scale.</span>
<span style="color: #BA2121; font-style: italic">The value z is the weighted sum of parameters involved in the learning algorithm.&quot;&quot;&quot;</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">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">import</span> <span style="color: #0000FF; font-weight: bold">math</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">mt</span>
z <span style="color: #666666">=</span> numpy<span style="color: #666666">.</span>arange(<span style="color: #666666">-5</span>, <span style="color: #666666">5</span>, <span style="color: #666666">.1</span>)
sigma_fn <span style="color: #666666">=</span> numpy<span style="color: #666666">.</span>vectorize(<span style="color: #008000; font-weight: bold">lambda</span> z: <span style="color: #666666">1/</span>(<span style="color: #666666">1+</span>numpy<span style="color: #666666">.</span>exp(<span style="color: #666666">-</span>z)))
sigma <span style="color: #666666">=</span> sigma_fn(z)
fig <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>figure()
ax <span style="color: #666666">=</span> fig<span style="color: #666666">.</span>add_subplot(<span style="color: #666666">111</span>)
ax<span style="color: #666666">.</span>plot(z, sigma)
ax<span style="color: #666666">.</span>set_ylim([<span style="color: #666666">-0.1</span>, <span style="color: #666666">1.1</span>])
ax<span style="color: #666666">.</span>set_xlim([<span style="color: #666666">-5</span>,<span style="color: #666666">5</span>])
ax<span style="color: #666666">.</span>grid(<span style="color: #008000">True</span>)
ax<span style="color: #666666">.</span>set_xlabel(<span style="color: #BA2121">&#39;z&#39;</span>)
ax<span style="color: #666666">.</span>set_title(<span style="color: #BA2121">&#39;sigmoid function&#39;</span>)
plt<span style="color: #666666">.</span>show()
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;Step Function&quot;&quot;&quot;</span>
z <span style="color: #666666">=</span> numpy<span style="color: #666666">.</span>arange(<span style="color: #666666">-5</span>, <span style="color: #666666">5</span>, <span style="color: #666666">.02</span>)
step_fn <span style="color: #666666">=</span> numpy<span style="color: #666666">.</span>vectorize(<span style="color: #008000; font-weight: bold">lambda</span> z: <span style="color: #666666">1.0</span> <span style="color: #008000; font-weight: bold">if</span> z <span style="color: #666666">&gt;=</span> <span style="color: #666666">0.0</span> <span style="color: #008000; font-weight: bold">else</span> <span style="color: #666666">0.0</span>)
step <span style="color: #666666">=</span> step_fn(z)
fig <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>figure()
ax <span style="color: #666666">=</span> fig<span style="color: #666666">.</span>add_subplot(<span style="color: #666666">111</span>)
ax<span style="color: #666666">.</span>plot(z, step)
ax<span style="color: #666666">.</span>set_ylim([<span style="color: #666666">-0.5</span>, <span style="color: #666666">1.5</span>])
ax<span style="color: #666666">.</span>set_xlim([<span style="color: #666666">-5</span>,<span style="color: #666666">5</span>])
ax<span style="color: #666666">.</span>grid(<span style="color: #008000">True</span>)
ax<span style="color: #666666">.</span>set_xlabel(<span style="color: #BA2121">&#39;z&#39;</span>)
ax<span style="color: #666666">.</span>set_title(<span style="color: #BA2121">&#39;step function&#39;</span>)
plt<span style="color: #666666">.</span>show()
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;tanh Function&quot;&quot;&quot;</span>
z <span style="color: #666666">=</span> numpy<span style="color: #666666">.</span>arange(<span style="color: #666666">-2*</span>mt<span style="color: #666666">.</span>pi, <span style="color: #666666">2*</span>mt<span style="color: #666666">.</span>pi, <span style="color: #666666">0.1</span>)
t <span style="color: #666666">=</span> numpy<span style="color: #666666">.</span>tanh(z)
fig <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>figure()
ax <span style="color: #666666">=</span> fig<span style="color: #666666">.</span>add_subplot(<span style="color: #666666">111</span>)
ax<span style="color: #666666">.</span>plot(z, t)
ax<span style="color: #666666">.</span>set_ylim([<span style="color: #666666">-1.0</span>, <span style="color: #666666">1.0</span>])
ax<span style="color: #666666">.</span>set_xlim([<span style="color: #666666">-2*</span>mt<span style="color: #666666">.</span>pi,<span style="color: #666666">2*</span>mt<span style="color: #666666">.</span>pi])
ax<span style="color: #666666">.</span>grid(<span style="color: #008000">True</span>)
ax<span style="color: #666666">.</span>set_xlabel(<span style="color: #BA2121">&#39;z&#39;</span>)
ax<span style="color: #666666">.</span>set_title(<span style="color: #BA2121">&#39;tanh function&#39;</span>)
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec28">Two parameters </h2>
<p>
We assume now that we have two classes with \( y_i \) either \( 0 \) or \( 1 \). Furthermore we assume also that we have only two parameters \( \beta \) in our fitting of the Sigmoid function, that is we define probabilities
$$
\begin{align*}
p(y_i=1|x_i,\hat{\beta}) &= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\
p(y_i=0|x_i,\hat{\beta}) &= 1 - p(y_i=1|x_i,\hat{\beta}),
\end{align*}
$$
where \( \hat{\beta} \) are the weights we wish to extract from data, in our case \( \beta_0 \) and \( \beta_1 \).
<p>
Note that we used
$$
p(y_i=0\vert x_i, \hat{\beta}) = 1-p(y_i=1\vert x_i, \hat{\beta}).
$$
<p>
<!-- !split -->
<h2 id="___sec29">Maximum likelihood </h2>
<p>
In order to define the total likelihood for all possible outcomes from a
dataset \( \mathcal{D}=\{(y_i,x_i)\} \), with the binary labels
\( y_i\in\{0,1\} \) and where the data points are drawn independently, we use the so-called <a href="https://en.wikipedia.org/wiki/Maximum_likelihood_estimation" target="_blank">Maximum Likelihood Estimation</a> (MLE) principle.
We aim thus at maximizing
the probability of seeing the observed data. We can then approximate the
likelihood in terms of the product of the individual probabilities of a specific outcome \( y_i \), that is
$$
\begin{align*}
P(\mathcal{D}|\hat{\beta})& = \prod_{i=1}^n \left[p(y_i=1|x_i,\hat{\beta})\right]^{y_i}\left[1-p(y_i=1|x_i,\hat{\beta}))\right]^{1-y_i}\nonumber \\
\end{align*}
$$
from which we obtain the log-likelihood and our <b>cost/loss</b> function
$$
\mathcal{C}(\hat{\beta}) = \sum_{i=1}^n \left( y_i\log{p(y_i=1|x_i,\hat{\beta})} + (1-y_i)\log\left[1-p(y_i=1|x_i,\hat{\beta}))\right]\right).
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec30">The cost function rewritten </h2>
<p>
Reordering the logarithms, we can rewrite the <b>cost/loss</b> function as
$$
\mathcal{C}(\hat{\beta}) = \sum_{i=1}^n \left(y_i(\beta_0+\beta_1x_i) -\log{(1+\exp{(\beta_0+\beta_1x_i)})}\right).
$$
<p>
The maximum likelihood estimator is defined as the set of parameters that maximize the log-likelihood where we maximize with respect to \( \beta \).
Since the cost (error) function is just the negative log-likelihood, for logistic regression we have that
$$
\mathcal{C}(\hat{\beta})=-\sum_{i=1}^n \left(y_i(\beta_0+\beta_1x_i) -\log{(1+\exp{(\beta_0+\beta_1x_i)})}\right).
$$
This equation is known in statistics as the <b>cross entropy</b>. Finally, we note that just as in linear regression,
in practice we often supplement the cross-entropy with additional regularization terms, usually \( L_1 \) and \( L_2 \) regularization as we did for Ridge and Lasso regression.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec31">Minimizing the cross entropy </h2>
<p>
The cross entropy is a convex function of the weights \( \hat{\beta} \) and,
therefore, any local minimizer is a global minimizer.
<p>
Minimizing this
cost function with respect to the two parameters \( \beta_0 \) and \( \beta_1 \) we obtain
$$
\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \beta_0} = -\sum_{i=1}^n \left(y_i -\frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}}\right),
$$
and
$$
\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \beta_1} = -\sum_{i=1}^n \left(y_ix_i -x_i\frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}}\right).
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec32">A more compact expression </h2>
<p>
Let us now define a vector \( \hat{y} \) with \( n \) elements \( y_i \), an
\( n\times p \) matrix \( \hat{X} \) which contains the \( x_i \) values and a
vector \( \hat{p} \) of fitted probabilities \( p(y_i\vert x_i,\hat{\beta}) \). We can rewrite in a more compact form the first
derivative of cost function as
$$
\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}} = -\hat{X}^T\left(\hat{y}-\hat{p}\right).
$$
<p>
If we in addition define a diagonal matrix \( \hat{W} \) with elements
\( p(y_i\vert x_i,\hat{\beta})(1-p(y_i\vert x_i,\hat{\beta}) \), we can obtain a compact expression of the second derivative as
$$
\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T} = \hat{X}^T\hat{W}\hat{X}.
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec33">Extending to more predictors </h2>
<p>
Within a binary classification problem, we can easily expand our model to include multiple predictors. Our ratio between likelihoods is then with \( p \) predictors
$$
\log{ \frac{p(\hat{\beta}\hat{x})}{1-p(\hat{\beta}\hat{x})}} = \beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p.
$$
Here we defined \( \hat{x}=[1,x_1,x_2,\dots,x_p] \) and \( \hat{\beta}=[\beta_0, \beta_1, \dots, \beta_p] \) leading to
$$
p(\hat{\beta}\hat{x})=\frac{ \exp{(\beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p)}}{1+\exp{(\beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p)}}.
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec34">Including more classes </h2>
<p>
Till now we have mainly focused on two classes, the so-called binary
system. Suppose we wish to extend to \( K \) classes. Let us for the sake
of simplicity assume we have only two predictors. We have then
following model
$$
\log{\frac{p(C=1\vert x)}{p(K\vert x)}} = \beta_{10}+\beta_{11}x_1,
$$
$$
\log{\frac{p(C=2\vert x)}{p(K\vert x)}} = \beta_{20}+\beta_{21}x_1,
$$
and so on till the class \( C=K-1 \) class
$$
\log{\frac{p(C=K-1\vert x)}{p(K\vert x)}} = \beta_{(K-1)0}+\beta_{(K-1)1}x_1,
$$
<p>
and the model is specified in term of \( K-1 \) so-called log-odds or
<b>logit</b> transformations.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec35">More classes </h2>
<p>
In our discussion of neural networks we will encounter the above again
in terms of a slightly modified function, the so-called <b>Softmax</b> function.
<p>
The softmax function is used in various multiclass classification
methods, such as multinomial logistic regression (also known as
softmax regression), multiclass linear discriminant analysis, naive
Bayes classifiers, and artificial neural networks. Specifically, in
multinomial logistic regression and linear discriminant analysis, the
input to the function is the result of \( K \) distinct linear functions,
and the predicted probability for the \( k \)-th class given a sample
vector \( \hat{x} \) and a weighting vector \( \hat{\beta} \) is (with two
predictors):
$$
p(C=k\vert \mathbf {x} )=\frac{\exp{(\beta_{k0}+\beta_{k1}x_1)}}{1+\sum_{l=1}^{K-1}\exp{(\beta_{l0}+\beta_{l1}x_1)}}.
$$
It is easy to extend to more predictors. The final class is
$$
p(C=K\vert \mathbf {x} )=\frac{1}{1+\sum_{l=1}^{K-1}\exp{(\beta_{l0}+\beta_{l1}x_1)}},
$$
<p>
and they sum to one. Our earlier discussions were all specialized to
the case with two classes only. It is easy to see from the above that
what we derived earlier is compatible with these equations.
<p>
To find the optimal parameters we would typically use a gradient
descent method. Newton's method and gradient descent methods are
discussed in the material on <a href="https://compphysics.github.io/MachineLearning/doc/pub/Splines/html/Splines-bs.html" target="_blank">optimization
methods</a>.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec36">A simple classification problem </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">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</span> <span style="color: #008000; font-weight: bold">import</span> datasets, linear_model
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">matplotlib.pyplot</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">plt</span>
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">generate_data</span>():
np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>seed(<span style="color: #666666">0</span>)
X, y <span style="color: #666666">=</span> datasets<span style="color: #666666">.</span>make_moons(<span style="color: #666666">200</span>, noise<span style="color: #666666">=0.20</span>)
<span style="color: #008000; font-weight: bold">return</span> X, y
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">visualize</span>(X, y, clf):
plot_decision_boundary(<span style="color: #008000; font-weight: bold">lambda</span> x: clf<span style="color: #666666">.</span>predict(x), X, y)
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">plot_decision_boundary</span>(pred_func, X, y):
<span style="color: #408080; font-style: italic"># Set min and max values and give it some padding</span>
x_min, x_max <span style="color: #666666">=</span> X[:, <span style="color: #666666">0</span>]<span style="color: #666666">.</span>min() <span style="color: #666666">-</span> <span style="color: #666666">.5</span>, X[:, <span style="color: #666666">0</span>]<span style="color: #666666">.</span>max() <span style="color: #666666">+</span> <span style="color: #666666">.5</span>
y_min, y_max <span style="color: #666666">=</span> X[:, <span style="color: #666666">1</span>]<span style="color: #666666">.</span>min() <span style="color: #666666">-</span> <span style="color: #666666">.5</span>, X[:, <span style="color: #666666">1</span>]<span style="color: #666666">.</span>max() <span style="color: #666666">+</span> <span style="color: #666666">.5</span>
h <span style="color: #666666">=</span> <span style="color: #666666">0.01</span>
<span style="color: #408080; font-style: italic"># Generate a grid of points with distance h between them</span>
xx, yy <span style="color: #666666">=</span> np<span style="color: #666666">.</span>meshgrid(np<span style="color: #666666">.</span>arange(x_min, x_max, h), np<span style="color: #666666">.</span>arange(y_min, y_max, h))
<span style="color: #408080; font-style: italic"># Predict the function value for the whole gid</span>
Z <span style="color: #666666">=</span> pred_func(np<span style="color: #666666">.</span>c_[xx<span style="color: #666666">.</span>ravel(), yy<span style="color: #666666">.</span>ravel()])
Z <span style="color: #666666">=</span> Z<span style="color: #666666">.</span>reshape(xx<span style="color: #666666">.</span>shape)
<span style="color: #408080; font-style: italic"># Plot the contour and training examples</span>
plt<span style="color: #666666">.</span>contourf(xx, yy, Z, cmap<span style="color: #666666">=</span>plt<span style="color: #666666">.</span>cm<span style="color: #666666">.</span>Spectral)
plt<span style="color: #666666">.</span>scatter(X[:, <span style="color: #666666">0</span>], X[:, <span style="color: #666666">1</span>], c<span style="color: #666666">=</span>y, cmap<span style="color: #666666">=</span>plt<span style="color: #666666">.</span>cm<span style="color: #666666">.</span>Spectral)
plt<span style="color: #666666">.</span>show()
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">classify</span>(X, y):
clf <span style="color: #666666">=</span> linear_model<span style="color: #666666">.</span>LogisticRegressionCV()
clf<span style="color: #666666">.</span>fit(X, y)
<span style="color: #008000; font-weight: bold">return</span> clf
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">main</span>():
X, y <span style="color: #666666">=</span> generate_data()
<span style="color: #408080; font-style: italic"># visualize(X, y)</span>
clf <span style="color: #666666">=</span> classify(X, y)
visualize(X, y, clf)
<span style="color: #008000; font-weight: bold">if</span> <span style="color: #19177C">__name__</span> <span style="color: #666666">==</span> <span style="color: #BA2121">&quot;__main__&quot;</span>:
main()
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec37">Cancer Data again now with Decision Trees and other Methods </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">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">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.datasets</span> <span style="color: #008000; font-weight: bold">import</span> load_breast_cancer
<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> LogisticRegression
<span style="color: #408080; font-style: italic"># Load the data</span>
cancer <span style="color: #666666">=</span> load_breast_cancer()
X_train, X_test, y_train, y_test <span style="color: #666666">=</span> train_test_split(cancer<span style="color: #666666">.</span>data,cancer<span style="color: #666666">.</span>target,random_state<span style="color: #666666">=0</span>)
<span style="color: #008000; font-weight: bold">print</span>(X_train<span style="color: #666666">.</span>shape)
<span style="color: #008000; font-weight: bold">print</span>(X_test<span style="color: #666666">.</span>shape)
<span style="color: #408080; font-style: italic"># Logistic Regression</span>
logreg <span style="color: #666666">=</span> LogisticRegression(solver<span style="color: #666666">=</span><span style="color: #BA2121">&#39;lbfgs&#39;</span>)
logreg<span style="color: #666666">.</span>fit(X_train, y_train)
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Test set accuracy with Logistic Regression: {:.2f}&quot;</span><span style="color: #666666">.</span>format(logreg<span style="color: #666666">.</span>score(X_test,y_test)))
<span style="color: #408080; font-style: italic">#now scale the data</span>
<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> StandardScaler
scaler <span style="color: #666666">=</span> StandardScaler()
scaler<span style="color: #666666">.</span>fit(X_train)
X_train_scaled <span style="color: #666666">=</span> scaler<span style="color: #666666">.</span>transform(X_train)
X_test_scaled <span style="color: #666666">=</span> scaler<span style="color: #666666">.</span>transform(X_test)
<span style="color: #408080; font-style: italic"># Logistic Regression</span>
logreg<span style="color: #666666">.</span>fit(X_train_scaled, y_train)
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Test set accuracy Logistic Regression with scaled data: {:.2f}&quot;</span><span style="color: #666666">.</span>format(logreg<span style="color: #666666">.</span>score(X_test_scaled,y_test)))
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec38">Other measures in classification studies: Cancer Data again </h2>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">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">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.datasets</span> <span style="color: #008000; font-weight: bold">import</span> load_breast_cancer
<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> LogisticRegression
<span style="color: #408080; font-style: italic"># Load the data</span>
cancer <span style="color: #666666">=</span> load_breast_cancer()
X_train, X_test, y_train, y_test <span style="color: #666666">=</span> train_test_split(cancer<span style="color: #666666">.</span>data,cancer<span style="color: #666666">.</span>target,random_state<span style="color: #666666">=0</span>)
<span style="color: #008000; font-weight: bold">print</span>(X_train<span style="color: #666666">.</span>shape)
<span style="color: #008000; font-weight: bold">print</span>(X_test<span style="color: #666666">.</span>shape)
<span style="color: #408080; font-style: italic"># Logistic Regression</span>
logreg <span style="color: #666666">=</span> LogisticRegression(solver<span style="color: #666666">=</span><span style="color: #BA2121">&#39;lbfgs&#39;</span>)
logreg<span style="color: #666666">.</span>fit(X_train, y_train)
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Test set accuracy with Logistic Regression: {:.2f}&quot;</span><span style="color: #666666">.</span>format(logreg<span style="color: #666666">.</span>score(X_test,y_test)))
<span style="color: #408080; font-style: italic">#now scale the data</span>
<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> StandardScaler
scaler <span style="color: #666666">=</span> StandardScaler()
scaler<span style="color: #666666">.</span>fit(X_train)
X_train_scaled <span style="color: #666666">=</span> scaler<span style="color: #666666">.</span>transform(X_train)
X_test_scaled <span style="color: #666666">=</span> scaler<span style="color: #666666">.</span>transform(X_test)
<span style="color: #408080; font-style: italic"># Logistic Regression</span>
logreg<span style="color: #666666">.</span>fit(X_train_scaled, y_train)
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Test set accuracy Logistic Regression with scaled data: {:.2f}&quot;</span><span style="color: #666666">.</span>format(logreg<span style="color: #666666">.</span>score(X_test_scaled,y_test)))
<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> LabelEncoder
<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_validate
<span style="color: #408080; font-style: italic">#Cross validation</span>
accuracy <span style="color: #666666">=</span> cross_validate(logreg,X_test_scaled,y_test,cv<span style="color: #666666">=10</span>)[<span style="color: #BA2121">&#39;test_score&#39;</span>]
<span style="color: #008000; font-weight: bold">print</span>(accuracy)
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Test set accuracy with Logistic Regression and scaled data: {:.2f}&quot;</span><span style="color: #666666">.</span>format(logreg<span style="color: #666666">.</span>score(X_test_scaled,y_test)))
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">scikitplot</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">skplt</span>
y_pred <span style="color: #666666">=</span> logreg<span style="color: #666666">.</span>predict(X_test_scaled)
skplt<span style="color: #666666">.</span>metrics<span style="color: #666666">.</span>plot_confusion_matrix(y_test, y_pred, normalize<span style="color: #666666">=</span><span style="color: #008000">True</span>)
plt<span style="color: #666666">.</span>show()
y_probas <span style="color: #666666">=</span> logreg<span style="color: #666666">.</span>predict_proba(X_test_scaled)
skplt<span style="color: #666666">.</span>metrics<span style="color: #666666">.</span>plot_roc(y_test, y_probas)
plt<span style="color: #666666">.</span>show()
skplt<span style="color: #666666">.</span>metrics<span style="color: #666666">.</span>plot_cumulative_gain(y_test, y_probas)
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
<!-- ------------------- end of main content --------------- -->
<center style="font-size:80%">
<!-- copyright --> &copy; 1999-2020, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license
</center>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink