Logistic Regression

In three lessons the target was a continuous number and the model said that number directly. The through-line problem of the course, though, is a classification: whether a subscriber carries illicit use or fault suspicion, and the baseline there is 0.7579. The same linear plane cannot be used directly on this problem. The output of the plane stretches from minus infinity to plus infinity; the label is only zero and one. Saying “-3.4 suspicious” of a subscriber is meaningless.

The bridge is an S-shaped function that squeezes the output of the plane between zero and one. The squeezed number is read as a probability, and in this form the model takes the name logistic regression — it is a regression procedure, because it predicts a continuous number; it becomes a classifier only after a threshold is applied to that number. The measure of this lesson lies exactly here: what the output carries before the threshold is applied. Choosing the threshold is a separate job and is left to the evaluation-and-tuning topic.

• LM27. The set, the split and the seed are as in the previous three lessons. The target is the suspicious label and the baseline is saying the most frequent class: 0.7579 on the test set.
• LM28. The feature set is the twelve columns left by the Data Preparation and Feature Engineering course and stays fixed throughout the course: mean consumption, volatility, household size, satisfaction, period count, zero-reading indicator, per-person consumption and five region indicators. The columns that the regression lessons left out are used here, because the target is no longer consumption but a label.
• LM29. The columns are scaled with the mean and spread learned from the train set. The Preparation Steps topic of the Data Preparation and Feature Engineering course had measured that this step is mandatory in a gradient-based model; here it is assumed.
• LM30. All five region indicators are used. In the regression lessons one had been left out as the reference; since the gradient-based solution does not solve a singular system, it is not needed here.
• LM31. The solution is gradient descent: the weights are updated step by step in the direction of the slope of the loss surface. The hyperparameters are the learning rate and the number of epochs.
• LM32. The number of candidates tried is 9: three learning rates, three epoch counts.
• LM33. The measure is accuracy and the test set is 252 subscribers; one subscriber is worth 0.0040. What accuracy hides and by which criterion it should be replaced are left to evaluation-and-tuning/01–02.
• LM34. The threshold sweep in this lesson does not make a choice; it is printed to show that a single model is independent of the threshold. Choosing the operating point by a cost calculation is the subject of evaluation-and-tuning/03.
• LM35. Whether the stated probability sits on the realized rate is a separate question and is left to evaluation-and-tuning/04.

From the Plane to a Probability

The setup generates the set, fixes the twelve columns and prints two rates.

# logistic.py — MODEL: the same FICTIONAL subscriber table, the same seed and the same split;
# in this lesson the target is the illicit use or fault suspicion label.
import math
import statistics

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
REGION = [("north", 0.28, 21, 0.00), ("south", 0.22, 17, -0.10),
          ("east", 0.18, 26, 0.30), ("west", 0.14, 14, -0.05),
          ("central", 0.18, 23, 0.05)]


def rng(t):
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

    def next_value():
        s[0] = (s[0] * 1664525 + 1013904223) & M32
        return s[0] / 4294967296
    return next_value


def pick(u, w):
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


SUBSCRIBERS, DATA = [], []
for i in range(RAW):
    r = rng(SEED + i)
    b = REGION[pick(r(), [x[1] for x in REGION])]
    household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05])
    rating = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    if r() >= 0.046:
        SUBSCRIBERS.append({"id": 10001 + i, "region": b, "household": household,
                            "satisfaction": rating})
for k in SUBSCRIBERS:
    r, v = rng(SEED + 7000 + k["id"]), []
    for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        v.append(0.0 if r() < 0.038 else math.floor(
            k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - n * 0.05) * 100 + 0.5) / 100)
        r()
    if not v:
        continue
    mean, vol, r = sum(v) / len(v), max(v) - min(v), rng(SEED + 51000 + k["id"])
    z = (0.052 * (mean - 20) + 0.85 * (min(v) == 0.0) + 0.026 * vol
         + 0.24 * (k["household"] >= 5) - 0.20 * (k["satisfaction"] >= 4) + k["region"][3]
         + 0.9 * (k["id"] - 10001) / 1399 + (r() + r() + r() - 1.5) * 1.30)
    x = {"household": k["household"], "satisfaction": k["satisfaction"], "period": len(v),
         "mean_consumption": round(mean, 2), "volatility": round(vol, 2),
         "zero_readings": int(min(v) == 0.0), "suspicious": int(z > 1.35)}
    x["per_person"] = round(x["mean_consumption"] / x["household"], 3)
    for b in REGION:
        x["area_" + b[0]] = int(k["region"][0] == b[0])
    DATA.append(x)


