Update with log regression

This commit is contained in:
mhjensen
2018-09-14 05:36:58 +02:00
parent 73474229be
commit df4a9d7f27
137 changed files with 79934 additions and 15967 deletions
+313 -80
View File
@@ -2196,26 +2196,26 @@ need for bootstrapping.
!split
===== Resampling methods: Jackknife =====
The Jackknife works by making many replicas of the estimator $\widehat{\vec{\theta}}$.
The jackknife is a resampling method, we explained that this happens by scrambling the data in some way. When using the jackknife, this is done by systematically leaving out one observation from the vector of observed values $\vec{X} = (X_1,X_2,\cdots,X_n)$.
Let $\vec{X}_i$ denote the vector
The Jackknife works by making many replicas of the estimator $\widehat{\theta}$.
The jackknife is a resampling method, we explained that this happens by scrambling the data in some way. When using the jackknife, this is done by systematically leaving out one observation from the vector of observed values $\hat{x} = (x_1,x_2,\cdots,X_n)$.
Let $\hat{x}_i$ denote the vector
!bt
\[
\vec{X}_i = (X_1,X_2,\cdots,X_{i-1},X_{i+1},\cdots,X_n),
\hat{x}_i = (x_1,x_2,\cdots,x_{i-1},x_{i+1},\cdots,x_n),
\]
!et
which equals the vector $\vec{X}$ with the exception that observation
which equals the vector $\hat{x}$ with the exception that observation
number $i$ is left out. Using this notation, define
$\widehat{\vec{\theta}}_i$ to be the estimator
$\widehat{\vec{\theta}}$ computed using $\vec{X}_i$.
$\widehat{\theta}_i$ to be the estimator
$\widehat{\theta}$ computed using $\vec{X}_i$.
!split
===== Resampling methods: Jackknife estimator =====
To get an estimate for the bias and
standard error of $\widehat{\vec{\theta}}$, use the following
estimators for each component of $\widehat{\vec{\theta}}$
standard error of $\widehat{\theta}$, use the following
estimators for each component of $\widehat{\theta}$
!bt
\[
@@ -2225,38 +2225,6 @@ estimators for each component of $\widehat{\vec{\theta}}$
!split
===== Resampling methods: Jackknife sample code =====
Sample code for the Jackknife method
!bc pycod
def jack(data, stat):
n = len(data); t = zeros(n); inds = arange(n); t0 = time()
# 'jackknifing' by leaving out an observation for each i
for i in range(n):
t[i] = stat(delete(data,i) )
return t
# define a function which returns your chosen estimator theta-hat
def stat(data):
theta-hat = mean(data)
return theta-hat
# Return the Jackknife sample
t = jack(X, stat)
!ec
Consider first the function _jack()_. This function repeatedly
estimates the function called _statistic()_ under the resampled
data by systematically leaving out one observation from the data. The
function _stat()_ is passed as an argument to
_jack()_. The array _t_ is eventually returned, which
contains all the estimates $\widehat{\vec{\theta}}$, and can be
plotted or analysed in other ways, such as by calling _std(t)_
from _numpy_ to estimate the standard error of
$\widehat{\vec{\theta}}$. The function _std(t)_ is just the
estimator $\widehat{\sigma}^2$.
!split
===== Resampling methods: Bootstrap =====
@@ -2275,32 +2243,32 @@ o It is relatively simple to apply the bootstrap to complex data-collection plan
!split
===== Resampling methods: Bootstrap background =====
Since $\widehat{\vec{\theta}} = \widehat{\vec{\theta}}(\vec{X})$ is a function of random variables,
$\widehat{\vec{\theta}}$ itself must be a random variable. Thus it has
Since $\widehat{\theta} = \widehat{\theta}(\hat{X})$ is a function of random variables,
$\widehat{\theta}$ itself must be a random variable. Thus it has
a pdf, call this function $p(\vec{t})$. The aim of the bootstrap is to
estimate $p(\vec{t})$ by the relative frequency of
$\widehat{\vec{\theta}}$. You can think of this as using a histogram
in the place of $p(\vec{t})$. If the relative frequency closely
estimate $p(\hat{t})$ by the relative frequency of
$\widehat{\theta}$. You can think of this as using a histogram
in the place of $p(\hat{t})$. If the relative frequency closely
resembles $p(\vec{t})$, then using numerics, it is straight forward to
estimate all the interesting parameters of $p(\vec{t})$ using point
estimate all the interesting parameters of $p(\hat{t})$ using point
estimators.
!split
===== Resampling methods: More Bootstrap background =====
In the case that $\widehat{\vec{\theta}}$ has
In the case that $\widehat{\theta}$ has
more than one component, and the components are independent, use the
same estimator on each component separately. If the probability
density function of $X_i$, $p(x)$, had been known, then it would have
been straight forward to do this by:
o Drawing lots of numbers from $p(x)$, suppose we call one such set of numbers $(X_1^*, X_2^*, \cdots, X_n^*)$.
o Then using these numbers, we could compute a replica of $\widehat{\vec{\theta}}$ called $\widehat{\vec{\theta}}^*$.
o Then using these numbers, we could compute a replica of $\widehat{\theta}$ called $\widehat{\theta}^*$.
By repeated use of (1) and (2), many
estimates of $\widehat{\vec{\theta}}$ could have been obtained. The
idea is to use the relative frequency of $\widehat{\vec{\theta}}^*$
(think of a histogram) as an estimate of $p(\vec{t})$.
estimates of $\widehat{\theta}$ could have been obtained. The
idea is to use the relative frequency of $\widehat{\theta}^*$
(think of a histogram) as an estimate of $p(\hat{t})$.
!split
===== Resampling methods: Bootstrap approach =====
@@ -2318,39 +2286,16 @@ result in some asymptotic sense? The answer is yes.
Instead of generating the histogram for the relative
frequency of the observation $X_i$, just draw the values
$(X_1^*,X_2^*,\cdots,X_n^*)$ with replacement from the vector
$\vec{X}$.
!split
===== Resampling methods: Bootstrap algorithm =====
!bc pycod
def boot(data, statistic, R):
t = zeros(R); n = len(data); inds = arange(n); t0 = time()
for i in range(R):
t[i] = statistic(data[randint(0,n,n)])
return t
# define a function which returns your chosen estimator theta-hat
def stat(data):
theta-hat = mean(data)
return theta-hat
t = boot(X, stat, 2**9)
!ec
Consider first the function _boot()_. In the _for_ loop, this function repeatedly estimates the function called _statistic()_ under the resampled data in _data(randint(0,n,n))_. The function _statistic()_ is passed as an argument to _boot()_. The array _t_ is eventually returned, which contains all the estimates $\widehat{\vec{\theta}}$, and can be plotted or analysed in other ways, such as by calling _std(t)_ from _numpy_ to estimate the standard error of $\widehat{\vec{\theta}}$. The function _std(t)_ is just the estimator $\widehat{\sigma}^2$.
$\hat{X}$.
!split
===== Resampling methods: Bootstrap steps =====
The independent bootstrap works like this:
o Draw with replacement $n$ numbers for the observed variables $\vec{x} = (x_1,x_2,\cdots,x_n)$.
o Define a vector $\vec{x}^*$ containing the values which were drawn from $\vec{x}$.
o Using the vector $\vec{x}^*$ compute $\widehat{\theta}^*$ by evaluating $\widehat \theta$ under the observations $\vec{x}^*$.
o Draw with replacement $n$ numbers for the observed variables $\hat{x} = (x_1,x_2,\cdots,x_n)$.
o Define a vector $\hat{x}^*$ containing the values which were drawn from $\hat{x}$.
o Using the vector $\hat{x}^*$ compute $\widehat{\theta}^*$ by evaluating $\widehat \theta$ under the observations $\hat{x}^*$.
o Repeat this process $k$ times.
When you are done, you can draw a histogram of the relative frequency of $\widehat \theta^*$. This is your estimate of the probability distribution $p(t)$. Using this probability distribution you can estimate any statistics thereof. In principle you never draw the histogram of the relative frequency of $\widehat{\theta}^*$. Instead you use the estimators corresponding to the statistic of interest. For example, if you are interested in estimating the variance of $\widehat \theta$, apply the esimator $\widehat \sigma^2$ to the values $\widehat \theta ^*$.
@@ -2369,7 +2314,7 @@ Assume $n = 2^d$ for some integer $d>1$ and $X_1,X_2,\cdots, X_n$ is a stationar
Moreover, assume that the time series is asymptotically uncorrelated. We switch to vector notation by arranging $X_1,X_2,\cdots,X_n$ in an $n$-tuple. Define:
!bt
\begin{align*}
\vec{X} = (X_1,X_2,\cdots,X_n).
\hat{X} = (X_1,X_2,\cdots,X_n).
\end{align*}
!et
@@ -2412,6 +2357,7 @@ $\vec{X}_i$. It follows by induction that $n_i = n/2^i$.
!split
===== Blocking Transformations =====
Using the
definition of the blocking transformation and the distributive
property of the covariance, it is clear that since $h =|i-j|$
@@ -2700,6 +2646,293 @@ dataAnalysis.printOutput()
!ec
!split
===== The bias-variance tradeoff =====
We begin with an unknown function $y=f(x)$ and fix a \emph{hypothesis set}
$\mathcal{H}$ consisting of all functions we are willing to consider,
defined also on the domain of $f$. This set may be uncountably
infinite (e.g.~if there are real-valued parameters to fit).
The
choice of which functions to include in $\mathcal{H}$ usually depends
on our intuition about the problem of interest. The function $f(x)$
produces a set of pairs $(x_i,y_i)$, $i=1\dots N$, which serve as the
observable data. Our goal is to select a function from the hypothesis
set $h\in\mathcal{H}$ which approximates $f(x)$ as best as possible,
namely, we would like to find $h\in\mathcal{H}$ such that $h\approx
f$ in some strict mathematical sense which we specify below. If this
is possible, we say that we \emph{learned} $f(x)$. But if the
function $f(x)$ can, in principle, take any value on
\emph{unobserved} inputs, how is it possible to learn in any
meaningful sense?
!split
===== Training and testing data =====
We will discuss the bias-variance tradeoff in the context of continuous predictions such as regression. However, many of the intuitions and ideas discussed here also carry over to classification tasks. Consider a dataset $\mathcal{L}$ consisting of the data $\mathbf{X}_\mathcal{L}=\{(y_j, \boldsymbol{x}_j), j=1\ldots N\}$. Let us assume that the true data is generated from a noisy model
!bt
\[
y=f(\boldsymbol{x}) + \epsilon
\]
!et
where $\epsilon$ is normally distributed with mean zero and standard deviation $\sigma_\epsilon$.
!split
===== Procedure to find a predictor =====
We have a statistical procedure (e.g. least-squares regression) for
forming a predictor $\hat{g}_{\mathcal{L}}(\boldsymbol{x})$ that gives the
prediction of our model for a new data point $\boldsymbol{x}$. This estimator
is chosen by minimizing a cost function which we take to be the
squared error
!bt
\[
\mathcal{C}( \boldsymbol{X}, \hat{g}(\boldsymbol{x})) = \sum_i (y_i - \hat{g}_\mathcal{L}(\boldsymbol{x}_i))^2.
\]
!et
!split
===== What we want =====
We are interested in the generalization error on all data drawn from
the true model, not just the error on the particular training dataset
$\mathcal{L}$ that we have in hand. This is just the expectation of
the cost function over many different data sets
$\{\mathcal{L}_j\}$. Denote this expectation value by
$E_{\mathcal{L}}$. In other words, we can view $\hat{g}_{\mathcal{L}}$
as a stochastic functional that depends on the dataset $\mathcal{L}$
and we can think of $E_{\mathcal{L}}$ as the expected value of the
functional if we drew an infinite number of datasets $\{\mathcal{L}_1,
\mathcal{L}_2, \ldots \}$.
!split
===== The expected generalization error =====
We would also like to average over different instances of the
``noise'' $\epsilon$ and we denote the expectation value over the
noise by $E_\epsilon$. Thus, we can decompose the expected
generalization error as
!bt
\begin{align}
E_\mathcal{L, \epsilon}[\mathcal{C}( \boldsymbol{X}, \hat{g}(\boldsymbol{x})) ]&= E_\mathcal{L,\epsilon}\left[ \sum_i ({y}_i - \hat{g}_\mathcal{L}(\boldsymbol{x}_i))^2 \right] \nonumber \\
&= E_\mathcal{L, \epsilon}\left[ \sum_{i}({y}_i -f(\boldsymbol{x}_i) +f(\boldsymbol{x}_i)- \hat{g}_\mathcal{L}(\boldsymbol{x}_i))^2\right] \nonumber \\
&= \sum_i E_\epsilon[ ({y}_i -f(\boldsymbol{x}_i))^2 ]+ E_\mathcal{L, \epsilon}[(f(\boldsymbol{x}_i)- \hat{g}_\mathcal{L}(\boldsymbol{x}_i))^2] + 2E_\epsilon[{y}_i -f(\boldsymbol{x}_i)]E_\mathcal{L}[f(\boldsymbol{x}_i)- \hat{g}_\mathcal{L}(\boldsymbol{x}_i)] \nonumber \\
&=\sum_i \sigma_\epsilon^2 + E_\mathcal{L}[(f(\boldsymbol{x}_i)- \hat{g}_\mathcal{L}(\boldsymbol{x}_i))^2],
\end{align}
!et
where in the last line we used the fact that our noise has zero mean
and variance $\sigma_\epsilon^2$ and the sum over $i$ applies to all
terms.
!split
===== Elaborating a little bit more =====
It is also helpful to further decompose the second term as
follows:
!bt
\begin{align}
E_\mathcal{L}[(f(\boldsymbol{x}_i)- \hat{g}_\mathcal{L}(\boldsymbol{x}_i))^2] &=E_\mathcal{L}[(f(\mathbf{x}_i)-E_\mathcal{L}[\hat{g}_\mathcal{L}(\boldsymbol{x}_i)]+ E_\mathcal{L}[\hat{g}_\mathcal{L}(\boldsymbol{x}_i)]- \hat{g}_\mathcal{L}(\boldsymbol{x}_i))^2] \nonumber \\
&=E_\mathcal{L}[(f(\boldsymbol{x}_i)-E_\mathcal{L}[\hat{g}_\mathcal{L}(\boldsymbol{x}_i)])^2] + E_\mathcal{L}[( \hat{g}_\mathcal{L}(\boldsymbol{x}_i)-E_\mathcal{L}[\hat{g}_\mathcal{L}(\boldsymbol{x}_i)])^2] \nonumber \\
&+2E_\mathcal{L}[(f(\boldsymbol{x}_i)-E_\mathcal{L}[\hat{g}_\mathcal{L}(\boldsymbol{x}_i)])( \hat{g}_\mathcal{L}(\boldsymbol{x}_i)-E_\mathcal{L}[\hat{g}_\mathcal{L}(\boldsymbol{x}_i)])] \nonumber \\
&=(f(\boldsymbol{x}_i)-E_\mathcal{L}[\hat{g}_\mathcal{L}(\boldsymbol{x}_i)])^2+E_\mathcal{L}[( \hat{g}_\mathcal{L}(\boldsymbol{x}_i)-E_\mathcal{L}[\hat{g}_\mathcal{L}(\boldsymbol{x}_i)])^2].
\end{align}
!et
!split
===== The bias =====
The first term is called the bias
!bt
\[
Bias^2= \sum_i (f(\boldsymbol{x}_i)-E_\mathcal{L}[\hat{g}_\mathcal{L}(\boldsymbol{x}_i)])^2
\]
!et
and measures the deviation of the expectation value of our estimator (i.e. the asymptotic value of our estimator in the infinite data limit) from the true value.
!split
===== The variance =====
The second term is called the variance
!bt
\[
Var=\sum_i E_\mathcal{L}[( \hat{g}_\mathcal{L}(\boldsymbol{x}_i)-E_\mathcal{L}[\hat{g}_\mathcal{L}(\boldsymbol{x}_i)])^2],
\]
!et
and measures how much our estimator fluctuates due to finite-sample effects. Combining these expressions, we see that the expected out-of-sample error of our model can be decomposed as
!bt
\[
E_\mathrm{out}=E_\mathcal{L, \epsilon}[\mathcal{C}( \boldsymbol{X}, \hat{g}(\boldsymbol{x})) ] = Bias^2 + Var + Noise.
\]
!et
The bias-variance tradeoff summarizes the fundamental tension in
machine learning, particularly supervised learning, between the
complexity of a model and the amount of training data needed to train
it. Since data is often limited, in practice it is often useful to
use a less-complex model with higher bias -- a model whose asymptotic
performance is worse than another model -- because it is easier to
train and less sensitive to sampling noise arising from having a
finite-sized training dataset (smaller variance).
!split
===== Summing up =====
The above equations tell us that in
order to minimize the expected test error, we need to select a
statistical learning method that simultaneously achieves low variance
and low bias. Note that variance is inherently a nonnegative quantity,
and squared bias is also nonnegative. Hence, we see that the expected
test MSE can never lie below $Var(\epsilon)$, the irreducible error.
What do we mean by the variance and bias of a statistical learning
method? The variance refers to the amount by which our model would change if we
estimated it using a different training data set. Since the training
data are used to fit the statistical learning method, different
training data sets will result in a different estimate. But ideally the
estimate for our model should not vary too much between training
sets. However, if a method has high variance then small changes in
the training data can result in large changes in the model. In general, more
flexible statistical methods have higher variance.
!split
===== Logistic Regression =====
So far we have focused on learning from datasets for which there is a
_continuous_ output. In linear regression we have been
concerned with learning the coefficients of a polynomial to predict
the response of a continuous variable $y_i$ on unseen data based on
its independent variables ${\bf x}_i$.
Classification problems,
however, are concerned with outcomes taking the form of discrete
variables (i.e. categories). For example, we may want to detect if
there's a cat or a dog in an image. Or given a specific system,
we'd like to identify its state, say whether it is an ordered or disordered system (typical situation in solid state physics).
(e.g. ordered/disordered).
_Logistic regression deals with binary, dichotomous outcomes (e.g. True or
False, Success or Failure, etc.). It is worth noting that logistic
regression is also commonly used in modern supervised Deep Learning
models_, as we will see later.
!split
===== Basics =====
We consider the case where the dependent variables $y_i\in\mathbb{Z}$
are discrete and only take values from $m=0,\dots,M-1$ (i.e. $M$
classes).
The goal is to predict the
output classes from the design matrix $X\in\mathbb{R}^{n\times p}$
made of $n$ samples, each of which bears $p$ features. Of cours e the
primary goal is to identify the classes to which new unseen samples
belong.
!split
===== Linear classifier =====
Let us start by considering a slightly simpler classifier: a linear classifier that categorizes examples using a weighted linear-combination of the features and an additive offset
!bt
\begin{equation}
s_i = \boldsymbol{x}_i^T\boldsymbol{w} + b_0 \equiv \mathbf{x}_i^T\mathbf{w},
\end{equation}
!et
where we use the short-hand notation
$\mathbf{x}_i = (1,\boldsymbol{x}_i)$ and $\mathbf{w}_i = (b_0,\boldsymbol{w}_i)$.
!split
===== Some selected properties =====
This function takes values on the entire real axis. In the case of logistic regression, however, the labels $y_i$ are discrete variables. One simple way to get a discrete output is to have sign functions that map the output of a linear regressor to $\{0,1\}$, $f(s_i)=$ sign$(s_i) = 1$ if $s_i\ge 0$ and 0 if otherwise. Indeed, this is commonly known as the ``perceptron" in the machine learning literature. This model is extremely simple, and it is favorable in many cases (e.g. noisy data) to have a ``soft" classifier that outputs the probability of a given category. For example, given $\mathbf{x}_i$, the classifier outputs the probability of being in category $m$. One such function is the logistic (or sigmoid) function:
!bt
\begin{equation}
f(s) = \frac{1}{1+\mathrm e^{-s}}.
\label{eq:log_fun}
\end{equation}
!et
Note that $1-f(s)= f(-s)$, which will be useful shortly.
!split
===== The cross-entropy as a cost function for logistic regression =====
The perceptron is an example of a ``hard classification'': 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 ``soft'' classifier that outputs the probability of a given category rather than a single value. For example, given $\mathbf{x}_i$, the classifier outputs the probability of being in category $m$.
Logistic regression is the most canonical example of a soft classifier. In logistic regression, the probability that a data point $\boldsymbol{x}_i$ belongs to a category $y_i=\{0,1\}$ is is given by
!bt
\begin{eqnarray}
P(y_i=1|\boldsymbol{x}_i,\boldsymbol{\theta)} &=& \frac{1}{1+\mathrm{e}^{-\mathbf{x}^T_i\mathbf{w}}},\nonumber\\
P(y_i=0|\boldsymbol{x}_i,\boldsymbol{\theta)} &=& 1 - P(y_i=1|\boldsymbol{x}_i,\boldsymbol{\theta)},
\end{eqnarray}
!et
where $\boldsymbol{\theta}=\mathbf{w}$ are the weights we wish to learn from the data.
Notice that in terms of the logistic function, we can write
!bt
\[
P(y_i=1) =f(\mathbf{x}_i^T\mathbf{w})=1-P(y_i=0).
\]
!et
!split
===== Maximum likelihood =====
We now define the cost function for logistic regression using Maximum
Likelihood Estimation (MLE). Recall, that in MLE we choose parameters
to maximize the probability of seeing the observed data. Consider a
dataset $\mathcal{D}=\{(y_i,\boldsymbol{x}_i)\}$ with binary labels
$y_i\in\{0,1\}$ where the data points are drawn independently. The
likelihood of the seeing the data under our model is just:
!bt
\begin{align}
P(\mathcal{D}|\mathbf{w})& = \prod_{i=1}^n \left[f(\mathbf{x}_i^T\mathbf{w})\right]^{y_i}\left[1-f(\mathbf{x}_i^T\mathbf{w})\right]^{1-y_i}\nonumber \\
\end{align}
!et
from which we can readily compute the log-likelihood:
!bt
\begin{equation}
l(\mathbf{w}) = \sum_{i=1}^n y_i\log f(\mathbf{x}_i^T\mathbf{w}) + (1-y_i)\log\left[1-f(\mathbf{x}_i^T\mathbf{w})\right].
\end{equation}
!et
!split
The maximum likelihood estimator is defined as the set of parameters that maximize the log-likelihood where we maximize with respect to $\theta$
!bt
\[
\hat{\mathbf{w}} = \sum_{i=1}^n y_i\log f(\mathbf{x}_i^T\mathbf{w}) + (1-y_i)\log\left[1-f(\mathbf{x}_i^T\mathbf{w})\right].
\]
!et
Since the cost (error) function is just the negative log-likelihood, for logistic regression we have that
!bt
\begin{eqnarray}
\mathcal{C}(\mathbf{w}) &=& - l(\mathbf{w}) \\
&=& \sum_{i=1}^n -y_i\log f(\mathbf{x}_i^T\mathbf{w}) - (1-y_i)\log\left[1-f(\mathbf{x}_i^T\mathbf{w})\right].\nonumber
\end{eqnarray}
!et
This equation is known in statistics as the \emph{cross entropy}. Finally, we note that just as in linear regression,
in practice we usually supplement the cross-entropy with additional regularization terms, usually $L_1$ and $L_2$ regularization as we did for Ridge and Lasso regression.
!split
===== Minimizing the cross entropy =====
The cross entropy is a convex function of the weights $\mathbf{w}$ and,
therefore, any local minimizer is a global minimizer. Minimizing this
cost function leads to the following equation
!bt
\begin{equation}
\boldsymbol{0}=\boldsymbol{\nabla} \mathcal{C}(\mathbf{w}) = \sum_{i=1}^n\left[f(\mathbf{x}_i^T\mathbf{w})-y_i\right]\mathbf{x}_i,
\end{equation}
!et
where we made use of the logistic function identity $\partial_z f(z) =
f(z)[1-f(z)]$. This equation defines a transcendental equation for
$\mathbf{w}$, the solution of which, unlike linear regression, cannot
be written in a closed form.
Here we need gradient descent methods!