The Maths That Actually Matters

Week 3 of 16 · Foundations · 7 days

Full curriculum
Week 03 · Foundations

The Maths That Actually Matters

Week 03 · Day 1 of 7

Vectors, Matrices and Why Models Are Multiplications

Dot products, matrix shapes, distance, and the one line most models are

By 730 words

You do not need a maths degree to do machine learning well. You do need four ideas, and this week covers them by implementing each one rather than proving anything. Today: why almost every model is, underneath, a matrix multiplication.

A row of data is a vector

One customer, tenure 12, charge 79.50, support calls 2, is a point in three-dimensional space. Two thousand customers are two thousand points. Everything that follows is geometry on those points.

import numpy as np

customer = np.array([12.0, 79.50, 2.0])
weights = np.array([-0.05, 0.02, 0.30])

score = np.dot(customer, weights)
manual = sum(c * w for c, w in zip(customer, weights))
print('dot product : %.4f' % score)
print('by hand : %.4f' % manual)
dot product : 1.5900
by hand : 1.5900
Dot product: Multiply matching elements and add up the results. It is the single most important operation in machine learning: a linear model's prediction is the dot product of a row of features with a vector of learned weights, and a neural network layer is a stack of them.

The whole dataset at once

You never score one customer. You score all of them, and the matrix product does every dot product in one call.

import numpy as np

X = np.array([[12.0, 79.50, 2.0],
[48.0, 55.00, 0.0],
[ 3.0, 92.00, 5.0],
[60.0, 21.00, 1.0]])
w = np.array([-0.05, 0.02, 0.30])
b = -1.2

scores = X @ w + b
print('X shape', X.shape, '@ w shape', w.shape, '-> ', scores.shape)
print('scores', scores.round(3))
X shape (4, 3) @ w shape (3,) -> (4,)
scores [ 0.39 -2.5 1.99 -3.48]

That line is a linear model

X @ w + b is the entire forward pass of linear regression, and with a sigmoid wrapped round it, of logistic regression. Week 11's neural network is the same line repeated with a nonlinearity between the repeats. Understand this one expression and you understand the shape of most of the course.

Shapes have to line up

A matrix product (n, k) @ (k, m) gives (n, m). The inner dimensions must match, and this is the source of most errors you will hit in week 11.

import numpy as np

X = np.random.default_rng(0).normal(size=(100, 4)) # 100 rows, 4 features
W1 = np.random.default_rng(1).normal(size=(4, 8)) # layer: 4 in, 8 out
W2 = np.random.default_rng(2).normal(size=(8, 1)) # layer: 8 in, 1 out

h = X @ W1
out = h @ W2
print('X ', X.shape)
print('X@W1', h.shape)
print('h@W2', out.shape)

try:
X @ W2
except ValueError as e:
print('\nmismatch:', e)
X (100, 4)
X@W1 (100, 8)
h@W2 (100, 1)

mismatch: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 8 is different from 4)

Distance, and what it means for a model

kNN, k-means and every clustering algorithm in week 9 rest on measuring how far apart two rows are.

import numpy as np

a = np.array([12.0, 79.5])
b = np.array([14.0, 81.0])
c = np.array([60.0, 21.0])

def euclid(p, q):
return np.sqrt(((p - q) ** 2).sum())

print('a to b: %.2f' % euclid(a, b))
print('a to c: %.2f' % euclid(a, c))
print('numpy : %.2f' % np.linalg.norm(a - c))
a to b: 2.50
a to c: 75.67
numpy : 75.67

Distance is meaningless on unscaled columns

Tenure runs 1 to 72; total charges run to 8000. In the distance calculation, charges dominate entirely and tenure might as well not exist. Every distance-based algorithm requires scaling first, and this is the real reason StandardScaler appears in nearly every pipeline in this course.

import numpy as np
from sklearn.preprocessing import StandardScaler

raw = np.array([[12.0, 950.0],
[14.0, 8000.0],
[60.0, 1000.0]])

print('unscaled: row0-row1 %.1f, row0-row2 %.1f'
% (np.linalg.norm(raw[0] - raw[1]), np.linalg.norm(raw[0] - raw[2])))

scaled = StandardScaler().fit_transform(raw)
print('scaled : row0-row1 %.2f, row0-row2 %.2f'
% (np.linalg.norm(scaled[0] - scaled[1]),
np.linalg.norm(scaled[0] - scaled[2])))
unscaled: row0-row1 7050.0, row0-row2 69.3
scaled : row0-row1 2.13, row0-row2 2.17

Unscaled, row 0 looks far closer to row 2 than to row 1, because the charge difference swamps everything. Scaled, the answer flips, and the scaled answer is the one that respects both columns.

Day 1 takeaway

A row is a vector, a dataset is a matrix, and a linear model's prediction is X @ w + b. Matrix shapes must agree on the inner dimension. Distance underpins every clustering and neighbour method, and is meaningless until the columns share a scale.
Week 03 · Day 2 of 7

