Outlier Management

The previous lesson measured the scarcity of suspicious subscribers in number. This lesson measures how far a few subscribers are in value. In the mean consumption and volatility columns there are subscribers outside the fence computed from the interquartile range of the training set. Some of them are real extremes — they come from a crowded household or from a connection with a high tariff. Some are recording errors produced by stuck meters. In the table the two look the same: a large number.

The interquartile range, the fence and the finding of extreme values were set up in the Data Collection and Preparation course of the Data Analytics curriculum; the measure there was how many rows stayed outside the fence. The question here is different: how much does the clipping decision move the number on the held-out set, and in which model.

• PS41. The stuck meter fault is fictional: the mean consumption of a small set of subscribers is recorded as thirty times higher. It is chosen independently of the label.
• PS42. This fault is different from the three defects in 01 and cannot be caught by their fixes: the value is in the right unit, written correctly and in a single row. It is only wrong.
• PS43. The fence is a fitted object: the quartiles are computed only from the training set, the fence is built with 1.5 times the interquartile range, and the same fence is applied to both sets.
• PS44. Clipping pulls a value outside the fence to the fence. No row is deleted; a prediction is produced for every subscriber in the test set.
• PS45. Two model families are compared: a gradient-based linear model and a depth 6 decision tree. The linear model is trained with standardized columns according to the discipline set up in the previous course; the tree is trained with raw columns, because the contribution of scaling on a tree is zero.

The Two Sources of an Extreme

The setup produces the table, puts the fault in and learns the fence from the training set.

# outlier.py — MODEL. The fictional subscriber table is regenerated with the same seed; the
# consumption of a small set of subscribers is recorded thirty times higher through a FICTIONAL
# stuck meter fault.
import math
import statistics

SEED, CANDIDATES, 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, READINGS = [], {}
for i in range(CANDIDATES):
    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])
    satisfaction = 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": satisfaction})
for k in SUBSCRIBERS:
    r = rng(SEED + 7000 + k["id"])
    for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        READINGS.setdefault(k["id"], []).append(0.0 if r() < 0.038 else math.floor(
            k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - d * 0.05) * 100 + 0.5) / 100)
        r()

DATA = []
for k in SUBSCRIBERS:
    v = READINGS.get(k["id"])
    if v is None:
        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)}
    for b in REGION:
        x["b_" + b[0]] = int(k["region"][0] == b[0])
    DATA.append(x)

fault_rng = rng(SEED + 81000)               # the stuck meter fault is FICTIONAL
for x in DATA:
    x["fault"] = (round(x["mean_consumption"] * 30, 2) if fault_rng() < 0.010
                  else x["mean_consumption"])


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, TEST = K[:756], K[1008:]
BASELINE = sum(x["suspicious"] == 0 for x in TEST) / len(TEST)
NUMERIC = ["mean_consumption", "volatility", "household", "satisfaction", "period"]
FIELDS = NUMERIC + ["zero_readings", "b_north", "b_south", "b_east", "b_west", "b_central"]
EXTREME = ["mean_consumption", "volatility"]


def quantile(v, q):
    s = sorted(v)
    i = q * (len(s) - 1)
    a = int(i)
    return s[a] if a + 1 >= len(s) else s[a] + (i - a) * (s[a + 1] - s[a])


def fit_fence(train, columns, k=1.5):       # FIT: interquartile range and fence
    d = {}
    for a in columns:
        q1, q3 = quantile([x[a] for x in train], 0.25), quantile([x[a] for x in train], 0.75)
        d[a] = (q1 - k * (q3 - q1), q3 + k * (q3 - q1))
    return d


def clip(c, s):                             # APPLY: no row is deleted
    return [dict(x, **{a: min(max(x[a], c[a][0]), c[a][1]) for a in c}) for x in s]


faulty = [x for x in TRAIN if x["fault"] != x["mean_consumption"]]
print(f"stuck meter: {len(faulty)} subscribers in train, "
      f"{sum(1 for x in TEST if x['fault'] != x['mean_consumption'])} in test; "
      f"{sum(x['suspicious'] for x in faulty)} of them are suspicious")
c = fit_fence(TRAIN, EXTREME)
print(f"{'column':<16} {'upper fence':>11} {'largest':>9} {'outside fence, train':>20} "
      f"{'test':>7}")
for a in EXTREME:
    print(f"{a:<16} {c[a][1]:>11.2f} {max(x[a] for x in TRAIN):>9.2f} "
          f"{sum(1 for x in TRAIN if x[a] > c[a][1]):>20} "
          f"{sum(1 for x in TEST if x[a] > c[a][1]):>7}")
u = [x for x in TRAIN if any(x[a] > c[a][1] for a in EXTREME)]
print(f"suspicious rate of the {len(u)} subscribers outside the fence "
      f"{sum(x['suspicious'] for x in u) / len(u):.4f}, baseline {BASELINE:.4f}")

