Logistic Regression — exercises

(a) $\sigma (-5.77+1.28\cdot 5)=\sigma (0.63)=1/(1+e^{-0.63})\approx 0.65$.

(b) $P=0.5\iff {\beta }_0+{\beta }_{1}b=0\iff b=-{\beta }_0/{\beta }_1=5.77/1.28\approx 4.5$. The model is undecided at biomarker level 4.5.

(c) The model says $\log \frac{P}{1-P}={\beta }_0+{\beta }_{1}b$. A unit increase in $b$ raises the log-odds by ${\beta }_1=1.28$, multiplying the odds by $e^{1.28}\approx 3.60$. So each additional unit of biomarker multiplies the disease odds by 3.6×. That's the natural-language sentence epidemiologists actually write into clinical papers — odds ratios are what logistic regression gives you for free.

import numpy as np

def sigmoid(z):
    return np.where(z >= 0,
                    1.0 / (1.0 + np.exp(-z)),
                    np.exp(z) / (1.0 + np.exp(z)))

# Ten labelled (age, impressions, clicked) rows.
data = np.array([
    [22,  2, 1],  [28,  5, 1],  [31,  9, 0],  [35,  3, 1],  [40, 15, 0],
    [45,  4, 1],  [52, 18, 0],  [55,  7, 1],  [60, 20, 0],  [62, 12, 0],
])
X_raw = data[:, :2].astype(float)
y     = data[:, 2].astype(float)
n     = len(y)

# Standardize features so GD has comparable scales.
X_raw = (X_raw - X_raw.mean(axis=0)) / X_raw.std(axis=0)
X     = np.column_stack([np.ones(n), X_raw])

beta = np.zeros(3)
for _ in range(20000):
    p = sigmoid(X @ beta)
    grad = X.T @ (p - y) / n
    beta -= 0.5 * grad

print("β        :", beta)
print("accuracy :", np.mean((sigmoid(X @ beta) >= 0.5) == y))

Fit returns approximately ${\beta }_0\approx -0.3$, ${\beta }_1\approx -0.7$ (on standardized age), ${\beta }_2\approx -2.0$ (on standardized imps). In-sample accuracy 9/10.

The number of impressions has the largest-magnitude negative coefficient — the more impressions of the ad the user has already seen, the less likely they are to click. (Banner blindness.) Age also leans negative but is dominated by imps.

Important caveat with $n=10$: these coefficients are noisy. The point of the exercise is the mechanics; in production you'd never trust a 10-sample fit.

(a) Boundary equation: $1+2x_1-3x_2=0$, or rearranged $x_2=(1+2x_1)/3$.