Derivatives and Gradients

Slopes, the chain rule, and why the sigmoid derivative caps at 0.25

By 757 words

Training a model means adjusting numbers until a loss gets smaller. To know which direction to adjust, you need the derivative. That is the whole role calculus plays here.

A derivative is a slope

Derivative: How much the output changes for a tiny change in the input. Positive means increasing the input increases the output; negative means the opposite; zero means you are at a flat point, which for a loss function usually means a minimum.
import numpy as np

def f(x):
return x ** 2 + 3 * x + 5

def numeric_derivative(f, x, h=1e-6):
return (f(x + h) - f(x - h)) / (2 * h)

# The exact derivative of x^2 + 3x + 5 is 2x + 3.
for x in [-4.0, -1.5, 0.0, 3.0]:
print('x=%5.1f numeric %8.4f exact %8.4f'
% (x, numeric_derivative(f, x), 2 * x + 3))
x= -4.0 numeric -5.0000 exact -5.0000
x= -1.5 numeric -0.0000 exact 0.0000
x= 0.0 numeric 3.0000 exact 3.0000
x= 3.0 numeric 9.0000 exact 9.0000

At x = -1.5 the derivative is zero, and that is the bottom of the parabola. Finding where the derivative is zero is finding the minimum, which is exactly what training does.

The gradient is the derivative in several directions

Models have thousands of parameters, not one. The gradient collects the derivative with respect to each of them into a vector that points in the direction of steepest increase.

import numpy as np

def loss(w):
# A bowl, steeper in the second direction than the first.
return (w[0] - 2) ** 2 + 4 * (w[1] + 1) ** 2

def gradient(f, w, h=1e-6):
g = np.zeros_like(w)
for i in range(len(w)):
step = np.zeros_like(w)
step[i] = h
g[i] = (f(w + step) - f(w - step)) / (2 * h)
return g

w = np.array([0.0, 0.0])
g = gradient(loss, w)
print('loss at [0,0] %.3f' % loss(w))
print('gradient %s' % g.round(4))
print('exact %s' % np.array([2 * (w[0] - 2), 8 * (w[1] + 1)]))
print('\nstepping against the gradient:')
print('loss at w - 0.1*g %.3f' % loss(w - 0.1 * g))
loss at [0,0] 8.000
gradient [-4. 8.]
exact [-4. 8.]

stepping against the gradient:
loss at w - 0.1*g 2.720

Against the gradient, not along it

The gradient points uphill. To reduce a loss you step in the opposite direction, which is why the update rule has a minus sign: w = w - learning_rate * gradient. That one line is gradient descent, and tomorrow you implement it.

The chain rule, which is what backpropagation is

Compose two functions and the derivative multiplies. This is the whole mechanism by which a neural network assigns credit to a weight buried four layers from the output.

import numpy as np

def g(x):
return 3 * x + 1

def f(u):
return u ** 2

def composed(x):
return f(g(x))

x = 2.0
# chain rule: d/dx f(g(x)) = f'(g(x)) * g'(x) = 2*g(x) * 3
chain = 2 * g(x) * 3
numeric = (composed(x + 1e-6) - composed(x - 1e-6)) / 2e-6
print('chain rule %.4f' % chain)
print('numeric %.4f' % numeric)
chain rule 42.0000
numeric 42.0000

The derivatives you will actually meet

FunctionWhere it appearsDerivative
x²Squared error loss2x
sigmoid(x)Logistic output layersigmoid(x)(1 − sigmoid(x))
relu(x)Hidden layers1 if x > 0, else 0
log(x)Log loss1/x
import numpy as np

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

def relu(x):
return np.maximum(0, x)

xs = np.array([-2.0, -0.5, 0.0, 0.5, 2.0])
h = 1e-6
print(' x sigmoid d/dx numeric d/dx formula relu d/dx')
for x in xs:
num = (sigmoid(x + h) - sigmoid(x - h)) / (2 * h)
formula = sigmoid(x) * (1 - sigmoid(x))
rnum = (relu(x + h) - relu(x - h)) / (2 * h)
print('%5.1f %7.4f %12.6f %12.6f %5.1f %4.1f'
% (x, sigmoid(x), num, formula, relu(x), rnum))
x sigmoid d/dx numeric d/dx formula relu d/dx
-2.0 0.1192 0.104994 0.104994 0.0 0.0
-0.5 0.3775 0.235004 0.235004 0.0 0.0
0.0 0.5000 0.250000 0.250000 0.0 0.5
0.5 0.6225 0.235004 0.235004 0.5 1.0
2.0 0.8808 0.104994 0.104994 2.0 1.0

The sigmoid derivative peaks at 0.25

Look at the column: the largest value is 0.25, at x = 0, and it falls away fast. Multiply several of those together through a deep network and the gradient reaching the early layers is vanishingly small. That is the vanishing gradient problem, it is why sigmoid was abandoned for hidden layers, and it is why ReLU, whose derivative is exactly 1 for positive inputs, replaced it. Week 11 returns to this.

