Chapter 3: Logistic Regression

Neurons can be categorized mainly into two types based on their output. If the output belongs to a discrete set, then the neuron is of a classification type. On the other hand, if the output is any real number in an interval, then the neuron is said to be of regression type. There is also a third type of neuron, whose output is a real number in the interval \((0, 1)\), which is interpreted as a class probability and used for probabilistic classification.

The perceptron, studied in the previous chapter, involves a hard classification rule, which is suitable for a linearly separable dataset. However, such a hard decision boundary does not measure confidence in the prediction when the dataset structure is unknown or the data involves random noise. In such cases, it is desirable to have a model that outputs a probability that the input belongs to each class, allowing both a probabilistic interpretation and more flexible decision rules at inference time.

In this chapter, we introduce a statistical model, called logistic regression, for probabilistic binary classification. Assuming that the log-likelihood ratio between two classes is affine in the feature vector, the model obtains the posterior probability of either class in the form of a logistic sigmoid composed with an affine function of the input. Thus, the model has a single neuron structure.

In Subsection 3.1.1, we derive logistic regression from first principles using Bayes' theorem. In Subsection 3.1.2, we observe that logistic regression is precisely a single artificial neuron with the sigmoid activation. The natural cost function for training a logistic regression model, derived from maximum likelihood estimation, is the cross-entropy cost function, discussed in Subsection 3.2.1. We also present a training algorithm based on the gradient descent method in Subsection 3.2.2 and end the chapter with a few real-world applications discussed in Section 3.3.

3.1 The Logistic Regression Model

Consider a dataset

\[ \mathcal{D} = \big(\mathcal{C}^-\times \{0\}\big) \cup \big(\mathcal{C}^+\times \{1\}\big).\]

If \(\mathcal{D}\) is linearly separable, then a perceptron can learn a separating hyperplane with weight \(\boldsymbol{w}\) and bias \(b\), and classify inputs by the rule

\[ \hat{y} = \begin{cases} 1, & \langle \boldsymbol{w}, \boldsymbol{x}\rangle + b > 0, \\ 0, & \text{otherwise}. \end{cases} \]

This is a hard classification rule that works only when the dataset is linearly separable. However, there are many instances where the data is noisy and consequently, the dataset is not linearly separable. In such cases, instead of making a hard decision, we may ask the following question:

What is the probability that an input \(\boldsymbol{x}\) belongs to \(\mathcal{C}^+\)?

This leads to a probabilistic model, which is preferable to a hard classification rule. We begin with a motivating example.

Example 3.1 (Hours studied versus exam outcome):

A class of \(N_\text{train} = 100\) students each study for a certain number of hours and subsequently sit an exam. The input is \(x_i \in [0, 10]\), the number of hours studied by student \(i\), and the label is \(y_i \in \{0,1\}\), where \(y_i = 1\) denotes a pass and \(y_i = 0\) a fail. The dataset \(\mathcal{D}\) is shown in Figure 3.1.

Since the data is one-dimensional, the separating hyperplane is just a point in this case. From the illustrating figure, we can see that the point cannot be placed anywhere in the interval \([0,10]\) so that the failed students are separated from the students who passed the exam. Therefore, a perceptron classification will induce a misclassification in this problem. An alternative is to seek a probabilistic approach, where the probability of passing as a function of hours studied is illustrated as a solid curve in Figure 3.1.

Figure 3.1: The hours-studied dataset of Example 3.1 is visually shown where plus (\(+\)) denote failures (\(y_i=0\)) and filled circles denote passes (\(y_i=1\)). The probability of passing as a function of hours studied, \(\hat{p}(x)\), is shown as a solid black curve.

3.1.1 Derivation using Bayes' Formula

In this subsection, we derive a probabilistic model, called the logistic regression model systematically using Bayes' theorem.

Using Bayes' theorem, the posterior probability for the class \(\mathcal{C}^+\) can be written as

\[ \mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x}) = \frac{p(\boldsymbol{x} \mid \mathcal{C}^+) \mathbb{P}(\mathcal{C}^+)} {p(\boldsymbol{x} \mid \mathcal{C}^+) \mathbb{P}(\mathcal{C}^+) + p(\boldsymbol{x} \mid \mathcal{C}^-) \mathbb{P}(\mathcal{C}^-)}, \]

where \(\mathbb{P}\) denotes a probability measure, and \(p\) is the probability density (or mass) function associated with \(\mathbb{P}\).

Note:

The notation \(\mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x})\) should be read as follows:

"Given that we have observed the input vector \(\boldsymbol{x}\), the probability that it belongs to the class \(\mathcal{C}^+\)''

This notation is standard in machine learning literature and is a shorthand notation for the conditional probability \(\mathbb{P}(Y = 1 \mid \boldsymbol{X} = \boldsymbol{x})\), where \(Y \in \{0,1\}\) is the class-label random variable and \(\boldsymbol{X}\) is the input random vector.

