Lab: logistic regression in NumPy

Time to put the whole section together. By the end of this lab you’ll have a working classifier written entirely by you: about thirty lines of NumPy that do what LogisticRegression in a big library does.

Run the cells in order: they share one Python session, like a notebook, so later cells use the variables earlier ones made. If something goes wrong, use Run all above to rebuild the state.

Step 1: the data

We’ll simulate the admissions problem from the decision-boundary lesson: 200 applicants with two exam scores out of 10. Behind the scenes, an applicant’s chance of admission rises with $x_1 + x_2$, and the outcome is a coin flip with that chance, so the classes overlap a little, the way real data does. We keep 150 applicants for training and hold 50 back for testing.

import numpy as np

rng = np.random.default_rng(42)
m_all = 200
X_all = rng.uniform(0, 10, size=(m_all, 2))                       # exam 1, exam 2
true_p = 1 / (1 + np.exp(-(1.2 * (X_all[:, 0] + X_all[:, 1]) - 12)))
y_all = (rng.uniform(size=m_all) < true_p).astype(int)            # 1 = admitted

X_train, y_train = X_all[:150], y_all[:150]
X_test, y_test = X_all[150:], y_all[150:]
print("training set:", X_train.shape, " admitted:", y_train.sum())
print("test set:    ", X_test.shape, " admitted:", y_test.sum())
print("first three applicants:\n", X_train[:3].round(2), y_train[:3])

Step 2: scale the features

Section 2 taught you to put features on a similar scale before gradient descent. We compute the mean and standard deviation on the training set only and apply the same numbers to the test set; using test data here would be leakage (Python for Data, section 3).

mu = X_train.mean(axis=0)
sigma = X_train.std(axis=0)

def scale(X):
    return (X - mu) / sigma

Xs_train, Xs_test = scale(X_train), scale(X_test)
print("scaled means:", Xs_train.mean(axis=0).round(3), " scaled stds:", Xs_train.std(axis=0).round(3))

Step 3: the model and the cost

These are the pieces you wrote in the exercises earlier in the section. Here they are together, vectorised over all examples.

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

def predict_proba(X, w, b):
    return sigmoid(X @ w + b)

def compute_cost(X, y, w, b, lambda_=0.0):
    f = predict_proba(X, w, b)
    m = len(y)
    log_loss = -np.mean(y * np.log(f) + (1 - y) * np.log(1 - f))
    return log_loss + lambda_ / (2 * m) * np.sum(w ** 2)

def compute_gradient(X, y, w, b, lambda_=0.0):
    m = len(y)
    err = predict_proba(X, w, b) - y
    dj_dw = X.T @ err / m + (lambda_ / m) * w
    dj_db = np.mean(err)
    return dj_dw, dj_db

w0, b0 = np.zeros(2), 0.0
print("cost at w = 0, b = 0:", round(compute_cost(Xs_train, y_train, w0, b0), 4), "(that's log 2)")

At $\vec{w} = 0$, $b = 0$ every prediction is 0.5 and the cost is $-\log 0.5 = \log 2 \approx 0.693$. That’s a useful sanity check: a model that knows nothing scores 0.693.

Step 4: train it

Now the heart of it. Write gradient descent: repeat the update iterations times and record the cost so you can check convergence.

import numpy as np

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

def compute_cost(X, y, w, b, lambda_=0.0):
    f = sigmoid(X @ w + b)
    return -np.mean(y * np.log(f) + (1 - y) * np.log(1 - f)) + lambda_ / (2 * len(y)) * np.sum(w ** 2)

def compute_gradient(X, y, w, b, lambda_=0.0):
    err = sigmoid(X @ w + b) - y
    return X.T @ err / len(y) + (lambda_ / len(y)) * w, np.mean(err)

def gradient_descent(X, y, w, b, alpha, iterations, lambda_=0.0):
    history = []
    for i in range(iterations):
        pass  # replace: gradient, update w and b, record the cost
    return w, b, history

