Anomaly Detection

The first three lessons of this topic did the same job in three different ways: the points were reduced to fewer numbers and each time what was lost was measured — the reconstruction error, the cost of rank, the neighborhood preservation rate. All three measures were read by taking all points together. What a point that stands alone is, or even whether such a point exists, was not asked.

Anomaly detection asks this question, and asks it as an unlabeled problem: no column says that a row is an anomaly, so the definition of “anomaly” comes from the method itself. In such a setup the only way to say which flagging is good is to be able to count how many of the flagged are real. This lesson does that by deliberately planting anomalies into the fictional dataset.

• DR31. The dataset is the fictional subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering, and Supervised Learning courses, seed 20260218, 1,260 rows and five numeric columns. The label column is not read in this course and is not computed here either.
• DR32. The anomalies are deliberately planted into the table. Which row was planted is known only for measurement; no detector sees this list. This is the place of the fiction: in a real dataset this list does not exist and the numbers of this lesson cannot be computed there.
• DR33. anomaly and outlier are kept apart. Outlier is the term of the Data Collection and Preparation course and describes a single value in one column standing at the extreme. Anomaly here is the row as a whole not fitting the dataset; a row with no value at the extreme may be an anomaly.
• DR34. The baseline is the interquartile range, fence and standard score procedure of the Data Collection and Preparation course. The procedure was built there and is not repeated here; because a fixed number of rows must be flagged, instead of the fence the continuous form of the same procedure, the magnitude of the standard score, is called as the ranking criterion.
• DR35. The contamination rate is a hyperparameter: how many rows will be counted as anomalies does not come out of the data. It is swept and the number of each rate is written separately.

The Planted Anomaly

# COMMON — FICTIONAL. The subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering and Supervised Learning courses is rebuilt with the same seed;
# the label column is never read in this course. The anomalies are PLANTED into the table ON PURPOSE.
import math

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


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 v:
        DATA.append({"id": k["id"], "household": k["household"], "satisfaction": k["satisfaction"],
                     "period": len(v), "mean_consumption": round(sum(v) / len(v), 2),
                     "volatility": round(max(v) - min(v), 2)})

# Planting. Three types: A and B are moved to the extreme in a single column, C breaks the
# relationship between columns without going to the extreme in any column (three readings, with almost no difference among them).
r, taken, PLANTED = rng(SEED + 31000), set(), {}


def choose(cond, count):
    s = [i for i in range(len(DATA)) if cond(DATA[i]) and i not in taken]
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    taken.update(s[:count])
    return s[:count]


for i in choose(lambda x: 12 <= x["mean_consumption"] <= 30 and x["period"] >= 2, 24):
    DATA[i]["mean_consumption"], PLANTED[i] = round(DATA[i]["mean_consumption"] * 5.0, 2), "A"
for i in choose(lambda x: x["period"] == 3 and 12 <= x["volatility"] <= 20, 15):
    DATA[i]["volatility"], PLANTED[i] = round(DATA[i]["volatility"] * 5.0, 2), "B"
for i in choose(lambda x: x["period"] == 3 and 22 <= x["mean_consumption"] <= 33, 24):
    DATA[i]["volatility"], PLANTED[i] = round(DATA[i]["mean_consumption"] * 0.02, 2), "C"

FIELDS = ["mean_consumption", "volatility", "household", "satisfaction", "period"]
n, p = len(DATA), len(FIELDS)
MEAN = [sum(x[a] for x in DATA) / n for a in FIELDS]
STD_DEV = [math.sqrt(sum((x[a] - MEAN[j]) ** 2 for x in DATA) / n)
           for j, a in enumerate(FIELDS)]
P = [[(x[a] - MEAN[j]) / STD_DEV[j] for j, a in enumerate(FIELDS)] for x in DATA]

print(f"fictional table {n} rows, {p} columns; planted anomalies {len(PLANTED)} "
      f"({len(PLANTED) / n:.4f})")
