Classical machine learning all 5 written

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56# Twenty points from a line, with noise added. The noise is the entire difficulty.
import random

TRUE_W, TRUE_B = 3.0, -2.0
rng = random.Random(42)             # seeded, so these numbers are the same every build

XS = [i * 0.5 for i in range(20)]
NOISE = [rng.gauss(0.0, 1.5) for _ in XS]
YS = [TRUE_W * x + TRUE_B + e for x, e in zip(XS, NOISE)]


def predict(w, b, x):
    return w * x + b


def mse(w, b):
    """Mean squared error. One number for how wrong the whole line is."""
    n = len(XS)
    return sum((predict(w, b, x) - y) ** 2 for x, y in zip(XS, YS)) / n


def gradient(w, b):
    """dL/dw and dL/db, by the chain rule, exactly as in lesson 7."""
    n = len(XS)
    errors = [predict(w, b, x) - y for x, y in zip(XS, YS)]
    dw = sum(2.0 * e * x for e, x in zip(errors, XS)) / n
    db = sum(2.0 * e for e in errors) / n
    return dw, db


print("     x         y      noiseless      noise")
for i in (0, 1, 2, 3, 18, 19):
    line = TRUE_W * XS[i] + TRUE_B
    print("  {:4.1f}  {:9.4f}   {:9.4f}   {:+8.4f}".format(XS[i], YS[i], line, YS[i] - line))

honest = mse(TRUE_W, TRUE_B)
print()
print("loss at the parameters that GENERATED the data: {:.6f}".format(honest))
print("the noise variance we asked for               : {:.6f}".format(1.5 ** 2))
print("Twenty draws is a small sample. This one came out quieter than requested.")

# Now the part worth sitting with. If w = 3 and b = -2 were the best answer, the
# gradient there would be zero. It is not.
dw, db = gradient(TRUE_W, TRUE_B)
print()
print("gradient there: dL/dw = {:+.6f}   dL/db = {:+.6f}".format(dw, db))
nudged = mse(TRUE_W - 0.01 * dw, TRUE_B - 0.01 * db)
print("one small step downhill from them: {:.6f}".format(nudged))

test("the noise moved every single point off the line",
     all(not close(y, TRUE_W * x + TRUE_B) for x, y in zip(XS, YS)), True)
test("so the loss at the generating parameters is not zero", honest > 0.0, True)
test("the loss there is exactly the mean squared noise",
     close(honest, sum(e * e for e in NOISE) / len(NOISE)), True)
test("the gradient there is not zero either", [close(dw, 0.0), close(db, 0.0)], [False, False])
test("one step downhill from the true parameters lowers the loss", nudged < honest, True)

     x         y      noiseless      noise
   0.0    -2.2161     -2.0000    -0.2161
   0.5    -0.7594     -0.5000    -0.2594
   1.0     0.8330      1.0000    -0.1670
   1.5     3.5530      2.5000    +1.0530
   9.0    25.0625     25.0000    +0.0625
   9.5    26.3405     26.5000    -0.1595

loss at the parameters that GENERATED the data: 0.930950
the noise variance we asked for               : 2.250000
Twenty draws is a small sample. This one came out quieter than requested.

gradient there: dL/dw = -1.273844   dL/db = -0.077312
one small step downhill from them: 0.919768
ok   the noise moved every single point off the line
ok   so the loss at the generating parameters is not zero
ok   the loss there is exactly the mean squared noise
ok   the gradient there is not zero either
ok   one step downhill from the true parameters lowers the loss

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62# Gradient descent on noisy data. Note the learning rate, and why it is small.
import random

TRUE_W, TRUE_B = 3.0, -2.0
rng = random.Random(42)
XS = [i * 0.5 for i in range(20)]
YS = [TRUE_W * x + TRUE_B + rng.gauss(0.0, 1.5) for x in XS]


def mse(w, b):
    n = len(XS)
    return sum((w * x + b - y) ** 2 for x, y in zip(XS, YS)) / n


def gradient(w, b):
    n = len(XS)
    errors = [w * x + b - y for x, y in zip(XS, YS)]
    return (sum(2.0 * e * x for e, x in zip(errors, XS)) / n,
            sum(2.0 * e for e in errors) / n)


def descend(lr, steps, w=0.0, b=0.0):
    trail = [(w, b, mse(w, b))]
    for _ in range(steps):
        dw, db = gradient(w, b)
        w -= lr * dw
        b -= lr * db
        trail.append((w, b, mse(w, b)))
    return trail


trail = descend(lr=0.01, steps=5000)
losses = [row[2] for row in trail]

print("  step            w            b           loss")
for step in (0, 1, 2, 5, 20, 100, 500, 1000, 2500, 5000):
    w, b, value = trail[step]
    print("  {:5}   {:10.6f}   {:10.6f}   {:12.8f}".format(step, w, b, value))

final_w, final_b, final_loss = trail[-1]
print()
print("finished at w = {:.6f}, b = {:.6f}".format(final_w, final_b))
print("true values were w = {}, b = {}".format(TRUE_W, TRUE_B))
print("loss at the true values: {:.8f}".format(mse(TRUE_W, TRUE_B)))
dw, db = gradient(final_w, final_b)
print("gradient at the end: {:+.3e}  {:+.3e}".format(dw, db))

# The slope is found in twenty steps. The intercept takes thousands, because only
# the leftmost rows say much about it, and it has to creep across the whole range.
strict = next(i for i in range(len(losses) - 1) if not losses[i + 1] < losses[i])
plateau = max(losses[3000:]) - min(losses[3000:])
print()
print("the loss fell strictly for {} steps, then only jittered by {:.1e}".format(strict, plateau))

test("the loss fell on every one of the first two thousand steps",
     all(losses[i + 1] < losses[i] for i in range(2000)), True)
test("after that it just jitters in the last digits of the float", plateau < 1e-12, True)
test("descent found a lower loss than the true parameters give", final_loss < mse(TRUE_W, TRUE_B), True)
test("the slope came back within a fifth of three", close(final_w, TRUE_W, 0.2), True)
test("the intercept came back within one of minus two", close(final_b, TRUE_B, 1.0), True)
test("the gradient at the end is flat to six decimal places",
     [close(dw, 0.0, 1e-6), close(db, 0.0, 1e-6)], [True, True])

  step            w            b           loss
      0     0.000000     0.000000   229.47285775
      1     1.675238     0.245773    33.41019432
      2     2.292669     0.327483     6.87665466
      5     2.639109     0.339397     2.68327603
     20     2.687148     0.147286     2.41429802
    100     2.813587    -0.667565     1.55404107
    500     3.025297    -2.031962     0.91429515
   1000     3.052439    -2.206883     0.90478459
   2500     3.054532    -2.220371     0.90473556
   5000     3.054533    -2.220376     0.90473556

finished at w = 3.054533, b = -2.220376
true values were w = 3.0, b = -2.0
loss at the true values: 0.93095004
gradient at the end: -7.645e-13  +4.909e-12

the loss fell strictly for 2912 steps, then only jittered by 3.5e-14
ok   the loss fell on every one of the first two thousand steps
ok   after that it just jitters in the last digits of the float
ok   descent found a lower loss than the true parameters give
ok   the slope came back within a fifth of three
ok   the intercept came back within one of minus two
ok   the gradient at the end is flat to six decimal places

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81# The exact answer, two ways, checked against five thousand steps of descent.
import random

TRUE_W, TRUE_B = 3.0, -2.0
rng = random.Random(42)
XS = [i * 0.5 for i in range(20)]
YS = [TRUE_W * x + TRUE_B + rng.gauss(0.0, 1.5) for x in XS]


def fit_covariance(xs, ys):
    """The textbook form: w is covariance over variance, b puts the line through the means."""
    n = len(xs)
    xbar = sum(xs) / n
    ybar = sum(ys) / n
    sxx = sum((x - xbar) ** 2 for x in xs)
    if sxx < 1e-12:
        raise ValueError("every x is identical, so no single line is determined")
    sxy = sum((x - xbar) * (y - ybar) for x, y in zip(xs, ys))
    w = sxy / sxx
    return w, ybar - w * xbar


def fit_normal_equations(xs, ys):
    """The same thing as a 2x2 linear system, solved by Cramer's rule.

        [ sum(x*x)  sum(x) ] [ w ]   [ sum(x*y) ]
        [ sum(x)         n ] [ b ] = [ sum(y)   ]
    """
    n = float(len(xs))
    sx = sum(xs)
    sxx = sum(x * x for x in xs)
    sy = sum(ys)
    sxy = sum(x * y for x, y in zip(xs, ys))
    det = sxx * n - sx * sx
    if abs(det) < 1e-12:
        raise ValueError("the normal equations are singular")
    return (sxy * n - sx * sy) / det, (sxx * sy - sx * sxy) / det


w_cov, b_cov = fit_covariance(XS, YS)
w_ne, b_ne = fit_normal_equations(XS, YS)

w_gd, b_gd = 0.0, 0.0
for _ in range(5000):
    n = len(XS)
    errors = [w_gd * x + b_gd - y for x, y in zip(XS, YS)]
    w_gd -= 0.01 * sum(2.0 * e * x for e, x in zip(errors, XS)) / n
    b_gd -= 0.01 * sum(2.0 * e for e in errors) / n

print("covariance form   : w = {:.12f}   b = {:.12f}".format(w_cov, b_cov))
print("normal equations  : w = {:.12f}   b = {:.12f}".format(w_ne, b_ne))
print("gradient descent  : w = {:.12f}   b = {:.12f}".format(w_gd, b_gd))
print()
print("closed form == gradient descent?  ==  says", w_cov == w_gd)
print("                                  close says", close(w_cov, w_gd, 1e-6))

# Two properties every least-squares fit with an intercept has, worth checking
# rather than believing: the residuals cancel, and the line goes through the means.
residuals = [y - (w_cov * x + b_cov) for x, y in zip(XS, YS)]
xbar = sum(XS) / len(XS)
ybar = sum(YS) / len(YS)
print()
print("sum of residuals at the fit: {:.3e}".format(sum(residuals)))
print("the fit at xbar = {:.3f} is {:.9f}, and ybar is {:.9f}".format(
    xbar, w_cov * xbar + b_cov, ybar))

try:
    fit_covariance([2.0, 2.0, 2.0], [1.0, 5.0, 9.0])
    refused = False
except ValueError as error:
    refused = True
    print()
    print("a vertical cloud of points is refused:", error)

test("the two closed forms agree", [close(w_cov, w_ne), close(b_cov, b_ne)], [True, True])
test("five thousand descent steps reached the same answer",
     [close(w_cov, w_gd, 1e-6), close(b_cov, b_gd, 1e-6)], [True, True])
test("but not bit for bit, which is why close() exists", w_cov == w_gd, False)
test("the residuals of a least-squares fit sum to zero", close(sum(residuals), 0.0, 1e-9), True)
test("the fitted line passes through the mean of the data", close(w_cov * xbar + b_cov, ybar), True)
test("no line can be fitted when every x is the same", refused, True)

covariance form   : w = 3.054533098252   b = -2.220376321641
normal equations  : w = 3.054533098252   b = -2.220376321641
gradient descent  : w = 3.054533098251   b = -2.220376321632

closed form == gradient descent?  ==  says False
                                  close says True

sum of residuals at the fit: -2.132e-14
the fit at xbar = 4.750 is 12.288655895, and ybar is 12.288655895

a vertical cloud of points is refused: every x is identical, so no single line is determined
ok   the two closed forms agree
ok   five thousand descent steps reached the same answer
ok   but not bit for bit, which is why close() exists
ok   the residuals of a least-squares fit sum to zero
ok   the fitted line passes through the mean of the data
ok   no line can be fitted when every x is the same

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59# R squared: the loss, rescaled against the laziest possible model.
import math
import random

TRUE_W, TRUE_B = 3.0, -2.0
rng = random.Random(42)
XS = [i * 0.5 for i in range(20)]
YS = [TRUE_W * x + TRUE_B + rng.gauss(0.0, 1.5) for x in XS]


def fit(xs, ys):
    n = len(xs)
    xbar = sum(xs) / n
    ybar = sum(ys) / n
    sxx = sum((x - xbar) ** 2 for x in xs)
    sxy = sum((x - xbar) * (y - ybar) for x, y in zip(xs, ys))
    w = sxy / sxx
    return w, ybar - w * xbar


def r_squared(xs, ys, w, b):
    """1 - (what the line got wrong) / (what predicting the mean gets wrong)."""
    ybar = sum(ys) / len(ys)
    rss = sum((y - (w * x + b)) ** 2 for x, y in zip(xs, ys))
    tss = sum((y - ybar) ** 2 for y in ys)
    return 1.0 - rss / tss


w, b = fit(XS, YS)
ybar = sum(YS) / len(YS)
clean_y = [TRUE_W * x + TRUE_B for x in XS]
cw, cb = fit(XS, clean_y)

print("R2 of the fitted line          : {:.6f}".format(r_squared(XS, YS, w, b)))
print("R2 of always predicting ybar   : {:.6f}   <- zero, by construction".format(
    r_squared(XS, YS, 0.0, ybar)))
print("R2 of a deliberately bad line  : {:.6f}   <- worse than the mean".format(
    r_squared(XS, YS, -1.0, 10.0)))
print("R2 on noiseless data           : {:.6f}".format(r_squared(XS, clean_y, cw, cb)))