stuck meter: 9 subscribers in train, 1 in test; 2 of them are suspicious
column           upper fence   largest outside fence, train    test
mean_consumption       35.80     42.90                   11       4
volatility             29.41     40.57                   23       7
suspicious rate of the 34 subscribers outside the fence 0.6765, baseline 0.7579

The last line states the cost of the clipping decision in advance. The suspicious rate of the 34 subscribers outside the fence is 0.6765; over the whole training set this rate is 0.2222. Two thirds of the subscribers at the extreme are suspicious; that is, the extreme values are not noise but the most informative region of the set.

The Same Clipping, Two Models

def gini(s):
    p = sum(x["suspicious"] for x in s) / len(s) if s else 0.0
    return 2 * p * (1 - p)


def tree(s, depth, fields, min_size=2):     # MODEL: decision tree, depth 6
    p = sum(x["suspicious"] for x in s) / len(s)
    best = None
    if depth and len(s) >= 2 * min_size and 0.0 < p < 1.0:
        for a in fields:
            d = sorted({x[a] for x in s})
            for v in (d[1:] if len(d) < 10 else
                      [d[int(i * len(d) / 10)] for i in range(1, 10)]):
                left = [x for x in s if x[a] < v]
                right = [x for x in s if x[a] >= v]
                if min(len(left), len(right)) < min_size:
                    continue
                k = gini(s) - (len(left) * gini(left) + len(right) * gini(right)) / len(s)
                if best is None or k > best[0]:
                    best = (k, a, v, left, right)
    if best is None or best[0] <= 1e-9:
        return {"prediction": int(p > 0.5)}
    return {"field": best[1], "threshold": best[2],
            "left": tree(best[3], depth - 1, fields, min_size),
            "right": tree(best[4], depth - 1, fields, min_size)}


def tree_predict(d, x):
    while "prediction" not in d:
        d = d["left"] if x[d["field"]] < d["threshold"] else d["right"]
    return d["prediction"]


def squash(z):
    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)


def linear_model(train, fields, step=0.5, rounds=400):  # MODEL: gradient-based linear model
    w, b, n = [0.0] * len(fields), 0.0, len(train)
    for _ in range(rounds):
        gw, gb = [0.0] * len(fields), 0.0
        for x in train:
            h = squash(sum(w[i] * x[a] for i, a in enumerate(fields)) + b) - x["suspicious"]
            for i, a in enumerate(fields):
                gw[i] += h * x[a]
            gb += h
        for i in range(len(fields)):
            w[i] -= step * gw[i] / n
        b -= step * gb / n
    return w, b


def standardize(train, sets):               # FIT + APPLY (previous lesson)
    p = {a: (statistics.fmean(x[a] for x in train),
             statistics.pstdev([x[a] for x in train])) for a in NUMERIC}
    return [[dict(x, **{a: (x[a] - p[a][0]) / p[a][1] for a in NUMERIC}) for x in k]
            for k in sets]


def build(source, clipping):
    e = [dict(x, mean_consumption=x[source]) for x in TRAIN]
    s = [dict(x, mean_consumption=x[source]) for x in TEST]
    if clipping:
        c = fit_fence(e, EXTREME)
        e, s = clip(c, e), clip(c, s)
    return e, s


RESULT = {}
print(f"{'case':<38} {'gradient-based':>14} {'tree':>8}")
print(f"{'baseline':<38} {BASELINE:>14.4f} {BASELINE:>8.4f}")
for name, source, clipping in (("real extremes only, no clipping", "mean_consumption", False),
                               ("real extremes only, clipped with fence", "mean_consumption", True),
                               ("stuck meter too, no clipping", "fault", False),
                               ("stuck meter too, clipped with fence", "fault", True)):
    e, s = build(source, clipping)
    es, ss = standardize(e, [e, s])
    m = linear_model(es, FIELDS)
    d = sum(int(squash(sum(m[0][i] * x[a] for i, a in enumerate(FIELDS)) + m[1]) > 0.5)
            == x["suspicious"] for x in ss) / len(ss)
    t = tree(e, 6, FIELDS)
    g = sum(tree_predict(t, x) == x["suspicious"] for x in s) / len(s)
    RESULT[name] = (d, g)
    print(f"{name:<38} {d:>14.4f} {g:>8.4f}")

case                                   gradient-based     tree
baseline                                       0.7579   0.7579
real extremes only, no clipping                0.8016   0.7778
real extremes only, clipped with fence         0.8056   0.7500
stuck meter too, no clipping                   0.7698   0.7857
stuck meter too, clipped with fence            0.8016   0.7500

The tree column gives three different numbers in four cases, and what is interesting is that the two clipped rows read the same number. Clipping brings both tables to the same place: when all the values above the fence are gathered at a single point, the faulty record and the real extreme become indistinguishable.

How Many Subscribers the Contribution Is Worth

print(f"{'clipping decision':<20} {'gradient-based':>14} {'subscribers':>11} "
      f"{'tree':>8} {'subscribers':>11}")