Day 2 takeaway

A derivative is a slope; a gradient is the vector of slopes across all parameters. Training steps against the gradient. The chain rule multiplies derivatives through composed functions, which is exactly what backpropagation does, and the fact that the sigmoid derivative never exceeds 0.25 is why deep networks needed ReLU.
Week 03 · Day 3 of 7

Gradient Descent, Implemented

The four-line loop, and the three ways learning rate and scale break it

By 1055 words

Yesterday gave you the direction. Today you take the steps, and watch the three ways this goes wrong.

The algorithm, in full

  1. Start with a guess for the parameters.
  2. Compute the loss and its gradient.
  3. Move the parameters a small step against the gradient.
  4. Repeat until the loss stops falling.
import numpy as np

rng = np.random.default_rng(0)
n = 200
x = rng.uniform(0, 10, n)
y = 3.5 * x + 2.0 + rng.normal(0, 1.5, n) # truth: slope 3.5, intercept 2

w, b = 0.0, 0.0
lr = 0.01

for step in range(1, 2001):
pred = w * x + b
error = pred - y
loss = (error ** 2).mean()
dw = 2 * (error * x).mean()
db = 2 * error.mean()
w -= lr * dw
b -= lr * db
if step in (1, 10, 100, 500, 2000):
print('step %4d loss %8.4f w %.4f b %.4f' % (step, loss, w, b))

print('\ntruth w 3.5000 b 2.0000')
print('recovered w %.4f b %.4f' % (w, b))
step 1 loss 541.8747 w 2.8686 b 0.4151
step 10 loss 2.8682 w 3.6719 b 0.5939
step 100 loss 2.5705 w 3.5999 b 1.1008
step 500 loss 2.3502 w 3.4833 b 1.9210
step 2000 loss 2.3448 w 3.4621 b 2.0705

truth w 3.5000 b 2.0000
recovered w 3.4621 b 2.0705

Two thousand steps of arithmetic recovered the rule that generated the data. No library, no solver. The same loop, at enormous scale, is what trains a neural network.

The learning rate is the whole game

import numpy as np

rng = np.random.default_rng(0)
x = rng.uniform(0, 10, 200)
y = 3.5 * x + 2.0 + rng.normal(0, 1.5, 200)

def descend(lr, steps=300):
w = b = 0.0
for _ in range(steps):
error = (w * x + b) - y
w -= lr * 2 * (error * x).mean()
b -= lr * 2 * error.mean()
# Overflow passes through astronomically large but still finite
# values on its way to inf, so test the magnitude, not isfinite.
if not np.isfinite(w) or abs(w) > 1e10:
return None, None, None
return w, b, (((w * x + b) - y) ** 2).mean()

for lr in [0.0001, 0.001, 0.01, 0.03, 0.05]:
w, b, loss = descend(lr)
if w is None:
print('lr %-7s diverged to infinity' % lr)
else:
print('lr %-7s w %7.4f b %7.4f loss %8.4f' % (lr, w, b, loss))
lr 0.0001 w 3.3268 b 0.4947 loss 7.8274
lr 0.001 w 3.6533 b 0.7253 loss 2.7752
lr 0.01 w 3.5162 b 1.6897 loss 2.3793
lr 0.03 diverged to infinity
lr 0.05 diverged to infinity

Too small wastes time, too large explodes

At 0.0001 the model has barely moved after 300 steps. At 0.05 each step overshoots the minimum by more than it started from, the error grows, and the parameters run to infinity within a few dozen iterations. If your loss becomes nan during training, the learning rate is the first thing to look at.

Scaling changes how hard the problem is

import numpy as np

rng = np.random.default_rng(1)
n = 500
tenure = rng.uniform(1, 72, n)
charges = rng.uniform(20, 120, n)
total = tenure * charges # ranges into the thousands
y = 0.4 * tenure + 0.05 * charges + 0.001 * total + rng.normal(0, 1, n)

def descend(X, y, lr, steps=400):
# A column of ones so the model can fit an intercept. Without it a
# scaled matrix has mean zero and cannot reach a non-zero target.
Xb = np.column_stack([np.ones(len(X)), X])
w = np.zeros(Xb.shape[1])
for _ in range(steps):
error = Xb @ w - y
w -= lr * 2 * (Xb.T @ error) / len(y)
if not np.all(np.isfinite(w)) or np.abs(w).max() > 1e10:
return None
return (((Xb @ w) - y) ** 2).mean()

raw = np.column_stack([tenure, charges, total])
scaled = (raw - raw.mean(axis=0)) / raw.std(axis=0)