Taking the ratio of posteriors, we get

\[ \frac{\mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x})}{\mathbb{P}(\mathcal{C}^- \mid \boldsymbol{x})} = \frac{p(\boldsymbol{x} \mid \mathcal{C}^+) \mathbb{P}(\mathcal{C}^+)}{p(\boldsymbol{x} \mid \mathcal{C}^-) \mathbb{P}(\mathcal{C}^-)}, \]

which is referred to as the odds. The logarithm of the odds is called the logit or log-odds, and from the above expression, we see that the log-odds can be written as

\[ \log \frac{\mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x})}{\mathbb{P}(\mathcal{C}^- \mid \boldsymbol{x})} = \log \Lambda(\boldsymbol{x}) + \log\left( \frac{\mathbb{P}(\mathcal{C}^+)}{\mathbb{P}(\mathcal{C}^-)}\right), \]

where

\[ \Lambda(\boldsymbol{x}) = \frac{p(\boldsymbol{x} \mid \mathcal{C}^+)}{p(\boldsymbol{x} \mid \mathcal{C}^-)} \]

is the likelihood ratio.

Assume that the log-likelihood ratio is linear in \(\boldsymbol{x}\), that is, there is an augmented weight vector \(\overline{\boldsymbol{w}}\) such that

\begin{eqnarray} \log \Lambda(\boldsymbol{x}) = \langle \overline{\boldsymbol{w}}, \overline{\boldsymbol{x}} \rangle, \end{eqnarray}(3.1)

for all input vectors \(\boldsymbol{x}\). Then the log-odds can be written as

\[ \log \frac{\mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x})}{\mathbb{P}(\mathcal{C}^- \mid \boldsymbol{x})} = \langle \boldsymbol{w}, \boldsymbol{x} \rangle + b, \]

where \(b = w_0 + \log\left( \tfrac{\mathbb{P}(\mathcal{C}^+)}{\mathbb{P}(\mathcal{C}^-)}\right)\) is treated as the bias.

Noting that \(\mathbb{P}(\mathcal{C}^- \mid \boldsymbol{x}) = 1 - \mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x})\) and solving the above equation for \(\mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x})\), we obtain

\begin{eqnarray} \mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x}) = \mathscr{S}_1\big(\langle \boldsymbol{w}, \boldsymbol{x} \rangle + b\big), \end{eqnarray}(3.2)

where \(\mathscr{S}_1\) is the logistic function defined in (1.5).

From the above expression, we see that estimating a weight vector \(\boldsymbol{w}\) and a bias \(b\) allows us to assign probabilities for each of the two classes \(\mathcal{C}^+\) and \(\mathcal{C}^-\). At the inference stage, by applying a suitable decision rule (for instance, thresholding at \(0.5\)), we obtain a binary (hard) classification. In the context of statistical learning, such a probabilistic classification model is called logistic regression.

Remark 3.1:

The name logistic regression is explained as follows:

The word logistic refers to the logistic activation function used to model the posterior probability.
The word regression reflects the assumption that the log-likelihood ratio is affine in \(\boldsymbol{x}\). This terminology is in contrast with the statistical convention in which "regression'' refers to function approximation, however, here the model is used for classification.

The following problem establishes that the inverse of the logistic function is precisely the logit function, which confirms the consistency of the terminology.

Problem 3.1:

Show that \(\mathscr{S}_1^{-1}\) is the logit or log-odds function.

Logistic regression assumes that the log-likelihood ratio between the two classes is an affine function of the input as given in (3.1). This assumption holds exactly for certain class-conditional distributions. One such distribution is given in the following problem.

Problem 3.2 (Gaussian Class-Conditionals Justify Linearity):

Consider a binary classification problem with classes \(\mathcal{C}^+\) and \(\mathcal{C}^-\). Let \(\boldsymbol{X} \in \mathbb{R}^{n_0}\) be a feature vector whose class-conditional distributions are multivariate Gaussian with class means \(\boldsymbol{\mu}^+, \boldsymbol{\mu}^- \in \mathbb{R}^{n_0}\) and a common positive-definite covariance matrix \(\Sigma \in \mathbb{R}^{n_0 \times n_0}\), that is,

\[ \boldsymbol{X} \mid \mathcal{C}^+ \sim \mathcal{N}(\boldsymbol{\mu}^+,\,\Sigma),\qquad \boldsymbol{X} \mid \mathcal{C}^- \sim \mathcal{N}(\boldsymbol{\mu}^-,\,\Sigma). \]

Find the augmented weight vector \(\overline{\boldsymbol{w}}\in \mathbb{R}\times \mathbb{R}^{n_0}\) in terms of \(\boldsymbol{\mu}^+\), \(\boldsymbol{\mu}^-\), and \(\Sigma\) such that the log-likelihood ratio is given by (3.1),

\[ \log\Lambda(\boldsymbol{x}) \;=\; \langle\boldsymbol{w},\boldsymbol{x}\rangle + w_0. \]