for t, name in (("A", "scale shift"), ("B", "volatility spike"),
                ("C", "steady high reading")):
    g = [i for i in PLANTED if PLANTED[i] == t]
    print(f"  type {t} {name:<21} {len(g):>3} rows  "
          f"mean_consumption {min(DATA[i]['mean_consumption'] for i in g):>6.2f}-"
          f"{max(DATA[i]['mean_consumption'] for i in g):<6.2f} "
          f"volatility {min(DATA[i]['volatility'] for i in g):>5.2f}-"
          f"{max(DATA[i]['volatility'] for i in g):.2f}")
print(f"  non-planted rows have mean_consumption at most "
      f"{max(DATA[i]['mean_consumption'] for i in range(n) if i not in PLANTED):.2f}, "
      f"volatility at most "
      f"{max(DATA[i]['volatility'] for i in range(n) if i not in PLANTED):.2f}")

fictional table 1260 rows, 5 columns; planted anomalies 63 (0.0500)
  type A scale shift            24 rows  mean_consumption  62.45-144.40 volatility  1.71-23.22
  type B volatility spike       15 rows  mean_consumption   9.50-35.29  volatility 61.10-98.55
  type C steady high reading    24 rows  mean_consumption  22.09-31.68  volatility  0.44-0.63
  non-planted rows have mean_consumption at most 43.48, volatility at most 51.15

The difference among the three types stands in the last line. The consumption of type A is between 62.45 and 144.40, while the largest among the non-planted is 43.48; the volatility of type B is between 61.10 and 98.55, while the largest among the non-planted is 51.15. Both stand at the extreme in one column and, in the term of the Data Collection and Preparation course, each carries an outlier in that column.

Type C goes to the extreme in no column. Its consumption is between 22.09 and 31.68, in the upper quarter of the dataset but well inside the limit. Its volatility is between 0.44 and 0.63, and there are 187 rows in the dataset with volatility below 1, because in a subscriber with a single reading the volatility is by definition zero. The anomaly is not in the values but in the togetherness of the values: a subscriber read over three periods, with high consumption and almost no difference among the three readings.

Five Detectors, Three Baselines

# Five detectors. The first three are baselines, the last two are multivariate.
def components(k, seed):
    # Covariance matrix and the first k eigenvectors: power iteration and deflation.
    # The component itself was built in the first two lessons of this topic; here it is an input.
    C = [[sum(P[i][a] * P[i][b] for i in range(n)) / n for b in range(p)]
         for a in range(p)]
    rr, found = rng(seed), []
    for _ in range(k):
        v = [rr() - 0.5 for _ in range(p)]
        for _ in range(300):
            w = [sum(C[a][b] * v[b] for b in range(p)) for a in range(p)]
            u = math.sqrt(sum(x * x for x in w)) or 1.0
            v = [x / u for x in w]
        lam = sum(v[a] * sum(C[a][b] * v[b] for b in range(p)) for a in range(p))
        found.append(v)
        C = [[C[a][b] - lam * v[a] * v[b] for b in range(p)] for a in range(p)]
    return found


def reconstruction_error(k):
    V = components(k, SEED + 31500)
    s = []
    for i in range(n):
        pr = [sum(P[i][a] * v[a] for a in range(p)) for v in V]
        rec = [sum(pr[j] * V[j][a] for j in range(k)) for a in range(p)]
        s.append(math.sqrt(sum((P[i][a] - rec[a]) ** 2 for a in range(p))))
    return s


FAR = []
for i in range(n):
    d = sorted(sum((P[i][a] - P[j][a]) ** 2 for a in range(p))
               for j in range(n) if j != i)
    FAR.append(math.sqrt(d[9]))          # distance of the 10th neighbor

SCORES = [("random", "random flagging", None),
          ("single column", "by single column: mean_consumption", [abs(P[i][0]) for i in range(n)]),
          ("column by column", "each column on its own: largest standard score",
           [max(abs(P[i][a]) for a in range(p)) for i in range(n)]),
          ("neighbor", "10th neighbor distance", FAR),
          ("reconstruct 2", "reconstruction error with two components", reconstruction_error(2)),
          ("reconstruct 4", "reconstruction error with four components", reconstruction_error(4))]