def shuffle(data, seed):
    r, s = rng(seed), list(range(len(data)))
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    return [data[i] for i in s]


K = shuffle(DATA, SEED + 90000)
TRAIN, VAL, TEST = K[:756], K[756:1008], K[1008:]
FIELDS = ["mean_consumption", "volatility", "household", "satisfaction", "period",
          "zero_readings", "per_person", "area_north", "area_south", "area_east",
          "area_west", "area_central"]
MEANS = {a: statistics.fmean(x[a] for x in TRAIN) for a in FIELDS}
STD_DEVS = {a: max(1e-9, statistics.pstdev([x[a] for x in TRAIN])) for a in FIELDS}
BASELINE = sum(x["suspicious"] == 0 for x in TEST) / len(TEST)


def scaled(x):                              # FIT from train, APPLY to every set
    return [1.0] + [(x[a] - MEANS[a]) / STD_DEVS[a] for a in FIELDS]


print(f"train {len(TRAIN)}, validation {len(VAL)}, test {len(TEST)}, columns {len(FIELDS)}")
print(f"suspicious rate on train {sum(x['suspicious'] for x in TRAIN) / len(TRAIN):.4f}, "
      f"on test {sum(x['suspicious'] for x in TEST) / len(TEST):.4f}")
print(f"baseline (most frequent class) on test {BASELINE:.4f}")

train 756, validation 252, test 252, columns 12
suspicious rate on train 0.2222, on test 0.2421
baseline (most frequent class) on test 0.7579

The squeeze and gradient descent stand together. The squeezed number is the same thing as reading the linear score as the log-odds: when the score is zero the probability is 0.5, as the score grows it goes toward one, as it shrinks toward zero, and it never touches the ends.

def squash(z):                              # S-shaped function, overflow-safe
    if z >= 0:
        return 1.0 / (1.0 + math.exp(-z)) if z < 700 else 1.0
    e = math.exp(z) if z > -700 else 0.0
    return e / (1.0 + e)


X = [scaled(x) for x in TRAIN]
Y = [x["suspicious"] for x in TRAIN]
P = len(X[0])


def logistic(rate, epochs):                 # MODEL: logistic regression
    w = [0.0] * P
    for _ in range(epochs):
        g = [0.0] * P
        for i in range(len(X)):
            h = squash(sum(w[j] * X[i][j] for j in range(P))) - Y[i]
            for j in range(P):
                g[j] += h * X[i][j]
        for j in range(P):
            w[j] -= rate * g[j] / len(X)
    return w


def probability(w, x):
    return squash(sum(a * b for a, b in zip(w, scaled(x))))


def accuracy(w, s, threshold=0.5):
    return sum((probability(w, x) > threshold) == x["suspicious"] for x in s) / len(s)


RESULTS = []
print(f"{'learning rate':>13} {'epochs':>6} {'train':>9} {'validation':>10} {'test':>8}")
for rate in (0.1, 0.3, 1.0):
    for epochs in (100, 300, 600):
        w = logistic(rate, epochs)
        RESULTS.append((rate, epochs, accuracy(w, VAL), accuracy(w, TEST), w))
        print(f"{rate:>13} {epochs:>6} {accuracy(w, TRAIN):>9.4f} "
              f"{RESULTS[-1][2]:>10.4f} {RESULTS[-1][3]:>8.4f}")