# For a single input there is an identity worth knowing: R2 is the correlation squared.
xbar = sum(XS) / len(XS)
sxx = sum((x - xbar) ** 2 for x in XS)
syy = sum((y - ybar) ** 2 for y in YS)
sxy = sum((x - xbar) * (y - ybar) for x, y in zip(XS, YS))
corr = sxy / math.sqrt(sxx * syy)
print()
print("correlation between x and y : {:.9f}".format(corr))
print("that squared                : {:.9f}".format(corr * corr))
print("R2 again                    : {:.9f}".format(r_squared(XS, YS, w, b)))

test("the fit explains most of the variance", 0.9 < r_squared(XS, YS, w, b) < 1.0, True)
test("predicting the mean every time scores exactly zero",
     close(r_squared(XS, YS, 0.0, ybar), 0.0), True)
test("a bad enough line scores below zero", r_squared(XS, YS, -1.0, 10.0) < 0.0, True)
test("a perfect fit on noiseless data scores one", close(r_squared(XS, clean_y, cw, cb), 1.0), True)
test("for one input, R2 is the squared correlation", close(corr * corr, r_squared(XS, YS, w, b)), True)
test("not even the true parameters beat the least-squares line on this sample",
     r_squared(XS, YS, w, b) > r_squared(XS, YS, TRUE_W, TRUE_B), True)

R2 of the fitted line          : 0.988469
R2 of always predicting ybar   : 0.000000   <- zero, by construction
R2 of a deliberately bad line  : -1.384582   <- worse than the mean
R2 on noiseless data           : 1.000000

correlation between x and y : 0.994217830
that squared                : 0.988469094
R2 again                    : 0.988469094
ok   the fit explains most of the variance
ok   predicting the mean every time scores exactly zero
ok   a bad enough line scores below zero
ok   a perfect fit on noiseless data scores one
ok   for one input, R2 is the squared correlation
ok   not even the true parameters beat the least-squares line on this sample

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59# One value changed out of twenty, and what it does to the line.
import random

TRUE_W, TRUE_B = 3.0, -2.0
rng = random.Random(42)
XS = [i * 0.5 for i in range(20)]
YS = [TRUE_W * x + TRUE_B + rng.gauss(0.0, 1.5) for x in XS]


def fit(xs, ys):
    n = len(xs)
    xbar = sum(xs) / n
    ybar = sum(ys) / n
    sxx = sum((x - xbar) ** 2 for x in xs)
    sxy = sum((x - xbar) * (y - ybar) for x, y in zip(xs, ys))
    w = sxy / sxx
    return w, ybar - w * xbar


def r_squared(xs, ys, w, b):
    ybar = sum(ys) / len(ys)
    rss = sum((y - (w * x + b)) ** 2 for x, y in zip(xs, ys))
    tss = sum((y - ybar) ** 2 for y in ys)
    return 1.0 - rss / tss


BROKEN = list(YS)
BROKEN[3] = YS[3] + 30.0            # one typo, one row, near the left of the range

w_ok, b_ok = fit(XS, YS)
w_bad, b_bad = fit(XS, BROKEN)

print("               slope     intercept        R2")
print("clean       {:8.4f}   {:9.4f}   {:8.4f}".format(w_ok, b_ok, r_squared(XS, YS, w_ok, b_ok)))
print("one typo    {:8.4f}   {:9.4f}   {:8.4f}".format(
    w_bad, b_bad, r_squared(XS, BROKEN, w_bad, b_bad)))
print()
print("One row in twenty, and the slope moved {:.4f} while R2 fell by half.".format(
    abs(w_bad - w_ok)))

# The residuals are how you find it. The offender is not subtle once you look.
residuals = [abs(y - (w_bad * x + b_bad)) for x, y in zip(XS, BROKEN)]
worst = residuals.index(max(residuals))
print()
print("largest residual is at row {}, size {:.3f}".format(worst, residuals[worst]))
print("second largest is           {:.3f}".format(sorted(residuals)[-2]))

keep = [i for i in range(len(XS)) if i != worst]
w_trim, b_trim = fit([XS[i] for i in keep], [BROKEN[i] for i in keep])
print("refit without that row: slope {:.4f}, against {:.4f} on the clean data".format(w_trim, w_ok))

test("one bad row in twenty moved the slope by more than a tenth", abs(w_bad - w_ok) > 0.1, True)
test("and dragged the intercept the other way to compensate", b_bad > b_ok, True)
test("R2 fell once the outlier was in",
     r_squared(XS, BROKEN, w_bad, b_bad) < r_squared(XS, YS, w_ok, b_ok), True)
test("the largest residual points straight at the damaged row", worst, 3)
test("it is more than twice the next worst, so it stands out",
     residuals[worst] > 2.0 * sorted(residuals)[-2], True)
test("dropping it recovers the original slope", close(w_trim, w_ok, 0.05), True)

               slope     intercept        R2
clean         3.0545     -2.2204     0.9885
one typo      2.4681      2.0653     0.5330

One row in twenty, and the slope moved 0.5865 while R2 fell by half.

largest residual is at row 3, size 27.786
second largest is           4.982
refit without that row: slope 3.0808, against 3.0545 on the clean data
ok   one bad row in twenty moved the slope by more than a tenth
ok   and dragged the intercept the other way to compensate
ok   R2 fell once the outlier was in
ok   the largest residual points straight at the damaged row
ok   it is more than twice the next worst, so it stands out
ok   dropping it recovers the original slope

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55# The sigmoid. One line is not enough.
import math


def sigmoid_naive(z):
    return 1.0 / (1.0 + math.exp(-z))


def sigmoid(z):
    """Two branches, so the exponent is never large and positive."""
    if z >= 0.0:
        return 1.0 / (1.0 + math.exp(-z))
    e = math.exp(z)
    return e / (1.0 + e)


print("      z     sigmoid(z)")
for z in (-8.0, -2.0, -0.5, 0.0, 0.5, 2.0, 8.0):
    print("  {:5.1f}   {:.9f}".format(z, sigmoid(z)))

print()
print("sigmoid(  40) =", sigmoid(40.0))
print("sigmoid( 800) =", sigmoid(800.0), " <- exactly 1.0, not nearly")
print("sigmoid(-800) =", sigmoid(-800.0), " <- exactly 0.0")

# A score of -800 is not exotic. Give a trained model an outlier, or let the
# weights grow for a few thousand steps, and exp(800) is what the naive line asks for.
try:
    sigmoid_naive(-800.0)
    overflowed = False
except OverflowError as error:
    overflowed = True
    print()
    print("the one-line version raises here: OverflowError:", error)
print("the two-branch version returns   :", sigmoid(-800.0))

# The derivative is unusually tidy, and worth checking rather than trusting.
d = sigmoid(1.3) * (1.0 - sigmoid(1.3))
numeric = (sigmoid(1.3 + 1e-5) - sigmoid(1.3 - 1e-5)) / 2e-5
print()
print("s'(1.3) claimed as s(1 - s)  : {:.12f}".format(d))
print("s'(1.3) by central difference : {:.12f}".format(numeric))

test("sigmoid of zero is a half", close(sigmoid(0.0), 0.5), True)
test("it is symmetric, so one minus s(z) is s(-z)",
     [close(1.0 - sigmoid(2.3), sigmoid(-2.3)), close(1.0 - sigmoid(-7.0), sigmoid(7.0))],
     [True, True])
test("it never leaves zero to one",
     all(0.0 <= sigmoid(z) <= 1.0 for z in range(-900, 900, 37)), True)
test("it only ever increases",
     all(sigmoid(z / 10.0) < sigmoid((z + 1) / 10.0) for z in range(-60, 60)), True)
test("the naive one-liner overflows on a large negative input", overflowed, True)
test("the two-branch version saturates to exactly zero and one instead",
     [sigmoid(-800.0), sigmoid(800.0)], [0.0, 1.0])
test("its derivative really is s times one minus s", close(d, numeric, 1e-9), True)

      z     sigmoid(z)
   -8.0   0.000335350
   -2.0   0.119202922
   -0.5   0.377540669
    0.0   0.500000000
    0.5   0.622459331
    2.0   0.880797078
    8.0   0.999664650

sigmoid(  40) = 1.0
sigmoid( 800) = 1.0  <- exactly 1.0, not nearly
sigmoid(-800) = 0.0  <- exactly 0.0

the one-line version raises here: OverflowError: math range error
the two-branch version returns   : 0.0

s'(1.3) claimed as s(1 - s)  : 0.168298362469
s'(1.3) by central difference : 0.168298362474
ok   sigmoid of zero is a half
ok   it is symmetric, so one minus s(z) is s(-z)
ok   it never leaves zero to one
ok   it only ever increases
ok   the naive one-liner overflows on a large negative input
ok   the two-branch version saturates to exactly zero and one instead
ok   its derivative really is s times one minus s

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58# Log loss, the clip it needs, and the reason squared error is not used here.
import math


def sigmoid(z):
    if z >= 0.0:
        return 1.0 / (1.0 + math.exp(-z))
    e = math.exp(z)
    return e / (1.0 + e)


def log_loss(p, y):
    """How surprised we were by the truth. Clipped, or log(0.0) raises."""
    p = min(max(p, 1e-15), 1.0 - 1e-15)
    return -(y * math.log(p) + (1.0 - y) * math.log(1.0 - p))


def log_loss_unclipped(p, y):
    return -(y * math.log(p) + (1.0 - y) * math.log(1.0 - p))


print("  label   predicted p     log loss")
for y, p in ((1, 0.99), (1, 0.75), (1, 0.5), (1, 0.25), (1, 0.01), (0, 0.01), (0, 0.99)):
    print("    {}       {:6.2f}      {:9.4f}".format(y, p, log_loss(p, y)))

saturated = sigmoid(-800.0)
print()
print("a very confident wrong prediction gives p =", saturated)
try:
    log_loss_unclipped(saturated, 1)
    raised = False
except ValueError as error:
    raised = True
    print("log loss without clipping raises: ValueError:", error)
print("with the 1e-15 clip it returns  : {:.6f}".format(log_loss(saturated, 1)))
print("which is -log(1e-15)            : {:.6f}".format(-math.log(1e-15)))

# Now the question people ask: why not just square the error, as in lesson 8?
# Because of what happens to the gradient when the model is badly wrong.
z = -20.0
p = sigmoid(z)
d_log = p - 1.0                                  # d/dz of log loss, with y = 1
d_squared = 2.0 * (p - 1.0) * p * (1.0 - p)      # d/dz of (p - y)^2, with y = 1
print()
print("at z = -20 with the true label 1, the model is confidently wrong.")
print("  p                     = {:.3e}".format(p))
print("  d(log loss)/dz        = {:+.9f}".format(d_log))
print("  d(squared error)/dz   = {:+.3e}".format(d_squared))
print("Squared error has no gradient left. Log loss still pushes at full strength.")

test("a confident correct prediction costs almost nothing", log_loss(0.999999, 1) < 1e-5, True)
test("a confident wrong prediction costs a great deal", log_loss(1e-9, 1) > 20.0, True)
test("without clipping, a saturated probability makes log() raise", raised, True)
test("the clip caps the worst possible loss at about 34.5",
     close(log_loss(saturated, 1), -math.log(1e-15), 1e-9), True)
test("log loss still has a full-strength slope when confidently wrong",
     close(d_log, -1.0, 1e-8), True)
test("squared error has none left, so it would stop learning", abs(d_squared) < 1e-8, True)

  label   predicted p     log loss
    1         0.99         0.0101
    1         0.75         0.2877
    1         0.50         0.6931
    1         0.25         1.3863
    1         0.01         4.6052
    0         0.01         0.0101
    0         0.99         4.6052

a very confident wrong prediction gives p = 0.0
log loss without clipping raises: ValueError: math domain error
with the 1e-15 clip it returns  : 34.538776
which is -log(1e-15)            : 34.538776

at z = -20 with the true label 1, the model is confidently wrong.
  p                     = 2.061e-09
  d(log loss)/dz        = -0.999999998
  d(squared error)/dz   = -4.122e-09
Squared error has no gradient left. Log loss still pushes at full strength.
ok   a confident correct prediction costs almost nothing
ok   a confident wrong prediction costs a great deal
ok   without clipping, a saturated probability makes log() raise
ok   the clip caps the worst possible loss at about 34.5
ok   log loss still has a full-strength slope when confidently wrong
ok   squared error has none left, so it would stop learning

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79# Logistic regression by gradient descent. The loop is unchanged from lesson 8.
import math
import random


def sigmoid(z):
    if z >= 0.0:
        return 1.0 / (1.0 + math.exp(-z))
    e = math.exp(z)
    return e / (1.0 + e)


def log_loss(p, y):
    p = min(max(p, 1e-15), 1.0 - 1e-15)
    return -(y * math.log(p) + (1.0 - y) * math.log(1.0 - p))


rng = random.Random(42)
ROWS = []
for _ in range(30):
    ROWS.append(([rng.gauss(2.0, 0.8), rng.gauss(2.0, 0.8)], 0))
for _ in range(30):
    ROWS.append(([rng.gauss(4.0, 0.8), rng.gauss(4.0, 0.8)], 1))


def probability(w, b, x):
    return sigmoid(sum(wi * xi for wi, xi in zip(w, x)) + b)


def mean_loss(w, b):
    return sum(log_loss(probability(w, b, x), y) for x, y in ROWS) / len(ROWS)


