Local Outlier Factor: Judging a Point by Its Neighbours (Part 3)

This is the third post in the series on model-based outlier detection. The first post covered the basics and robust statistics, and the second introduced Isolation Forest, which is strongest at global anomalies. In this post I shall discuss the Local Outlier Factor (LOF), which judges each point against the density of its own neighbourhood, and build it from scratch before using scikit-learn.

Why "local" matters

Suppose a dataset has a tight cluster and a loose cluster. Points in the loose cluster are naturally far from their neighbours, so a method that only looks at absolute distances (or at how easily a point is separated from the rest) treats the loose cluster's fringe as suspicious. Meanwhile, a point sitting a short distance outside the tight cluster looks harmless in absolute terms, even though it is clearly out of place relative to its surroundings.

LOF (Breunig, Kriegel, Ng and Sander, 2000) fixes this by asking a relative question: is this point's neighbourhood much sparser than the neighbourhoods of its own neighbours?

The definitions

Fix an integer \(k\) (called MinPts in the original paper) and a distance \(d\).

1. k-distance. The distance from a point \(A\) to its \(k\)-th nearest neighbour:

$$k\text{-dist}(A) = d(A, o_k)$$

where \(o_k\) is the \(k\)-th nearest neighbour of \(A\). The k-neighbourhood \(N_k(A)\) contains every other point within that distance, which is \(k\) points unless there are ties.

2. Reachability distance. For two points \(A\) and \(B\):

$$\text{reach-dist}_k(A, B) = \max\{\, k\text{-dist}(B),\; d(A, B) \,\}$$

If \(A\) is far from \(B\), this is just the actual distance. If \(A\) lies inside \(B\)'s own neighbourhood, the distance is rounded up to \(k\text{-dist}(B)\). This smooths out the random fluctuations of very small distances and makes the density estimates more stable.

3. Local reachability density (lrd). The inverse of the average reachability distance from \(A\) to its neighbours:

$$\text{lrd}_k(A) = \left( \frac{1}{|N_k(A)|} \sum_{B \in N_k(A)} \text{reach-dist}_k(A, B) \right)^{-1}$$

A point surrounded by close neighbours has a high density, and an isolated point has a low one.

4. Local Outlier Factor. The average ratio of the neighbours' density to the point's own density:

$$\text{LOF}_k(A) = \frac{1}{|N_k(A)|} \sum_{B \in N_k(A)} \frac{\text{lrd}_k(B)}{\text{lrd}_k(A)}$$

The interpretation is simple.

• \(\text{LOF} \approx 1\): the point is about as dense as its neighbours, so it is a regular point.
• \(\text{LOF} > 1\): the point is sparser than its neighbours, which makes it an outlier candidate. The larger the value, the stronger the evidence.
• \(\text{LOF} < 1\): the point is denser than its neighbours, which is typical of a cluster's core.

Implementing it from scratch

The code below is a direct translation of the four definitions. It uses a full distance matrix, so it is only suitable for small datasets. scikit-learn uses tree-based neighbour search instead.

import numpy as np
from scipy.spatial.distance import cdist

def lof_scores(X, k=20):
    """Local Outlier Factor from scratch (O(n^2) memory, fine for small data)."""
    n = len(X)
    D = cdist(X, X)
    np.fill_diagonal(D, np.inf)                       # a point is not its own neighbour
    k_dist = np.sort(D, axis=1)[:, k - 1]             # k-distance: distance to the k-th neighbour
    nbrs = [np.where(D[i] <= k_dist[i])[0] for i in range(n)]   # N_k(i), ties included

    lrd = np.empty(n)                                 # local reachability density
    for i in range(n):
        reach = np.maximum(k_dist[nbrs[i]], D[i, nbrs[i]])      # reach-dist_k(i, o)
        lrd[i] = 1.0 / reach.mean()

    return np.array([lrd[nbrs[i]].mean() / lrd[i] for i in range(n)])

A dataset where "local" matters

Now build a dataset to show the difference. It has a tight cluster (300 points), a loose cluster (100 points), 4 planted local outliers close to the tight cluster, and 3 planted global outliers far from everything. I constructed this dataset deliberately to illustrate the point, so on real data the gap between methods will usually be smaller. The code also compares our implementation with scikit-learn's LocalOutlierFactor.

from scipy.stats import spearmanr
from sklearn.neighbors import LocalOutlierFactor

rng = np.random.RandomState(7)
dense = rng.normal([0, 0], 0.4, size=(300, 2))        # tight cluster
sparse = rng.normal([6, 6], 1.6, size=(100, 2))       # loose cluster
local_out = np.array([[1.8, 0.2], [-0.3, 1.9], [-1.6, -1.2], [0.6, -1.9]])   # near the dense cluster
global_out = np.array([[-5.0, 8.0], [12.0, -3.0], [10.0, 12.0]])             # far from everything

X = np.vstack([dense, sparse, local_out, global_out])
n_reg = len(dense) + len(sparse)
y = np.r_[np.zeros(n_reg), np.ones(len(local_out) + len(global_out))]        # 1 = planted outlier
is_local = np.r_[np.zeros(n_reg), np.ones(len(local_out)), np.zeros(len(global_out))].astype(bool)
is_global = np.r_[np.zeros(n_reg + len(local_out)), np.ones(len(global_out))].astype(bool)

k = 20
lof_mine = lof_scores(X, k)
lof_sk = LocalOutlierFactor(n_neighbors=k).fit(X)
lof_sklearn = -lof_sk.negative_outlier_factor_        # sklearn stores the NEGATIVE of LOF

print("max |ours - sklearn| = %.2e" % np.abs(lof_mine - lof_sklearn).max())
print("Spearman correlation = %.4f" % spearmanr(lof_mine, lof_sklearn).correlation)
print("LOF, regular points: median = %.2f, 99th percentile = %.2f, max = %.2f"
      % (np.median(lof_mine[y == 0]), np.percentile(lof_mine[y == 0], 99), lof_mine[y == 0].max()))
print("LOF, local outliers :", np.round(lof_mine[is_local], 2))
print("LOF, global outliers:", np.round(lof_mine[is_global], 2))

flagged = LocalOutlierFactor(n_neighbors=k).fit_predict(X) == -1
print("offset_ ('auto') = %.1f; flagged = %d, planted caught = %d/%d, false alarms = %d"
      % (lof_sk.offset_, flagged.sum(), (flagged & (y == 1)).sum(), y.sum(), (flagged & (y == 0)).sum()))

The output shows the following.

• Our implementation matches scikit-learn. The largest difference between the two sets of scores is about \(6\times10^{-9}\), and the rank correlation is 1.0.
• Regular points score near 1, with a median LOF of 1.02, though a few reach 2.73 (the 99th percentile is 2.04).
• Planted outliers score much higher. The four local outliers have LOF values of 3.70, 4.87, 4.26 and 3.95, and the global outliers have 21.32, 5.84 and 3.56.
• The default threshold is crude. With contamination='auto', scikit-learn flags points with LOF above 1.5 (offset_ is −1.5). Here that catches all 7 planted outliers but also flags 23 regular points. Treat the score as a ranking and choose the cutoff from the cost of false alarms, as in the Isolation Forest post.

Note also that the global outlier at (−5, 8) has a lower LOF (3.56) than the local outliers. LOF measures relative density, not distance, so it should not be read as "how far away a point is."

One scikit-learn detail: negative_outlier_factor_ stores the negative of LOF, so higher means more normal. This is why the code uses -lof_sk.negative_outlier_factor_.

LOF against a global method and Isolation Forest

Now compare LOF with two methods that look at the data more globally: the distance to the \(k\)-th nearest neighbour (the raw quantity that LOF normalises) and Isolation Forest. The code ranks all points by each method's score, checks how many of the 7 planted outliers appear in the top 7, and computes the ROC AUC (the probability that a random planted outlier scores higher than a random regular point, so 1.0 is perfect). It then varies \(k\) for LOF.

from sklearn.ensemble import IsolationForest
from sklearn.metrics import roc_auc_score
from sklearn.neighbors import NearestNeighbors