def flag(score, m, seed=0):
    if score is None:
        rr, s = rng(seed), list(range(n))
        for i in range(n - 1, 0, -1):
            j = int(rr() * (i + 1))
            s[i], s[j] = s[j], s[i]
        return s[:m]
    return sorted(range(n), key=lambda i: (-score[i], i))[:m]


M = 63
print(f"contamination rate 0.0500 -> {M} rows are flagged; planted 63\n")
print(f"{'detector':<46} {'real':>7} {'precision':>9} "
      f"{'A':>4} {'B':>4} {'C':>4}")
for _, name, score in SCORES:
    if score is None:
        d = [sum(1 for i in flag(None, M, SEED + 32000 + t) if i in PLANTED)
             for t in range(200)]
        g = sum(d) / len(d)
        print(f"{name:<46} {g:>7.2f} {g / M:>9.4f} {'':>4} {'':>4} {'':>4}")
        continue
    s = flag(score, M)
    counts = {t: sum(1 for i in s if PLANTED.get(i) == t) for t in "ABC"}
    print(f"{name:<46} {sum(counts.values()):>7} {sum(counts.values()) / M:>9.4f} "
          f"{counts['A']:>4} {counts['B']:>4} {counts['C']:>4}")

contamination rate 0.0500 -> 63 rows are flagged; planted 63

detector                                          real precision    A    B    C
random flagging                                   3.08    0.0488               
by single column: mean_consumption                  25    0.3968   24    1    0
each column on its own: largest standard score      39    0.6190   24   15    0
10th neighbor distance                              39    0.6190   24   15    0
reconstruction error with two components            40    0.6349   23   15    2
reconstruction error with four components           31    0.4921    1   15   15

The baseline is in the first row: when 63 rows are flagged at random, the mean of 200 draws is 3.08 real anomalies, a precision of 0.0488 — the contamination rate itself. Next to this number no detector counts as having found anything.

The second baseline looks at a single column and gives 0.3968 — eight times the random one. What it found is in the third column: it caught all 24 of type A, only one of type B, none of type C. The third baseline applies the same procedure to the five columns one by one and ranks by the largest standard score; the precision rises to 0.6190 and the missing type B is completed too. Type C is still zero.

The last three rows are multivariate. The distance of the tenth neighbor gives the same number as the column-by-column baseline: 39 real, the same distribution. The reconstruction error with two components is slightly better, 0.6349. The reconstruction error with four components, however, falls behind all of them in precision, 0.4921 — and by itself finds fifteen of type C.

The precision ranking is misleading in this table. The three detectors with the highest precision see none of the 24 anomalies that go to the extreme in no column; the one with the lowest precision sees fifteen. The difference comes not from the sensitivity of the method but from what it looks at. Type A has gone far along the strongest component of the dataset, and four components reconstruct it almost exactly; its residual is close to zero. Type C, on the other hand, is far along no component; it stands in the direction left outside the components, because what it breaks is not a value but the relationship between the number of periods and the volatility and consumption. For the same reason the reconstruction with four components misses type A: the residual it looks for is not a residual A ever produces.

When the Contamination Rate Is Swept

# The contamination rate is a hyperparameter: how many rows will be counted as anomalies
# comes not from the data but from a decision. Five rates are swept.
RATES = [0.02, 0.05, 0.10, 0.15, 0.20]
HEADER = f"{'rate':>5} {'flagged':>7} " + " ".join(f"{a:>16}" for a, _, _ in SCORES)
print("precision: how many of the flagged are planted anomalies\n" + HEADER)
for o in RATES:
    m, row = int(round(o * n)), []
    for _, _, score in SCORES:
        if score is None:
            d = [sum(1 for i in flag(None, m, SEED + 32000 + t) if i in PLANTED)
                 for t in range(200)]
            row.append(sum(d) / len(d) / m)
        else:
            row.append(sum(1 for i in flag(score, m) if i in PLANTED) / m)
    print(f"{o:>5.2f} {m:>7} " + " ".join(f"{x:>16.4f}" for x in row))