(b) The gradient of ${\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2$ in the $(x_1,x_2)$-plane is $(2,-3)$. This is the direction of increasing log-odds, hence of increasing probability. Normalize: unit vector $=(2,-3)/\sqrt{4+9}=(2/\sqrt{13},-3/\sqrt{13})\approx (0.555,-0.832)$.

(c) The signed distance from the origin to the hyperplane $a{x}_1+b{x}_2+c=0$ is $-c/\sqrt{a^2+b^2}$. Here $a=2,b=-3,c=1$, so distance $=-1/\sqrt{13}\approx -0.277$. The origin is on the $y=0$ side (negative).

(a) The classes are separated by the curve $x_1^2+x_2^2=1.5^2$ (a circle), which is not a straight line. A linear-in-$(x_1,x_2)$ model can only carve out a half-plane, so any choice of $\beta$ miscalls at least half the rings.

(b) Code:

import numpy as np

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

rng = np.random.default_rng(2)
n = 50
# Class 0: small radius ring.  Class 1: large radius ring.
r0 = rng.uniform(0, 1, n); t0 = rng.uniform(0, 2*np.pi, n)
r1 = rng.uniform(2, 3, n); t1 = rng.uniform(0, 2*np.pi, n)
inner = np.column_stack([r0 * np.cos(t0), r0 * np.sin(t0)])
outer = np.column_stack([r1 * np.cos(t1), r1 * np.sin(t1)])
X_raw = np.vstack([inner, outer])
y     = np.concatenate([np.zeros(n), np.ones(n)])

def fit(X, lr=0.5, n_iter=20000):
    beta = np.zeros(X.shape[1])
    for _ in range(n_iter):
        p = sigmoid(X @ beta)
        beta -= lr * X.T @ (p - y) / len(y)
    return beta

# Linear features only.
X_lin = np.column_stack([np.ones(2*n), X_raw])
beta_lin = fit(X_lin)
acc_lin = np.mean((sigmoid(X_lin @ beta_lin) >= 0.5) == y)

# Add x², y² features.
X_poly = np.column_stack([np.ones(2*n), X_raw, X_raw[:,0]**2, X_raw[:,1]**2])
beta_poly = fit(X_poly)
acc_poly = np.mean((sigmoid(X_poly @ beta_poly) >= 0.5) == y)

print("linear features  accuracy:", acc_lin)   # ~ 0.6, near chance
print("polynomial-feat. accuracy:", acc_poly)  # ~ 1.0

Linear-feature accuracy ≈ 0.6 (essentially chance). Polynomial-feature accuracy ≈ 1.0.

Why it works: the boundary is now ${\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_1^2+{\beta }_4{x}_2^2=0$, which is a conic section. For roughly circular data the fit converges with ${\beta }_3\approx {\beta }_4>0$ and ${\beta }_1,{\beta }_2\approx 0$, recovering ${\beta }_3(x_1^2+x_2^2)=-{\beta }_0$ — a circle of radius $\sqrt{-{\beta }_0/{\beta }_3}$.

(c) The model is still linear in parameters (it's the same OLS-style design matrix, just with extra columns). The hidden trick is that polynomial features are a poor person's kernel — neural networks generalize this by learning the nonlinear features rather than hand-engineering them.

Weighted loss:

${\mathcal{L}}_w(\beta )=-\frac{1}{n}{\sum }_i[w_{+}{y}_i\log {\hat{p}}_i+w_{-}(1-y_i)\log (1-{\hat{p}}_i)]$ where ${\hat{p}}_i=\sigma (x_i^{\top }\beta )$.

Differentiating and using ${\sigma }^{\prime }(z)=\sigma (z)(1-\sigma (z))$, the contribution from row $i$ is $w_{+}{y}_i({\hat{p}}_i-1)+w_{-}(1-y_i){\hat{p}}_i=(w_{+}{y}_i+w_{-}(1-y_i))({\hat{p}}_i-y_i)$.

So the gradient is $\nabla {\mathcal{L}}_w(\beta )=\frac{1}{n}{X}^{\top }W(\hat{p}-y)$ with $W=\mathrm{diag}(w_{+}{y}_i+w_{-}(1-y_i))$ — the unweighted form times a row-rescaling.

Setting $w_{+}=1/p_{+}$ and $w_{-}=1/p_{-}$ (inverse class frequencies) makes positives and negatives contribute equal total weight, so the loss "sees" a balanced dataset. Practical equivalent: oversample the minority class to 50%. Both fixes solve the same problem with the same algebra under the hood.

(a) Perfect calibration is the 45° diagonal: bucket $b$ has mean predicted probability $\hat{p}\approx 0.05+0.1b$ and empirical rate equal to that same number. A well-calibrated model lies along $y=x$.

(b) Over-confidence at the extremes looks like a sigmoid: shallow at the ends (predicted 0.05 maps to actual 0.15, so the curve sits above the diagonal for low $\hat{p}$), steep in the middle, shallow again at the high end (predicted 0.95 maps to actual 0.85, below the diagonal). The curve bows above the diagonal on the left and below it on the right.

(c) Standard fixes: (i) Platt scaling — fit a second 1D logistic regression that maps the model's raw logits to calibrated probabilities. (ii) Isotonic regression — fit a monotone step function from raw to calibrated probabilities. (iii) Temperature scaling for neural networks — divide logits by a learned scalar $T>1$ to soften over-confident outputs. (iv) Add regularization during training; over-confidence often comes from unconstrained $\Vert \beta \Vert$.

(a) Negative log-likelihood for one point: $\mathcal{L}=-\log P(y\mid x)=-{\sum }_k{y}_k\log {\hat{p}}_k$ with ${\hat{p}}_k=e^{x^{\top }{\beta }_k}/Z$ and $Z={\sum }_j{e}^{x^{\top }{\beta }_j}$.

(b) Differentiate $\log {\hat{p}}_k=x^{\top }{\beta }_k-\log Z$ with respect to ${\beta }_{k^{\prime }}$:

${\partial }_{{\beta }_{k^{\prime }}}\log {\hat{p}}_k=1[k=k^{\prime }]x-\frac{\partial \log Z}{\partial {\beta }_{k^{\prime }}}=(1[k=k^{\prime }]-{\hat{p}}_{k^{\prime }})x$

Since only one $y_k$ is 1 (call it $y_{\ast }$), ${\partial }_{{\beta }_{k^{\prime }}}\mathcal{L}=-(1[y_{\ast }=k^{\prime }]-{\hat{p}}_{k^{\prime }})x=({\hat{p}}_{k^{\prime }}-y_{k^{\prime }})x$ ✓

(c) Stack across the dataset: with $\hat{P}\in {\mathbb{R}}^{n\times K}$ holding predicted probabilities and $Y\in {\mathbb{R}}^{n\times K}$ the one-hot label matrix, the gradient with respect to the weight matrix $B\in {\mathbb{R}}^{p\times K}$ is ${\nabla }_B\mathcal{L}=X^{\top }(\hat{P}-Y)/n$. Exact same shape as binary; softmax is binary LogR generalized cleanly. This is also the gradient at the output layer of every classification neural network.