knn_dist = NearestNeighbors(n_neighbors=k + 1).fit(X).kneighbors(X)[0][:, k]   # column 0 is the point itself
iso_score = -IsolationForest(random_state=0).fit(X).score_samples(X)
scores = {"k-NN distance": knn_dist, "Isolation Forest": iso_score, "LOF": lof_mine}

n_planted = int(y.sum())
for name, s in scores.items():
    top = np.argsort(s)[::-1][:n_planted]             # the n_planted highest-scoring points
    print(f"{name:17s} ROC AUC = {roc_auc_score(y, s):.3f} | top {n_planted}: "
          f"{int(is_local[top].sum())}/{is_local.sum()} local, "
          f"{int(is_global[top].sum())}/{is_global.sum()} global")

print()
for kk in [5, 10, 20, 40, 80]:
    s = lof_scores(X, kk)
    top = np.argsort(s)[::-1][:n_planted]
    print(f"LOF k={kk:2d}: ROC AUC = {roc_auc_score(y, s):.3f}, planted in top {n_planted} = {int(y[top].sum())}")

Method	ROC AUC	Local outliers in top 7 (of 4)	Global outliers in top 7 (of 3)
k-NN distance	0.925	0	3
Isolation Forest	0.978	1	3
LOF (k = 20)	1.000	4	3

All three methods find the global outliers. Only LOF finds the local ones. The figure shows why.

import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 2, figsize=(11, 9), sharex=True, sharey=True)

def draw(ax, title, top=None):
    ax.scatter(X[y == 0, 0], X[y == 0, 1], s=10, c="lightgrey", edgecolors="none")
    if top is None:                                    # panel 1: show the planted outliers
        ax.scatter(X[is_local, 0], X[is_local, 1], s=60, c="tab:orange", edgecolors="k", label="planted local outliers")
        ax.scatter(X[is_global, 0], X[is_global, 1], s=60, c="tab:red", edgecolors="k", label="planted global outliers")
    else:                                              # other panels: show the top-scoring points
        hit = np.isin(np.arange(len(X)), top) & (y == 1)
        miss = np.isin(np.arange(len(X)), top) & (y == 0)
        ax.scatter(X[hit, 0], X[hit, 1], s=60, c="tab:green", edgecolors="k", label="top-scored, planted outlier")
        ax.scatter(X[miss, 0], X[miss, 1], s=60, c="tab:blue", marker="X", edgecolors="k", label="top-scored, regular point")
    ax.set_title(title)
    ax.legend(loc="upper left", fontsize=8)

draw(axes[0, 0], "The data and the planted outliers")
for ax, (name, s) in zip([axes[0, 1], axes[1, 0], axes[1, 1]], scores.items()):
    top = np.argsort(s)[::-1][:n_planted]
    draw(ax, f"{name}: top {n_planted} points", top)
for ax in axes[1]:
    ax.set_xlabel("x1")
for ax in axes[:, 0]:
    ax.set_ylabel("x2")
plt.tight_layout()
plt.show()

Figure 1: Top-left: the data, with planted local (orange) and global (red) outliers. The other panels show the 7 highest-scoring points for each method: green points are planted outliers, and blue crosses are regular points.

Both the k-NN distance and Isolation Forest spend part of their top 7 on regular points at the fringe of the loose cluster (the blue crosses), because those points really are far from their neighbours or easy to isolate in absolute terms. They miss the planted points beside the tight cluster for the opposite reason: those points are close to something, so they look unremarkable. LOF compares each point with the density around it, so it ranks all 7 planted points on top.

Choosing k

The output of the last loop in the previous block shows how LOF depends on \(k\).

k	ROC AUC	Planted outliers in top 7 (of 7)
5	0.999	5
10	1.000	7
20	1.000	7
40	1.000	7
80	0.995	5

A very small \(k\) makes the density estimates noisy, and a very large \(k\) blurs the neighbourhood across clusters, so the sweet spot is in the middle. The original paper recommends using at least 10 and, if you can, evaluating a range of values and taking each point's maximum LOF over that range. In practice, try several values of \(k\) and favour points that stay high across the range. Your own data's cluster sizes matter here: \(k\) should be smaller than the smallest group you want to treat as a cluster.

Novelty detection with LOF

By default, LocalOutlierFactor is for outlier detection on the training data itself, and you use fit_predict or negative_outlier_factor_. To score new observations against clean training data, set novelty=True. Then predict, decision_function and score_samples become available, and fit_predict no longer is.

X_clean = np.vstack([dense, sparse])                  # training data assumed clean
lof_nov = LocalOutlierFactor(n_neighbors=20, novelty=True).fit(X_clean)

X_new = np.array([[0.1, -0.1],     # typical point in the dense cluster
                  [6.5, 5.0],      # typical point in the sparse cluster
                  [1.7, 0.1],      # just outside the dense cluster
                  [9.0, -4.0]])    # far from everything
print("predict:", lof_nov.predict(X_new))                       # +1 inlier, -1 outlier
print("LOF    :", np.round(-lof_nov.score_samples(X_new), 2))   # score_samples = -LOF

The typical points get +1 (LOF about 1.00 and 0.97), the point just outside the tight cluster gets −1 (LOF 3.31), and the distant point gets −1 (LOF 5.87). Remember that score_samples again returns the negative of LOF.

Practical notes

• Scale your features. LOF is built on distances, so a feature with a large range dominates. Standardise first (I skipped it here only because both coordinates are on the same scale).
• Distances degrade in high dimensions, so consider PCA or feature selection before LOF on wide data.
• Duplicates can cause trouble. If more than \(k\) points are identical, the reachability distances become zero and the density blows up. Deduplicate or add tiny jitter.
• LOF is not a probability. A value of 1.5 means something different in every dataset, as the 2.73 maximum among our regular points shows.
• Cost: the neighbour search dominates. scikit-learn's tree-based search handles moderate datasets comfortably, but for very large data Isolation Forest is much cheaper.
• No single winner. Isolation Forest is fast and good at global anomalies, and LOF is better when density varies across the data. Using both and looking at where they disagree is often informative.

Conclusion

In this post I discussed how LOF compares each point's density with that of its neighbours, built it from scratch, checked it against scikit-learn, and showed on a deliberately constructed dataset how it catches local outliers that a global method misses. In the next post I shall discuss the One-Class SVM, which learns a boundary around the normal data in a kernel-induced feature space.

Reference: M. M. Breunig, H.-P. Kriegel, R. T. Ng and J. Sander, "LOF: Identifying Density-Based Local Outliers", Proceedings of the ACM SIGMOD International Conference on Management of Data, 2000.

The code in this post was tested with Python 3 and scikit-learn 1.8.0.

Isolation Forest: Finding Outliers by Isolating Them (Part 2)

This is the second post in the series on model-based outlier detection. The first post discussed what outliers are, why z-scores fail through masking, and how robust statistics help. In this post I shall discuss Isolation Forest, one of the most widely used methods for detecting anomalies, and build it from scratch before using scikit-learn.

The key idea: anomalies are easy to isolate

Most methods profile what "normal" looks like and then flag whatever does not fit. Isolation Forest (Liu, Ting and Zhou, 2008) turns this around. Anomalies are few and different, so if we keep cutting the feature space at random, an anomaly is separated from the rest after only a few cuts, while a point deep inside a dense cluster needs many more.

The number of cuts needed to isolate a point is its path length. Short paths suggest anomalies and long paths suggest regular points. No distance, density or distributional assumption is needed.

Building an isolation tree

An isolation tree (iTree) is built on a sample of \(\psi\) points as follows.

1. Pick a feature \(q\) at random.
2. Pick a split value \(p\) uniformly at random between the minimum and maximum of that feature in the current node.
3. Send points with \(x_q < p\) to the left child and the rest to the right child.
4. Repeat on each child until a node has one point, all its points are identical, or the tree reaches the height limit \(l = \lceil \log_2 \psi \rceil\).

The height limit is deliberate. We only care about points that are isolated early, so there is no need to grow the tree fully.

From path length to anomaly score

For a point \(x\), let \(h(x)\) be the number of edges from the root to the leaf where \(x\) ends up. If that leaf still holds \(n\) unseparated points (because the height limit was reached), we add an adjustment \(c(n)\) for the splits the tree would still have needed. Averaging over all trees gives \(E[h(x)]\).

The adjustment is the average path length of an unsuccessful search in a binary search tree built on \(n\) points:

$$c(n) = 2H(n-1) - \frac{2(n-1)}{n}, \qquad H(i) \approx \ln(i) + \gamma$$

where \(\gamma \approx 0.5772\) is Euler's constant. It also serves as a normaliser, so that scores are comparable across sample sizes. The anomaly score is

$$s(x, \psi) = 2^{-\frac{E[h(x)]}{c(\psi)}}$$

The score behaves as follows.

• As \(E[h(x)] \to 0\), \(s \to 1\). The point is almost certainly an anomaly.
• When \(E[h(x)] = c(\psi)\), \(s = 0.5\). The point is average.
• As \(E[h(x)] \to \psi - 1\), \(s \to 0\). The point is very safely regular.

If every point scores close to 0.5, the data has no distinct anomalies. For the default \(\psi = 256\), \(c(\psi) \approx 10.24\).

Why subsampling helps

Each tree is built on a small random subsample (256 points by default in scikit-learn), not on the whole dataset. This sounds like it should hurt, but it helps. With the full data, anomalies sitting close to normal points (swamping) or clustered together (masking) become hard to isolate. A small sample thins out both effects. It also makes the method fast, since training cost depends on \(\psi\) and the number of trees rather than the dataset size, which is why Isolation Forest scales to large datasets.

Implementing it from scratch

The code below follows the description above. The three classes together are about 50 lines.

import numpy as np

def c(n):
    """Average path length of an unsuccessful search in a binary search tree."""
    if n <= 1:
        return 0.0
    if n == 2:
        return 1.0
    return 2.0 * (np.log(n - 1) + np.euler_gamma) - 2.0 * (n - 1) / n

class IsoTree:
    def __init__(self, X, height, max_height, rng):
        n = X.shape[0]
        self.size = n
        self.leaf = (height >= max_height) or (n <= 1)
        if not self.leaf:
            q = rng.randint(X.shape[1])                  # random feature
            lo, hi = X[:, q].min(), X[:, q].max()
            if lo == hi:                                 # cannot split further
                self.leaf = True
            else:
                self.q, self.p = q, rng.uniform(lo, hi)  # random split value
                mask = X[:, q] < self.p
                self.left = IsoTree(X[mask], height + 1, max_height, rng)
                self.right = IsoTree(X[~mask], height + 1, max_height, rng)

    def path_length(self, x, height=0):
        if self.leaf:
            return height + c(self.size)                 # adjustment for unexpanded leaf
        child = self.left if x[self.q] < self.p else self.right
        return child.path_length(x, height + 1)

class MyIsolationForest:
    def __init__(self, n_trees=100, psi=256, seed=0):
        self.n_trees, self.psi, self.seed = n_trees, psi, seed

    def fit(self, X):
        rng = np.random.RandomState(self.seed)
        self.psi_ = min(self.psi, len(X))
        max_h = int(np.ceil(np.log2(self.psi_)))         # height limit
        self.trees = [
            IsoTree(X[rng.choice(len(X), self.psi_, replace=False)], 0, max_h, rng)
            for _ in range(self.n_trees)
        ]
        return self

    def avg_path_length(self, X):
        return np.array([[t.path_length(x) for t in self.trees] for x in X]).mean(axis=1)

    def anomaly_score(self, X):
        return 2.0 ** (-self.avg_path_length(X) / c(self.psi_))

Comparing with scikit-learn

Now generate a dataset with two clusters of different densities plus 30 uniformly scattered outliers, and run both our forest and scikit-learn's IsolationForest on it. (Run the code above first, since this block reuses c and MyIsolationForest.)

from scipy.stats import spearmanr
from sklearn.datasets import make_blobs
from sklearn.ensemble import IsolationForest
from sklearn.preprocessing import StandardScaler

rng = np.random.RandomState(42)
X_in, _ = make_blobs(n_samples=[300, 300], centers=[[-2, -2], [3, 3]],
                     cluster_std=[0.5, 1.2], random_state=42)
X_out = rng.uniform(-6, 8, size=(30, 2))
X = StandardScaler().fit_transform(np.vstack([X_in, X_out]))
is_out = np.r_[np.zeros(len(X_in), bool), np.ones(len(X_out), bool)]

# 1. Our own forest
mine = MyIsolationForest(n_trees=100, psi=256, seed=0).fit(X)
eh = mine.avg_path_length(X)          # E[h(x)]
s_mine = mine.anomaly_score(X)        # s(x)
print("c(256) = %.2f" % c(256))
print("mean E[h]: regular = %.2f, planted outliers = %.2f"
      % (eh[~is_out].mean(), eh[is_out].mean()))
print("mean s(x): regular = %.3f, planted outliers = %.3f"
      % (s_mine[~is_out].mean(), s_mine[is_out].mean()))

# 2. scikit-learn
iso = IsolationForest(n_estimators=100, random_state=0).fit(X)
s_sk = -iso.score_samples(X)          # score_samples returns the NEGATIVE of s(x)
print("max_samples_ = %d, offset_ = %.1f" % (iso.max_samples_, iso.offset_))
print("Spearman correlation (ours vs sklearn): %.3f" % spearmanr(s_mine, s_sk).correlation)

# 3. Effect of the contamination parameter
def report(name, model):
    flagged = model.predict(X) == -1
    tp = (flagged & is_out).sum()
    fp = (flagged & ~is_out).sum()
    print(f"{name}: flagged={flagged.sum()}, caught={tp}/{is_out.sum()}, false alarms={fp}")

report("contamination='auto'", iso)
for cont in [0.02, 0.05, 0.10, 0.20]:
    report(f"contamination={cont}",
           IsolationForest(n_estimators=100, contamination=cont, random_state=0).fit(X))

# 4. Sensitivity to the random seed (contamination fixed at 0.05)
for seed in range(3):
    report(f"seed={seed}, contamination=0.05",
           IsolationForest(n_estimators=100, contamination=0.05, random_state=seed).fit(X))

The output is as follows.

• The average path length \(E[h(x)]\) is 12.16 for regular points and 6.72 for the planted outliers. The mean anomaly score \(s(x)\) is 0.441 for regular points and 0.641 for the outliers.
• scikit-learn uses 256 samples per tree by default, and with contamination='auto' its decision threshold corresponds to \(s = 0.5\) (offset_ is −0.5).
• The Spearman rank correlation between our scores and scikit-learn's is 0.975. The two forests are built from different random trees and differ in some implementation details, yet they rank the points almost identically.

One detail to remember: scikit-learn's score_samples returns the negative of \(s(x)\), so higher means more normal. This is why the code above uses -iso.score_samples(X). Likewise decision_function is positive for inliers and negative for outliers, and predict returns +1 or −1.

The contamination parameter matters a lot

The last loop in the code above varies the contamination parameter, which is the fraction of points you expect to be anomalies. The table shows what happened with 600 regular points and 30 planted outliers.

contamination	flagged	outliers caught (of 30)	false alarms
'auto' (threshold s = 0.5)	97	27	70
0.02	13	13	0
0.05	32	22	10
0.10	63	27	36
0.20	126	27	99

Three lessons follow from this.

• The default is not a safe choice. With 'auto', the broad cluster's fringe scores above 0.5, so 70 regular points are flagged. A threshold of 0.5 is not a universal rule for "anomalous."
• Raising contamination trades precision for recall. Going from 0.02 to 0.10 catches more outliers, but false alarms jump from 0 to 36.
• The seed matters only slightly. With contamination=0.05, three different seeds all flag 32 points, and catch 22, 23 and 22 of the outliers.

In practice, treat the score as a ranking, look at the distribution of scores, and choose the threshold from the cost of a missed anomaly versus a false alarm. Reviewing the top-ranked cases by hand is often the best calibration you can get without labels.

Visualising the result

import matplotlib.pyplot as plt

m5 = IsolationForest(n_estimators=100, contamination=0.05, random_state=0).fit(X)
flagged = m5.predict(X) == -1

xx, yy = np.meshgrid(np.linspace(X[:, 0].min() - 0.5, X[:, 0].max() + 0.5, 300),
                     np.linspace(X[:, 1].min() - 0.5, X[:, 1].max() + 0.5, 300))
Z = -m5.score_samples(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8))