best_val = max(RESULTS, key=lambda t: t[2])
best_test = max(RESULTS, key=lambda t: t[3])
print(f"\nbaseline {BASELINE:.4f}, candidates tried {len(RESULTS)}")
print(f"setting chosen on validation ({best_val[0]}, {best_val[1]}): validation {best_val[2]:.4f}, "
      f"test {best_val[3]:.4f}")
print(f"setting if chosen on test ({best_test[0]}, {best_test[1]}): validation {best_test[2]:.4f}, "
      f"test {best_test[3]:.4f}")
print(f"difference between the two choices {best_test[3] - best_val[3]:.4f}, that is "
      f"{round((best_test[3] - best_val[3]) * len(TEST))} subscribers")
print(f"candidates with the same highest number on validation "
      f"{sum(1 for t in RESULTS if abs(t[2] - best_val[2]) < 1e-12)}")
W = best_val[4]

learning rate epochs     train validation     test
          0.1    100    0.8175     0.8056   0.8016
          0.1    300    0.8201     0.8373   0.7937
          0.1    600    0.8188     0.8373   0.8016
          0.3    100    0.8201     0.8373   0.7937
          0.3    300    0.8175     0.8333   0.8016
          0.3    600    0.8148     0.8333   0.8016
          1.0    100    0.8161     0.8333   0.8016
          1.0    300    0.8148     0.8373   0.8016
          1.0    600    0.8148     0.8373   0.8016

baseline 0.7579, candidates tried 9
setting chosen on validation (0.1, 300): validation 0.8373, test 0.7937
setting if chosen on test (0.1, 100): validation 0.8056, test 0.8016
difference between the two choices 0.0079, that is 2 subscribers
candidates with the same highest number on validation 5

The triple is complete: baseline 0.7579, the family’s number 0.7937, candidates tried 9. The difference is 0.0357, that is, nine subscribers on the test set. The test column of the nine candidates takes only two distinct values — 0.7937 and 0.8016 — and the gap between them is two subscribers. The learning rate and the epoch count substitute for each other in this problem: six hundred epochs with a small rate and three hundred with a large one arrive at the same place.

The last line shows the rule of the course from an unexpected place. Five of the nine candidates give exactly the same number on the validation set, 0.8373. The validation set is 252 rows and measures in steps of 0.0040; it is ordinary for five different settings to fall on the same digit. But those five candidates do not give the same number on the test set. The rule that breaks the tie — taking the first place in the table — is not a property of the model, and it is worth two subscribers.

The Output Is Continuous Before the Threshold

The chosen model is a single setting, but it does not give a single number. What it gives is 252 probabilities, and extracting a class from those probabilities is a separate decision.

tied = [t for t in RESULTS if abs(t[2] - best_val[2]) < 1e-12]
print(f"the {len(tied)} candidates tied on validation, test numbers: "
      f"{sorted({round(t[3], 4) for t in tied})}")
p = sorted((probability(W, x), x["suspicious"]) for x in TEST)
print(f"distinct probability values on test {len({round(q, 9) for q, _ in p})}/{len(TEST)}, "
      f"lowest {p[0][0]:.4f}, highest {p[-1][0]:.4f}\n")
print(f"{'probability range':<18} {'subscribers':>11} {'suspicious':>10} {'suspicious rate':>15}")
for k in range(5):
    a, b = k / 5, (k + 1) / 5
    d = [q for q in p if a <= q[0] < b or (k == 4 and q[0] >= b)]
    g = sum(t[1] for t in d)
    print(f"[{a:.1f}, {b:.1f})".ljust(18)
          + f" {len(d):>11} {g:>10} {g / len(d) if d else 0.0:>15.4f}")
print(f"\n{'threshold':>9} {'accuracy':>9} {'flagged':>12} {'caught':>10} "
      f"{'missed':>10}")
for e in (0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80):
    flagged = sum(1 for q, _ in p if q > e)
    caught = sum(1 for q, y in p if q > e and y)
    print(f"{e:>9.2f} {accuracy(W, TEST, e):>9.4f} {flagged:>12} {caught:>10} "
          f"{sum(t[1] for t in p) - caught:>10}")