print("\ncaught: how many of the 63 anomalies / how many of the 24 of type C\n" + HEADER)
for o in RATES:
    m, row = int(round(o * n)), []
    for _, _, score in SCORES:
        s = flag(score, m, SEED + 32000) if score is None else flag(score, m)
        row.append(f"{sum(1 for i in s if i in PLANTED)}/"
                   f"{sum(1 for i in s if PLANTED.get(i) == 'C')}")
    print(f"{o:>5.2f} {m:>7} " + " ".join(f"{x:>16}" for x in row))

precision: how many of the flagged are planted anomalies
 rate flagged           random    single column column by column         neighbor    reconstruct 2    reconstruct 4
 0.02      25           0.0506           0.9600           1.0000           1.0000           1.0000           0.6000
 0.05      63           0.0488           0.3968           0.6190           0.6190           0.6349           0.4921
 0.10     126           0.0496           0.2143           0.3175           0.3175           0.3333           0.3333
 0.15     189           0.0501           0.1534           0.2116           0.2169           0.2275           0.2328
 0.20     252           0.0501           0.1230           0.1587           0.1706           0.1706           0.1865

caught: how many of the 63 anomalies / how many of the 24 of type C
 rate flagged           random    single column column by column         neighbor    reconstruct 2    reconstruct 4
 0.02      25              1/1             24/0             25/0             25/0             25/0             15/0
 0.05      63              2/1             25/0             39/0             39/0             40/2            31/15
 0.10     126              3/1             27/0             40/1             40/1             42/3            42/24
 0.15     189              6/3             29/1             40/1             41/2             43/4            44/24
 0.20     252              7/3             31/2             40/1             43/4             43/4            47/24

The first row shows how little it takes to raise the precision. When the rate is pulled to 2%, three detectors give 1.0000: all 25 of the 25 flagged rows are real. In the second table the same row says the opposite — only 25 of the 63 anomalies have been found, none of them of type C. Perfect precision is the product of a narrow threshold and has been bought by never seeing 38 anomalies.

As the rate grows, the precision of every detector falls, because the number of flagged rows exceeds the number of planted anomalies and the excess is necessarily wrong. The only thing that stays constant is random flagging: 0.0506, 0.0488, 0.0496, 0.0501, 0.0501. The baseline is independent of the rate, and precision can be read only next to this number.

The second table counts the limit of the univariate procedures. Flagging by a single column carries the caught anomalies only from 24 to 31 even at a rate of 20%; the column-by-column procedure stops at 40 and does not progress at all after 10%. In these two procedures raising the rate does not show what is unseen, it only multiplies the wrongly flagged rows. The reconstruction error with four components, however, has found all 24 of type C at 10%.

At a rate of 20% the best precision is 0.1865: 252 rows flagged, 47 of them real. The remaining 205 rows have been set aside for inspection and none of them holds a planted anomaly. Choosing the flag rate is accepting the cost of these 205 rows; the criterion does not choose the rate, the rate determines the number of the criterion.

Summary

• Of the 63 anomalies deliberately planted into the fictional table, 24 go to the extreme in no column: the anomaly is not in the values but in the togetherness of the values.
• At a contamination rate of 5% the baseline is 0.0488; flagging by a single column gives 0.3968, the procedure that looks at each column on its own 0.6190, the distance of the tenth neighbor 0.6190, and the reconstruction error with four components 0.4921.
• The precision ranking is misleading: the three detectors with the highest precision find none of the 24 multivariate anomalies, the one with the lowest precision finds fifteen; at a rate of 10% it finds all twenty-four.
• When the rate is pulled to 2%, three detectors give a precision of 1.0000 and at the same time catch only 25 of the 63 anomalies; perfect precision is the product of a narrow threshold.
• The precision of random flagging stays between 0.0488 and 0.0506 at all rates; the baseline is independent of the rate and the other numbers are read only next to it.

Next Step

In this lesson the row was handled with five numeric columns and whether a point stands alone was measured. Subscriber records have one more side that is not a number: which events occur together. Is the joint occurrence of two events in the same subscriber a rule, or can it already be explained by how frequent the two are separately. The next lesson builds association rules, counts how many of the high-confidence rules come only from a frequent item, and reads the baseline from the support that random matching produces.