def gradient(w, b):
    """dL/dw is mean((p - y) * x), dL/db is mean(p - y). The sigmoid cancels out."""
    n = len(ROWS)
    dw = [0.0] * len(w)
    db = 0.0
    for x, y in ROWS:
        e = probability(w, b, x) - y
        for j in range(len(w)):
            dw[j] += e * x[j]
        db += e
    return [g / n for g in dw], db / n


def accuracy(w, b):
    right = sum(1 for x, y in ROWS if (1 if probability(w, b, x) >= 0.5 else 0) == y)
    return right / len(ROWS)


w = [0.0, 0.0]
b = 0.0
history = [mean_loss(w, b)]
print("  step        w1        w2         b       loss    accuracy")
for step in range(1, 4001):
    dw, db = gradient(w, b)
    w = [wi - 0.2 * g for wi, g in zip(w, dw)]
    b -= 0.2 * db
    history.append(mean_loss(w, b))
    if step in (1, 2, 10, 50, 200, 1000, 4000):
        print("  {:5}  {:8.4f}  {:8.4f}  {:8.4f}   {:8.5f}    {:.4f}".format(
            step, w[0], w[1], b, history[-1], accuracy(w, b)))

print()
print("starting loss was {:.6f}, which is log(2) = {:.6f}".format(history[0], math.log(2.0)))
print("final weights: w =", [round(v, 4) for v in w], " b = {:.4f}".format(b))
print("final loss {:.6f}, accuracy {:.4f}".format(history[-1], accuracy(w, b)))

test("with all parameters at zero every prediction is a coin flip",
     close(history[0], math.log(2.0)), True)
test("the loss fell on every one of the four thousand steps",
     all(history[i + 1] < history[i] for i in range(len(history) - 1)), True)
test("both weights came out positive, since both features rise with the class",
     [w[0] > 0.0, w[1] > 0.0], [True, True])
test("the bias came out negative to balance them", b < 0.0, True)
test("it separates the training data almost perfectly", accuracy(w, b) >= 0.95, True)
test("every predicted probability is a real probability",
     all(0.0 <= probability(w, b, x) <= 1.0 for x, _ in ROWS), True)

  step        w1        w2         b       loss    accuracy
      1    0.1063    0.1011    0.0000    0.64122    0.5000
      2    0.1063    0.0994   -0.0302    0.63678    0.5000
     10    0.1454    0.1258   -0.2543    0.60410    0.5167
     50    0.3129    0.2521   -1.2346    0.47830    0.8000
    200    0.7213    0.5811   -3.6053    0.27556    0.9833
   1000    1.5548    1.2529   -8.2008    0.11718    0.9833
   4000    2.5242    2.0416  -13.3321    0.06527    0.9833

starting loss was 0.693147, which is log(2) = 0.693147
final weights: w = [2.5242, 2.0416]  b = -13.3321
final loss 0.065266, accuracy 0.9833
ok   with all parameters at zero every prediction is a coin flip
ok   the loss fell on every one of the four thousand steps
ok   both weights came out positive, since both features rise with the class
ok   the bias came out negative to balance them
ok   it separates the training data almost perfectly
ok   every predicted probability is a real probability

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86# The boundary is a line. Here it is, with the probability checked along it.
import math
import random


def sigmoid(z):
    if z >= 0.0:
        return 1.0 / (1.0 + math.exp(-z))
    e = math.exp(z)
    return e / (1.0 + e)


rng = random.Random(42)
ROWS = []
for _ in range(30):
    ROWS.append(([rng.gauss(2.0, 0.8), rng.gauss(2.0, 0.8)], 0))
for _ in range(30):
    ROWS.append(([rng.gauss(4.0, 0.8), rng.gauss(4.0, 0.8)], 1))


def fit(rows, lr=0.2, steps=4000):
    w = [0.0, 0.0]
    b = 0.0
    n = len(rows)
    for _ in range(steps):
        dw = [0.0, 0.0]
        db = 0.0
        for x, y in rows:
            e = sigmoid(w[0] * x[0] + w[1] * x[1] + b) - y
            dw[0] += e * x[0]
            dw[1] += e * x[1]
            db += e
        w = [w[0] - lr * dw[0] / n, w[1] - lr * dw[1] / n]
        b -= lr * db / n
    return w, b


w, b = fit(ROWS)


def score(x):
    """The number before the squashing. Positive means class 1."""
    return w[0] * x[0] + w[1] * x[1] + b


print("the model is p = sigmoid({:+.4f} * x1 {:+.4f} * x2 {:+.4f})".format(w[0], w[1], b))
slope = -w[0] / w[1]
intercept = -b / w[1]
print()
print("set that inner expression to zero and solve for x2:")
print("    x2 = {:+.6f} * x1 {:+.6f}".format(slope, intercept))
print()
print("      x1        x2      score          p")
for x1 in (1.0, 2.0, 3.0, 4.0, 5.0):
    x2 = slope * x1 + intercept
    print("  {:6.2f}  {:8.4f}   {:+8.1e}   {:.9f}".format(
        x1, x2, score([x1, x2]), sigmoid(score([x1, x2]))))

mean0 = [sum(x[j] for x, y in ROWS if y == 0) / 30 for j in (0, 1)]
mean1 = [sum(x[j] for x, y in ROWS if y == 1) / 30 for j in (0, 1)]
print()
print("class 0 mean {} scores {:+.4f}".format([round(v, 3) for v in mean0], score(mean0)))
print("class 1 mean {} scores {:+.4f}".format([round(v, 3) for v in mean1], score(mean1)))

# Dividing the score by the length of w turns it into a signed distance to the line,
# which is the useful way to read "how far into the wrong side was this row".
length = math.sqrt(w[0] ** 2 + w[1] ** 2)
wrong = [(x, y, sigmoid(score(x))) for x, y in ROWS
         if (1 if sigmoid(score(x)) >= 0.5 else 0) != y]
print()
print("misclassified rows:", len(wrong))
for x, y, p in wrong:
    print("   label {} at {} got p = {:.4f}, distance to the boundary {:+.4f}".format(
        y, [round(v, 3) for v in x], p, score(x) / length))

test("every point on the boundary line has probability one half",
     all(close(sigmoid(score([x1, slope * x1 + intercept])), 0.5, 1e-12)
         for x1 in (1.0, 2.0, 3.0, 4.0, 5.0)), True)
test("the boundary really is straight: the slope is the same everywhere",
     close((slope * 5.0 + intercept - (slope * 1.0 + intercept)) / 4.0, slope), True)
test("the two class means fall on opposite sides of it",
     [score(mean0) < 0.0, score(mean1) > 0.0], [True, True])
test("a point far into class 1 territory is called with confidence",
     sigmoid(score([8.0, 8.0])) > 0.999, True)
test("only one row is misclassified, and it sits nearest the boundary",
     [len(wrong), all(abs(p - 0.5) < 0.3 for _, _, p in wrong)], [1, True])

the model is p = sigmoid(+2.5242 * x1 +2.0416 * x2 -13.3321)

set that inner expression to zero and solve for x2:
    x2 = -1.236398 * x1 +6.530318

      x1        x2      score          p
    1.00    5.2939   +1.8e-15   0.500000000
    2.00    4.0575   +0.0e+00   0.500000000
    3.00    2.8211   +0.0e+00   0.500000000
    4.00    1.5847   +0.0e+00   0.500000000
    5.00    0.3483   +0.0e+00   0.500000000

class 0 mean [2.037, 1.985] scores -4.1375
class 1 mean [4.163, 4.006] scores +5.3546

misclassified rows: 1
   label 1 at [2.527, 3.134] got p = 0.3647, distance to the boundary -0.1710
ok   every point on the boundary line has probability one half
ok   the boundary really is straight: the slope is the same everywhere
ok   the two class means fall on opposite sides of it
ok   a point far into class 1 territory is called with confidence
ok   only one row is misclassified, and it sits nearest the boundary

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93# Exclusive-or: four points, and a wall.
import math


def sigmoid(z):
    if z >= 0.0:
        return 1.0 / (1.0 + math.exp(-z))
    e = math.exp(z)
    return e / (1.0 + e)


def log_loss(p, y):
    p = min(max(p, 1e-15), 1.0 - 1e-15)
    return -(y * math.log(p) + (1.0 - y) * math.log(1.0 - p))


XOR = [([-1.0, -1.0], 0), ([-1.0, 1.0], 1), ([1.0, -1.0], 1), ([1.0, 1.0], 0)]


def plain(x):
    return list(x)


def with_product(x):
    """The same two inputs, plus the one feature that makes the problem linear."""
    return [x[0], x[1], x[0] * x[1]]


def fit(rows, features, lr=0.5, steps=5000):
    width = len(features(rows[0][0]))
    w = [0.0] * width
    b = 0.0
    n = len(rows)
    for _ in range(steps):
        dw = [0.0] * width
        db = 0.0
        for x, y in rows:
            f = features(x)
            e = sigmoid(sum(wi * fi for wi, fi in zip(w, f)) + b) - y
            for j in range(width):
                dw[j] += e * f[j]
            db += e
        w = [wi - lr * g / n for wi, g in zip(w, dw)]
        b -= lr * db / n
    return w, b


def report(rows, features, w, b):
    loss = 0.0
    right = 0
    for x, y in rows:
        f = features(x)
        p = sigmoid(sum(wi * fi for wi, fi in zip(w, f)) + b)
        loss += log_loss(p, y)
        right += 1 if (1 if p >= 0.5 else 0) == y else 0
    return loss / len(rows), right / len(rows)


w_plain, b_plain = fit(XOR, plain)
loss_plain, acc_plain = report(XOR, plain, w_plain, b_plain)
w_prod, b_prod = fit(XOR, with_product)
loss_prod, acc_prod = report(XOR, with_product, w_prod, b_prod)

print("the label is 1 when the two inputs disagree")
print("    x1     x2   label")
for x, y in XOR:
    print("  {:5.1f}  {:5.1f}      {}".format(x[0], x[1], y))

print()
print("fitted on x1 and x2 alone")
print("  weights {}  bias {:.6f}".format([round(v, 6) for v in w_plain], b_plain))
print("  loss {:.6f}   accuracy {:.2f}".format(loss_plain, acc_plain))
print("  log(2) for comparison: {:.6f}".format(math.log(2.0)))
print()
print("fitted on x1, x2 and their product")
print("  weights {}  bias {:.6f}".format([round(v, 4) for v in w_prod], b_prod))
print("  loss {:.6f}   accuracy {:.2f}".format(loss_prod, acc_prod))
print()
print("    x1     x2   label       p")
for x, y in XOR:
    f = with_product(x)
    print("  {:5.1f}  {:5.1f}      {}   {:.6f}".format(
        x[0], x[1], y, sigmoid(sum(wi * fi for wi, fi in zip(w_prod, f)) + b_prod)))

test("on the raw inputs, five thousand steps learn nothing at all",
     all(close(v, 0.0, 1e-9) for v in w_plain + [b_plain]), True)
test("the loss is stuck at log 2, the cost of guessing", close(loss_plain, math.log(2.0), 1e-9), True)
test("so it gets half of four points right, which is chance", close(acc_plain, 0.5), True)
test("adding the product of the two inputs fixes it completely", close(acc_prod, 1.0), True)
test("and the product is the feature carrying all the weight",
     abs(w_prod[2]) > 10.0 * max(abs(w_prod[0]), abs(w_prod[1]), abs(b_prod)), True)
test("its weight is negative, since agreement means label zero", w_prod[2] < 0.0, True)
test("the loss on the fixed model is far below log 2", loss_prod < 0.05, True)

the label is 1 when the two inputs disagree
    x1     x2   label
   -1.0   -1.0      0
   -1.0    1.0      1
    1.0   -1.0      1
    1.0    1.0      0

fitted on x1 and x2 alone
  weights [0.0, 0.0]  bias 0.000000
  loss 0.693147   accuracy 0.50
  log(2) for comparison: 0.693147

fitted on x1, x2 and their product
  weights [0.0, 0.0, -7.822]  bias -0.000000
  loss 0.000401   accuracy 1.00

    x1     x2   label       p
   -1.0   -1.0      0   0.000401
   -1.0    1.0      1   0.999599
    1.0   -1.0      1   0.999599
    1.0    1.0      0   0.000401
ok   on the raw inputs, five thousand steps learn nothing at all
ok   the loss is stuck at log 2, the cost of guessing
ok   so it gets half of four points right, which is chance
ok   adding the product of the two inputs fixes it completely
ok   and the product is the feature carrying all the weight
ok   its weight is negative, since agreement means label zero
ok   the loss on the fixed model is far below log 2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55# Four counts. Every metric in this lesson is arithmetic on these.
TRUTH     = [1] * 8 + [0] * 12
PREDICTED = [1, 1, 1, 1, 1, 1, 0, 0] + [1, 1, 1] + [0] * 9


def confusion(truth, predicted):
    if len(truth) != len(predicted):
        raise ValueError("{} labels against {} predictions".format(len(truth), len(predicted)))
    counts = {"tp": 0, "fp": 0, "fn": 0, "tn": 0}
    for y, p in zip(truth, predicted):
        if p == 1:
            counts["tp" if y == 1 else "fp"] += 1
        else:
            counts["fn" if y == 1 else "tn"] += 1
    return counts


m = confusion(TRUTH, PREDICTED)