ax = axes[0]
cs = ax.contourf(xx, yy, Z, levels=20, cmap="viridis_r", alpha=0.85)
plt.colorbar(cs, ax=ax, label="anomaly score s(x)")
ax.scatter(X[~flagged, 0], X[~flagged, 1], s=12, c="white",
           edgecolors="k", linewidths=0.4, label="not flagged")
ax.scatter(X[flagged, 0], X[flagged, 1], s=30, c="red",
           edgecolors="k", linewidths=0.5, label="flagged")
ax.set_title("Isolation Forest anomaly-score map")
ax.set_xlabel("x1 (standardised)")
ax.set_ylabel("x2 (standardised)")
ax.legend(loc="lower right", fontsize=8, framealpha=0.9)

ax = axes[1]
bins = np.linspace(eh.min(), eh.max(), 30)
ax.hist(eh[~is_out], bins=bins, alpha=0.7, density=True, label="regular points")
ax.hist(eh[is_out], bins=bins, alpha=0.7, density=True, color="tab:red",
        label="planted outliers")
ax.axvline(c(mine.psi_), color="k", ls="--", lw=1, label=r"$c(\psi)$ (score = 0.5)")
ax.set_title("Average path length E[h(x)] (our own forest)")
ax.set_xlabel("E[h(x)]")
ax.set_ylabel("density")
ax.legend(fontsize=8)

plt.tight_layout()
plt.show()

 

Figure 1: Left: anomaly-score map from scikit-learn (contamination = 0.05), with flagged points in red. Right: average path lengths from our own forest.

In the left panel, yellow-green regions have a low score (regular) and purple regions a high score (anomalous). In the right panel, most planted outliers have an average path length between about 4 and 9, while most regular points are above 10. The dashed line marks \(c(\psi)\), where the score is 0.5.

Notice the blocky, axis-aligned bands in the score map, especially in the corners. They come from the axis-parallel splits, which can give regions along the axes slightly different scores from regions in other directions. Variants such as the Extended Isolation Forest use randomly oriented splits to reduce this artefact.

Why were three outliers never caught?

In the table, no setting catches more than 27 of the 30 planted outliers. The reason is in the data. The outliers were drawn uniformly from a wide box, so a few of them landed inside the clusters, where they are not really different from the regular points. The code below checks this.

from sklearn.neighbors import NearestNeighbors

m20 = IsolationForest(n_estimators=100, contamination=0.20, random_state=0).fit(X)
missed = np.where(is_out & (m20.predict(X) == 1))[0]       # planted outliers not flagged

nn = NearestNeighbors(n_neighbors=11).fit(X[~is_out])
d_missed = nn.kneighbors(X[missed])[0][:, 10]              # distance to 10th regular neighbour
d_regular = nn.kneighbors(X[~is_out])[0][:, 10]
print("missed outliers:", missed)
print("their 10th-neighbour distances:", np.round(d_missed, 3))
print("regular points: median = %.3f, 90th percentile = %.3f"
      % (np.median(d_regular), np.percentile(d_regular, 90)))

The three missed points have a 10th-nearest-neighbour distance of 0.14 to 0.17, which is within the range of regular points (median 0.115, 90th percentile 0.248). They are not outliers in any useful sense, and no detector should flag them without producing many false alarms. This is a good reminder that "planted outlier" and "detectable outlier" are not the same thing.

Practical notes

• Feature scaling is not required. Each split is drawn within a feature's own range, so rescaling a feature does not change the trees. I standardised the data here only to make plotting convenient.
• Irrelevant features hurt. Random features dilute the splits, so drop noisy features or select them first.
• Local anomalies are a weakness. Isolation Forest is strongest at global anomalies. A point that is unusual only relative to its immediate neighbourhood, especially next to a dense cluster, may score as normal. That is where density-based methods do better.
• Discrete and categorical features need care. Splits are on numeric ranges, so encode categories thoughtfully.
• Novelty detection: you can fit on clean data and call predict or decision_function on new observations.

Conclusion

In this post I discussed how Isolation Forest isolates anomalies with random splits, how path length becomes an anomaly score, and how to build the method from scratch and compare it with scikit-learn. I also showed that the choice of threshold matters at least as much as the algorithm. In the next post I shall discuss the Local Outlier Factor, which compares each point's density with that of its neighbours and so handles clusters of different densities.

Reference: F. T. Liu, K. M. Ting and Z.-H. Zhou, "Isolation Forest", Proceedings of the 8th IEEE International Conference on Data Mining (ICDM), 2008.

The code in this post was tested with Python 3 and scikit-learn 1.8.0.

Outliers, Anomalies and Novelties: Where to Begin (Part 1)

This is the first post in a series on model-based outlier detection. In the coming posts I shall discuss Isolation Forest, Local Outlier Factor and One-Class SVM. Before that, I want to build the foundation: what exactly an outlier is, why the simple rules we learned first fail, and how scikit-learn thinks about the problem.

What is an outlier?

An outlier is an observation that deviates so much from the rest of the data that it seems to come from a different mechanism. Sometimes it is an error (a mistyped price). Sometimes it is the most valuable record in the dataset (a fraudulent transaction, a failing machine, a market crash day). The method does not decide which, so a domain expert should review whatever is flagged.

Three types are worth separating:

• Point anomalies: a single observation is unusual on its own, such as a transaction of ₹50 lakh in an account that normally spends ₹5,000.
• Contextual anomalies: a value is normal in general but unusual in context. A temperature of 35°C is fine in May but strange in January.
• Collective anomalies: each point looks fine, but together they are abnormal, such as a flat line in a sensor signal that should oscillate.

Why the z-score rule fails

The usual rule flags a point when \(|z_i| > 3\), where

$$z_i = \frac{x_i - \bar{x}}{s}$$

There are two problems.

Masking. The mean \(\bar{x}\) and standard deviation \(s\) are themselves computed from the data, including the outliers. Outliers pull the mean toward them and inflate \(s\), so they hide each other.

A hard ceiling. For a sample of size \(n\), no point can have a z-score larger than

$$|z_i| \le \frac{n-1}{\sqrt{n}}$$

(using \(s\) with the \(n-1\) divisor). For \(n = 9\) this is about 2.67, so the "3-sigma rule" can never flag anything in such a sample.

Consider this tiny dataset: 10, 11, 9, 10, 12, 11, 10, 45, 50. Two values are clearly far from the rest, yet their z-scores are only about 1.6 and 1.9.

A robust alternative: median and MAD

Replace the mean with the median and the standard deviation with the median absolute deviation, \(\text{MAD} = \text{median}(|x_i - \tilde{x}|)\). The modified z-score is

$$M_i = \frac{0.6745\,(x_i - \tilde{x})}{\text{MAD}}$$

A common rule (Iglewicz and Hoaglin) flags points with \(|M_i| > 3.5\). The median and MAD barely move when a few outliers are present, so masking largely disappears.

import numpy as np

x = np.array([10, 11, 9, 10, 12, 11, 10, 45, 50])

z = (x - x.mean()) / x.std(ddof=1)
med = np.median(x)
mad = np.median(np.abs(x - med))
mz = 0.6745 * (x - med) / mad

print("z-score :", np.round(z, 2))
print("modified:", np.round(mz, 2))

The z-scores of 45 and 50 come out at 1.60 and 1.91, while their modified z-scores are 22.93 and 26.31. The robust version separates them from the rest without any difficulty.

The multivariate case: Mahalanobis distance

With several features, an observation can be unremarkable on every single axis yet unusual in combination (for example, a very tall and very light person). The Mahalanobis distance accounts for the correlation between features:

$$D^2(\mathbf{x}) = (\mathbf{x}-{\mu})^\top {\Sigma}^{-1} (\mathbf{x}-{\mu})$$

If the data are approximately multivariate normal in \(p\) dimensions, \(D^2\) follows a \(\chi^2_p\) distribution, so a natural cutoff is the 97.5th percentile \(\chi^2_{p,\,0.975}\).

The same masking problem returns, because \({\mu}\) and \({\Sigma}\) are estimated from contaminated data. The fix is the Minimum Covariance Determinant (MCD) estimator. It searches for the subset of \(h\) observations whose sample covariance matrix has the smallest determinant, and computes the mean and covariance from that subset alone. It can tolerate up to nearly 50% contamination, and scikit-learn implements it in MinCovDet.