X = np.array([[-1.2, -0.8], [-0.5, -1.0], [-0.9, 0.1], [0.4, 1.1], [1.0, 0.3], [0.6, 0.9]])
y = np.array([0, 0, 0, 1, 1, 1])
w, b, history = gradient_descent(X, y, np.zeros(2), 0.0, alpha=0.5, iterations=200)
print("final cost:", history[-1] if history else "no history yet")

Once your function passes, run the cell below. It trains on the scaled training data and plots the cost history. (If you skipped the exercise, the cell uses a reference version.)

if "gradient_descent" not in globals():
    def gradient_descent(X, y, w, b, alpha, iterations, lambda_=0.0):
        history = []
        for _ in range(iterations):
            dj_dw, dj_db = compute_gradient(X, y, w, b, lambda_)
            w = w - alpha * dj_dw
            b = b - alpha * dj_db
            history.append(compute_cost(X, y, w, b, lambda_))
        return w, b, history

w, b, history = gradient_descent(Xs_train, y_train, np.zeros(2), 0.0, alpha=0.5, iterations=500)
print(f"cost: start {history[0]:.4f} -> end {history[-1]:.4f}")
print("w =", w.round(3), " b =", round(b, 3))
bpl.plot(np.arange(len(history)), np.array(history), kind="line", title="Training cost per iteration")

The curve should drop steeply and then flatten: gradient descent has converged. If yours rises or zig-zags, the learning rate is too large.

Step 5: how good is it?

Accuracy is the fraction of examples classified correctly. The number that matters is on the test set, the 50 applicants the model never saw.

def predict(X, w, b, threshold=0.5):
    return (predict_proba(X, w, b) >= threshold).astype(int)

train_acc = np.mean(predict(Xs_train, w, b) == y_train)
test_acc = np.mean(predict(Xs_test, w, b) == y_test)
print(f"training accuracy: {train_acc:.1%}")
print(f"test accuracy:     {test_acc:.1%}")

new = np.array([[6.5, 7.0], [4.0, 4.5]])
for scores, p in zip(new, predict_proba(scale(new), w, b)):
    print(f"applicant with scores {scores}: chance of admission {p:.0%}")

Training and test accuracy are close, so the model isn’t overfitting; with only three parameters it hardly can. It isn’t 100 % because the data really is noisy: some applicants with high scores were unlucky. Notice the model gives real probabilities, not just yes or no.

Step 6: a curved boundary

Now a harder shape. In this dataset, class 1 lives inside a circle, like the Ring dataset in the widget below. Set its degree to 1 and to 2 and press Run to see why a straight line can’t work and a degree-2 model can.

To do the same in NumPy you need polynomial features: build the columns $x_1$, $x_2$, $x_1^2$, $x_1 x_2$, $x_2^2$ from the two inputs, then run exactly the same logistic regression on the wider matrix. Write the feature map, then train with a little regularisation.

import numpy as np

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

def train(X, y, alpha=0.5, iterations=2000, lambda_=0.1):
    w, b, m = np.zeros(X.shape[1]), 0.0, len(y)
    for _ in range(iterations):
        err = sigmoid(X @ w + b) - y
        w = w - alpha * (X.T @ err / m + (lambda_ / m) * w)
        b = b - alpha * np.mean(err)
    return w, b

def map_features(X):
    x1, x2 = X[:, 0], X[:, 1]
    # return the columns x1, x2, x1^2, x1*x2, x2^2 side by side
    return X  # replace this line

rng = np.random.default_rng(7)
X_ring = rng.uniform(-2, 2, size=(160, 2))
y_ring = (X_ring[:, 0] ** 2 + X_ring[:, 1] ** 2 < 1.6).astype(int)

F = map_features(X_ring)
w, b = train(F, y_ring)
accuracy = np.mean((sigmoid(F @ w + b) >= 0.5) == y_ring)
print("features per example:", F.shape[1], " training accuracy:", round(accuracy, 3))

That’s logistic regression, built from nothing but NumPy: a linear part, a sigmoid, the log loss, its gradient and a loop. Neural networks in the next section stack many of these units together, and everything you’ve learned here, from the cost to gradient descent to overfitting and regularisation, carries straight over.

Section quiz

Ten questions on the whole section: classification, the sigmoid, decision boundaries, the log loss, gradient descent, overfitting and regularisation.