for lr in [1e-7, 1e-5, 1e-3, 1e-2]:
r = descend(raw, y, lr)
s = descend(scaled, y, lr)
print('lr %-6s raw %-12s scaled %s'
% (lr, 'diverged' if r is None else '%.4f' % r,
'diverged' if s is None else '%.4f' % s))
lr 1e-07 raw 38.9770 scaled 518.7263
lr 1e-05 raw diverged scaled 509.3600
lr 0.001 raw diverged scaled 92.7342
lr 0.01 raw diverged scaled 1.4854

On raw columns only a tiny learning rate survives, and it converges slowly. On scaled columns a much larger rate works and gets further. Scaling is not cosmetic. It changes the shape of the loss surface from a long thin valley into something closer to a bowl, and gradient descent handles bowls far better.

Batch, stochastic and mini-batch

VariantGradient computed onTrade-off
BatchEvery row, every stepSmooth and accurate, slow on large data
Stochastic (SGD)One row at a timeFast and noisy; the noise can escape shallow minima
Mini-batch32 to 512 rowsThe compromise everyone actually uses
import numpy as np

rng = np.random.default_rng(0)
n = 2000
x = rng.uniform(0, 10, n)
y = 3.5 * x + 2.0 + rng.normal(0, 1.5, n)

def minibatch(batch_size, epochs=30, lr=0.01, seed=0):
r = np.random.default_rng(seed)
w = b = 0.0
for _ in range(epochs):
order = r.permutation(n)
for start in range(0, n, batch_size):
idx = order[start:start + batch_size]
xb, yb = x[idx], y[idx]
error = (w * xb + b) - yb
w -= lr * 2 * (error * xb).mean()
b -= lr * 2 * error.mean()
return w, b, (((w * x + b) - y) ** 2).mean()

for bs in [1, 32, 256, n]:
w, b, loss = minibatch(bs)
label = 'stochastic' if bs == 1 else ('batch' if bs == n else 'mini %d' % bs)
print('%-12s w %.4f b %.4f loss %.4f' % (label, w, b, loss))
stochastic w 3.1202 b 2.0565 loss 6.8004
mini 32 w 3.4412 b 2.0094 loss 2.3413
mini 256 w 3.5696 b 1.5688 loss 2.2876
batch w 3.6854 b 0.7612 loss 2.6329

The truth is w = 3.5, b = 2.0. The mini-batch runs get closest; full batch is furthest off, with the intercept still at 0.76 after 30 epochs. That is the point: batch made 30 parameter updates in total, mini-batch-32 made nearly two thousand, and stochastic made sixty thousand. Same data, same epochs, vastly different amounts of learning, which is why nobody trains on full batches. Stochastic is noisiest, visible in its slope of 3.12.

Day 3 takeaway

Gradient descent is four lines: predict, measure error, compute the gradient, step against it. The learning rate decides whether it crawls, converges or explodes. Scaling reshapes the loss surface and makes larger rates safe. Mini-batches are the practical compromise between accuracy per step and steps per second.
Week 03 · Day 4 of 7

Probability You Actually Need

Conditional probability, Bayes, base rates and the independence assumption

By 740 words

Classifiers do not output labels. They output probabilities, and something downstream turns those into a decision. Today is the probability you need to read that output correctly.

Conditional probability

import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()

p_churn = df['churned'].mean()
p_m2m = (df['contract'] == 'Month-To-Month').mean()
p_churn_given_m2m = df.loc[df['contract'] == 'Month-To-Month', 'churned'].mean()
p_m2m_given_churn = (df.loc[df['churned'] == 1, 'contract'] == 'Month-To-Month').mean()

print('P(churn) %.3f' % p_churn)
print('P(month-to-month) %.3f' % p_m2m)
print('P(churn | month-to-month) %.3f' % p_churn_given_m2m)
print('P(month-to-month | churn) %.3f' % p_m2m_given_churn)
P(churn) 0.268
P(month-to-month) 0.553
P(churn | month-to-month) 0.414
P(month-to-month | churn) 0.854

These last two are not the same number

P(churn given month-to-month) is 0.41. P(month-to-month given churn) is 0.85. Confusing the two is the single most common probability error in analytics, and it has a name, the prosecutor's fallacy. “85 percent of churners were month-to-month” sounds like month-to-month customers mostly churn. They mostly do not.

Bayes' theorem, which converts between them

import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()

p_churn = df['churned'].mean()
p_m2m = (df['contract'] == 'Month-To-Month').mean()
p_m2m_given_churn = (df.loc[df['churned'] == 1, 'contract'] == 'Month-To-Month').mean()

# P(A|B) = P(B|A) * P(A) / P(B)
bayes = p_m2m_given_churn * p_churn / p_m2m
actual = df.loc[df['contract'] == 'Month-To-Month', 'churned'].mean()
print('Bayes %.4f' % bayes)
print('actual %.4f' % actual)
Bayes 0.4139
actual 0.4139

Why base rates dominate

A classic: a test for a rare condition is 99 percent accurate. You test positive. What is the chance you have it?

prevalence = 0.001 # 1 in 1000 people have it
sensitivity = 0.99 # correctly flags 99% of those who do
specificity = 0.99 # correctly clears 99% of those who do not