print("twenty emails, 1 means spam.")
print("truth     ", "".join(str(v) for v in TRUTH))
print("predicted ", "".join(str(v) for v in PREDICTED))
print()
print("                     predicted spam   predicted not")
print("  actually spam      {:14}   {:13}".format(m["tp"], m["fn"]))
print("  actually not       {:14}   {:13}".format(m["fp"], m["tn"]))
print()
print("true positives  {:3}   caught spam".format(m["tp"]))
print("false positives {:3}   a real email sent to the spam folder".format(m["fp"]))
print("false negatives {:3}   spam that landed in the inbox".format(m["fn"]))
print("true negatives  {:3}   a real email left alone".format(m["tn"]))

perfect = confusion(TRUTH, TRUTH)
print()
print("the same matrix for a perfect predictor:", perfect)

# The length check is not decoration. Predictions and labels get shuffled, filtered
# and re-indexed on the way to a report, and zip() would quietly score the shorter one.
try:
    confusion([1, 0, 1], [1, 0])
    refused = False
except ValueError as error:
    refused = True
    print("mismatched lengths are refused:", error)

test("the four counts add up to every row", sum(m.values()), len(TRUTH))
test("true positives plus false negatives is the number of actual positives",
     m["tp"] + m["fn"], sum(TRUTH))
test("true positives plus false positives is the number we flagged",
     m["tp"] + m["fp"], sum(PREDICTED))
test("this filter caught six of eight and flagged three innocents",
     [m["tp"], m["fn"], m["fp"]], [6, 2, 3])
test("a perfect predictor has nothing off the diagonal",
     [perfect["fp"], perfect["fn"]], [0, 0])
test("labels and predictions of different lengths are refused", refused, True)

twenty emails, 1 means spam.
truth      11111111000000000000
predicted  11111100111000000000

                     predicted spam   predicted not
  actually spam                   6               2
  actually not                    3               9

true positives    6   caught spam
false positives   3   a real email sent to the spam folder
false negatives   2   spam that landed in the inbox
true negatives    9   a real email left alone

the same matrix for a perfect predictor: {'tp': 8, 'fp': 0, 'fn': 0, 'tn': 12}
mismatched lengths are refused: 3 labels against 2 predictions
ok   the four counts add up to every row
ok   true positives plus false negatives is the number of actual positives
ok   true positives plus false positives is the number we flagged
ok   this filter caught six of eight and flagged three innocents
ok   a perfect predictor has nothing off the diagonal
ok   labels and predictions of different lengths are refused

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77# The four metrics, from the four counts, with the zero denominators handled.
TRUTH     = [1] * 8 + [0] * 12
PREDICTED = [1, 1, 1, 1, 1, 1, 0, 0] + [1, 1, 1] + [0] * 9


def confusion(truth, predicted):
    counts = {"tp": 0, "fp": 0, "fn": 0, "tn": 0}
    for y, p in zip(truth, predicted):
        if p == 1:
            counts["tp" if y == 1 else "fp"] += 1
        else:
            counts["fn" if y == 1 else "tn"] += 1
    return counts


def ratio(top, bottom):
    """Zero when there is nothing to divide by. The alternative is a crash mid-report."""
    return top / bottom if bottom else 0.0


def accuracy(m):
    return ratio(m["tp"] + m["tn"], m["tp"] + m["fp"] + m["fn"] + m["tn"])


def precision(m):
    """Of the things we flagged, how many should have been."""
    return ratio(m["tp"], m["tp"] + m["fp"])


def recall(m):
    """Of the things that should have been flagged, how many were."""
    return ratio(m["tp"], m["tp"] + m["fn"])


def f1(m):
    """The harmonic mean of precision and recall. It follows the smaller one."""
    p, r = precision(m), recall(m)
    return ratio(2.0 * p * r, p + r)


m = confusion(TRUTH, PREDICTED)
p, r = precision(m), recall(m)

print("accuracy  {:.6f}   (6 + 9) / 20".format(accuracy(m)))
print("precision {:.6f}   6 / (6 + 3)".format(p))
print("recall    {:.6f}   6 / (6 + 2)".format(r))
print("F1        {:.6f}".format(f1(m)))
print()
print("arithmetic mean of precision and recall: {:.6f}".format((p + r) / 2.0))
print("F1, the harmonic mean                  : {:.6f}   <- always the smaller".format(f1(m)))

# Why the harmonic mean and not the ordinary one. A model that flags everything has
# perfect recall for free, and the arithmetic mean lets it look half competent.
flags_everything = confusion(TRUTH, [1] * len(TRUTH))
print()
print("a model that flags every email:")
print("  precision {:.6f}  recall {:.6f}".format(
    precision(flags_everything), recall(flags_everything)))
print("  arithmetic mean {:.6f}  <- looks respectable".format(
    (precision(flags_everything) + recall(flags_everything)) / 2.0))
print("  F1              {:.6f}  <- does not".format(f1(flags_everything)))

flags_nothing = confusion(TRUTH, [0] * len(TRUTH))
print()
print("a model that flags nothing has no precision to speak of:", precision(flags_nothing))
print("without the guard that is 0 / 0, and the report dies instead of the model.")

test("accuracy is fifteen right out of twenty", close(accuracy(m), 0.75), True)
test("precision is six of the nine we flagged", close(p, 6.0 / 9.0), True)
test("recall is six of the eight that were spam", close(r, 0.75), True)
test("F1 sits between the two and below their average",
     [min(p, r) <= f1(m) <= max(p, r), f1(m) < (p + r) / 2.0], [True, True])
test("flagging everything buys perfect recall and terrible precision",
     [close(recall(flags_everything), 1.0), close(precision(flags_everything), 0.4)], [True, True])
test("and F1 refuses to be impressed by it", f1(flags_everything) < 0.6, True)
test("flagging nothing returns zero precision rather than raising",
     close(precision(flags_nothing), 0.0), True)

accuracy  0.750000   (6 + 9) / 20
precision 0.666667   6 / (6 + 3)
recall    0.750000   6 / (6 + 2)
F1        0.705882

arithmetic mean of precision and recall: 0.708333
F1, the harmonic mean                  : 0.705882   <- always the smaller

a model that flags every email:
  precision 0.400000  recall 1.000000
  arithmetic mean 0.700000  <- looks respectable
  F1              0.571429  <- does not

a model that flags nothing has no precision to speak of: 0.0
without the guard that is 0 / 0, and the report dies instead of the model.
ok   accuracy is fifteen right out of twenty
ok   precision is six of the nine we flagged
ok   recall is six of the eight that were spam
ok   F1 sits between the two and below their average
ok   flagging everything buys perfect recall and terrible precision
ok   and F1 refuses to be impressed by it
ok   flagging nothing returns zero precision rather than raising

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66# The lesson. A model with no logic in it, scoring 99%.
def confusion(truth, predicted):
    counts = {"tp": 0, "fp": 0, "fn": 0, "tn": 0}
    for y, p in zip(truth, predicted):
        if p == 1:
            counts["tp" if y == 1 else "fp"] += 1
        else:
            counts["fn" if y == 1 else "tn"] += 1
    return counts


def ratio(top, bottom):
    return top / bottom if bottom else 0.0


def accuracy(m):
    return ratio(m["tp"] + m["tn"], sum(m.values()))


def precision(m):
    return ratio(m["tp"], m["tp"] + m["fp"])


def recall(m):
    return ratio(m["tp"], m["tp"] + m["fn"])


def f1(m):
    p, r = precision(m), recall(m)
    return ratio(2.0 * p * r, p + r)


TRUTH = [1] * 10 + [0] * 990          # ten ill patients in a thousand

# Model A contains no logic at all. It says "healthy" and stops.
lazy = confusion(TRUTH, [0] * 1000)

# Model B is a real screening test: eight of the ten found, forty false alarms
# that a second and more expensive test would clear.
useful = confusion(TRUTH, ([1] * 8 + [0] * 2) + ([1] * 40 + [0] * 950))

print("1000 patients, 10 of them ill. The base rate is 1%.")
print()
print("                          accuracy   precision   recall       F1")
for name, mm in (("predicts healthy always", lazy), ("a real screening test  ", useful)):
    print("  {}  {:9.4f}   {:9.4f}   {:6.4f}   {:6.4f}".format(
        name, accuracy(mm), precision(mm), recall(mm), f1(mm)))

print()
print("the useless model: ", lazy)
print("the useful model:  ", useful)
print()
print("Accuracy prefers the model that has never found a single patient.")
print("It misses {} cases out of {}, and is right 99% of the time.".format(
    lazy["fn"], sum(TRUTH)))

test("the do-nothing model scores ninety nine percent accuracy",
     close(accuracy(lazy), 0.99), True)
test("it finds nobody, so its recall is zero", close(recall(lazy), 0.0), True)
test("and its F1 is zero as well", close(f1(lazy), 0.0), True)
test("the useful model has lower accuracy than the useless one",
     accuracy(useful) < accuracy(lazy), True)
test("but it finds eight of the ten patients",
     [useful["tp"], close(recall(useful), 0.8)], [8, True])
test("and its F1 is not zero, which is the only score that noticed",
     [f1(useful) > 0.25, f1(lazy) == 0.0], [True, True])

1000 patients, 10 of them ill. The base rate is 1%.

                          accuracy   precision   recall       F1
  predicts healthy always     0.9900      0.0000   0.0000   0.0000
  a real screening test       0.9580      0.1667   0.8000   0.2759

the useless model:  {'tp': 0, 'fp': 0, 'fn': 10, 'tn': 990}
the useful model:   {'tp': 8, 'fp': 40, 'fn': 2, 'tn': 950}

Accuracy prefers the model that has never found a single patient.
It misses 10 cases out of 10, and is right 99% of the time.
ok   the do-nothing model scores ninety nine percent accuracy
ok   it finds nobody, so its recall is zero
ok   and its F1 is zero as well
ok   the useful model has lower accuracy than the useless one
ok   but it finds eight of the ten patients
ok   and its F1 is not zero, which is the only score that noticed

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74# Fifty scored examples, and the same model at eleven different thresholds.
import random


def counts_at(scored, threshold):
    m = {"tp": 0, "fp": 0, "fn": 0, "tn": 0}
    for score, y in scored:
        if score >= threshold:
            m["tp" if y == 1 else "fp"] += 1
        else:
            m["fn" if y == 1 else "tn"] += 1
    return m


def ratio(top, bottom):
    return top / bottom if bottom else 0.0


def precision(m):
    return ratio(m["tp"], m["tp"] + m["fp"])


def recall(m):
    return ratio(m["tp"], m["tp"] + m["fn"])


def f1(m):
    p, r = precision(m), recall(m)
    return ratio(2.0 * p * r, p + r)


rng = random.Random(42)
SCORED = []
for _ in range(40):
    SCORED.append((min(max(rng.gauss(0.35, 0.18), 0.0), 1.0), 0))
for _ in range(10):
    SCORED.append((min(max(rng.gauss(0.70, 0.18), 0.0), 1.0), 1))

BASE_RATE = sum(y for _, y in SCORED) / len(SCORED)
THRESHOLDS = [i / 10.0 for i in range(11)]

print("fifty scored examples, ten of them positive. base rate {:.2f}".format(BASE_RATE))
print()
print("  threshold   flagged   tp   fp   fn   precision   recall       F1")
table = {}
for t in THRESHOLDS:
    m = counts_at(SCORED, t)
    table[t] = m
    print("  {:9.1f}   {:7}   {:2}   {:2}   {:2}   {:9.4f}   {:6.4f}   {:6.4f}".format(
        t, m["tp"] + m["fp"], m["tp"], m["fp"], m["fn"], precision(m), recall(m), f1(m)))

recalls = [recall(table[t]) for t in THRESHOLDS]
best = max(THRESHOLDS, key=lambda t: f1(table[t]))
print()
print("best F1 is {:.4f}, at a threshold of {:.1f}".format(f1(table[best]), best))
print("at 0.3 : precision {:.4f}  recall {:.4f}   flag nearly everything".format(
    precision(table[0.3]), recall(table[0.3])))
print("at 0.7 : precision {:.4f}  recall {:.4f}   flag only the certain ones".format(
    precision(table[0.7]), recall(table[0.7])))
print()
print("One model, eleven behaviours, no retraining. Neither number means")
print("anything without the other one printed beside it.")

test("at a threshold of zero everything is flagged, so recall is perfect",
     close(recalls[0], 1.0), True)
test("and precision is then just the base rate", close(precision(table[0.0]), BASE_RATE), True)
test("recall never rises as the threshold rises",
     all(recalls[i + 1] <= recalls[i] for i in range(len(recalls) - 1)), True)
test("moving from 0.3 to 0.7 buys precision and pays in recall",
     [precision(table[0.7]) > precision(table[0.3]), recall(table[0.7]) < recall(table[0.3])],
     [True, True])
test("the best F1 is at neither extreme of the sweep", best not in (0.0, 1.0), True)
test("above every score nothing is flagged, and both numbers collapse to zero",
     [precision(table[1.0]), recall(table[1.0])], [0.0, 0.0])

fifty scored examples, ten of them positive. base rate 0.20

  threshold   flagged   tp   fp   fn   precision   recall       F1
        0.0        50   10   40    0      0.2000   1.0000   0.3333
        0.1        46   10   36    0      0.2174   1.0000   0.3571
        0.2        41   10   31    0      0.2439   1.0000   0.3922
        0.3        38   10   28    0      0.2632   1.0000   0.4167
        0.4        24   10   14    0      0.4167   1.0000   0.5882
        0.5        15   10    5    0      0.6667   1.0000   0.8000
        0.6         8    8    0    2      1.0000   0.8000   0.8889
        0.7         6    6    0    4      1.0000   0.6000   0.7500
        0.8         3    3    0    7      1.0000   0.3000   0.4615
        0.9         0    0    0   10      0.0000   0.0000   0.0000
        1.0         0    0    0   10      0.0000   0.0000   0.0000