Explain what happens when the two Gaussian class-conditional distributions have different covariance matrices.

Remark 3.2:

From the above problem, we see that, in the Gaussian equal-covariance setting, logistic regression is exact (not an approximation), and its weight vector \(\boldsymbol{w} = \Sigma^{-1}(\boldsymbol{\mu}^+-\boldsymbol{\mu}^-)\) is the same direction used in Fisher's Linear Discriminant Analysis.

3.1.2 Single Neuron Interpretation

The logistic regression model derived in the previous subsection has a natural interpretation in the language of artificial neurons with the neuron function

\[ \text{𝕗}^{(1)}(\boldsymbol{x};\overline{\boldsymbol{w}}^{(1)}) = \mathscr{S}_1 \big( \text{𝕒}(\overline{\boldsymbol{x}};\overline{\boldsymbol{w}}^{(1)}) \big), \]

for any given input vector \(\boldsymbol{x}\in \mathbb{R}^{n_0}\), where

\[\text{𝕒}(\overline{\boldsymbol{x}};\overline{\boldsymbol{w}}^{(1)}) = \langle \boldsymbol{w}, \boldsymbol{x}\rangle + b\]

is the affine pre-activation. Thus, logistic regression can be trained as an artificial neuron to approximate the posterior probability of the positive class of a given dataset, so that the output is

\[ \hat{p} = \mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x}). \]

Hard classification is then obtained during inference by applying a decision rule

\[ \hat{y} = \begin{cases} 1, & \text{if} ~ \hat{p} \ge \hat{p}^*, \\ 0, & \text{otherwise}, \end{cases} \]

where \(\hat{p}^*\in (0,1)\) is chosen based on the application. For instance, \(\hat{p}^*=0.5\) is a default choice.

Note:

Logistic regression finds the best linear boundary in a probabilistic sense that provides a measure of confidence. Hard classification is obtained only at the inference stage by applying a decision rule (e.g., thresholding). Because of this probabilistic formulation, logistic regression can also be applied to datasets that are not linearly separable.

3.2 Training Logistic Regression

Recall from Subsection 1.4.2 that to learn parameters \((b, \boldsymbol{w})\) of a model, we have to choose a suitable cost function. There is a natural choice of the cost function to train logistic regression, which we discuss in the following subsection.

3.2.1 Cross-Entropy Cost Function

The natural training principle for logistic regression is maximum likelihood estimation (MLE), where the parameters \((b,\boldsymbol{w})\) are chosen to maximize the likelihood of the observed labels given the inputs.

For a single training example \((\boldsymbol{x}_i,y_i)\), the model predicts a probability \(\hat{p}_i = \mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x}_i)\). The likelihood of observing the true label \(y_i \in \{0,1\}\) under this prediction is

\[ L_i = \hat{p}_i^{\,y_i}\,(1-\hat{p}_i)^{\,1-y_i}. \]

For an entire training dataset \(\mathcal{D}_\text{train} = \big\{(\boldsymbol{x}_i,y_i) \mid i=1,2,\ldots, N_\text{train}\big\}\) having independent samples (examples) from the underlying joint distribution, the joint likelihood is the product

\[ L(b, \boldsymbol{w}) = \prod_{i=1}^{N_\text{train}} \hat{p}_i^{y_i} (1 - \hat{p}_i)^{1-y_i}. \]

Rather than maximizing this product directly, it is more convenient to work with the log-likelihood

\[ \log L(b, \boldsymbol{w}) = \sum_{i=1}^{N_\text{train}} \Big[ y_i \log \hat{p}_i + (1-y_i)\log(1-\hat{p}_i) \Big]. \]

Thus, the MLE for logistic regression is to obtain the parameters \((b,\boldsymbol{w})\) that maximise the log-likelihood. This is equivalent to minimizing the negative log-likelihood, normalised by \(N_\text{train}\), given by

\begin{eqnarray} \mathcal{L}(b, \boldsymbol{w}) = - \frac{1}{N_\text{train}}\sum_{i=1}^{N_\text{train}} \Big[ y_i \log \hat{p}_i + (1-y_i)\log(1-\hat{p}_i) \Big]. \end{eqnarray}(3.3)

The above expression is the cross-entropy cost function for binary labels \(y_i\in\{0,1\}\) under the Bernoulli model. The following definition generalizes this concept to arbitrary probability distributions.

Definition 3.1 (Cross-Entropy and Entropy):

Let \(p\) and \(q\) be two densities (or mass functions) on the same space \(\mathcal{X}\) and let \(\mathbb{E}_p\) be the expectation with respect to \(p\). The cross-entropy of \(p\) with respect to \(q\) is denoted by \(\mathfrak{C}(p,q)\) and is defined as

\[ \mathfrak{C}(p,q) = \mathbb{E}_p(-\log q). \]