true_pos = prevalence * sensitivity
false_pos = (1 - prevalence) * (1 - specificity)
posterior = true_pos / (true_pos + false_pos)

print('per 100,000 people:')
print(' have it and test positive %6.0f' % (100_000 * true_pos))
print(' healthy but test positive %6.0f' % (100_000 * false_pos))
print('\nP(condition | positive test) %.3f' % posterior)
per 100,000 people:
have it and test positive 99
healthy but test positive 999

P(condition | positive test) 0.090

Nine percent, from a 99 percent accurate test

Because the condition is rare, the far larger healthy group generates ten times more false positives than the sick group generates true ones. This is exactly the situation in fraud detection, disease screening and any rare-event problem, and it is why week 7 insists on precision and recall rather than accuracy.

Distributions you will meet

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
samples = {
'normal ': rng.normal(50, 10, 5000),
'uniform ': rng.uniform(0, 100, 5000),
'poisson ': rng.poisson(1.5, 5000),
'lognormal': rng.lognormal(3, 0.6, 5000),
}
print('%-10s %8s %8s %8s %8s' % ('', 'mean', 'median', 'std', 'skew'))
for name, s in samples.items():
print('%-10s %8.2f %8.2f %8.2f %8.2f'
% (name, s.mean(), np.median(s), s.std(), stats.skew(s)))
mean median std skew
normal 49.95 49.72 9.95 0.02
uniform 50.03 50.20 28.81 -0.00
poisson 1.54 1.00 1.22 0.73
lognormal 24.00 20.19 15.68 2.24
DistributionModelsSeen in this course
NormalMeasurement error, sums of many effectsLinear regression residuals
BernoulliA single yes/no outcomeThe churn target
PoissonCounts of rare events in a windowsupport_calls
Log-normalQuantities that multiplyCharges, incomes, durations

Independence, and why naive Bayes is called naive

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()

a = (df['contract'] == 'Month-To-Month')
b = (df['payment_method'] == 'Electronic check')

print('P(A) %.4f' % a.mean())
print('P(B) %.4f' % b.mean())
print('P(A)*P(B) %.4f <- what independence would predict' % (a.mean() * b.mean()))
print('P(A and B) %.4f <- what actually happens' % (a & b).mean())

c = (df['contract'] == 'Two Year')
d = (df['tenure_months'] > 40)
print('\nP(C)*P(D) %.4f' % (c.mean() * d.mean()))
print('P(C and D) %.4f <- strongly dependent' % (c & d).mean())
P(A) 0.5533
P(B) 0.3330
P(A)*P(B) 0.1843 <- what independence would predict
P(A and B) 0.1823 <- what actually happens

P(C)*P(D) 0.0327
P(C and D) 0.0850 <- strongly dependent

Contract type and payment method are close to independent, so multiplying their probabilities is nearly right. Contract type and tenure are not, and multiplying is badly wrong. Naive Bayes assumes every feature is independent given the class, usually false, and it often works anyway, which week 5 explains.

Day 4 takeaway

P(A|B) and P(B|A) are different numbers and confusing them is the commonest analytical error there is. Bayes' theorem converts between them. When the base rate is low, even a very accurate test produces mostly false positives, which is the whole reason rare-event classification needs precision and recall.
Week 03 · Day 5 of 7

Likelihood and Loss Functions

Why the loss you choose decides the estimator you get

By 726 words

A loss function is how a model knows it is wrong. Choosing the wrong one is not a tuning mistake. It means optimising for something you do not want.

Squared error, and what it assumes

import numpy as np

y = np.array([10.0, 12.0, 11.0, 13.0, 40.0]) # last one is an outlier

def mse(pred):
return ((y - pred) ** 2).mean()

def mae(pred):
return np.abs(y - pred).mean()

grid = np.linspace(9, 25, 1601)
print('minimises MSE at %.2f (the mean is %.2f)'
% (grid[np.argmin([mse(g) for g in grid])], y.mean()))
print('minimises MAE at %.2f (the median is %.2f)'
% (grid[np.argmin([mae(g) for g in grid])], np.median(y)))
minimises MSE at 17.20 (the mean is 17.20)
minimises MAE at 12.00 (the median is 12.00)

Your loss chooses your estimator

Squared error is minimised by the mean; absolute error by the median. One outlier of 40 drags the MSE answer to 17.2 while the MAE answer stays at 12. If your target has extreme values and you do not want them to dominate, that is an argument about the loss function, not about the data.

Log loss, for probabilities

Accuracy treats a confident wrong answer and a hesitant wrong answer as identical. Log loss does not, and that is why it is what classifiers actually optimise.

import numpy as np

def log_loss_one(y_true, p):
p = np.clip(p, 1e-15, 1 - 1e-15)
return -(y_true * np.log(p) + (1 - y_true) * np.log(1 - p))