The code below plants a small cluster of outliers in a correlated dataset and compares the classical and robust estimates.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2
from sklearn.covariance import EmpiricalCovariance, MinCovDet

rng = np.random.RandomState(42)
n_in, n_out = 200, 30
cov = [[1.0, 0.8], [0.8, 1.0]]
X_in = rng.multivariate_normal([0, 0], cov, size=n_in)
X_out = rng.multivariate_normal([4, -3], [[0.3, 0], [0, 0.3]], size=n_out)
X = np.vstack([X_in, X_out])
is_out = np.r_[np.zeros(n_in, bool), np.ones(n_out, bool)]

classical = EmpiricalCovariance().fit(X)
robust = MinCovDet(random_state=0).fit(X)

d2_classical = classical.mahalanobis(X)   # squared distances
d2_robust = robust.mahalanobis(X)
cutoff = chi2.ppf(0.975, df=X.shape[1])

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharex=True, sharey=True)
for ax, (name, d2) in zip(axes, [("Classical", d2_classical),
                                 ("Robust (MCD)", d2_robust)]):
    flagged = d2 > cutoff
    tp = (flagged & is_out).sum()
    fp = (flagged & ~is_out).sum()
    print(f"{name}: flagged={flagged.sum()}, "
          f"true outliers caught={tp}/{n_out}, false alarms={fp}")
    ax.scatter(X[~flagged, 0], X[~flagged, 1], s=18, alpha=0.7, label="not flagged")
    ax.scatter(X[flagged, 0], X[flagged, 1], s=28, c="red", label="flagged")
    ax.set_title(name)
    ax.set_xlabel("x1")
axes[0].set_ylabel("x2")
axes[0].legend()
plt.tight_layout()
plt.show()

Running this gives the following counts.

• Classical: 13 points flagged, 11 of the 30 planted outliers caught, 2 false alarms.
• Robust (MCD): 36 points flagged, all 30 planted outliers caught, 6 false alarms.

The classical estimate misses most of the outliers, because they drag the mean and covariance toward themselves. The robust version catches all of them, at the cost of a few false alarms. This is the usual trade-off with a 97.5% cutoff, and it is worth keeping in mind when choosing a threshold. The exact counts can change slightly with the random seed and library version.

Figure 1: Points flagged (red) by the classical Mahalanobis distance (left) and the robust MCD-based distance (right).

 

Outlier detection versus novelty detection

scikit-learn distinguishes two tasks, and it is important to choose the right one:

• Outlier detection: the training data itself contains outliers, and the goal is to find them. The methods are unsupervised, and you typically call fit_predict.
• Novelty detection: the training data is assumed clean, and the goal is to decide whether a new observation looks like the training data. You call fit on clean data and then predict on new data.

In scikit-learn, predictions use the convention +1 for inliers and −1 for outliers. The estimators I shall cover in this series are IsolationForest, LocalOutlierFactor (which needs novelty=True for the novelty setting) and OneClassSVM (designed for novelty detection and sensitive to contamination in the training set).

Conclusion

In this post I discussed what outliers are, why mean-and-standard-deviation rules fail through masking, how median/MAD and the MCD-based Mahalanobis distance give robust alternatives, and how outlier detection differs from novelty detection. In the next post I shall discuss Isolation Forest, an elegant method that isolates anomalies with random splits and works well without any distributional assumption.

The code in this post was tested with Python 3 and scikit-learn 1.8.0.

Let us understand Logistic Regression (Part 2)

This is the second blog post on Logistic Regression. The first post on this topic is available here. The first post discussed the formulation of Logistic Regression and the estimation of its parameters using the concept of Maximum Likelihood Estimation. In this post I shall discuss some other aspects associated with Logistic Regression. 

Handling Multi-class problems

Logistic regression, by default, is a binary classifier. But, with some modifications, the same can be adopted for a multiclass classification problem also. Most commonly used methods are:

• One versus Rest method
• Multinomial Logistic Regression (with deviance loss)

One versus Rest method

This method is very much intuitive in nature. As the name suggests, one of the classes is considered the positive class and the remaining all are put into the negative class to make it binomial in nature. Hence, if there are \(K\) classes, there will be \(K\) logistic regression models. When an input is given, the input is passed through all the models. Each model predicts the probability of occurrence of the corresponding class. Finally, that class is chosen which has the highest probability among the \(K\) probabilities. Mathematically, if for the \(i^{th}\) data, the class is \(k\), then

$$y^{(k)} = \begin{cases} 1 & \text{if } y = k \\ 0 & \text{otherwise}\end{cases}$$

Correspondingly, the associated probability is calculated as per the equation given below.

$$P(y = k \mid \mathbf{x}) = \sigma(\mathbf{w}_k^\top \mathbf{x} + b_k) = \frac{1}{1 + e^{-(\mathbf{w}_k^\top \mathbf{x} + b_k)}}$$ where \(\mathbf{x}\) is the input, \(\mathbf{w}_k\) is the weight associated with the \(k^{th}\) model and \(b_k\) is the bias (or intercept). The final class is predicted based on the class having the highest probability. $$\hat{y} = \arg\max_{k} P(y = k \mid \mathbf{x})$$

Multinomial Logistic Regression (with deviance loss)

Another method of performing multiclass classification using logistic regression is to use numerical optimisation algorithms on a loss function. In this process, for \(K\) classes, not \(K\) models are built. Instead, all classes are modelled using the softmax function. Suppose \( W\in R^{K\times d}\) is the weight matrix (\(d\) being the dimension of the data) and \(\mathbf{b} \in R^K\) is the bias. In that case, a score function is used for every class as given below $$\mathbf{z}=W\mathbf{x}+\mathbf{b}$$ If \(w_k\) denotes the parameters associated with the \(k^{th}\) class, the score function \(z_k\) will be \(z_k=w_k^T\mathbf{x}+b_k\). This scoring function is converted to probability using the softmax function as shown below.

$$P(y = k \mid \mathbf{x}) = \frac{e^{\mathbf{w}_k^\top \mathbf{x} + b_k}}{\sum_{j=1}^{K} e^{\mathbf{w}_j^\top \mathbf{x} + b_j}}$$ \(P(y=k|\mathbf{x})\) is the probability that the input \(\mathbf{x}\) belong to class \(k\). With these probabilities, negative log loss is calculated as given below. $$L(W, b | \mathbf{X}) = -\frac{1}{n}\sum_{i=1}^n\log P(y_i \mid \mathbf{x_i}) = -\frac{1}{n}\sum_{i=1}^n\log\left( \frac{e^{z_y}}{\sum_{k=1}^{K} e^{z_k}} \right) = -\frac{1}{n}\sum_{i=1}^n \left [z_y - \log\left( \sum_{k=1}^{K} e^{z_k} \right) \right ]$$ where \(n\) is the number of data points and \(z_k=w_k^T\mathbf{x_i}+b_k\). Once the loss function is defined, numerical methods such as L-BFGS, Genetic Algorithm, Conjugate Descent or Newton's Method can be employed easily to estimate the parameters. Sometimes it is required to regularize the parameters to make the model more generalized. Hence, the loss function is modified to accommodate the regularisation term (\(L_2\) regularisation) as given below. $$L(W,b |\mathbf{X})=\frac{1}{n}\sum_{i=1}^nL(W,b|\mathbf{x_i})+\frac{\lambda}{2}||W||^2$$ where \(\lambda\) is the regularisation parameter. Afterward, a numerical method can be applied to this loss function to estimate the regularised parameter. Because it is a multiclass classification task, to get the feature importance for a particular class, the respective weight vector needs to be looked at. The python code below generates a synthetic multiclass classification dataset and applies multinomial logistic regression, and plots feature importance plots for each class as heatmaps.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler

# Generate synthetic dataset
X, y = make_classification(n_samples=1000, n_features=5, n_informative=5,
                           n_redundant=0, n_classes=4, n_clusters_per_class=1,
                           random_state=42)

# Standardize features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Fit One-vs-Rest Logistic Regression
ovr_model = LogisticRegression(multi_class='ovr', solver='lbfgs', max_iter=1000)
ovr_model.fit(X_scaled, y)