In particular, if \(p=q\), we call \(\mathfrak{C}(p,p) =: \mathfrak{S}(p)\), the entropy of \(p\).

Entropy is a concept in information theory which measures the uncertainty of a random variable \(X\). It quantifies the average amount of information we gain while observing the outcome of \(X\). When \(X\) is more predictable, the entropy is low and it takes the value zero if \(X\) is deterministic. On the other hand, if \(X\) is uniformly distributed over many possible outcomes, then entropy is high.

Example 3.2 (Entropy of the Bernoulli distribution):

Let \(Y \sim \mathrm{Bernoulli}(r)\), \(r \in (0,1)\), so that \(\mathbb{P}(Y=1)=r\) and \(\mathbb{P}(Y=0)=1-r\). Applying the discrete case of Definition 3.1 with \(\mathcal{X}=\{0,1\}\) and \(p=q\),

\[ \mathfrak{S}(r) = -r\log r - (1-r)\log(1-r). \]

Three values are instructive:

\(r=0\) or \(r=1\): the outcome is deterministic, so \(\mathfrak{S}=0\) (using the convention \(0\log 0 = 0\)).
\(r=\tfrac{1}{2}\): both outcomes are equally likely, giving the maximum \(\mathfrak{S}(\tfrac{1}{2}) = \ln 2 \approx 0.693\) nats.
For all other \(r \in (0,1)\): the entropy is strictly between \(0\) and \(\ln 2\).

In the logistic regression context, \(\hat{p}_i = \mathscr{S}_1(\langle\boldsymbol{w},\boldsymbol{x}_i\rangle - b)\) is the model's predicted probability of the positive class. The entropy \(\mathfrak{S}(\hat{p}_i)\) measures the model's uncertainty about prediction \(i\): it is zero when the model is fully confident (\(\hat{p}_i\) near \(0\) or \(1\)) and maximal when the model is completely uncertain (\(\hat{p}_i = \tfrac{1}{2}\)). The cross-entropy cost function (3.3) goes further: it penalises not just uncertainty but assigning low probability to the correct label.

Example 3.3 (Entropy of fair and biased dice):

Consider two six-faced dice with outcomes \(\{1,2,3,4,5,6\}\), where Die-1 is fair and Die-2 is biased:

Die-1: \(\mathbb{P}(X=k) = \frac{1}{6},~k=1,2,\ldots,6.\)
Die-2: \(\mathbb{P}(X=1) = 0.9,~\mathbb{P}(X=k) = 0.02,~k=2,\ldots,6.\)

Since any outcome in Die-1 is equally likely, there is high randomness, and we expect a larger entropy value. In Die-2, face-1 dominates with very high probability, so there is less uncertainty. Let us verify this by computing the entropy of the distribution of Die-k, denoted by \(p_i\), for \(i=1,2\) to get

\[ \mathfrak{S}(p_1) = -\sum_{k=1}^6 \frac{1}{6}\ln\!\left(\frac{1}{6}\right) = \ln 6 \approx 1.792 \text{ nats}, \]

\[ \mathfrak{S}(p_2) = -\!\left(0.9\ln(0.9) + \sum_{k=2}^6 0.02\ln(0.02)\right) \approx 0.486 \text{ nats.} \]

From the above calculations, we see that Die-1 has higher entropy, which confirms a greater uncertainty.

Problem 3.3:

Let \(p\) and \(q\) be two probability densities on \(\mathbb{R}\) such that the integral

\[\mathfrak{C}(p,q) = -\int\limits_{\mathbb{R}} p(x)\log q(x)\,dx\]

is well-defined. Show that \(\mathfrak{C}(p,p) \leq \mathfrak{C}(p,q)\).

Observe that the entropy \(\mathfrak{S}(p)\) measures the inherent uncertainty of the true distribution \(p\), whereas the cross-entropy \(\mathfrak{C}(p,q)\) measures the average penalty incurred when using \(q\) to model outcomes drawn from \(p\). From the above problem, we see that minimizing cross-entropy is equivalent to making \(q\) as close as possible to the true distribution \(p\).

In logistic regression, \(p_i\) is the point mass at the true label \(y_i \in \{0,1\}\) and \(q_i = \mathrm{Bernoulli}(\hat{p}_i)\) is the model's predicted distribution for sample \(i\). Computing the discrete cross-entropy over \(\mathcal{Y}=\{0,1\}\) and averaging over the training set recovers exactly (3.3) as

\[ \mathcal{L}(b,\boldsymbol{w}) = \frac{1}{N_\text{train}}\sum_{i=1}^{N_\text{train}}\mathfrak{C}(p_i,q_i). \]

There is another important concept in information theory, closely related to entropy and cross-entropy, which we define below.

Definition 3.2 (Kullback--Leibler Divergence):

Assume that \(p(x)\) is the true distribution and \(q(x)\) is a model approximating it. The Kullback--Leibler (KL) divergence of \(q\) from \(p\) is defined as