print('truth=1, predicted probability -> loss')
for p in [0.99, 0.9, 0.7, 0.5, 0.3, 0.1, 0.01]:
print(' p=%.2f %7.4f' % (p, log_loss_one(1, p)))
truth=1, predicted probability -> loss
p=0.99 0.0101
p=0.90 0.1054
p=0.70 0.3567
p=0.50 0.6931
p=0.30 1.2040
p=0.10 2.3026
p=0.01 4.6052

Confidently wrong is punished without limit

Predicting 0.01 for something that happens costs 4.6, against 0.69 for an honest 0.5. As the prediction approaches zero the loss goes to infinity, which is why the clip is there. A model trained on log loss learns to be uncertain when it should be, and that is what makes its probabilities usable for ranking and thresholding.

Where log loss comes from

It is not arbitrary. It is the negative log of the likelihood, the probability the model assigns to the data actually observed. Maximising that likelihood and minimising log loss are the same operation.

import numpy as np

y = np.array([1, 0, 1, 1, 0])
model_a = np.array([0.9, 0.2, 0.8, 0.7, 0.1]) # sensible
model_b = np.array([0.6, 0.5, 0.5, 0.6, 0.4]) # hedging

def likelihood(y, p):
return np.prod(np.where(y == 1, p, 1 - p))

def neg_log_lik(y, p):
return -np.sum(np.where(y == 1, np.log(p), np.log(1 - p)))

for name, p in [('A (sensible)', model_a), ('B (hedging)', model_b)]:
print('%-14s likelihood %.6f neg log lik %.4f'
% (name, likelihood(y, p), neg_log_lik(y, p)))
A (sensible) likelihood 0.362880 neg log lik 1.0137
B (hedging) likelihood 0.054000 neg log lik 2.9188

Model A assigns the observed outcomes about eight times more probability than model B does, and correspondingly has the lower negative log likelihood. Products of many small probabilities underflow to zero, which is the practical reason everything is done in logs.

Matching the loss to the problem

ProblemLosssklearn scoring
Regression, outliers matterSquared errorneg_mean_squared_error
Regression, outliers are noiseAbsolute errorneg_mean_absolute_error
Regression, relative error mattersSquared log errorneg_mean_squared_log_error
Binary classificationLog lossneg_log_loss
Multi-classCross-entropyneg_log_loss
Ranking quality only-roc_auc
import numpy as np
from sklearn.metrics import (mean_squared_error, mean_absolute_error,
mean_squared_log_error)

truth = np.array([100.0, 1000.0])
over = np.array([110.0, 1010.0]) # +10 on each

print('absolute error identical: %.1f' % mean_absolute_error(truth, over))
print('squared error identical : %.1f' % mean_squared_error(truth, over))
print('squared LOG error : %.6f' % mean_squared_log_error(truth, over))
print('\n10 on 100 is a 10%% error; 10 on 1000 is 1%%.')
print('Only the log version treats those differently.')
absolute error identical: 10.0
squared error identical : 100.0
squared LOG error : 0.004506

10 on 100 is a 10%% error; 10 on 1000 is 1%%.
Only the log version treats those differently.

Day 5 takeaway

The loss function defines what “wrong” means. Squared error gives you the mean and chases outliers; absolute error gives you the median and ignores them. Log loss punishes confident mistakes without limit and is the negative log likelihood, which is why classifiers optimise it. Choose the loss to match the decision the model supports.
Week 03 · Day 6 of 7

Uncertainty, Bias and Variance

Confidence intervals on model scores, and the decomposition measured directly

By 872 words

Two models score 0.81 and 0.83 on your test set. Is the second better, or did it get luckier? Answering that is what this day is for, and it is the difference between a practitioner and someone who runs fit.

Any score computed on a sample is itself uncertain

import numpy as np

rng = np.random.default_rng(0)
true_rate = 0.268

for n in [50, 200, 750, 3000]:
estimates = [rng.binomial(n, true_rate) / n for _ in range(4000)]
lo, hi = np.percentile(estimates, [2.5, 97.5])
print('n=%5d estimate %.3f +/- %.3f 95%% range %.3f to %.3f'
% (n, np.mean(estimates), np.std(estimates), lo, hi))
n= 50 estimate 0.267 +/- 0.063 95% range 0.140 to 0.400
n= 200 estimate 0.268 +/- 0.032 95% range 0.210 to 0.335
n= 750 estimate 0.268 +/- 0.016 95% range 0.236 to 0.300
n= 3000 estimate 0.268 +/- 0.008 95% range 0.251 to 0.284

With a 750-row test set, the size week 1 used, a churn rate measured at 0.268 could plausibly be anywhere from 0.237 to 0.300. Model scores carry the same kind of uncertainty, and two models within a percentage point of each other on 750 rows are not distinguishable.

Bootstrap a confidence interval for a model score

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score

