Files
FYS-STK4155/doc/pub/LogReg/html/LogReg.html
T
2018-09-26 08:27:07 +02:00

477 lines
23 KiB
HTML

<!--
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="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': [('Logistic Regression', 2, None, '___sec0'),
('Optimization and Deep learning', 2, None, '___sec1'),
('Basics', 2, None, '___sec2'),
('Linear classifier', 2, None, '___sec3'),
('Some selected properties', 2, None, '___sec4'),
('The logistic function', 2, None, '___sec5'),
('Two parameters', 2, None, '___sec6'),
('Maximum likelihood', 2, None, '___sec7'),
('The cost function rewritten', 2, None, '___sec8'),
('Minimizing the cross entropy', 2, None, '___sec9'),
('A more compact expression', 2, None, '___sec10'),
('Extending to more predictors', 2, None, '___sec11'),
('Including more classes', 2, None, '___sec12'),
('Optimizing the cost function', 2, None, '___sec13'),
('A _scikit-learn_ example', 2, None, '___sec14'),
('A simple classification problem', 2, None, '___sec15')]}
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 26, 2018</h4></center> <!-- date -->
<br>
<p>
<!-- !split -->
<h2 id="___sec0">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 (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>
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="___sec1">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 optmization 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
regression are also commonly used in modern supervised Deep Learning
models, as we will see later.
<p>
<!-- !split -->
<h2 id="___sec2">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 who 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="___sec3">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{_auto1}
\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="___sec4">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.
<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>
The code for plotting the perceptron can be seen here. This si nothing but the standard <a href="https://en.wikipedia.org/wiki/Heaviside_step_function" target="_blank">Heaviside step function</a>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec5">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) \).
The following code plots the logistic function.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec6">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="___sec7">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="___sec8">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="___sec9">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="___sec10">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="___sec11">Extending to more predictors </h2>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec12">Including more classes </h2>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec13">Optimizing the cost function </h2>
<p>
Newton's method and gradient descent methods
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec14">A <b>scikit-learn</b> example </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</span> <span style="color: #008000; font-weight: bold">import</span> datasets
iris <span style="color: #666666">=</span> datasets<span style="color: #666666">.</span>load_iris()
<span style="color: #008000">list</span>(iris<span style="color: #666666">.</span>keys())
[<span style="color: #BA2121">&#39;data&#39;</span>, <span style="color: #BA2121">&#39;target_names&#39;</span>, <span style="color: #BA2121">&#39;feature_names&#39;</span>, <span style="color: #BA2121">&#39;target&#39;</span>, <span style="color: #BA2121">&#39;DESCR&#39;</span>]
X <span style="color: #666666">=</span> iris[<span style="color: #BA2121">&quot;data&quot;</span>][:, <span style="color: #666666">3</span>:] <span style="color: #408080; font-style: italic"># petal width</span>
y <span style="color: #666666">=</span> (iris[<span style="color: #BA2121">&quot;target&quot;</span>] <span style="color: #666666">==</span> <span style="color: #666666">2</span>)<span style="color: #666666">.</span>astype(np<span style="color: #666666">.</span>int) <span style="color: #408080; font-style: italic"># 1 if Iris-Virginica, else 0</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> LogisticRegression
log_reg <span style="color: #666666">=</span> LogisticRegression()
log_reg<span style="color: #666666">.</span>fit(X, y)
X_new <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linspace(<span style="color: #666666">0</span>, <span style="color: #666666">3</span>, <span style="color: #666666">1000</span>)<span style="color: #666666">.</span>reshape(<span style="color: #666666">-1</span>, <span style="color: #666666">1</span>)
y_proba <span style="color: #666666">=</span> log_reg<span style="color: #666666">.</span>predict_proba(X_new)
plt<span style="color: #666666">.</span>plot(X_new, y_proba[:, <span style="color: #666666">1</span>], <span style="color: #BA2121">&quot;g-&quot;</span>, label<span style="color: #666666">=</span><span style="color: #BA2121">&quot;Iris-Virginica&quot;</span>)
plt<span style="color: #666666">.</span>plot(X_new, y_proba[:, <span style="color: #666666">0</span>], <span style="color: #BA2121">&quot;b--&quot;</span>, label<span style="color: #666666">=</span><span style="color: #BA2121">&quot;Not Iris-Virginica&quot;</span>)
plt<span style="color: #666666">.</span>show()
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec15">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):
<span style="color: #408080; font-style: italic"># plt.scatter(X[:, 0], X[:, 1], s=40, c=y, cmap=plt.cm.Spectral)</span>
<span style="color: #408080; font-style: italic"># plt.show()</span>
plot_decision_boundary(<span style="color: #008000; font-weight: bold">lambda</span> x: clf<span style="color: #666666">.</span>predict(x), X, y)
plt<span style="color: #666666">.</span>title(<span style="color: #BA2121">&quot;Logistic Regression&quot;</span>)
<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>
<!-- ------------------- end of main content --------------- -->
<center style="font-size:80%">
<!-- copyright --> &copy; 1999-2018, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license
</center>
</body>
</html>