best F1 is 0.8889, at a threshold of 0.6
at 0.3 : precision 0.2632  recall 1.0000   flag nearly everything
at 0.7 : precision 1.0000  recall 0.6000   flag only the certain ones

One model, eleven behaviours, no retraining. Neither number means
anything without the other one printed beside it.
ok   at a threshold of zero everything is flagged, so recall is perfect
ok   and precision is then just the base rate
ok   recall never rises as the threshold rises
ok   moving from 0.3 to 0.7 buys precision and pays in recall
ok   the best F1 is at neither extreme of the sweep
ok   above every score nothing is flagged, and both numbers collapse to zero

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79# Macro against micro, on 900 common cases, 90 uncommon and 10 rare.
def ratio(top, bottom):
    return top / bottom if bottom else 0.0


CLASSES = ("common", "uncommon", "rare")

# rows are the truth, columns are what was predicted
ALWAYS_COMMON = {
    "common":   {"common": 900, "uncommon": 0, "rare": 0},
    "uncommon": {"common":  90, "uncommon": 0, "rare": 0},
    "rare":     {"common":  10, "uncommon": 0, "rare": 0},
}
BALANCED = {
    "common":   {"common": 800, "uncommon": 80, "rare": 20},
    "uncommon": {"common":  20, "uncommon": 65, "rare":  5},
    "rare":     {"common":   3, "uncommon":  1, "rare":  6},
}


def one_class(matrix, label):
    """Collapse a multi-class matrix into this class against everything else."""
    tp = matrix[label][label]
    fp = sum(matrix[other][label] for other in CLASSES if other != label)
    fn = sum(matrix[label][other] for other in CLASSES if other != label)
    return {"tp": tp, "fp": fp, "fn": fn}


def f1_of(m):
    p = ratio(m["tp"], m["tp"] + m["fp"])
    r = ratio(m["tp"], m["tp"] + m["fn"])
    return ratio(2.0 * p * r, p + r)


def accuracy(matrix):
    right = sum(matrix[c][c] for c in CLASSES)
    total = sum(matrix[a][b] for a in CLASSES for b in CLASSES)
    return right / total


def macro_f1(matrix):
    """Average the per-class F1 scores. Every class counts the same."""
    return sum(f1_of(one_class(matrix, c)) for c in CLASSES) / len(CLASSES)


def micro_f1(matrix):
    """Pool every count first, then compute one F1. Large classes dominate."""
    pooled = {"tp": 0, "fp": 0, "fn": 0}
    for c in CLASSES:
        for key, value in one_class(matrix, c).items():
            pooled[key] += value
    return f1_of(pooled)


print("900 common cases, 90 uncommon, 10 rare.")
for name, matrix in (("always says common", ALWAYS_COMMON), ("spreads its bets ", BALANCED)):
    print()
    print(name)
    print("  accuracy {:.6f}   micro F1 {:.6f}   macro F1 {:.6f}".format(
        accuracy(matrix), micro_f1(matrix), macro_f1(matrix)))
    for c in CLASSES:
        print("     {:9} F1 {:.6f}   {}".format(
            c, f1_of(one_class(matrix, c)), one_class(matrix, c)))

print()
print("The second model is WORSE on accuracy and far better at the job.")
print("Micro F1 agrees with accuracy and hides that. Macro F1 does not.")

test("always guessing the common class scores ninety percent accuracy",
     close(accuracy(ALWAYS_COMMON), 0.9), True)
test("micro F1 is exactly accuracy, so it adds nothing here",
     close(micro_f1(ALWAYS_COMMON), accuracy(ALWAYS_COMMON)), True)
test("macro F1 lands near a third, because two classes score nothing",
     close(macro_f1(ALWAYS_COMMON), 0.315789, 1e-5), True)
test("both minority classes score exactly zero F1",
     [f1_of(one_class(ALWAYS_COMMON, "uncommon")), f1_of(one_class(ALWAYS_COMMON, "rare"))],
     [0.0, 0.0])
test("the better model has lower accuracy", accuracy(BALANCED) < accuracy(ALWAYS_COMMON), True)
test("and nearly double the macro F1", macro_f1(BALANCED) > 1.8 * macro_f1(ALWAYS_COMMON), True)

900 common cases, 90 uncommon, 10 rare.

always says common
  accuracy 0.900000   micro F1 0.900000   macro F1 0.315789
     common    F1 0.947368   {'tp': 900, 'fp': 100, 'fn': 0}
     uncommon  F1 0.000000   {'tp': 0, 'fp': 0, 'fn': 90}
     rare      F1 0.000000   {'tp': 0, 'fp': 0, 'fn': 10}

spreads its bets 
  accuracy 0.871000   micro F1 0.871000   macro F1 0.590714
     common    F1 0.928613   {'tp': 800, 'fp': 23, 'fn': 100}
     uncommon  F1 0.550847   {'tp': 65, 'fp': 81, 'fn': 25}
     rare      F1 0.292683   {'tp': 6, 'fp': 25, 'fn': 4}

The second model is WORSE on accuracy and far better at the job.
Micro F1 agrees with accuracy and hides that. Macro F1 does not.
ok   always guessing the common class scores ninety percent accuracy
ok   micro F1 is exactly accuracy, so it adds nothing here
ok   macro F1 lands near a third, because two classes score nothing
ok   both minority classes score exactly zero F1
ok   the better model has lower accuracy
ok   and nearly double the macro F1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99# Sixteen noisy points, and a polynomial of any degree fitted to them.
import math
import random

rng = random.Random(42)                 # seeded, so this page is the same every build

XS = [-1.0 + i * (2.0 / 15.0) for i in range(16)]
YS = [math.sin(3.0 * x) + rng.gauss(0.0, 0.10) for x in XS]


def solve(matrix, rhs):
    """Gaussian elimination with partial pivoting. Small, and enough for this."""
    n = len(rhs)
    a = [list(row) + [rhs[i]] for i, row in enumerate(matrix)]
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(a[r][col]))
        if abs(a[pivot][col]) < 1e-14:
            raise ValueError("the system is singular at column {}".format(col))
        a[col], a[pivot] = a[pivot], a[col]
        for r in range(col + 1, n):
            factor = a[r][col] / a[col][col]
            for c in range(col, n + 1):
                a[r][c] -= factor * a[col][c]
    out = [0.0] * n
    for r in range(n - 1, -1, -1):
        out[r] = (a[r][n] - sum(a[r][c] * out[c] for c in range(r + 1, n))) / a[r][r]
    return out


def polyfit(xs, ys, degree):
    """Least squares for c0 + c1*x + c2*x^2 + ..., by the normal equations of lesson 8."""
    width = degree + 1
    if len(xs) < width:
        raise ValueError("{} points cannot pin down {} coefficients".format(len(xs), width))
    rows = [[x ** k for k in range(width)] for x in xs]
    ata = [[sum(r[i] * r[j] for r in rows) for j in range(width)] for i in range(width)]
    aty = [sum(r[i] * y for r, y in zip(rows, ys)) for i in range(width)]
    return solve(ata, aty)


def poly(coeffs, x):
    total = 0.0
    for c in reversed(coeffs):          # Horner's method, evaluated from the inside out
        total = total * x + c
    return total


def rmse(coeffs, xs, ys):
    """Root mean squared error, in the same units as y so it can be read."""
    return math.sqrt(sum((poly(coeffs, x) - y) ** 2 for x, y in zip(xs, ys)) / len(xs))


print("sixteen points from sin(3x), with gaussian noise of 0.10 added")
print("      x        y     sin(3x)")
for i in (0, 1, 7, 8, 14, 15):
    print("  {:5.2f}  {:7.4f}   {:8.4f}".format(XS[i], YS[i], math.sin(3.0 * XS[i])))

print()
print("  degree   parameters   error on the data it was fitted to")
errors = []
for degree in range(10):
    errors.append(rmse(polyfit(XS, YS, degree), XS, YS))
    print("  {:5}   {:9}   {:.6f}".format(degree, degree + 1, errors[-1]))

print()
print("Every extra parameter lowered it. Not one of those numbers is evidence.")

# A polynomial with as many coefficients as there are points has no work left to
# do: it passes through every one of them exactly, whatever they happen to be.
FOUR_X = [-1.0, 0.0, 1.0, 2.0]
FOUR_Y = [1.0, 2.0, 0.0, 5.0]
cubic = polyfit(FOUR_X, FOUR_Y, 3)
print()
print("four arbitrary points, one cubic:", [round(v, 6) for v in cubic])
print("what it predicts at each of them:", [round(poly(cubic, x), 9) for x in FOUR_X])
print("error: {:.3e}   <- a flawless fit, and it means nothing".format(rmse(cubic, FOUR_X, FOUR_Y)))

try:
    polyfit(FOUR_X, FOUR_Y, 7)
    refused = False
except ValueError as error:
    refused = True
    print()
    print("asking for more coefficients than points is refused:", error)

drops = [errors[i] - errors[i + 1] for i in range(len(errors) - 1)]

test("the noise moved every point off the true curve",
     all(not close(y, math.sin(3.0 * x)) for x, y in zip(XS, YS)), True)
test("a flat line does badly, as it should", errors[0] > 0.6, True)
test("the training error falls at every single degree",
     all(errors[i + 1] < errors[i] for i in range(len(errors) - 1)), True)
test("the largest single improvement is the step up to degree three",
     drops.index(max(drops)), 2)
test("the last four parameters bought less than a fifth of the error",
     (errors[5] - errors[9]) / errors[5] < 0.2, True)
test("a cubic through four points hits all four exactly",
     close(rmse(cubic, FOUR_X, FOUR_Y), 0.0, 1e-12), True)
test("more coefficients than data points is refused", refused, True)

sixteen points from sin(3x), with gaussian noise of 0.10 added
      x        y     sin(3x)
  -1.00  -0.1555    -0.1411
  -0.87  -0.5328    -0.5155
  -0.07  -0.2254    -0.1987
   0.07   0.1770     0.1987
   0.87   0.4417     0.5155
   1.00   0.0397     0.1411

  degree   parameters   error on the data it was fitted to
      0           1   0.723787
      1           2   0.488223
      2           3   0.487946
      3           4   0.093185
      4           5   0.089890
      5           6   0.044409
      6           7   0.041709
      7           8   0.041566
      8           9   0.039136
      9          10   0.038182

Every extra parameter lowered it. Not one of those numbers is evidence.

four arbitrary points, one cubic: [2.0, -2.166667, -1.5, 1.666667]
what it predicts at each of them: [1.0, 2.0, 0.0, 5.0]
error: 2.483e-15   <- a flawless fit, and it means nothing

asking for more coefficients than points is refused: 4 points cannot pin down 8 coefficients
ok   the noise moved every point off the true curve
ok   a flat line does badly, as it should
ok   the training error falls at every single degree
ok   the largest single improvement is the step up to degree three
ok   the last four parameters bought less than a fifth of the error
ok   a cubic through four points hits all four exactly
ok   more coefficients than data points is refused

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99# The same sixteen points, split in two before anything is fitted.
import math
import random


def solve(matrix, rhs):
    n = len(rhs)
    a = [list(row) + [rhs[i]] for i, row in enumerate(matrix)]
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(a[r][col]))
        if abs(a[pivot][col]) < 1e-14:
            raise ValueError("singular")
        a[col], a[pivot] = a[pivot], a[col]
        for r in range(col + 1, n):
            factor = a[r][col] / a[col][col]
            for c in range(col, n + 1):
                a[r][c] -= factor * a[col][c]
    out = [0.0] * n
    for r in range(n - 1, -1, -1):
        out[r] = (a[r][n] - sum(a[r][c] * out[c] for c in range(r + 1, n))) / a[r][r]
    return out


def polyfit(xs, ys, degree):
    width = degree + 1
    rows = [[x ** k for k in range(width)] for x in xs]
    ata = [[sum(r[i] * r[j] for r in rows) for j in range(width)] for i in range(width)]
    aty = [sum(r[i] * y for r, y in zip(rows, ys)) for i in range(width)]
    return solve(ata, aty)


def poly(coeffs, x):
    total = 0.0
    for c in reversed(coeffs):
        total = total * x + c
    return total


def rmse(coeffs, xs, ys):
    return math.sqrt(sum((poly(coeffs, x) - y) ** 2 for x, y in zip(xs, ys)) / len(xs))


def train_test_split(xs, ys, n_train, seed):
    """Shuffle the row numbers with a named seed, then cut. The seed is the point:
    an unseeded split scores differently every run, so nothing is comparable."""
    order = list(range(len(xs)))
    random.Random(seed).shuffle(order)
    train, held = sorted(order[:n_train]), sorted(order[n_train:])
    return (([xs[i] for i in train], [ys[i] for i in train]),
            ([xs[i] for i in held], [ys[i] for i in held]))


rng = random.Random(42)
XS = [-1.0 + i * (2.0 / 15.0) for i in range(16)]
YS = [math.sin(3.0 * x) + rng.gauss(0.0, 0.10) for x in XS]

(xtr, ytr), (xte, yte) = train_test_split(XS, YS, 10, seed=42)