for name, a, b in (("real extremes only", "real extremes only, no clipping",
                    "real extremes only, clipped with fence"),
                   ("stuck meter too", "stuck meter too, no clipping",
                    "stuck meter too, clipped with fence")):
    kd, kg = RESULT[b][0] - RESULT[a][0], RESULT[b][1] - RESULT[a][1]
    print(f"{name:<20} {kd:>+14.4f} {round(kd * 252):>11} "
          f"{kg:>+8.4f} {round(kg * 252):>11}")
e, s = build("fault", False)
c = fit_fence(e, EXTREME)
print(f"\nupper fence learned on the faulty table: mean_consumption "
      f"{c['mean_consumption'][1]:.2f}, volatility {c['volatility'][1]:.2f}")
print(f"values clipped: train {sum(1 for x in e for a in EXTREME if x[a] > c[a][1])}, "
      f"test {sum(1 for x in s for a in EXTREME if x[a] > c[a][1])}; "
      f"rows deleted 0")
print(f"cost of the fault to the gradient-based model "
      f"{RESULT['stuck meter too, no clipping'][0] - RESULT['real extremes only, no clipping'][0]:+.4f}, "
      f"to the tree {RESULT['stuck meter too, no clipping'][1] - RESULT['real extremes only, no clipping'][1]:+.4f}")
print(f"features used in six lessons {len(FIELDS)}, new quantities produced 0")

clipping decision    gradient-based subscribers     tree subscribers
real extremes only          +0.0040           1  -0.0278          -7
stuck meter too             +0.0317           8  -0.0357          -9

upper fence learned on the faulty table: mean_consumption 36.20, volatility 29.41
values clipped: train 41, test 12; rows deleted 0
cost of the fault to the gradient-based model -0.0317, to the tree +0.0079
features used in six lessons 11, new quantities produced 0

The signs are opposite in both rows. Clipping gives a positive contribution in the gradient-based model and a negative one in the tree, and when the stuck meter is present both directions rise to a readable size: eight subscribers and nine subscribers. The mechanism is in the way the two models handle an extreme value. The gradient-based model sums the error of all rows in every round and shifts the weights according to that sum; a row thirty times away from the mean determines that sum alone and pulls all eleven weights at once. The tree, on the other hand, can set an extreme row aside with a single split and leaves the rest of the set untouched — but clipping removes that split, because the ordering above the fence is erased, and the suspicious rate of 0.6765 of the 34 subscribers remaining there can no longer be turned into a threshold.

The fault row confirms this from the second side. A recording error of nine rows makes the gradient-based model lose 0.0317 and makes the tree gain 0.0079 — that is, for the tree an unmeasured difference. When clipping is applied the gradient-based model returns exactly to 0.8016, the number it would have read if the fault had never existed. The same step on the same data fully undoes the damage in one model, while it lowers the other down to 0.0079 above the baseline.

The procedural side of clipping is also worth noting. The fence was learned only from the training set, the upper fence moved to 36.20 on the faulty table, 41 values were clipped and zero rows were deleted. For an extreme of unknown source there is no single right answer; the only thing that can be measured is how many subscribers that answer is worth in which model.

Summary

• The suspicious rate of the 34 subscribers outside the fence is 0.6765, over the whole set 0.2222; in this table the extreme values are not noise but the most informative region.
• The same clipping decision gives a positive contribution in the gradient-based linear model and a negative one in the decision tree: +0.0040 and −0.0278 on real extremes only, +0.0317 and −0.0357 when the stuck meter record is also present.
• The gradient-based model updates its weights by summing the error of all rows, and a row thirty times away from the mean determines that sum alone; while the tree can separate an extreme row with a single split, clipping removes that split.
• A recording error of nine rows makes the gradient-based model lose 0.0317 and leaves an unmeasured difference of 0.0079 on the tree; clipping brings the gradient-based model back exactly to 0.8016 of the fault-free state.
• The fence was learned only from the training set, 41 values were clipped and no row was deleted; a prediction was produced for every subscriber in the test set.

Next Step

The last line of the last output states the common limit of the six lessons: the number of features used is eleven and the number of new quantities produced is zero. Everything done in the six lessons was changing the way an existing field is read — fixing the unit, learning the scale from training, turning a category into a number, writing the gap itself into a column, changing the weight of rows, pulling the extreme to the fence. In no lesson was the ratio, difference or product of two columns taken; no seasonal indicator was extracted from a date and no length from a text. And nothing was done in the opposite direction either: no column was removed from the table on the grounds that “its contribution was not measured”, and the eleven features stayed as they were throughout six lessons. The next topic takes on these two jobs. First, new columns are derived from domain knowledge and the contribution of each is measured one by one — those that come out at zero are written down too. Then the opposite question is asked: which columns should be removed from the table, how many features do three selection methods choose on the same set, and how much does the selection itself inflate the number on the held-out set.