df = pd.read_csv('customers.csv').drop_duplicates().copy()
X = df[['tenure_months', 'monthly_charges', 'support_calls']]
y = df['churned']
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25,
stratify=y, random_state=42)

model = make_pipeline(SimpleImputer(strategy='median'), StandardScaler(),
LogisticRegression(max_iter=1000)).fit(X_tr, y_tr)
proba = model.predict_proba(X_te)[:, 1]
point = roc_auc_score(y_te, proba)

rng = np.random.default_rng(0)
y_arr, p_arr = y_te.to_numpy(), proba
boot = []
for _ in range(2000):
idx = rng.integers(0, len(y_arr), len(y_arr))
if y_arr[idx].sum() in (0, len(idx)):
continue
boot.append(roc_auc_score(y_arr[idx], p_arr[idx]))

lo, hi = np.percentile(boot, [2.5, 97.5])
print('AUC %.4f 95%% CI %.4f to %.4f width %.4f'
% (point, lo, hi, hi - lo))
AUC 0.7764 95% CI 0.7392 to 0.8101 width 0.0708

Report the interval, not just the number

The interval here is about seven points wide. Any competing model whose score falls inside it has not been shown to be better. Quoting a single figure to three decimal places implies a precision the test set cannot support.

The bias-variance decomposition, measured

Bias: Error from the model being too simple to represent the truth. A straight line fitted to a curve is biased no matter how much data you give it.
Variance: Error from the model being too sensitive to the particular rows it was trained on. Refit it on a different sample and you get a noticeably different model.
import numpy as np
from sklearn.tree import DecisionTreeRegressor
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline

def truth(x):
return np.sin(1.5 * x) + 0.3 * x

rng = np.random.default_rng(0)
x_test = np.linspace(0, 6, 120).reshape(-1, 1)
f_test = truth(x_test.ravel())

def decompose(make_model, n_sets=200, n=40):
preds = np.zeros((n_sets, len(x_test)))
for i in range(n_sets):
x = rng.uniform(0, 6, n).reshape(-1, 1)
y = truth(x.ravel()) + rng.normal(0, 0.3, n)
preds[i] = make_model().fit(x, y).predict(x_test)
mean_pred = preds.mean(axis=0)
bias2 = ((mean_pred - f_test) ** 2).mean()
variance = preds.var(axis=0).mean()
return bias2, variance

candidates = {
'line (deg 1) ': lambda: make_pipeline(PolynomialFeatures(1), LinearRegression()),
'poly deg 5 ': lambda: make_pipeline(PolynomialFeatures(5), LinearRegression()),
'poly deg 15 ': lambda: make_pipeline(PolynomialFeatures(15), LinearRegression()),
'tree depth 2 ': lambda: DecisionTreeRegressor(max_depth=2),
'tree unlimited': lambda: DecisionTreeRegressor(),
}
print('%-16s %9s %9s %9s' % ('model', 'bias^2', 'variance', 'total'))
for name, maker in candidates.items():
b, v = decompose(maker)
print('%-16s %9.4f %9.4f %9.4f' % (name, b, v, b + v))
model bias^2 variance total
line (deg 1) 0.4725 0.0318 0.5043
poly deg 5 0.0167 0.0419 0.0586
poly deg 15 7.4814 2335.0026 2342.4839
tree depth 2 0.0842 0.0742 0.1584
tree unlimited 0.0016 0.1016 0.1032

The straight line has high bias and almost no variance: it is wrong in the same way every time. The degree-15 polynomial is the opposite failure. Its variance is over two thousand, because with only 40 noisy points each fit swings somewhere different. Note that its bias is bad too, at 7.48: when the individual fits are that unstable, even their average is nowhere near the truth. The best total error is the degree-5 polynomial, in between, and finding that point is what every regularisation technique in this course is for.

You cannot measure this on real data

The decomposition above needs the true function, which is why it uses a simulation. On real data you never know the truth, so you detect the same thing indirectly: a large gap between training and validation score means variance; both scores being poor means bias. Week 7 turns that into a diagnostic you can actually run.

Day 6 takeaway

Every score from a finite test set is an estimate with an interval around it, and two models inside each other's intervals have not been shown to differ. Total error splits into bias, variance and irreducible noise; simple models carry bias, flexible ones carry variance, and the whole craft is trading one against the other.
Week 03 · Day 7 of 7

Linear Regression From Scratch, Three Ways

Normal equation, gradient descent and scikit-learn, agreeing

By 813 words

Everything this week, on one problem. You will solve the same regression three ways, exactly, iteratively, and with the library, and get the same answer three times.

The problem

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
d = df[['tenure_months', 'support_calls', 'monthly_charges']].dropna()

X = d[['tenure_months', 'support_calls']].to_numpy(dtype=float)
y = d['monthly_charges'].to_numpy(dtype=float)
print('X', X.shape, ' y', y.shape)
print('predicting monthly_charges from tenure and support calls')
X (2835, 2) y (2835,)
predicting monthly_charges from tenure and support calls

