Logistic Regression From Scratch

Question

You want to model the probability that a limit order gets filled within one second, as a function of features such as queue position and spread.

First, derive the gradient of the log-likelihood for logistic regression. Then implement the training loop with gradient descent using numpy only (no scikit-learn).
def fit_logistic_regression(X: np.ndarray, y: np.ndarray, lr: float, n_iters: int) -> tuple:
    pass
Input
X is an (n_samples, n_features) numpy array, y is an (n_samples,) array of 0/1 labels, lr is the learning rate and n_iters the number of iterations.

Output
Return the tuple (w, b) of learned parameters.

Also answer
Why is logistic regression trained on the log-likelihood instead of the mean squared error?

Accepted Answer

Logistic regression models the fill probability as \begin{equation} p_i = P(y_i = 1 | x_i) = \sigma(z_i), \qquad z_i = x_i^T w + b \end{equation} where $\sigma$ is the sigmoid function \begin{equation} \sigma(z) = \frac{1}{1 + e^{-z}} \end{equation} which squashes any real number into a valid probability in (0, 1).

Step 1: a property of the sigmoid we'll needDifferentiate $\sigma(z) = (1 + e^{-z})^{-1}$ with the chain rule: \begin{equation} \sigma'(z) = -(1 + e^{-z})^{-2} \cdot (-e^{-z}) = \frac{e^{-z}}{(1 + e^{-z})^2} \end{equation} Split the fraction into two factors: \begin{equation} \sigma'(z) = \frac{1}{1 + e^{-z}} \cdot \frac{e^{-z}}{1 + e^{-z}} = \sigma(z) \big(1 - \sigma(z)\big) \end{equation} where the last step uses $\frac{e^{-z}}{1+e^{-z}} = \frac{1 + e^{-z} - 1}{1+e^{-z}} = 1 - \sigma(z)$.

Step 2: the likelihoodEach label is a Bernoulli draw with success probability $p_i$, which we can write in one line as $p_i^{y_i} (1-p_i)^{1-y_i}$ (check: it equals $p_i$ when $y_i = 1$ and $1 - p_i$ when $y_i = 0$). Assuming independent observations, the likelihood of the data is the product over all $i$, and taking the log turns the product into a sum: \begin{equation} \ell(w, b) = \sum_{i=1}^{n} \Big[ y_i \log p_i + (1 - y_i) \log (1 - p_i) \Big] \end{equation} Maximising $\ell$ is the maximum likelihood principle; the negative of $\ell$ is exactly the cross-entropy loss used in deep learning.

Step 3: the gradientWe differentiate one term of $\ell$ with respect to $z_i$ first, then push through to $w$. Using $\frac{\partial p_i}{\partial z_i} = p_i(1-p_i)$ from Step 1: \begin{equation} \frac{\partial}{\partial z_i} \Big[ y_i \log p_i + (1-y_i) \log(1-p_i) \Big] = y_i \frac{p_i(1-p_i)}{p_i} - (1-y_i)\frac{p_i(1-p_i)}{1-p_i} \end{equation} \begin{equation} = y_i (1 - p_i) - (1 - y_i) p_i = y_i - y_i p_i - p_i + y_i p_i = y_i - p_i \end{equation} A remarkably clean result: the sensitivity of the log-likelihood to the score $z_i$ is just the prediction error $y_i - p_i$. Since $z_i = x_i^T w + b$, the chain rule gives $\frac{\partial z_i}{\partial w} = x_i$ and $\frac{\partial z_i}{\partial b} = 1$, so \begin{equation} 
abla_w \ell = \sum_{i=1}^{n} (y_i - p_i) \, x_i, \qquad \frac{\partial \ell}{\partial b} = \sum_{i=1}^{n} (y_i - p_i) \end{equation} To maximise $\ell$ we do gradient ascent — or equivalently, gradient descent on the negative log-likelihood with gradient $X^T(p - y)$, which is what the code below implements (scaled by $1/n$ so the learning rate is sample-size independent).

Step 4: implementation
import numpy as np

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

def fit_logistic_regression(X: np.ndarray, y: np.ndarray, lr: float = 0.5, n_iters: int = 20000) -> tuple:
    n_samples, n_features = X.shape
    w = np.zeros(n_features)
    b = 0.0
    for _ in range(n_iters):
        p = sigmoid(X @ w + b)              # current fill probabilities
        grad_w = (1 / n_samples) * (X.T @ (p - y))   # gradient of the negative log-likelihood
        grad_b = (1 / n_samples) * np.sum(p - y)
        w -= lr * grad_w
        b -= lr * grad_b
    return w, b

# Test on synthetic data with known parameters
rng = np.random.default_rng(7)
true_w, true_b = np.array([1.5, -2.0]), 0.5
X = rng.normal(size=(5000, 2))
y = (rng.uniform(size=5000) < sigmoid(X @ true_w + true_b)).astype(float)

w, b = fit_logistic_regression(X, y)
print(w, b)  # approximately [1.5, -2.0] and 0.5

On this synthetic data the procedure recovers $w \approx (1.45, -1.98)$ and $b \approx 0.49$ — the true parameters up to sampling noise.

Step 5: why not mean squared error?Two reasons. First, plugging the sigmoid into the MSE loss creates a non-convex objective with flat plateaus: when the model is confidently wrong ($p_i$ near 0 or 1 on the wrong side), the factor $p_i(1-p_i)$ that appears in the MSE gradient is nearly zero, so learning stalls precisely where the model most needs correcting. The log-likelihood gradient $y_i - p_i$ has no such factor — a confidently wrong prediction produces a large gradient. Second, the negative log-likelihood is convex in $(w, b)$, so gradient descent finds the global optimum and the result is reproducible. As a bonus, maximum likelihood produces calibrated probabilities, which is what you actually need to make expected-value decisions on order placement. In production you would add an L2 penalty (which also handles separable data, where unregularised weights diverge to infinity) and validate the calibration out of sample.

Title	Category	Subcategory	Difficulty	Status
Naive Bayes From First Principles	Machine Learning	Classification	Medium	Unlock
The SVM Margin and the Kernel Trick	Machine Learning	Classification	Hard	Unlock

The Interview Guide for Quants and Traders

Logistic Regression From Scratch

What do you think of this page?