# Fit Softmax (Multinomial) Logistic Regression
softmax_model = LogisticRegression(multi_class='multinomial', solver='lbfgs', max_iter=1000)
softmax_model.fit(X_scaled, y)

# Plot heatmap of feature importance for OvR
plt.figure(figsize=(10, 4))
sns.heatmap(ovr_model.coef_, annot=True, cmap='coolwarm', center=0,
            xticklabels=[f'Feature {i}' for i in range(X.shape[1])],
            yticklabels=[f'Class {i}' for i in range(ovr_model.coef_.shape[0])])
plt.title("Feature Importance - One-vs-Rest Logistic Regression")
plt.xlabel("Features")
plt.ylabel("Classes")
plt.tight_layout()
plt.show()

# Plot heatmap of feature importance for Softmax
plt.figure(figsize=(10, 4))
sns.heatmap(softmax_model.coef_, annot=True, cmap='coolwarm', center=0,
            xticklabels=[f'Feature {i}' for i in range(X.shape[1])],
            yticklabels=[f'Class {i}' for i in range(softmax_model.coef_.shape[0])])
plt.title("Feature Importance - Softmax (Multinomial) Logistic Regression")
plt.xlabel("Features")
plt.ylabel("Classes")
plt.tight_layout()
plt.show()

The corresponding heat maps are shown below. 

As per the above plots, if feature 0 and feature 4 are increased by one unit, the log odds for class 2 will increase by 0.92 and 0.77, respectively, if softmax logistic regression is used. In a similar way, other parameters can be explained for different classes as well for both OvR and Softmax Logistic regression.

Conclusion

This post discussed the theoretical aspects associated with multiclass classification using logistic regression. The Python code is also provided for better understanding with feature importance plots. In the next post, I shall discuss the metrics associated with measuring the model's performance. 

Let us understand Logistic Regression (Part 1)

Introduction

Logistic regression, as the name suggests, deals with 1) the logit function and 2) regression. To many people, this algorithm's name is rather confusing because the name contains the term regression, whereas the algorithm actually does classification. In this post, I will provide some details related to this algorithm and explain how the name came into existence.

Logistic regression is essentially binomial in nature, i.e., it can be used for binary classification only. However, with some smart tweaks, the same algorithm can also be used for multi-class classification. I am not going to discuss the multi-class classification in this post. It will be discussed in another post. So let us start with the basics. Below you can see the regression line for two separate cases, i.e., regression and classification. In part (a), there is no problem, but in the case of part (b) (the classification case), there are quite a few problems. 

Problem 1: Assumption of normality of residuals

    For linear regression, an important assumption is the normality of the residuals. Without this assumption, the estimated parameters will become biased. For the classification task, the target variable has only 2 values, i.e., either 0 or 1. Hence, the residuals are more likely to deviate from the normal distribution.

Problem 2: Assumption of homoskedasticity of error variance

    Another important assumption is the homoskedasticity of error variance for the levels of the predictor variables. If this assumption is not valid, the model cannot be generalized for different levels of the input variables. This will also result in the estimation of an incorrect confidence interval. As the target variable has only two values, the errors are also guaranteed to be heteroskedastic in nature.

Problem 3: Predictions are unbounded

    A linear regression line is unbounded and hence not restricted to be bounded between 0 and 1. This problem is, however solvable by adopting a clipping mechanism. That is, if any value goes beyond 1, then it will be brought back to 1 and if any value goes below 0, it will be clipped to 0. This will not solve the problems associated with the assumptions of linear regression.

A typical error analysis plot for the classification data is given below for getting a better understanding. The residuals look normally distributed, but at both ends they deviated. The error variance is also very different from white noise (homoskedastic nature). Hence, it is for sure that linear regression is highly unsuitable for classification tasks.

Finding solution

Probability of success (corresponding to target value 1) can range from 0 to 1, and there can be, theoretically, an infinite number of probabilities between these two values. This is good for a regression task. However, the probability is bound between 0 and 1. As stated in Problem 3, regression outcomes are unbounded, and hence there is a need to extend this binding nature of probability to some other entity which is unbounded. One option is to look for a ratio, i.e., \(\frac{P(success)}{P(failure)}=\frac{P(success)}{1-P(success)}\).  So, as \(P(success) \to 0\), \(\frac{P(success)}{P(failure)} \to \infty\). This is a partial solution because this ratio can never go beyond 0. The solution is to take the logarithm of this ratio so that when \(P(success) <  0.5\), \(log(\frac{P(success)}{P(failure)}) < 0\) and this value can stretch from \(-\infty\) (when \(P(success) = 0\)) to \(\infty\) (when \(P(success) = 1\)). Not only this, the new entity is purely continuous in nature within \((-\infty, \infty)\). This new entity is the \(log\ of\ odd\ ratio\) and the entity \(\frac{P(success)}{P(failure)}\) is called the odd ratio. So, suppose the log of odd ratio is modelled as a linear combination of the features. What we have in our hand is a model: $$log\left (\frac{P(success)}{P(failure)}\right ) = \beta_0 + \beta_1X_1 + \beta_2X_2 + ... +\beta_mX_m$$ The equation looks like a linear regression model which tries to predict log of odd ratio using linear combination of the features. This is why the model is named logistic regression.

The next challenge is to estimate the \(\beta\) values. The problem is, \(P(success)\) is not available in the dataset. The target values cannot be considered as \(P(success)\) because in that case, \(log\left (\frac{P(success)}{P(failure)}\right )\) would be either \(\infty\) or \(-\infty\). This is where some more calculations are required, and depending on the method chosen, \(\beta\) values are estimated. In this post, I shall discuss the process of Maximum Likelihood Estimation (MLE). The reader can get a detailed coverage of this process in my post on MLE.

Maximum Likelihood Estimation of \(\beta\)

Since \(P(success)\) is unknown, a different route is required to find the \(\beta\) values. Let us assume \(Z=\beta_0 + \beta_1X_1 + \beta_2X_2 + ... +\beta_mX_m\). This would mean, \(log\left (\frac{P(success)}{P(failure)}\right )=Z\). For simplicity, let us assume \(P(success) = p\). Then, $$log\left (\frac{p}{1-p}\right ) = \beta_0 + \beta_1X_1 + \beta_2X_2 + ... +\beta_mX_m=Z$$. After a simple calculation, the expression of \(p\) becomes: \(p=\frac{e^z}{1 + e^z}\). \(\frac{e^z}{1 + e^z}\) is the logit function and thus, the probability of success is modelled as the logit function. The next thing is to use the concept of MLE, and for this, an assumption of distribution is needed. To use MLE, it is assumed that the outcomes of the target variable follow a Bernoulli distribution. If \(y_i\) is the outcome of the \(i^{th}\) data point, \(y_i \in \{0,1\}\). Thus, probability mass function associated with the \(i^{th}\) data point is: $$P(X_i|\beta_0,\beta_1,...\beta_m)=p_i^{y_i}(1-p_i)^{1-y_i}$$Assuming that the data points are all following i.i.d condition, the likelihood function is defined as: $$L=\prod_{i=1}^np_i^{y_i}(1-p_i)^{1-y_i}$$For practical purposes and mathematical convenience, we work with log of this likelihood function. The log-likelihood funtion is given as: $$l=log(L)=\sum_{i=1}^n\left [y_ilog(p_i)+(1-y_i)log(1-p_i)\right ]$$The expression within the square brackets is also known as the binary cross entropy, and, it is a popular loss function which can be used for parameter optimisation using different iterative optimisation methods such as Gradient Descent, L-BFGS or Newton's Method (just to name a few!). Clearly, \(l\) is a function of \(p\), which is a function of \(\beta\)s. Hence, to optimise (say) \(\beta_k\) (the parameter associated with the \(k^{th}\) variable) , we need to find the expression of \(\frac{\partial l}{\partial \beta_k}\) and equate it to \(0\). Mathematically, $$\frac{\partial l}{\partial \beta_k}=\sum_{i=1}^{n}y_i \frac{\partial log(p_i)}{\partial \beta_k}+\sum_{i=1}^{n}[(1-y_i) \frac{\partial log(1-p_i)}{\partial \beta_k}]=0$$After putting the value of \(p_i = \frac{e^{z_i}}{1+e^{z_i}}\) in the above expression and doing some further calculations it will be seen that $$\frac{\partial l}{\partial \beta_k} = \sum_{i=1}^n [y_ix_{ik} - p_ix_{ik}]=0$$\(x_{ik}\) is the \(k^{th}\) variable in the \(i^{th}\) data point. There will be \(m+1\) number of such equations where \(m\) is the number of features in the dataset. This final equation is not a closed form equation for \(\beta_k\) and hence, we need to apply numerical methods to estimate the value of \(\beta_k\). Newton Raphson's method is quite popular to find the values of \(\beta_k\). However, other methods are also there for this purpose. Another important thing to note here is that the Hessian matrix need to be -ve definite to ensure that the likelihood value is maximized at the estimated values of \(\beta_k\). If we differentiate \(\frac{\partial l}{\partial \beta_k}\) again with respect to \(\beta_{r}\) (\(\beta\) associated with the \(r^{th}\) variable), $$\frac{\partial l}{\partial \beta_k \partial \beta_{r}}=-\sum_{i=1}^nx_{ik}p_i(1-p_i)x_{ir}$$ In the matrix format the Hessian matrix will take the form $$H = -X^T\Sigma X$$ where \(\Sigma\) is a diagonal matrix having diagonal elements \(p_i(1-p_i) \in [0, 0.25]\). This matrix is positive definite. Hence, as long as the rank of the matrix \(X\) is not rank deficient, \(-X^T\Sigma X\) is -ve definite. Thus, parameters estimated using MLE for logistic regression is going to maximize the likelihood of observing the data. The below graph shows the logit function after fitting the data using logistic regression.