\[ \mathcal{C}_{\texttt{KL}}(p\|q) = \mathfrak{C}(p,q) - \mathfrak{S}(p), \]

which is the extra cost incurred by encoding data with \(q\) instead of \(p\).

Note:

If \(\mathcal{X}\) is continuous, then

\[ \mathcal{C}_{\texttt{KL}}(p\|q) = -\int\limits_{\mathbb{R}} p(x)\ln\frac{q(x)}{p(x)}\,dx. \]

From Definition 3.1, we see that \(\mathcal{C}_{\texttt{KL}}(p\|q) = 0\) if and only if \(q = p\). Since \(\mathfrak{S}(p)\) does not involve \(b\) and \(\boldsymbol{w}\), minimizing the cross-entropy cost is the same as minimizing the KL divergence between the true and the predicted distribution.

As mentioned in the definition, KL divergence measures the discrepancy between two probability distributions. It can be used as a cost function in neural networks whenever the network output represents a probability distribution, because minimizing KL encourages the predicted distribution to match the target distribution.

Problem 3.4:

In binary classification (logistic regression), let \(p_i\) be the point mass at \(y_i\in\{0,1\}\) and \(q_i=\mathrm{Bernoulli}(\hat{p}_i)\). Show that minimizing \(\dfrac{1}{N_\text{train}}\displaystyle\sum_{i=1}^{N_\text{train}}\mathcal{C}_{\texttt{KL}}(p_i\|q_i)\) over \((b,\boldsymbol{w})\) is equivalent to minimizing the cross-entropy cost function (3.3).

Problem 3.5:

\(~\)

Show that \(\mathcal{C}_{\texttt{KL}}(p\|q) \geq 0\).
Is \(\mathcal{C}_{\texttt{KL}}(p\|q)\) symmetric? Justify your answer.
Does the following triangle inequality hold? Justify your answer.
\[ \mathcal{C}_{\texttt{KL}}(p\|q) \leq \mathcal{C}_{\texttt{KL}}(p\|r) + \mathcal{C}_{\texttt{KL}}(r\|q). \]

3.2.2 Gradient Descent Algorithm

The cost function \(\mathcal{L}(b, \boldsymbol{w})\) given by (3.3) is smooth and convex, but it has no closed-form minimizer. Hence, we must rely on a computational technique to approximate

\[ (b^*, \boldsymbol{w}^*) = \text{argmin}_{(b, \boldsymbol{w})} \mathcal{L}(b, \boldsymbol{w}), \]

if it exists. A commonly used optimization method is gradient descent, which we present here in its simplest form and postpone a detailed discussion to a later chapter.

For a given initial guess \((b^{(0)}, \boldsymbol{w}^{(0)})\), the iterative formula for the gradient descent method is given by

\begin{eqnarray} b^{(k+1)} &=& b^{(k)} - \eta \frac{\partial \mathcal{L}}{\partial b}(b^{(k)}, \boldsymbol{w}^{(k)}),\\ \boldsymbol{w}^{(k+1)} &=& \boldsymbol{w}^{(k)} - \eta \nabla_{\boldsymbol{w}} \mathcal{L}(b^{(k)}, \boldsymbol{w}^{(k)}), ~~k=0,1,\ldots. \end{eqnarray}(3.4)

Since \(\nabla_{\boldsymbol{w}} s_i = \boldsymbol{x}_i\) and \(\partial s_i/\partial b = 1\), the chain rule gives

\begin{eqnarray} \nabla_{\boldsymbol{w}}\,\mathcal{L}(b, \boldsymbol{w}) = \frac{1}{N_\text{train}}\sum_{i=1}^{N_\text{train}} \big(\hat{p}_i - y_i\big)\boldsymbol{x}_i, \end{eqnarray}(3.5)

\begin{eqnarray} \frac{\partial \mathcal{L}}{\partial b}(b, \boldsymbol{w}) = \frac{1}{N_\text{train}}\sum_{i=1}^{N_\text{train}} \big(\hat{p}_i - y_i\big). \end{eqnarray}(3.6)

Starting from an initial guess \((b^{(0)}, \boldsymbol{w}^{(0)})\), the parameters are updated iteratively as described in the following algorithm.

Algorithm 3.1 (Gradient Descent for Logistic Regression):

Input:

Training dataset \(\mathcal{D}_\text{train} = \big\{(\boldsymbol{x}_i, y_i) \mid i = 1,2,\ldots, N_\text{train}\big\}\)
Initial parameters \((b^{(0)}, \boldsymbol{w}^{(0)})\),
learning rate \(\eta > 0\), and
tolerance \(\varepsilon > 0\)

Processing:
For each \(k = 0, 1, 2, \ldots\), perform the following steps:

Step 1: Compute:
\(\hat{p}_i^{(k)} = \mathscr{S}_1\!\big(\langle \boldsymbol{w}^{(k)}, \boldsymbol{x}_i \rangle + b^{(k)}\big)\), for each \(i = 1, 2, \ldots, N_\text{train}\).
Also, the gradients \(\nabla_{\boldsymbol{w}}\mathcal{L}\) and \(\frac{\partial\mathcal{L}}{\partial b}\).
Step 2: Check:
If the condition
\[\displaystyle\sqrt{\|\nabla_{\boldsymbol{w}}\mathcal{L}\|^2 + \left(\frac{\partial\mathcal{L}}{\partial b}\right)^2} < \varepsilon\]
is satisfied, then go to the Output step below.
Else, go to Step 3.
Step 3: Update:
\[ \boldsymbol{w}^{(k+1)} = \boldsymbol{w}^{(k)} - \frac{\eta}{N_\text{train}}\sum_{i=1}^{N_\text{train}} \big(\hat{p}_i^{(k)} - y_i\big)\boldsymbol{x}_i, \]

\[ b^{(k+1)} = b^{(k)} - \frac{\eta}{N_\text{train}}\sum_{i=1}^{N_\text{train}} \big(\hat{p}_i^{(k)} - y_i\big). \]
Increment \(k\) and go to Step 1.

Output: Learned parameters \((\hat{b}^*,\, \hat{\boldsymbol{w}}^*) = (b^{(k)}, \boldsymbol{w}^{(k)})\).

Note:

Since \(\mathcal{L}(b, \boldsymbol{w})\) is smooth and convex, gradient descent converges to the global minimizer \((b^*, \boldsymbol{w}^*)\) for any sufficiently small learning rate \(\eta > 0\).

Code 3.1:

Develop a Python code to train a logistic regression model for the problem illustrated in Example 3.1 by going through the following steps:

Generate a synthetic dataset \(\mathcal{D}=\{(x_i, y_i)~\mid~i=1,2,\ldots,N\}\) by considering two classes
\[x \mid \mathcal{C}^- \sim \mathcal{N}(2.5,\, 1.5),\quad x \mid \mathcal{C}^+ \sim \mathcal{N}(7.5,\, 1.5),\]
where \(\mathcal{N}(\mu, \sigma^2)\) denotes the Gaussian distribution with mean \(\mu\) and variance \(\sigma^2\). Note that numpy.random.normal(mu, sigma) in Python takes the standard deviation \(\sigma\) as its second argument; use \(\sigma = \sqrt{1.5}\). Take \(N = \#(\mathcal{C}^-) + \#(\mathcal{C}^+) = 50 + 50.\)
Choose 70 examples uniformly at random and store them as \(\mathcal{D}_{\rm train}\). Store the remaining as \(\mathcal{D}_{\rm test}.\)
Train a logistic regression on the training dataset \(\mathcal{D}_{\rm train}\) using Algorithm 3.1.
Evaluate the trained model on \(\mathcal{D}_{\rm test}\) by computing
- the misclassification rate:
  \[\text{Error} = \frac{1}{N_{\rm test}}\sum_{i \in \mathcal{D}_{\rm test}} \mathbf{1}[\hat{y}_i \neq y_i],\quad \hat{y}_i = \begin{cases} 1 & \text{if}~\hat{p}_i \geq 0.5,\\ 0 & \text{if}~\hat{p}_i < 0.5, \end{cases}\]
- the test cross-entropy loss:
  \[\mathcal{L}_{\rm test} = -\frac{1}{N_{\rm test}}\sum_{i \in \mathcal{D}_{\rm test}} \big[y_i \log \hat{p}_i + (1-y_i)\log(1-\hat{p}_i)\big].\]
Compare \(\mathcal{L}_{\rm test}\) with the training loss \(\mathcal{L}_{\rm train}\) at convergence.

3.3 Applications of Logistic Regression

Logistic regression is well-suited for binary classification problems in which the log-odds of the outcome are approximately linear in the input features, and the output probability has a direct, interpretable meaning in the application domain. We illustrate this with five canonical examples drawn from finance, medicine, natural language processing, business analytics, and fraud detection.

3.3.1 Finance: Credit Default Prediction

A bank maintains a labelled dataset of past borrowers,

\[ \mathcal{D} = \big\{(\boldsymbol{x}_i, y_i) \mid i = 1,2,\ldots, N\big\} \subset \mathbb{R}^n \times \{0,1\}, \]

where the feature vector \(\boldsymbol{x}_i = (x_{i1}, x_{i2}, \ldots, x_{in})^\top \in \mathbb{R}^n\) records financial attributes of borrower \(i\) such as annual income, debt-to-income ratio, credit score, requested loan amount, and employment tenure, and the label is

\[ y_i = \begin{cases} 1, & \text{borrower} i \text{ defaulted on the loan,}\\ 0, & \text{otherwise.}\end{cases} \]

Given a new applicant with profile \(\boldsymbol{x}\), the aim is to estimate the probability of default

\[ \hat{p}(\boldsymbol{x}) = \mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x}), \quad \text{where}\ \mathcal{C}^+ = \{\text{default}\}. \]