print("ten rows to fit on, six held back and not looked at until the end")
print("train x:", [round(v, 3) for v in xtr])
print("test  x:", [round(v, 3) for v in xte])
print()
print("  degree        train           test     verdict")
train_err = []
test_err = []
for degree in range(10):
    coeffs = polyfit(xtr, ytr, degree)
    train_err.append(rmse(coeffs, xtr, ytr))
    test_err.append(rmse(coeffs, xte, yte))
    if degree == 0:
        verdict = ""
    elif test_err[-1] < min(test_err[:-1]):
        verdict = "better"
    else:
        verdict = "WORSE, and getting worse"
    print("  {:5}   {:12.8f}   {:11.6f}   {}".format(
        degree, train_err[-1], test_err[-1], verdict))

best = min(range(10), key=lambda d: test_err[d])
print()
print("lowest test error at degree {}: {:.6f}".format(best, test_err[best]))
print("degree 9 fits the training rows to {:.1e} and is {:.0f}x worse on the held-out six".format(
    train_err[9], test_err[9] / test_err[best]))
print()
print("The two columns move in opposite directions after degree {}. That gap is".format(best))
print("the whole of overfitting: the model is memorising ten points, noise included.")

test("training error falls all the way down, as it always will",
     all(train_err[i + 1] < train_err[i] for i in range(9)), True)
test("with ten rows and ten coefficients the training error reaches zero",
     close(train_err[9], 0.0, 1e-8), True)
test("test error is worse than training error at every degree",
     all(test_err[d] > train_err[d] for d in range(10)), True)
test("the best degree on held-out data is neither the smallest nor the largest",
     best not in (0, 9), True)
test("past that degree the test error rises every time",
     all(test_err[d + 1] > test_err[d] for d in range(best, 9)), True)
test("the flawless fit is two hundred times worse than the best one",
     test_err[9] / test_err[best] > 200.0, True)

ten rows to fit on, six held back and not looked at until the end
train x: [-0.867, -0.733, -0.333, -0.2, -0.067, 0.067, 0.2, 0.333, 0.6, 0.867]
test  x: [-1.0, -0.6, -0.467, 0.467, 0.733, 1.0]

  degree        train           test     verdict
      0     0.68250093      0.787863   
      1     0.38416891      0.661124   better
      2     0.38416119      0.661557   WORSE, and getting worse
      3     0.05224970      0.239638   better
      4     0.05166971      0.233405   better
      5     0.03995617      0.072441   better
      6     0.03369870      0.298805   WORSE, and getting worse
      7     0.02382779      1.091486   WORSE, and getting worse
      8     0.00738631      6.332269   WORSE, and getting worse
      9     0.00000000     15.920658   WORSE, and getting worse

lowest test error at degree 5: 0.072441
degree 9 fits the training rows to 3.8e-11 and is 220x worse on the held-out six

The two columns move in opposite directions after degree 5. That gap is
the whole of overfitting: the model is memorising ten points, noise included.
ok   training error falls all the way down, as it always will
ok   with ten rows and ten coefficients the training error reaches zero
ok   test error is worse than training error at every degree
ok   the best degree on held-out data is neither the smallest nor the largest
ok   past that degree the test error rises every time
ok   the flawless fit is two hundred times worse than the best one

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110# Ridge regression: the same solve, with one number added to the diagonal.
import math
import random


def solve(matrix, rhs):
    n = len(rhs)
    a = [list(row) + [rhs[i]] for i, row in enumerate(matrix)]
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(a[r][col]))
        if abs(a[pivot][col]) < 1e-14:
            raise ValueError("singular")
        a[col], a[pivot] = a[pivot], a[col]
        for r in range(col + 1, n):
            factor = a[r][col] / a[col][col]
            for c in range(col, n + 1):
                a[r][c] -= factor * a[col][c]
    out = [0.0] * n
    for r in range(n - 1, -1, -1):
        out[r] = (a[r][n] - sum(a[r][c] * out[c] for c in range(r + 1, n))) / a[r][r]
    return out


def polyfit(xs, ys, degree, penalty=0.0):
    """Minimise squared error PLUS penalty times the sum of the squared coefficients.

    Setting the derivative of that to zero adds the penalty to the diagonal of the
    normal equations, and changes nothing else. Index 0 is left alone: penalising
    the intercept would only drag every prediction towards zero."""
    width = degree + 1
    rows = [[x ** k for k in range(width)] for x in xs]
    ata = [[sum(r[i] * r[j] for r in rows) for j in range(width)] for i in range(width)]
    aty = [sum(r[i] * y for r, y in zip(rows, ys)) for i in range(width)]
    for i in range(1, width):
        ata[i][i] += penalty
    return solve(ata, aty)


def poly(coeffs, x):
    total = 0.0
    for c in reversed(coeffs):
        total = total * x + c
    return total


def rmse(coeffs, xs, ys):
    return math.sqrt(sum((poly(coeffs, x) - y) ** 2 for x, y in zip(xs, ys)) / len(xs))


def train_test_split(xs, ys, n_train, seed):
    order = list(range(len(xs)))
    random.Random(seed).shuffle(order)
    train, held = sorted(order[:n_train]), sorted(order[n_train:])
    return (([xs[i] for i in train], [ys[i] for i in train]),
            ([xs[i] for i in held], [ys[i] for i in held]))


rng = random.Random(42)
XS = [-1.0 + i * (2.0 / 15.0) for i in range(16)]
YS = [math.sin(3.0 * x) + rng.gauss(0.0, 0.10) for x in XS]
(xtr, ytr), (xte, yte) = train_test_split(XS, YS, 10, seed=42)

# First the symptom. The degree 9 fit is not slightly wrong, it is violent.
wild = polyfit(xtr, ytr, 9)
tame = polyfit(xtr, ytr, 5)
print("degree 9, no penalty:", [round(v, 1) for v in wild])
print("degree 5, no penalty:", [round(v, 3) for v in tame])
print()
print("largest coefficient: {:.1f} against {:.1f}".format(
    max(abs(v) for v in wild), max(abs(v) for v in tame)))
print("at x = -1 the truth is {:.4f} and degree 9 predicts {:.2f}".format(
    math.sin(-3.0), poly(wild, -1.0)))
print("Huge terms nearly cancelling. That is what a curve which threads ten points")
print("and then leaves the page between them is made of.")

print()
print("    penalty       train        test    largest coefficient")
PENALTIES = (0.0, 1e-6, 1e-4, 1e-3, 3e-3, 1e-2, 1e-1, 1.0, 10.0)
result = {}
for penalty in PENALTIES:
    coeffs = polyfit(xtr, ytr, 9, penalty)
    result[penalty] = (rmse(coeffs, xtr, ytr), rmse(coeffs, xte, yte),
                       max(abs(v) for v in coeffs))
    print("  {:9.0e}   {:9.5f}   {:9.5f}   {:12.4f}".format(penalty, *result[penalty]))

best = min(PENALTIES, key=lambda p: result[p][1])
print()
print("best penalty {:.0e}: test error {:.5f}, down from {:.5f}".format(
    best, result[best][1], result[0.0][1]))
print("plain degree 5, no penalty at all:  {:.5f}".format(rmse(tame, xte, yte)))
print()
print("The penalty rescued an oversized model and still did not beat the right-sized")
print("one. Regularisation buys tolerance for a bad guess at how big to go.")

test("the unpenalised degree nine coefficients run into the hundreds",
     max(abs(v) for v in wild) > 300.0, True)
test("and it predicts minus thirty one where the answer is near zero",
     poly(wild, -1.0) < -30.0, True)
test("a small penalty cuts the largest coefficient a hundredfold",
     result[best][2] < result[0.0][2] / 100.0, True)
test("and cuts the test error a hundredfold with it",
     result[best][1] < result[0.0][1] / 100.0, True)
test("the training error only ever rises as the penalty rises",
     all(result[PENALTIES[i + 1]][0] > result[PENALTIES[i]][0]
         for i in range(len(PENALTIES) - 1)), True)
test("so the best penalty is in the middle, because too much of it underfits",
     [best not in (PENALTIES[0], PENALTIES[-1]), result[10.0][1] > result[best][1]],
     [True, True])
test("even penalised, degree nine does not beat plain degree five",
     result[best][1] > rmse(tame, xte, yte), True)

degree 9, no penalty: [-0.0, 3.1, 2.7, -11.4, -35.5, 111.2, 83.2, -309.6, -54.4, 234.0]
degree 5, no penalty: [-0.023, 3.17, 0.114, -4.545, -0.18, 1.42]

largest coefficient: 309.6 against 4.5
at x = -1 the truth is -0.1411 and degree 9 predicts -31.32
Huge terms nearly cancelling. That is what a curve which threads ten points
and then leaves the page between them is made of.

    penalty       train        test    largest coefficient
      0e+00     0.00000    15.92066       309.5691
      1e-06     0.01941     1.79600        21.4250
      1e-04     0.02922     0.46044         3.0344
      1e-03     0.03479     0.21319         3.1667
      3e-03     0.04578     0.15307         2.8997
      1e-02     0.07429     0.32293         2.6517
      1e-01     0.17796     0.34037         1.9012
      1e+00     0.40255     0.53508         0.8157
      1e+01     0.57601     0.69381         0.2154

best penalty 3e-03: test error 0.15307, down from 15.92066
plain degree 5, no penalty at all:  0.07244

The penalty rescued an oversized model and still did not beat the right-sized
one. Regularisation buys tolerance for a bad guess at how big to go.
ok   the unpenalised degree nine coefficients run into the hundreds
ok   and it predicts minus thirty one where the answer is near zero
ok   a small penalty cuts the largest coefficient a hundredfold
ok   and cuts the test error a hundredfold with it
ok   the training error only ever rises as the penalty rises
ok   so the best penalty is in the middle, because too much of it underfits
ok   even penalised, degree nine does not beat plain degree five

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110# k-fold cross validation: choose the degree without ever touching the test set.
import math
import random


def solve(matrix, rhs):
    n = len(rhs)
    a = [list(row) + [rhs[i]] for i, row in enumerate(matrix)]
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(a[r][col]))
        if abs(a[pivot][col]) < 1e-14:
            raise ValueError("singular")
        a[col], a[pivot] = a[pivot], a[col]
        for r in range(col + 1, n):
            factor = a[r][col] / a[col][col]
            for c in range(col, n + 1):
                a[r][c] -= factor * a[col][c]
    out = [0.0] * n
    for r in range(n - 1, -1, -1):
        out[r] = (a[r][n] - sum(a[r][c] * out[c] for c in range(r + 1, n))) / a[r][r]
    return out


def polyfit(xs, ys, degree):
    width = degree + 1
    if len(xs) < width:
        raise ValueError("{} points cannot pin down {} coefficients".format(len(xs), width))
    rows = [[x ** k for k in range(width)] for x in xs]
    ata = [[sum(r[i] * r[j] for r in rows) for j in range(width)] for i in range(width)]
    aty = [sum(r[i] * y for r, y in zip(rows, ys)) for i in range(width)]
    return solve(ata, aty)


def poly(coeffs, x):
    total = 0.0
    for c in reversed(coeffs):
        total = total * x + c
    return total


def rmse(coeffs, xs, ys):
    return math.sqrt(sum((poly(coeffs, x) - y) ** 2 for x, y in zip(xs, ys)) / len(xs))


def train_test_split(xs, ys, n_train, seed):
    order = list(range(len(xs)))
    random.Random(seed).shuffle(order)
    train, held = sorted(order[:n_train]), sorted(order[n_train:])
    return (([xs[i] for i in train], [ys[i] for i in train]),
            ([xs[i] for i in held], [ys[i] for i in held]))


def folds_of(n, k, seed):
    """Shuffle the rows, then deal them round like cards into k piles."""
    order = list(range(n))
    random.Random(seed).shuffle(order)
    return [sorted(order[i::k]) for i in range(k)]


def cross_validate(xs, ys, degree, k, seed):
    """Fit on k-1 folds, score the one left out, rotate, pool the squared errors."""
    total = 0.0
    for held in folds_of(len(xs), k, seed):
        keep = [i for i in range(len(xs)) if i not in held]
        coeffs = polyfit([xs[i] for i in keep], [ys[i] for i in keep], degree)
        total += sum((poly(coeffs, xs[i]) - ys[i]) ** 2 for i in held)
    return math.sqrt(total / len(xs))


rng = random.Random(42)
XS = [-1.0 + i * (2.0 / 15.0) for i in range(16)]
YS = [math.sin(3.0 * x) + rng.gauss(0.0, 0.10) for x in XS]
(xtr, ytr), (xte, yte) = train_test_split(XS, YS, 10, seed=42)

FOLDS = folds_of(len(xtr), 5, seed=42)
print("ten training rows dealt into five folds:", FOLDS)
print("each fold is the scoring set once and training data the other four times")
print()
print("  degree   5-fold cv    held-out test")
cv = {}
for degree in range(8):
    cv[degree] = cross_validate(xtr, ytr, degree, 5, seed=42)
    print("  {:5}   {:9.5f}   {:13.5f}".format(
        degree, cv[degree], rmse(polyfit(xtr, ytr, degree), xte, yte)))

chosen = min(cv, key=lambda d: cv[d])
honest = min(range(8), key=lambda d: rmse(polyfit(xtr, ytr, d), xte, yte))
print()
print("cross validation chooses degree {}, from training rows only".format(chosen))
print("the held-out set would have chosen degree {}".format(honest))
print()
print("Degrees 8 and 9 are missing from the table on purpose. A fold leaves eight")
print("rows to fit on, and eight rows cannot pin down nine or ten coefficients:")
try:
    cross_validate(xtr, ytr, 9, 5, seed=42)
    refused = False