Conclusion

In this post I have discussed why logistic regression got its name and the very working of MLE to estimate the parameters of this model. Since \(Z\) is a linear combination of predictor variables, and the predictor variables are not raised to any other power (other than 1), the discission boundary of logistic regression is linear. Logistic regression falls under linear model. However, with the adaptation of kernels, it is possible to make logistic regression map nonlinearity. I shall discuss that in a separate post. 

EM Algorithm and its usage (Part 2)

EM algorithm is discussed in the previous post related to the tossing of coins. The same algorithm is quite helpful in many other Machine Learning (ML) applications. EM algorithm can be used in missing value imputation, topic modelling (variational EM), anomaly detection, clustering etc. In this post, I am going to discuss how the EM algorithm can be used in clustering.

EM Clustering (Gaussian Mixture Model)

We all know that clustering is an unsupervised machine learning method. Hard clustering methods such as KMeans, Hierarchical and DBSCAN etc., attach a single data point to a single cluster only. EM clustering, on the other hand, associates each data point to all clusters with varying degrees of association. That is why, it is also called a soft clustering technique. EM clustering is essentially a Gaussian Mixture Model which assumes that each cluster is generated from some Gaussian. Thus, EM clustering actually generates a model with likelihood scores. To understand the details of this process, we need to understand the math behind it. 

Gaussian Mixture

As mentioned above, EM Clustering assumes that clusters are essentially samples drawn from some Gaussian and, if there are K numbers of clusters available, the distribution of the data is essentially a mixture of K Gaussian. This is where the concept of multivariate normal distribution comes into play. The generalized equation of a d-dimensional normal distribution is: $$f_i(x|\mu_i, \Sigma_i)=\frac{1}{(2\pi)^{d/2}|\Sigma_i|^{d/2}}exp\left[-\frac{(x-\mu_i)^T\Sigma_i^{-1}(x-\mu_i)}{2}\right]$$ where \(\mu_i\) and \(\Sigma_i\) are the mean vector and covariance matrix respectively. As per the mixture model, the probability density function (pdf) of \(x\) is the summation of pdf of all the K normals. That is, $$f(x) = \Sigma_{i=1}^Kf_i(x)P(C_i)=\Sigma_{i=1}^Kf_i(x|\mu_i, \Sigma_i)P(C_i)$$ \(P(C_i)\) is called the mixture parameter which is a prior probability (probability of occurrence of the \(i^{th}\) cluster itself) and it must satisfy the condition $$\Sigma_{i=1}^KP(C_i)=1$$ Thus, the number of parameters to estimate is a set $$\theta = \{\mu_1, \Sigma_1, P(C_1), ... \mu_K, \Sigma_K, P(C_K)\}$$ It is quite easy to understand that the number of parameters to estimate is quite large compared to (say) KMeans Clustering algorithm. 

Maximum Likelihood Estimation