3.3.1.1 Model and training.

From (3.2), the logistic regression model suggests

\[ \mathbb{P}(\mathcal{C}^+ \mid \boldsymbol{x}) = \mathscr{S}_1\!\big(\langle \boldsymbol{w}, \boldsymbol{x} \rangle + b\big), \]

under the assumption that the log-likelihood ratio \(\log \Lambda(\boldsymbol{x})\) is linear in \(\boldsymbol{x}\). The parameters \((b, \boldsymbol{w})\) are estimated by minimizing the cross-entropy cost function (3.3).

We can use Algorithm 1 to obtain the estimate of \((b, \boldsymbol{w})\), which is the augmented weight vector \((b^*, \boldsymbol{w}^*)\) to which the algorithm converges.

Remark 3.3 (Interpretation):

From (3.1), we see that the \(j^{\rm th}\) component \(w^{*}_{j}\) of the fitted weight \(\boldsymbol{w}^{*}\) represents the change in log-odds of default per unit increase in the corresponding feature \(x_j\), with all other features held fixed. Equivalently, \(e^{w_j}\) is the multiplicative factor by which the odds of default change per unit increase in \(x_j\). This is analogous to an odds-ratio interpretation used in classical credit scorecards.

Code 3.2 (Credit Default Prediction from Real Data):

(a) Downloading the Dataset: Download the Default of Credit Card Clients dataset from the

UCI Credit Card Default Dataset

The dataset includes \(N = 30000\) credit card holders of a Taiwanese bank, each described by \(n = 23\) numerical attributes \(\boldsymbol{x}_i\) (credit limit, age, six monthly repayment statuses, six bill amounts, six payment amounts, and three coded demographic variables), with \(y_i = 1\) if the holder defaulted in the following month and \(y_i = 0\) otherwise.

Compute the empirical default rate

\[\hat{\pi} = \frac{1}{N}\sum_{i=1}^{N} y_i,\]

which estimates the prior \(\mathbb{P}(\mathcal{C}^+)\).

(b) Dataset Preparation:

Split the index set \(\{1, 2, \ldots, N\}\) in the ratio \(70\!:\!30\) by drawing the training indices uniformly at random without replacement, initialising the random number generator with a fixed seed. State the seed you used. Denote the resulting training and test sets by \(\mathcal{D}_{\text{train}}\) and \(\mathcal{D}_{\text{test}}\), of sizes \(N_{\text{train}}\) and \(N_{\text{test}}\).

Next standardise each feature,

\[\tilde{x}_{ij} = \frac{x_{ij}-\mu_j}{\sigma_j},\]

where \(\mu_j\) and \(\sigma_j\) are the mean and the standard deviation of the \(j\)-th feature computed over \(\mathcal{D}_{\text{train}}\) alone. Apply this same transformation to the test points.

Use the standardised features in training the model.

Also compute the default rate within \(\mathcal{D}_{\text{train}}\),

\[ \hat{\pi}_{\text{train}} = \frac{1}{N_{\text{train}}} \sum_{i \in \mathcal{D}_{\text{train}}} y_i , \]

and the corresponding rate \(\hat{\pi}_{\text{test}}\) within \(\mathcal{D}_{\text{test}}\). Verify that the three rates \(\hat{\pi}\), \(\hat{\pi}_{\text{train}}\) and \(\hat{\pi}_{\text{test}}\) agree to within sampling error.

(c) Training Algorithm: Train a logistic regression model on the training dataset \(\mathcal{D}_{\rm train}\) using Algorithm 3.1, with \((b^{(0)},\boldsymbol{w}^{(0)}) = (0,\boldsymbol{0})\), \(\eta = 0.1\), \(\varepsilon = 10^{-5}\), and a cap of \(10^4\) iterations.

Verify the identity

\[ \frac{1}{N_{\text{train}}}\sum_{i \in \mathcal{D}_{\text{train}}} \hat{p}_i = \hat{\pi}_{\text{train}}, \]

which follows from \(\partial \mathcal{L}/\partial b = 0\) at the minimum and holds up to the tolerance \(\varepsilon\).

(d) Observing Convergence: Plot \(\mathcal{L}(b^{(k)},\boldsymbol{w}^{(k)})\) against \(k\). Repeat with \(\eta = 0.001\) and \(\eta = 10\) and plot all three curves together. Describe in one or two sentences what each learning rate does.

Compare your observations with the sufficient condition for descent,

\[ L \le \frac{\|\tilde{X}\|_2^2}{4N_{\text{train}}} \]

bounds the Lipschitz constant of \(\nabla \mathcal{L}\); note that this condition is sufficient but not necessary.