except ValueError as error:
    refused = True
    print("   ", error)

test("the folds cover every training row exactly once",
     sorted(i for fold in FOLDS for i in fold), list(range(len(xtr))))
test("five folds of ten rows means two rows scored per fit",
     [len(f) for f in FOLDS], [2, 2, 2, 2, 2])
test("cross validation rejects the flexible models too", cv[7] > 10.0 * cv[chosen], True)
test("it picks a small degree having never seen the test set", chosen <= 5, True)
test("it does not pick exactly what the test set would have picked", chosen == honest, False)
test("but both of them reject everything above degree five",
     [cv[6] > cv[chosen], cv[7] > cv[6]], [True, True])
test("a degree too large to fit inside a fold is refused, not fudged", refused, True)

ten training rows dealt into five folds: [[6, 7], [3, 9], [2, 4], [0, 8], [1, 5]]
each fold is the scoring set once and training data the other four times

  degree   5-fold cv    held-out test
      0     0.76478         0.78786
      1     0.52406         0.66112
      2     0.77012         0.66156
      3     0.11620         0.23964
      4     0.27002         0.23340
      5     0.51500         0.07244
      6     1.96450         0.29880
      7     6.46259         1.09149

cross validation chooses degree 3, from training rows only
the held-out set would have chosen degree 5

Degrees 8 and 9 are missing from the table on purpose. A fold leaves eight
rows to fit on, and eight rows cannot pin down nine or ten coefficients:
    8 points cannot pin down 10 coefficients
ok   the folds cover every training row exactly once
ok   five folds of ten rows means two rows scored per fit
ok   cross validation rejects the flexible models too
ok   it picks a small degree having never seen the test set
ok   it does not pick exactly what the test set would have picked
ok   but both of them reject everything above degree five
ok   a degree too large to fit inside a fold is refused, not fudged

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85# k-nearest neighbours. Sort by distance, count votes, done.
import math

# Four class 0 points around the origin, eight class 1 points far away, and one
# stray class 1 sitting almost on top of the query. The stray is a mislabelled row.
LABELLED = [
    ((0.10, 0.05), 1),          # the stray: nearest of all, and wrong
    ((0.40, 0.30), 0),
    ((-0.35, 0.45), 0),
    ((0.50, -0.40), 0),
    ((-0.50, -0.50), 0),
    ((2.30, 2.60), 1),
    ((2.70, 2.40), 1),
    ((2.50, 2.90), 1),
    ((2.90, 2.50), 1),
    ((2.10, 2.30), 1),
    ((2.60, 2.10), 1),
    ((3.00, 2.80), 1),
    ((2.40, 3.10), 1),
]
QUERY = (0.0, 0.0)


def distance(a, b):
    """Straight-line distance, from lesson 5. Everything here rests on it."""
    return math.sqrt(sum((p - q) ** 2 for p, q in zip(a, b)))


def neighbours(rows, query, k):
    """The k closest rows, nearest first. The sort is the whole algorithm."""
    if k < 1 or k > len(rows):
        raise ValueError("k must be between 1 and {}, got {}".format(len(rows), k))
    return sorted(rows, key=lambda row: distance(row[0], query))[:k]


def classify(rows, query, k):
    """Majority vote. A tie goes to whichever tied class has the nearer member."""
    near = neighbours(rows, query, k)
    votes = {}
    for point, label in near:
        votes[label] = votes.get(label, 0) + 1
    most = max(votes.values())
    tied = [label for label, n in votes.items() if n == most]
    if len(tied) == 1:
        return tied[0]
    for point, label in near:               # nearest first, so the first tied one wins
        if label in tied:
            return label


print("thirteen labelled points, sorted by distance from a query at the origin")
print("  distance   label   point")
for point, label in sorted(LABELLED, key=lambda row: distance(row[0], QUERY)):
    print("  {:8.4f}   {:5}   {}".format(distance(point, QUERY), label, point))

print()
print("    k   labels of the k nearest    votes 0/1   answer")
answer = {}
for k in (1, 2, 3, 5, 7, 9, 13):
    near = [label for _, label in neighbours(LABELLED, QUERY, k)]
    answer[k] = classify(LABELLED, QUERY, k)
    print("  {:3}   {:24}    {:3} /{:3}     {}".format(
        k, "".join(str(v) for v in near), near.count(0), near.count(1), answer[k]))

print()
print("k = 1   believes the one mislabelled row sitting next to the query.")
print("k = 3   outvotes it, and gives the answer the neighbourhood supports.")
print("k = 9   reaches into the far cluster, and the vote flips back.")
print("k = 13  is every row in the dataset, so it just reports the majority class.")

try:
    classify(LABELLED, QUERY, 20)
    refused = False
except ValueError as error:
    print()
    refused = True
    print("asking for more neighbours than there are rows is refused:", error)

test("the single nearest neighbour is the mislabelled stray", answer[1], 1)
test("three neighbours outvote it", answer[3], 0)
test("five and seven agree with three", [answer[5], answer[7]], [0, 0])
test("nine reaches too far and flips the answer back", answer[9], 1)
test("so one point, one dataset, and two different answers", answer[1] == answer[3], False)
test("an even k can tie, and then the tie-break decides it, not the data", answer[2], 1)
test("k outside one to the number of rows is refused", refused, True)

thirteen labelled points, sorted by distance from a query at the origin
  distance   label   point
    0.1118       1   (0.1, 0.05)
    0.5000       0   (0.4, 0.3)
    0.5701       0   (-0.35, 0.45)
    0.6403       0   (0.5, -0.4)
    0.7071       0   (-0.5, -0.5)
    3.1145       1   (2.1, 2.3)
    3.3422       1   (2.6, 2.1)
    3.4713       1   (2.3, 2.6)
    3.6125       1   (2.7, 2.4)
    3.8288       1   (2.5, 2.9)
    3.8288       1   (2.9, 2.5)
    3.9205       1   (2.4, 3.1)
    4.1037       1   (3.0, 2.8)

    k   labels of the k nearest    votes 0/1   answer
    1   1                             0 /  1     1
    2   10                            1 /  1     1
    3   100                           2 /  1     0
    5   10000                         4 /  1     0
    7   1000011                       4 /  3     0
    9   100001111                     4 /  5     1
   13   1000011111111                 4 /  9     1

k = 1   believes the one mislabelled row sitting next to the query.
k = 3   outvotes it, and gives the answer the neighbourhood supports.
k = 9   reaches into the far cluster, and the vote flips back.
k = 13  is every row in the dataset, so it just reports the majority class.

asking for more neighbours than there are rows is refused: k must be between 1 and 13, got 20
ok   the single nearest neighbour is the mislabelled stray
ok   three neighbours outvote it
ok   five and seven agree with three
ok   nine reaches too far and flips the answer back
ok   so one point, one dataset, and two different answers
ok   an even k can tie, and then the tie-break decides it, not the data
ok   k outside one to the number of rows is refused

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76# k = 1 has zero training error and has learned nothing. Lesson 11, again.
import math
import random


def distance(a, b):
    return math.sqrt(sum((p - q) ** 2 for p, q in zip(a, b)))


def classify(rows, query, k):
    near = sorted(rows, key=lambda row: distance(row[0], query))[:k]
    votes = {}
    for point, label in near:
        votes[label] = votes.get(label, 0) + 1
    most = max(votes.values())
    tied = [label for label, n in votes.items() if n == most]
    if len(tied) == 1:
        return tied[0]
    for point, label in near:
        if label in tied:
            return label


rng = random.Random(42)
ROWS = []
for _ in range(50):
    ROWS.append(((rng.gauss(0.0, 1.3), rng.gauss(0.0, 1.3)), 0))
for _ in range(50):
    ROWS.append(((rng.gauss(2.2, 1.3), rng.gauss(2.2, 1.3)), 1))

order = list(range(len(ROWS)))
random.Random(42).shuffle(order)
TRAIN = [ROWS[i] for i in order[:60]]
HELD = [ROWS[i] for i in order[60:]]


def right(data, k):
    """A count, not a ratio, so an assertion can compare whole numbers exactly."""
    return sum(1 for point, label in data if classify(TRAIN, point, k) == label)


majority = max(sum(1 for _, lab in TRAIN if lab == 0), sum(1 for _, lab in TRAIN if lab == 1))
print("two overlapping clusters. 60 rows to look things up in, 40 held back.")
print("the majority class holds {} of the {} training rows".format(majority, len(TRAIN)))
print()
print("     k    train      test       gap")
KS = (1, 3, 5, 7, 11, 17, 25, 39, 60)
train_acc = {}
test_acc = {}
for k in KS:
    train_acc[k] = right(TRAIN, k) / len(TRAIN)
    test_acc[k] = right(HELD, k) / len(HELD)
    print("   {:3}   {:.4f}   {:.4f}   {:+.4f}".format(
        k, train_acc[k], test_acc[k], train_acc[k] - test_acc[k]))

best = max(KS, key=lambda k: test_acc[k])
usable = [k for k in KS if k != 60]
print()
print("k = 1 is right about all sixty training rows, because every training row IS")
print("its own nearest neighbour. There is no fit, no parameter and no learning.")
print("best on held-out rows: k = {} at {:.4f}".format(best, test_acc[best]))
print()
print("k = 60 lets every row vote, so every query gets the same answer.")
print("  distinct predictions across the held-out set:",
      len(set(classify(TRAIN, point, 60) for point, _ in HELD)))
print("  which scores exactly the majority rate, {}/{}".format(majority, len(TRAIN)))

test("at k = 1 every single training row is classified correctly", right(TRAIN, 1), len(TRAIN))
test("yet it is the worst usable setting on rows it has not seen",
     max(usable, key=lambda k: -test_acc[k]), 1)
test("k = 1 also has the widest gap between the two columns",
     max(usable, key=lambda k: train_acc[k] - test_acc[k]), 1)
test("the best k is neither the smallest nor the largest", best not in (1, 60), True)
test("with every row voting there is only one answer left",
     len(set(classify(TRAIN, point, 60) for point, _ in HELD)), 1)
test("and that scores exactly the majority class rate", right(TRAIN, 60), majority)

two overlapping clusters. 60 rows to look things up in, 40 held back.
the majority class holds 31 of the 60 training rows

     k    train      test       gap
     1   1.0000   0.8500   +0.1500
     3   0.9333   0.9250   +0.0083
     5   0.9500   0.9500   +0.0000
     7   0.9500   0.9250   +0.0250
    11   0.9333   0.9000   +0.0333
    17   0.9500   0.9250   +0.0250
    25   0.9500   0.8750   +0.0750
    39   0.9500   0.9250   +0.0250
    60   0.5167   0.4750   +0.0417

k = 1 is right about all sixty training rows, because every training row IS
its own nearest neighbour. There is no fit, no parameter and no learning.
best on held-out rows: k = 5 at 0.9500

k = 60 lets every row vote, so every query gets the same answer.
  distinct predictions across the held-out set: 1
  which scores exactly the majority rate, 31/60
ok   at k = 1 every single training row is classified correctly
ok   yet it is the worst usable setting on rows it has not seen
ok   k = 1 also has the widest gap between the two columns
ok   the best k is neither the smallest nor the largest
ok   with every row voting there is only one answer left
ok   and that scores exactly the majority class rate

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100# The trap. Distance has units, and the largest unit wins.
import math
import random


def distance(a, b):
    return math.sqrt(sum((p - q) ** 2 for p, q in zip(a, b)))


def classify(rows, query, k):
    near = sorted(rows, key=lambda row: distance(row[0], query))[:k]
    votes = {}
    for point, label in near:
        votes[label] = votes.get(label, 0) + 1
    most = max(votes.values())
    tied = [label for label, n in votes.items() if n == most]
    if len(tied) == 1:
        return tied[0]
    for point, label in near:
        if label in tied:
            return label


# Age separates the two classes. Income does not: both classes are drawn from the
# same distribution. But income is measured in pounds.
rng = random.Random(42)
ROWS = []
for _ in range(50):
    ROWS.append(((rng.gauss(30.0, 5.0), rng.gauss(50000.0, 15000.0)), 0))
for _ in range(50):
    ROWS.append(((rng.gauss(45.0, 5.0), rng.gauss(50000.0, 15000.0)), 1))

order = list(range(len(ROWS)))
random.Random(42).shuffle(order)
TRAIN = [ROWS[i] for i in order[:60]]
HELD = [ROWS[i] for i in order[60:]]


def standardiser(rows):
    """Mean and standard deviation per column, from the TRAINING rows only."""
    n = len(rows)
    width = len(rows[0][0])
    means = [sum(p[j] for p, _ in rows) / n for j in range(width)]
    sds = []
    for j in range(width):
        variance = sum((p[j] - means[j]) ** 2 for p, _ in rows) / n
        sds.append(math.sqrt(variance) if variance > 1e-12 else 1.0)
    return means, sds


def rescale(rows, means, sds):
    """Subtract the mean, divide by the deviation. Every column now spans about the same."""
    return [(tuple((p[j] - means[j]) / sds[j] for j in range(len(p))), lab) for p, lab in rows]


MEANS, SDS = standardiser(TRAIN)
S_TRAIN = rescale(TRAIN, MEANS, SDS)
S_HELD = rescale(HELD, MEANS, SDS)