If the dataset is denoted by \(D\) having \(n\) number of data points, the likelihood of \(\theta\) is given by $$L = \Pi_{j=1}^nf(x_j)$$ It is to be kept in mind that the assumption of i.i.d of \(x_j\) is important for this expression to be valid. Usually, for clustering tasks, we assume that the data points are following i.i.d condition. For mathematical convenience, we work with log-likelihood function and not likelihood function. Log being monotonic in nature, optimum result with log function would result in the optimum value for the original function also. Thus, optimizing log-likelihood is equivalent to optimizing likelihood. The log-likelihood function is given by $$ln(L)=\Sigma_{j=1}^nln(f(x_j)=\Sigma_{j=1}^nln\left(\Sigma_{i=1}^Kf_i(x_j|\mu_i, \Sigma_i)P(C_i)\right)$$ 

Expectation Maximization to rescue

The above equation shows that optimizing the parameters using MLE directly is a prohibitively difficult task. In fact, MLE in such cases becomes computationally intractable when it is to optimize the number of clusters and the best set of parameters. Hence, to solve this problem, the number of clusters is fixed and the parameters are estimated iteratively. The main objective of the clustering process then turns out to be finding out the conditional probability of cluster membership when the data is given. Mathematically, \(P(C_i|x_j)\) is required to be found out. Since this is a conditional probability, Baye's Theorem can be used easily. Applying Baye's theorem, $$P(C_i|x_j)=\frac{P(x_j|C_i)P(C_i)}{\Sigma_{l=1}^KP(x_j|C_l)P(C_l)}$$ If pdf is known, the probability is found out by calculating the area under the curve. If a small interval \(\epsilon>0\) is considered centered at \(x_j\), we can calculate the conditional probability \(p(x_j|C_i)\) as $$p(x_j|C_i) \approx 2\epsilon .f_i(x_j|\mu_i, \Sigma_i)$$ Based on the above equation, the posterior probability \(P(C_i|x_j)\) becomes $$P(C_i|x_j)=\frac{f_i(x_j|\mu_i, \Sigma_i)P(C_i)}{\Sigma_{l=1}^Kf_l(x_j|\mu_l, \Sigma_l)P(C_l)}=w_{ij}$$ \(P(C_i|x_j)\) is considered as the weight or contribution of the point \(x_j\) to cluster \(C_i\). With \(w_{ij}\) extracted, the Expectation Step in EM algorithm is over. In the Maximization step, the same \(w_{ij}\) values are used to re-estimate the parameters as shown below. $$\mu_i=\frac{\Sigma_{j=1}^nw_{ij}x_j}{\Sigma_{j=1}^nw_{ij}}$$ $$\Sigma_i=\frac{\Sigma_{j=1}^nw_{ij}(x-\mu_i)(x-\mu_i)^T}{\Sigma_{j=1}^nw_{ij}}$$ $$P(C_i)=\frac{\Sigma_{j=1}^nw_{ij}}{n}$$

The parameters estimated in the Maximization step are fed back to the Expectation step and the process continues till convergence. With more features, the number of parameters to estimate becomes larger and the process becomes slower. However, EM clustering process gives more insight into the clusters than KMeans clustering process. One such information is the uncertainty measure. Since each data point is associated with each cluster with varying degrees of association (given by the probabilities), uncertainty is defined as $$uncertainty=1- max(probability)$$. This uncertainty is quite helpful to identify the data points which are lying near the boundaries of two clusters. The higher the uncertainty, the higher the chance that with a small push, the data point can be pushed into its neighbouring cluster. This aspect may be quite helpful in customer segmentation analysis with clustering.

EM Clustering with Python

This is the practical part of EM clustering using Python. For the said task, the "Wholesale customer data.csv" is downloaded from this Kaggle webpage. It is a small dataset (440 rows) with 8 features whose description is given below.

• Channel: Channel used to purchase the goods (categorical) 
• Region: The region of the buyer (categorical) 
• Fresh: Amount spent to buy fresh items 
• Milk: Amount spent to buy milk items 
• Grocery: Amount spent to buy grocery items 
• Frozen: Amount spent to buy frozen foods 
• Detergent Papers: Amount spent on detergent papers 
• Delicassen: Amount spent on delicatessen.

The first two features are not used in this analysis as they are categorical in nature and their values are nominal. The remaining six features are used for clustering. The features are standardized with scaling and are made ready for clustering. The Python codes are given below:

import pandas as pd
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler
import numpy as np

data = pd.read_csv('Wholesale customers data.csv')
data.head()

Since Gaussian Mixture Model (GMM) actually builds a model with probabilities associated with each data point, it is possible to extract the likelihood scores for every model and hence the Bayesian Information Criteria (BIC) scores can be calculated for each model. This BIC score can be used effectively to arrive at the optimum number of cluster. The model with lowest BIC score is the best model and the number of components in the model is the optimum number of clusters. The python code and the associated plot shows that the optimum number of clusters is 5

data_scaled = StandardScaler().fit_transform(data.drop(['Channel', 'Region'], axis=1))

# Optimize best cluster solution with BIC Scores
bic = []
for i in range(1, 21):
    gmm = GaussianMixture(n_components=i, n_init=5, random_state=42)
    gmm.fit(data_scaled)
    bic.append(gmm.bic(data_scaled))

plt.plot(range(1, 21), bic)
plt.xlabel('Number of components')
plt.ylabel('BIC')
plt.show()

For visualization purpose, the contour plots of the clusters are created along with uncertainty measures of data. As mentioned above uncertainty shows the level of uncertainty associated with a data point when it associated with a cluster having highest probability of association. The higher the uncertainty, the higher the possibility of the data point getting miss-clustered. This is important for customer segmentation. The below code and the plot shows the same. TSNE is used to plot the high dimensional data into 2 dimensions.

best_gmm = GaussianMixture(n_components=5, n_init=5, random_state=42)
best_gmm.fit(data_scaled)

clusters = best_gmm.predict(data_scaled)

cluster_prob = best_gmm.predict_proba(data_scaled)
uncertainty = 1 - np.max(cluster_prob, axis=1)

tsne = TSNE(n_components=2, random_state=42)
tsne_data = tsne.fit_transform(data.loc[:,'Fresh':'Delicassen'])

# Countour plot creation
xx, yy = np.meshgrid(np.linspace(-20, 30, 100), np.linspace(-20, 30, 100))
grid_data = np.array([xx.ravel(), yy.ravel()]).T
gmm_tsne = GaussianMixture(n_components=5, n_init=5, random_state=42)
gmm_tsne.fit(tsne_data)
zz = gmm_tsne.score_samples(grid_data).reshape(xx.shape)

plt.contourf(xx, yy, zz, alpha=0.3, cmap='tab20', levels=100)

# Plot scatter plot on the contour plot
plt.scatter(tsne_data[:, 0], tsne_data[:, 1], c=clusters, cmap='tab20', 
            s=uncertainty*100)
plt.xlabel('t-SNE 1')
plt.ylabel('t-SNE 2')
plt.show()

In the above plot, the larger the size of the bubble, larger the uncertainty. Data points with higher uncertainty need to be analyzed separately.

Conclusion

EM clustering is essentially a Gaussian Mixture Model which models each cluster as a sample from a specific Gaussian with its own set of parameters. KMeans clustering is, in fact, a special case of EM clustering. EM clustering gives many important insights about the clusters which KMeans clustering does not give directly. However, EM clustering is slower compared to KMeans clustering and it is not so scalable as KMeans clustering. But, a careful marketing analyst would, probably, want to see various aspects of the clusters formed before deciding next marketing moves and in such situations, EM clustering can be quite helpful.

EM Algorithm and its usage (Part 1)

Maximum likelihood estimate is a very powerful and convenient method of estimating parameters of a model. A detail of the same is given in an earlier article. This method is used in several statistical models to estimate the model parameters in an unbiased way. However, a key requirement of this method is to have complete information related to the data, e.g., the distribution to follow and other aspects related to the distribution. But, in real life cases, it is not unusual to have datasets for which some of the information are hidden and yet it is required to estimate parameters of a model to extract the pattern. A simple scenario is mentioned below. 

Two coins problem:

Assume that there are two coins each of which is biased either towards head or tail. Let us assume that coin A is biased towards head and coin B is biased towards tail. At random, one coin is chosen and it is tossed for 10 times. This experiment is done 5 times. The below image shows how easy it is to calculate the level of biasness when all information are known, that is which coin was chosen and how many heads and tails are observed. 

The estimates are nothing but Maximum Likelihood Estimates based on the available data and information. The same problem becomes quite complicated if one information is kept hidden, i.e., which coin was chosen. In such a situation, working with MLE is not advisable as MLE is bad at dealing with such situation. However, with Expectation Maximization (EM) algorithm, it is possible to estimate the parameters in an iterative way. As the name suggests, EM algorithm has an expectation step and a maximization step. The primary objective of the expectation step is calculating probabilities associated with the data points so that using the probabilities, an expression of likelihood can be found out. The likelihood function is then maximized by optimizing the parameters. The parameters then take part is recalculating the probabilities for each data point and the likelihood value is again maximized by optimizing the parameters. Thus, in an iterative manner the parameters keep updating until a breaking point (or threshold limit) is reached.  

The problem of two coins can be solved with EM algorithm when the critical information is mission, i.e., which coin was chosen for tossing. For applying the algorithm, an initialization is needed which is mostly done in a random manner. Afterward, the iterative process commences. The process is shown in the diagram below.

In the above diagram, how probabilities are calculated is not mentioned explicitly. This is mentioned in the below diagram for the first data point where there were 5 heads and 5 tails. $\theta_A^{(0)}=0.6$ and $\theta_B^{(0)}=0.5$. 

After 10 iterations, the values converged to 0.8 for coin A and 0.52 for coin B. Obviously, the number of experiments was less; hence, the iterative algorithm shows both coins as biased towards heads. The same experiment can be replicated with a Python routine. The code is given below:

import numpy as np
import scipy
import scipy

# Generate Data
'''
    Each coin is tossed 100 times and the experiment is done 120 times, 60 times for each coin
'''
coinA = np.random.binomial(n=1, p=0.8, size=(60,100)) # probability of head = 0.8
coinB = np.random.binomial(n=1, p=0.35, size=(60,100)) # probability of head = 0.35
data = np.vstack([coinA, coinB])

index = np.array(range(data.shape[0]))
np.random.shuffle(index)
data = data[index]

# Create frequency table for heads and tails counts
n_heads = data.sum(axis=1)
n_tails = data.shape[1] - n_heads
freq_table = np.hstack((n_heads.reshape(-1,1), n_tails.reshape(-1,1)))

# EM Algorithm starts
# --------------------

# Initialization
thetaA = 0.6
thetaB = 0.3
n_iter = 20

for it in range(n_iter):

    # E-Step starts:
    p_headsA = np.array([scipy.stats.binom.pmf(r, data.shape[1], thetaA) for r in freq_table[:,0]])
    p_headsB = np.array([scipy.stats.binom.pmf(r, data.shape[1], thetaB) for r in freq_table[:,0]])

    P_A = p_headsA/(p_headsA + p_headsB)
    P_B = p_headsB/(p_headsA + p_headsB)

    expected_count_A = P_A.reshape(-1,1) * freq_table
    expected_count_B = P_B.reshape(-1,1) * freq_table

    # E-Step ends

    # M-Step starts
    expected_total_A = expected_count_A.sum(axis=0)
    expected_total_B = expected_count_B.sum(axis=0)

    thetaA = expected_total_A[0]/ expected_total_A.sum()
    thetaB = expected_total_B[0]/ expected_total_B.sum()

    # M-Step ends

print(f"Probability of head in coin A: {thetaA}")
print(f"Probability of head in coin A: {thetaB}")

It will be seen that the values are reaching closer to the input parameters for creating the dataset. This is the EM algorithm in Python for solving the two coin problems.

I hope that the explanations given are understandable.

Comments are welcom.