print(f"\nsubscribers who change class between 0.40 and 0.60 "
      f"{sum(1 for q, _ in p if 0.40 < q <= 0.60)}")
print(f"threshold fixed at 0.5: accuracy {accuracy(W, TEST):.4f}, "
      f"baseline {BASELINE:.4f}, best threshold at 0.60: {accuracy(W, TEST, 0.60):.4f}")

the 5 candidates tied on validation, test numbers: [0.7937, 0.8016]
distinct probability values on test 252/252, lowest 0.0196, highest 0.9838

probability range  subscribers suspicious suspicious rate
[0.0, 0.2)                 150         15          0.1000
[0.2, 0.4)                  52         18          0.3462
[0.4, 0.6)                  27         11          0.4074
[0.6, 0.8)                  18         12          0.6667
[0.8, 1.0)                   5          5          1.0000

threshold  accuracy      flagged     caught     missed
     0.20    0.7183          102         46         15
     0.30    0.7659           68         35         26
     0.40    0.7817           50         28         33
     0.50    0.7937           31         20         41
     0.60    0.8016           23         17         44
     0.70    0.7897           16         12         49
     0.80    0.7778            5          5         56

subscribers who change class between 0.40 and 0.60 27
threshold fixed at 0.5: accuracy 0.7937, baseline 0.7579, best threshold at 0.60: 0.8016

The second line is the plainest finding of the lesson: all 252 of the 252 subscribers receive a different probability. The model does not produce a label, it produces a ranking. The lowest probability is 0.0196 and the highest 0.9838; the squeeze does not touch the ends, because an S-shaped function reaches zero and one only at infinity.

The range table shows that the ranking is not empty. Of the 150 subscribers in the lowest range, 15, that is one in ten, are suspicious; all five subscribers in the highest range are suspicious. The ranking goes in the same direction as the label. Whether the stated probability sits on the realized rate — for example, the rate coming out as 0.4074 in the 0.4–0.6 range — is a separate question and is left to evaluation-and-tuning/04.

The threshold table gives the number that carries the title of the lesson. A single model, a single weight set, a single training — and seven different accuracy numbers: from 0.7183 to 0.8016. At a threshold of 0.20 the model flags 102 subscribers and catches 46 of the 61 suspicious ones; at a threshold of 0.80 only five subscribers are flagged and all five are truly suspicious, but 56 suspicious ones escape. None of these is a different model.

For this reason the sentence “the accuracy of the model is 0.7937” is incomplete. The correct sentence is: it is 0.7937 at a 0.5 threshold. The threshold is not a property of the model but a decision, and if the decision is moved between 0.40 and 0.60, 27 subscribers change class. Which threshold to choose depends on the cost of the illicit use caught against the cost of the wasted field visit; that calculation is the subject of evaluation-and-tuning/03 and is not done here.

Summary

• Logistic regression is not a classifier but a regression procedure that predicts a probability; it becomes a classifier only after a threshold is applied.
• With twelve fixed columns the model gives 0.7937 on the test set; the baseline is 0.7579, the difference 0.0357 and the number of candidates tried 9.
• Five of the nine candidates tie at exactly 0.8373 on the validation set but split on test between 0.7937 and 0.8016; the rule that breaks the tie is not a property of the model and is worth two subscribers.
• All 252 of the 252 subscribers in the test set receive a different probability, the lowest 0.0196, the highest 0.9838. The model produces a ranking, not a label, and the ranking goes in the same direction as the label.
• When the threshold is swept from 0.20 to 0.80, the same model gives seven different accuracies between 0.7183 and 0.8016; the number of flagged subscribers falls from 102 to 5 and 27 subscribers change class between 0.40 and 0.60.

Next Step

The model of this lesson carries twelve coefficients, and those coefficients are readable: each belongs to a column and has a sign and a magnitude. This is the most often cited advantage of linear models. The next lesson treats that advantage not as praise but as a number: how many subscribers is the difference between the held-out number of a model with readable coefficients and the number of the best family in this topic, and what can and cannot be said by looking at a coefficient. The close of the topic is there too.