a, b = TRAIN[0][0], TRAIN[1][0]
age_part = (a[0] - b[0]) ** 2
income_part = (a[1] - b[1]) ** 2
print("two rows: age {:.1f}, income {:.0f}   against   age {:.1f}, income {:.0f}".format(
    a[0], a[1], b[0], b[1]))
print("  the age gap contributes    {:14.2f} to the squared distance".format(age_part))
print("  the income gap contributes {:14.2f}".format(income_part))
print("  so age is {:.6f}% of the distance between these two people".format(
    100.0 * age_part / (age_part + income_part)))
print()
print("training means {}, standard deviations {}".format(
    [round(v, 2) for v in MEANS], [round(v, 2) for v in SDS]))
print()
print("     k   raw features   standardised")
raw = {}
scaled = {}
for k in (1, 3, 5, 7, 11):
    raw[k] = sum(1 for p, lab in HELD if classify(TRAIN, p, k) == lab) / len(HELD)
    scaled[k] = sum(1 for p, lab in S_HELD if classify(S_TRAIN, p, k) == lab) / len(S_HELD)
    print("   {:3}   {:12.4f}   {:12.4f}".format(k, raw[k], scaled[k]))

print()
print("Same rows, same labels, same algorithm, two lines of arithmetic between")
print("chance and useful. Nothing raised to tell you which one you had.")

mean_age = sum(p[0] for p, _ in S_TRAIN) / len(S_TRAIN)
sd_age = math.sqrt(sum((p[0] - mean_age) ** 2 for p, _ in S_TRAIN) / len(S_TRAIN))
held_mean_age = sum(p[0] for p, _ in S_HELD) / len(S_HELD)

test("income is measured in a unit a thousand times larger than age",
     SDS[1] / SDS[0] > 1000.0, True)
test("so age contributes under a millionth of the raw squared distance",
     age_part / (age_part + income_part) < 1e-6, True)
test("on raw features every k is no better than a coin toss",
     all(raw[k] < 0.6 for k in raw), True)
test("standardising the same rows puts every k above ninety percent",
     all(scaled[k] >= 0.9 for k in scaled), True)
test("the standardised training column has mean zero and deviation one",
     [close(mean_age, 0.0, 1e-12), close(sd_age, 1.0, 1e-12)], [True, True])
test("the held-out rows were shifted by the training statistics, not their own",
     close(held_mean_age, 0.0, 1e-6), False)

two rows: age 27.2, income 56236   against   age 20.8, income 33763
  the age gap contributes             40.55 to the squared distance
  the income gap contributes   505062432.37
  so age is 0.000008% of the distance between these two people

training means [38.24, 49947.61], standard deviations [8.86, 13650.03]

     k   raw features   standardised
     1         0.5750         0.9250
     3         0.5000         0.9250
     5         0.4750         0.9250
     7         0.4500         0.9000
    11         0.4250         0.9250

Same rows, same labels, same algorithm, two lines of arithmetic between
chance and useful. Nothing raised to tell you which one you had.
ok   income is measured in a unit a thousand times larger than age
ok   so age contributes under a millionth of the raw squared distance
ok   on raw features every k is no better than a coin toss
ok   standardising the same rows puts every k above ninety percent
ok   the standardised training column has mean zero and deviation one
ok   the held-out rows were shifted by the training statistics, not their own

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93# k-means: assign, update, repeat. Both halves lower the score, so it must stop.
import math
import random

rng = random.Random(42)
POINTS = []
TRUE = []
for c, (cx, cy) in enumerate(((0.0, 0.0), (6.0, 0.5), (12.0, 0.0))):
    for _ in range(12):
        POINTS.append((rng.gauss(cx, 0.8), rng.gauss(cy, 0.8)))
        TRUE.append(c)               # kept only to grade the answer at the end


def distance(a, b):
    return math.sqrt(sum((p - q) ** 2 for p, q in zip(a, b)))


def assign(points, centres):
    """Every point joins its nearest centre. No label is consulted anywhere."""
    return [min(range(len(centres)), key=lambda c: distance(p, centres[c])) for p in points]


def update(points, labels, k):
    """Every centre moves to the mean of its own members."""
    fresh = []
    for c in range(k):
        members = [p for p, lab in zip(points, labels) if lab == c]
        if not members:
            raise ValueError("cluster {} lost every member".format(c))
        fresh.append(tuple(sum(p[j] for p in members) / len(members)
                           for j in range(len(points[0]))))
    return fresh


def inertia(points, labels, centres):
    """The only score k-means has: total squared distance from each point to its centre."""
    return sum(distance(p, centres[lab]) ** 2 for p, lab in zip(points, labels))


def kmeans(points, k, seed, max_steps=50):
    centres = random.Random(seed).sample(points, k)     # k of the points, chosen at random
    labels = assign(points, centres)
    trail = [("start ", inertia(points, labels, centres))]
    for step in range(max_steps):
        centres = update(points, labels, k)
        trail.append(("update", inertia(points, labels, centres)))
        fresh = assign(points, centres)
        trail.append(("assign", inertia(points, fresh, centres)))
        settled = fresh == labels
        labels = fresh
        if settled:
            return centres, labels, trail, step + 1
    raise ValueError("did not settle in {} passes".format(max_steps))


centres, labels, trail, passes = kmeans(POINTS, 3, seed=1)

print("36 unlabelled points. k-means is told k = 3 and nothing else.")
print()
print("   phase    inertia")
for phase, value in trail:
    print("   {}   {:9.4f}".format(phase, value))
print()
print("settled after {} full passes, inertia {:.4f}".format(passes, trail[-1][1]))
print("centres found :", [tuple(round(v, 3) for v in c) for c in centres])
print("centres used to make the data: (0, 0), (6, 0.5), (12, 0)")
print()
print("   found cluster   size   which real clusters its members came from")
for c in range(3):
    members = [TRUE[i] for i in range(len(POINTS)) if labels[i] == c]
    print("   {:13}   {:4}   {}".format(c, len(members), sorted(set(members))))

values = [v for _, v in trail]
try:
    update(POINTS, [0] * len(POINTS), 3)
    refused = False
except ValueError as error:
    refused = True
    print()
    print("if a centre ends up with no members at all:", error)

test("the score never rises, at either half of any step",
     all(values[i + 1] <= values[i] for i in range(len(values) - 1)), True)
test("it falls strictly until the final pass, which changes nothing",
     [all(values[i + 1] < values[i] for i in range(len(values) - 2)),
      close(values[-1], values[-2])], [True, True])
test("the last assignment moved no point, which is why it stopped", assign(POINTS, centres), labels)
test("it found three clusters of twelve, having never seen a label",
     sorted(labels.count(c) for c in range(3)), [12, 12, 12])
test("and each one is exactly one of the real clusters",
     [len(set(TRUE[i] for i in range(len(POINTS)) if labels[i] == c)) for c in range(3)],
     [1, 1, 1])
test("an empty cluster is refused rather than dividing by zero", refused, True)

36 unlabelled points. k-means is told k = 3 and nothing else.

   phase    inertia
   start     885.8303
   update    240.7501
   assign    228.2672
   update    223.5904
   assign    223.1046
   update    222.5889
   assign    220.3537
   update    210.8169
   assign    195.5938
   update    161.7129
   assign     89.6280
   update     34.6485
   assign     34.6485

settled after 6 full passes, inertia 34.6485
centres found : [(6.126, 0.285), (0.026, -0.041), (12.024, 0.3)]
centres used to make the data: (0, 0), (6, 0.5), (12, 0)

   found cluster   size   which real clusters its members came from
               0     12   [1]
               1     12   [0]
               2     12   [2]

if a centre ends up with no members at all: cluster 1 lost every member
ok   the score never rises, at either half of any step
ok   it falls strictly until the final pass, which changes nothing
ok   the last assignment moved no point, which is why it stopped
ok   it found three clusters of twelve, having never seen a label
ok   and each one is exactly one of the real clusters
ok   an empty cluster is refused rather than dividing by zero

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100# The real lesson: same data, same code, worse answer, no complaint.
import math
import random

rng = random.Random(42)
POINTS = []
TRUE = []
for c, (cx, cy) in enumerate(((0.0, 0.0), (6.0, 0.5), (12.0, 0.0))):
    for _ in range(12):
        POINTS.append((rng.gauss(cx, 0.8), rng.gauss(cy, 0.8)))
        TRUE.append(c)


def distance(a, b):
    return math.sqrt(sum((p - q) ** 2 for p, q in zip(a, b)))


def assign(points, centres):
    return [min(range(len(centres)), key=lambda c: distance(p, centres[c])) for p in points]


def update(points, labels, k):
    fresh = []
    for c in range(k):
        members = [p for p, lab in zip(points, labels) if lab == c]
        if not members:
            raise ValueError("cluster {} lost every member".format(c))
        fresh.append(tuple(sum(p[j] for p in members) / len(members)
                           for j in range(len(points[0]))))
    return fresh


def inertia(points, labels, centres):
    return sum(distance(p, centres[lab]) ** 2 for p, lab in zip(points, labels))


def kmeans(points, k, seed, max_steps=50):
    centres = random.Random(seed).sample(points, k)
    labels = assign(points, centres)
    for step in range(max_steps):
        centres = update(points, labels, k)
        fresh = assign(points, centres)
        settled = fresh == labels
        labels = fresh
        if settled:
            return centres, labels, inertia(points, labels, centres), step + 1
    raise ValueError("did not settle")


runs = dict((seed, kmeans(POINTS, 3, seed)) for seed in list(range(1, 21)) + [42])

print("the same 36 points, k = 3, and only the starting centres different")
print()
print("   seed   passes    inertia   cluster sizes")
for seed in (1, 2, 11, 12, 17, 42):
    _, labels, score, passes = runs[seed]
    print("   {:4}   {:6}   {:9.4f}   {}".format(
        seed, passes, score, sorted(labels.count(c) for c in range(3))))

good = runs[1]
unlucky = runs[42]
print()
print("seed 42 settled in {} passes and stopped, {:.2f} times worse than seed 1.".format(
    unlucky[3], unlucky[2] / good[2]))
print("its centres:", [tuple(round(v, 2) for v in c) for c in unlucky[0]])
print()
print("   found cluster   size   real clusters mixed into it")
for c in range(3):
    members = [TRUE[i] for i in range(len(POINTS)) if unlucky[1][i] == c]
    print("   {:13}   {:4}   {}".format(c, len(members), sorted(set(members))))
print()
print("Two of its centres landed inside the left-hand cluster and carved it up")
print("between them. The third was left covering both of the others. Every point")
print("is already in its nearest cluster, so no point wants to move, so the")
print("algorithm is finished. It is wrong, it is stable, and nothing raised.")

# The fix is not cleverness, it is repetition: run it from several starts and keep
# the run with the lowest inertia.
scores = dict((seed, runs[seed][2]) for seed in range(1, 21))
best_seed = min(scores, key=lambda s: scores[s])
bad = [s for s in scores if scores[s] > scores[best_seed] + 1e-9]
print()
print("twenty runs, each from a different start:")
print("   distinct answers reached:", sorted(set(round(v, 4) for v in scores.values())))
print("   {} of the 20 landed above the best one".format(len(bad)))
print("   keeping the lowest inertia gives {:.4f}".format(scores[best_seed]))

still = assign(POINTS, update(POINTS, unlucky[1], 3))

test("the unlucky start reached a genuine stopping point, not a half-finished one",
     still, unlucky[1])
test("its score is six times worse than the good run", unlucky[2] / good[2] > 6.0, True)
test("it merged two real clusters together and split a third",
     sorted(unlucky[1].count(c) for c in range(3)) == [12, 12, 12], False)
test("it also stopped sooner, so a pass count says nothing about quality",
     unlucky[3] < good[3], True)
test("inertia is the only signal that anything went wrong", unlucky[2] > good[2], True)
test("keeping the best of twenty restarts finds the good answer",
     close(scores[best_seed], good[2]), True)
test("but four of the twenty were bad, so a single run is a gamble", len(bad) >= 4, True)

the same 36 points, k = 3, and only the starting centres different

   seed   passes    inertia   cluster sizes
      1        6     34.6485   [12, 12, 12]
      2        3    224.9202   [5, 8, 23]
     11        3    251.2943   [3, 9, 24]
     12        1    254.4851   [2, 10, 24]
     17        1    251.2943   [3, 9, 24]
     42        2    224.5397   [3, 10, 23]

seed 42 settled in 2 passes and stopped, 6.48 times worse than seed 1.
its centres: [(-0.09, -1.06), (0.45, 0.27), (9.3, 0.3)]

   found cluster   size   real clusters mixed into it
               0      3   [0]
               1     10   [0, 1]
               2     23   [1, 2]

Two of its centres landed inside the left-hand cluster and carved it up
between them. The third was left covering both of the others. Every point
is already in its nearest cluster, so no point wants to move, so the
algorithm is finished. It is wrong, it is stable, and nothing raised.

twenty runs, each from a different start:
   distinct answers reached: [34.6485, 224.9202, 251.2943, 254.4851]
   4 of the 20 landed above the best one
   keeping the lowest inertia gives 34.6485
ok   the unlucky start reached a genuine stopping point, not a half-finished one
ok   its score is six times worse than the good run
ok   it merged two real clusters together and split a third
ok   it also stopped sooner, so a pass count says nothing about quality
ok   inertia is the only signal that anything went wrong
ok   keeping the best of twenty restarts finds the good answer
ok   but four of the twenty were bad, so a single run is a gamble