External Metrics: Adjusted Rand Index and Normalized Mutual Information

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

import numpy as np
from scipy.special import comb
from sklearn.metrics import rand_score
 
def contingency_table(labels_true, labels_pred):
    """
    Construct the contingency table from two label vectors.
    """
    classes = np.unique(labels_true)
    clusters = np.unique(labels_pred)
    
    n_classes = len(classes)
    n_clusters = len(clusters)
    
    # Map labels to indices
    class_to_idx = {c: i for i, c in enumerate(classes)}
    cluster_to_idx = {c: i for i, c in enumerate(clusters)}
    
    # Build contingency table
    table = np.zeros((n_classes, n_clusters), dtype=int)
    for true_label, pred_label in zip(labels_true, labels_pred):
        i = class_to_idx[true_label]
        j = cluster_to_idx[pred_label]
        table[i, j] += 1
    
    return table, classes, clusters
 
 
def rand_index_from_contingency(table):
    """
    Compute Rand Index from contingency table.
    """
    n = table.sum()
    
    # Row and column sums
    a = table.sum(axis=1)  # Class sizes
    b = table.sum(axis=0)  # Cluster sizes
    
    # TP: pairs together in both
    tp = sum(comb(table[i, j], 2) for i in range(table.shape[0]) for j in range(table.shape[1]))
    
    # Sum of C(a_i, 2) for all classes
    sum_a = sum(comb(ai, 2) for ai in a)
    
    # Sum of C(b_j, 2) for all clusters
    sum_b = sum(comb(bj, 2) for bj in b)
    
    # FP and FN
    fp = sum_b - tp
    fn = sum_a - tp
    
    # Total pairs
    total_pairs = comb(n, 2)
    
    # TN
    tn = total_pairs - tp - fp - fn
    
    # Rand Index
    ri = (tp + tn) / total_pairs
    
    return ri, {'TP': tp, 'FP': fp, 'FN': fn, 'TN': tn}
 
 
# Example
labels_true = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2])
labels_pred = np.array([0, 0, 1, 1, 1, 1, 2, 2, 2])
 
table, classes, clusters = contingency_table(labels_true, labels_pred)
print("Contingency Table:")
print(table)
 
