Logistic regression from scratch: fit probabilities and inspect errors

Logistic regression maps a linear score to a probability using the sigmoid function. Fitting adjusts the weights to reduce binary cross entropy. Implementing that small calculation yourself makes the target encoding, numerical stability and prediction threshold visible.

We will generate an original two-feature classification task, fit on a training partition and evaluate once on a held-out partition. The labels contain randomness, so even the data-generating probabilities cannot guarantee perfect classification. This is a synthetic experiment, not a result for a real application.

Derive the loss and gradient

For score z = Xw + b, the probability of class 1 is sigmoid(z). Binary cross entropy can be written stably as logaddexp(0, z) - y*z. Averaging gives the training objective. We avoid calculating log(sigmoid(z)) directly, which can encounter numerical extremes.

The derivative with respect to a score is p - y. The weight gradient is X.T @ (p-y) / n, and the bias gradient is the mean residual. The implementation adds a constant column so the last weight acts as the bias.

The sigmoid helper uses logaddexp too. It produces the same mathematical probability while avoiding an unguarded exponential of a large positive number. Stable arithmetic prevents one class of bug; it does not repair a bad label definition or invalid evaluation split.

Run the fitted classifier

import numpy as np
rng = np.random.default_rng(26)
X = rng.normal(size=(600, 2))
probability = np.exp(-np.logaddexp(0, -(1.4*X[:, 0] - .8*X[:, 1] + .2)))
y = rng.binomial(1, probability)
A = np.column_stack([X, np.ones(len(X))])
train, test = A[:400], A[400:]
y_train, y_test = y[:400], y[400:]
w = np.zeros(3)
def loss(matrix, labels, weights):
    score = matrix @ weights
    return np.mean(np.logaddexp(0, score) - labels*score)
start = loss(train, y_train, w)
for _ in range(1500):
    p = np.exp(-np.logaddexp(0, -(train @ w)))
    w -= .1 * (train.T @ (p-y_train) / len(train))
p_test = np.exp(-np.logaddexp(0, -(test @ w)))
prediction = p_test >= .5
baseline = np.full(len(test), y_train.mean())
baseline_loss = -np.mean(y_test*np.log(baseline)+(1-y_test)*np.log1p(-baseline))
print(f'train loss: {start:.4f} -> {loss(train, y_train, w):.4f}')
print(f'test loss={loss(test, y_test, w):.4f} constant baseline={baseline_loss:.4f}')
print('confusion [TN FP; FN TP]:', [[int(np.sum((y_test==a)&(prediction==b)))
                                    for b in [0, 1]] for a in [0, 1]])
print('weights and bias:', np.round(w, 3).tolist())
assert loss(train, y_train, w) < start
assert np.isfinite(p_test).all()

Checked output:

train loss: 0.6931 -> 0.4600
test loss=0.5308 constant baseline=0.7069
confusion [TN FP; FN TP]: [[67, 33], [16, 84]]
weights and bias: [1.693, -0.811, 0.403]

Read the result beyond a single score

The learned weights should have the same broad direction as the generating rule, but finite samples and noisy labels prevent an exact match. The held-out loss compares probability assignments with a constant predictor using only the training class frequency. That baseline is fitted without seeing held-out labels.

The confusion matrix uses threshold 0.5. Changing the threshold changes false positives and false negatives without changing the trained weights. A threshold should reflect the decision and be chosen on development data, not repeatedly optimized against this final test set.

Training loss decreasing is a necessary implementation check here, not proof of useful generalization. Add an irrelevant feature or reduce the training sample and compare the held-out result. A more complicated linear score can fit noise while retaining an apparently successful training curve.

Know what the implementation leaves out

This model has no regularization or automatic stopping criterion. Its fixed learning rate works for this declared scale; larger feature magnitudes can change optimization behavior. On perfectly separable data, unregularized logistic regression may keep increasing weight magnitudes rather than settling at finite maximum-likelihood parameters.

A linear decision boundary cannot represent every relationship. Try labeling points by whether they lie inside a circle. The limitation then belongs to the representation, not necessarily the optimizer. Compare richer features using a valid development process rather than calling the original implementation broken.

Extend the classifier

Check one gradient coordinate using central differences, then add an L2 penalty to the feature weights while leaving the bias unpenalized. State the changed objective and derive the changed gradient before running it.

The Machine Learning course connects classification to baselines and evaluation. Review squared-error gradient checks for the numerical-check method and ROC versus PR summaries for ranking decisions beyond threshold accuracy.

References

scikit-learn logistic regression guide — The probabilistic classifier, regularization and solver context.

NumPy logaddexp — Stable computation of the logarithm of a sum of exponentials.