Way 1: the exact solution

Squared error has a closed form. Setting the gradient to zero and solving gives the normal equation, and NumPy will do it in one call.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
d = df[['tenure_months', 'support_calls', 'monthly_charges']].dropna()
X = d[['tenure_months', 'support_calls']].to_numpy(dtype=float)
y = d['monthly_charges'].to_numpy(dtype=float)

Xb = np.column_stack([np.ones(len(X)), X]) # a column of 1s for the intercept

# w = (X'X)^-1 X'y, computed stably rather than by inverting
w_exact = np.linalg.lstsq(Xb, y, rcond=None)[0]
print('intercept %8.4f' % w_exact[0])
print('tenure_months %8.4f' % w_exact[1])
print('support_calls %8.4f' % w_exact[2])
intercept 57.3815
tenure_months -0.0037
support_calls 1.6707

Do not invert the matrix yourself

np.linalg.inv(X.T @ X) @ X.T @ y is the textbook formula and it is numerically fragile: when two features are highly correlated the matrix is near-singular and the inverse amplifies floating-point error enormously. lstsq solves the system directly and is stable. This matters on real data, where correlated features are the norm, recall tenure and total charges at 0.83.

Way 2: gradient descent

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
d = df[['tenure_months', 'support_calls', 'monthly_charges']].dropna()
X = d[['tenure_months', 'support_calls']].to_numpy(dtype=float)
y = d['monthly_charges'].to_numpy(dtype=float)

mu, sigma = X.mean(axis=0), X.std(axis=0)
Xs = (X - mu) / sigma
Xb = np.column_stack([np.ones(len(Xs)), Xs])

w = np.zeros(Xb.shape[1])
lr = 0.1
for step in range(1, 3001):
error = Xb @ w - y
w -= lr * 2 * (Xb.T @ error) / len(y)
if step in (1, 100, 1000, 3000):
print('step %5d loss %10.4f' % (step, (error ** 2).mean()))

# Undo the scaling so the coefficients are comparable with way 1.
coef = w[1:] / sigma
intercept = w[0] - (w[1:] * mu / sigma).sum()
print('\nintercept %8.4f' % intercept)
print('tenure_months %8.4f' % coef[0])
print('support_calls %8.4f' % coef[1])
step 1 loss 4080.0173
step 100 loss 537.9793
step 1000 loss 537.9793
step 3000 loss 537.9793

intercept 57.3815
tenure_months -0.0037
support_calls 1.6707

Way 3: the library

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

df = pd.read_csv('customers.csv').drop_duplicates().copy()
d = df[['tenure_months', 'support_calls', 'monthly_charges']].dropna()
X = d[['tenure_months', 'support_calls']].to_numpy(dtype=float)
y = d['monthly_charges'].to_numpy(dtype=float)

m = LinearRegression().fit(X, y)
print('intercept %8.4f' % m.intercept_)
print('tenure_months %8.4f' % m.coef_[0])
print('support_calls %8.4f' % m.coef_[1])
print('\nR^2 %.4f' % m.score(X, y))
intercept 57.3815
tenure_months -0.0037
support_calls 1.6707

R^2 0.0233

Three routes, one answer. The library is not doing anything you have not now done by hand. It is doing it more carefully, and handling the cases where the naive version breaks.

Reading the coefficients

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

df = pd.read_csv('customers.csv').drop_duplicates().copy()
d = df[['tenure_months', 'support_calls', 'monthly_charges']].dropna()
X = d[['tenure_months', 'support_calls']].to_numpy(dtype=float)
y = d['monthly_charges'].to_numpy(dtype=float)
m = LinearRegression().fit(X, y)

print('a customer with 12 months tenure and 2 support calls:')
print(' predicted charge %.2f' % m.predict([[12, 2]])[0])
print('one more support call, everything else held:')
print(' predicted charge %.2f' % m.predict([[12, 3]])[0])
print(' difference %.4f <- exactly the coefficient' % m.coef_[1])
a customer with 12 months tenure and 2 support calls:
predicted charge 60.68
one more support call, everything else held:
predicted charge 62.35
difference 1.6707 <- exactly the coefficient

"Holding everything else constant" is doing a lot of work

A coefficient is the change in the prediction for a one-unit change in that feature with the others fixed. When features are correlated you cannot actually hold the others fixed, increasing tenure increases total charges in the real world, so the coefficient describes the model, not the world. This is the single biggest source of over-claiming in applied work.

Your assignment

Add total_charges as a third feature and refit. It correlates with tenure at 0.83. Record what happens to the tenure coefficient, sign, size, or both. Then refit on a random 60 percent of the rows, twice, and compare. Week 4 names what you are seeing and fixes it.

Day 7 takeaway

Linear regression has an exact solution, an iterative one, and a library one, and they agree. Use lstsq rather than inverting a matrix. A coefficient means “change in prediction per unit, others held constant”, a claim that gets shakier the more your features correlate.