(e) Testing the Trained Model: For every \((\tilde{\boldsymbol{x}}_i, y_i)\in \mathcal{D}_{\text{test}}\), compute \(\hat{p}(\tilde{\boldsymbol{x}}_i)\) and predict default \(\hat{y}_i\) by defining the hard classification rule as

\[ \hat{y}_i = \begin{cases} 1, & \text{if} ~ \hat{p}(\tilde{\boldsymbol{x}}_i) \ge 0.5, \\ 0, & \text{otherwise}, \end{cases} \]

Let us define the terminology of classifier performance:

A test point is a true positive if \(\hat{y}_i = y_i = 1\), denoted by \(N_{\mathrm{TP}}\);
A test point a false positive if \(\hat{y}_i = 1\) but \(y_i = 0\), denoted by \(N_{\mathrm{FP}}\);
A test point a true negative if \(\hat{y}_i = y_i = 0\), denoted by \(N_{\mathrm{TN}}\); and
A test point a false negative if \(\hat{y}_i = 0\) but \(y_i = 1\), denoted by \(N_{\mathrm{FN}}\).

We have

\[ N_{\text{test}} = N_{\mathrm{TP}} + N_{\mathrm{FP}} + N_{\mathrm{TN}} + N_{\mathrm{FN}}, \]

False positives are also called type-1 errors and false negatives type-2 errors.

The confusion matrix is defined as

\[ \begin{array}{c|cc} & \hat{y} = 0 & \hat{y} = 1 \\ \hline y = 0 & N_{\mathrm{TN}} & N_{\mathrm{FP}} \\ y = 1 & N_{\mathrm{FN}} & N_{\mathrm{TP}} \end{array} \]

Compute the confusion matrix, together with
\[ \text{Accuracy} = \frac{N_{\mathrm{TP}} + N_{\mathrm{TN}}}{N_{\text{test}}}, \qquad \text{Precision} = \frac{N_{\mathrm{TP}}}{N_{\mathrm{TP}} + N_{\mathrm{FP}}}, \]

\[ \text{Recall} = \frac{N_{\mathrm{TP}}}{N_{\mathrm{TP}} + N_{\mathrm{FN}}}, \text{False positive rate} = \frac{N_{\mathrm{FP}}}{N_{\mathrm{FP}} + N_{\mathrm{TN}}}, \]
Compute the test error, which is the empirical risk under the \(0\)--\(1\) loss \(\ell_{01}(\hat{y}, y) = \mathbb{1}\{\hat{y} \neq y\}\), defined as
\[ \mathcal{E}_{\text{test}} = \frac{1}{N_{\text{test}}}\sum_{i=1}^{N_{\text{test}}} \mathbb{1}\{\hat{y}_i \neq y_i\} \left(= \frac{N_{\mathrm{FP}} + N_{\mathrm{FN}}}{N_{\text{test}}} = 1 - \text{Accuracy}\right). \]
Consider the trivial classifier \(\hat{y}_i^{\,\text{triv}} = 0\) for all \(i\), which grants every loan and uses none of the features. Write down its confusion matrix and verify that its accuracy equals \(N_{\mathrm{TN}}/N_{\text{test}}\), which is \(1-\hat{\pi}\) up to sampling error, and that its recall is \(0\). Compute the difference
\[ \text{gain} = \text{Accuracy}_{\text{model}} - (1-\hat{\pi}), \]
and comment on its size.
Examine the recall. Of the borrowers in \(\mathcal{D}_{\text{test}}\) who actually defaulted, what fraction did the model identify? In which cell of the confusion matrix do most of the model's correct answers lie, and does the trivial rule of (iii) obtain them as well?
Conclude: explain why accuracy is misleading when the two classes are strongly imbalanced, and state which of the quantities in (i) you would report to a credit risk committee instead.

(f) Further Experiments:

Lower the threshold from \(1/2\) to \(0.3\) and recompute the confusion matrix. Which errors increase and which decrease? Which of the two thresholds should a bank prefer, and why?
Which of \(N_{\mathrm{FP}}\) and \(N_{\mathrm{FN}}\) increases and which decreases? What happens to the accuracy, and what happens to the recall? Which of the two thresholds should a bank prefer, given that a type 2 error costs it the unrecovered principal whereas a type 1 error costs it only the interest margin it did not earn?
Rerun (c) without the standardisation of (b), keeping \(\eta = 0.1\). What happens, and why? Use the bound \(L \le \|\tilde{X}\|_2^2 / (4N_{\text{train}})\) of part (d) to explain the observation quantitatively, and find by trial a learning rate for which the unstandardised problem does converge.

A single logistic regression neuron is typically sufficient for smaller feature dimension, say \(n \le 100\). When the feature space is of high-dimansional or when there are significant dependencies between features, then the linearity assumption on the log-odds (through the log-likelihood ratio) is violated. In such situations, a fully connected feedforward network (MLP) with one or two hidden layers provides a flexible nonlinear extension. However, the interpretability (like the one given in Remark 3.3) of the model is difficult and often impossible.