ri, counts = rand_index_from_contingency(table)
print(f"
Pair Counts: {counts}")
print(f"Rand Index: {ri:.4f}")
 
# Verify with sklearn
ri_sklearn = rand_score(labels_true, labels_pred)
print(f"Sklearn Rand Index: {ri_sklearn:.4f}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

import numpy as np
from sklearn.metrics import rand_score, adjusted_rand_score
 
def demonstrate_chance_problem():
    """
    Demonstrate that random clusterings get high Rand Index.
    """
    np.random.seed(42)
    n = 1000
    n_classes = 10
    
    # Random ground truth
    labels_true = np.random.randint(0, n_classes, n)
    
    # Multiple random clusterings
    random_scores = []
    for _ in range(100):
        labels_random = np.random.randint(0, n_classes, n)
        ri = rand_score(labels_true, labels_random)
        ari = adjusted_rand_score(labels_true, labels_random)
        random_scores.append((ri, ari))
    
    ri_mean = np.mean([s[0] for s in random_scores])
    ri_std = np.std([s[0] for s in random_scores])
    ari_mean = np.mean([s[1] for s in random_scores])
    ari_std = np.std([s[1] for s in random_scores])
    
    print("Random Clustering vs Random Ground Truth")
    print(f"  Rand Index:          {ri_mean:.4f} ± {ri_std:.4f}")
    print(f"  Adjusted Rand Index: {ari_mean:.4f} ± {ari_std:.4f}")
    print()
    print("Expected values:")
    print(f"  RI for random: ~{1 - 2*n_classes*(n_classes-1)/(n-1):.4f} (theoretical)")
    print(f"  ARI for random: ~0.0 (by design)")
 
demonstrate_chance_problem()

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

import numpy as np
from scipy.special import comb
from sklearn.metrics import adjusted_rand_score
 
def adjusted_rand_index_manual(labels_true, labels_pred):
    """
    Manual implementation of ARI for understanding the formula.
    """
    # Build contingency table
    classes = np.unique(labels_true)
    clusters = np.unique(labels_pred)
    
    class_to_idx = {c: i for i, c in enumerate(classes)}
    cluster_to_idx = {c: i for i, c in enumerate(clusters)}
    
    n = len(labels_true)
    table = np.zeros((len(classes), len(clusters)), dtype=int)
    
    for true_label, pred_label in zip(labels_true, labels_pred):
        i = class_to_idx[true_label]
        j = cluster_to_idx[pred_label]
        table[i, j] += 1
    
    # Marginal sums
    a = table.sum(axis=1)  # Row sums (class sizes)
    b = table.sum(axis=0)  # Column sums (cluster sizes)
    
    # Index: sum of C(n_ij, 2)
    index = sum(comb(table[i, j], 2) 
                for i in range(table.shape[0]) 
                for j in range(table.shape[1]))
    
    # Sum of C(a_i, 2)
    sum_a = sum(comb(ai, 2) for ai in a)
    
    # Sum of C(b_j, 2)
    sum_b = sum(comb(bj, 2) for bj in b)
    
    # Total pairs
    n_pairs = comb(n, 2)
    
    # Expected index under null hypothesis (random permutation)
    expected_index = (sum_a * sum_b) / n_pairs
    
    # Maximum index
    max_index = 0.5 * (sum_a + sum_b)
    
    # ARI
    if max_index == expected_index:
        return 1.0  # Perfect agreement when only one cluster/class
    
    ari = (index - expected_index) / (max_index - expected_index)
    
    return ari, {
        'index': index,
        'expected_index': expected_index,
        'max_index': max_index,
        'sum_a': sum_a,
        'sum_b': sum_b,
        'n_pairs': n_pairs
    }
 
 
# Example with varying agreement levels
examples = [
    ("Perfect match", [0, 0, 0, 1, 1, 1, 2, 2, 2], [0, 0, 0, 1, 1, 1, 2, 2, 2]),
    ("Good match", [0, 0, 0, 1, 1, 1, 2, 2, 2], [0, 0, 1, 1, 1, 1, 2, 2, 2]),
    ("Poor match", [0, 0, 0, 1, 1, 1, 2, 2, 2], [0, 1, 2, 0, 1, 2, 0, 1, 2]),
    ("Permuted labels", [0, 0, 0, 1, 1, 1, 2, 2, 2], [2, 2, 2, 0, 0, 0, 1, 1, 1]),
]
 
print("ARI Examples:
")
for name, true, pred in examples:
    true = np.array(true)
    pred = np.array(pred)
    
    ari_manual, details = adjusted_rand_index_manual(true, pred)
    ari_sklearn = adjusted_rand_score(true, pred)
    
    print(f"{name}:")
    print(f"  True:    {true}")
    print(f"  Pred:    {pred}")
    print(f"  ARI:     {ari_sklearn:.4f}")
    print()
 
# Key property: ARI is invariant to label permutation!
print("Note: 'Perfect match' and 'Permuted labels' both get ARI = 1.0")
print("This label invariance is crucial - we care about structure, not label names.")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

import numpy as np
from sklearn.metrics import mutual_info_score, normalized_mutual_info_score
from sklearn.metrics import adjusted_mutual_info_score
 
def compute_entropy(labels):
    """Compute entropy of a label distribution."""
    n = len(labels)
    unique, counts = np.unique(labels, return_counts=True)
    probs = counts / n
    # Filter out zero probabilities to avoid log(0)
    probs = probs[probs > 0]
    return -np.sum(probs * np.log(probs))
 
 
def compute_mi_manual(labels_true, labels_pred):
    """Manual Mutual Information computation."""
    n = len(labels_true)
    
    # Joint distribution
    classes = np.unique(labels_true)
    clusters = np.unique(labels_pred)
    
    mi = 0
    for c in classes:
        for k in clusters:
            # Joint probability
            p_ij = np.sum((labels_true == c) & (labels_pred == k)) / n
            if p_ij == 0:
                continue
            
            # Marginal probabilities
            p_i = np.sum(labels_true == c) / n
            p_j = np.sum(labels_pred == k) / n
            
            mi += p_ij * np.log(p_ij / (p_i * p_j))
    
    return mi
 
 
# Example
labels_true = np.array([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2])
labels_pred = np.array([0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 2, 2])
 
h_true = compute_entropy(labels_true)
h_pred = compute_entropy(labels_pred)
mi_manual = compute_mi_manual(labels_true, labels_pred)
mi_sklearn = mutual_info_score(labels_true, labels_pred)
 
print("Information-Theoretic Analysis:")
print(f"  H(true):  {h_true:.4f} nats")
print(f"  H(pred):  {h_pred:.4f} nats")
print(f"  MI:       {mi_sklearn:.4f} nats")
print(f"  MI (manual): {mi_manual:.4f} nats")
print()
print(f"  MI / H(true): {mi_sklearn / h_true:.4f}")
print(f"  MI / H(pred): {mi_sklearn / h_pred:.4f}")
print()
 
# Show why normalization is needed
print("Effect of cluster count on MI:")
for n_clusters in [3, 5, 10, 20]:
    random_labels = np.random.randint(0, n_clusters, len(labels_true))
    mi = mutual_info_score(labels_true, random_labels)
    nmi = normalized_mutual_info_score(labels_true, random_labels)
    print(f"  k={n_clusters}: MI={mi:.4f}, NMI={nmi:.4f}")