Reading a Confusion Matrix - The Most-Confused Facies Pair

Try solving it yourself first — the hints walk you through it. The solution below is one valid approach; yours may differ and still be correct.

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import confusion_matrix


# ── Verified OML-58 facies generator (do not edit) ───────────────────────
FACIES = {0: "Shale", 1: "Brine Sand", 2: "Gas Sand", 3: "Siltstone"}

# Each facies' centre in log space (GR, RHOB, NPHI, ln RT) and its scatter (the *_sd
# terms). Resistivity is centred on a LOG scale because it ranges over decades. Defined
# once and shared by the training-set and new-well generators below, which guarantees the
# rock physics is identical in training and inference.
CENTROIDS = {0: (108, 2.51, 0.33, np.log(3.5),  16, 0.080, 0.045, 0.40),
             1: (46,  2.28, 0.27, np.log(9.0),   18, 0.075, 0.045, 0.45),
             2: (43,  2.21, 0.21, np.log(15.0),  18, 0.085, 0.050, 0.55),
             3: (74,  2.40, 0.29, np.log(7.0),   21, 0.090, 0.050, 0.48)}

def make_facies_dataset(n=1500, seed=11):
    """Self-contained labelled facies set for OML 58: four logs -> rock type.

    Each facies has a centroid in log space (GR, RHOB, NPHI, log-RT) and realistic
    scatter. The gas sand deliberately overlaps the brine sandstone - same low GR
    and high porosity - separated only by its lower density and neutron and higher
    resistivity. That overlap, plus its rarity (7% of the section), is what makes it
    the hard class, exactly as in the field.
    """
    rng = np.random.default_rng(seed)
    facies = rng.choice([0, 1, 2, 3], size=n, p=[0.42, 0.33, 0.07, 0.18])  # imbalanced
    GR = np.zeros(n); RHOB = np.zeros(n); NPHI = np.zeros(n); RT = np.zeros(n)
    for f, (gr, rhob, nphi, lnrt, grs, rhs, nps, rts) in CENTROIDS.items():
        m = facies == f; k = int(m.sum())
        GR[m]   = gr + rng.normal(0, grs, k)
        RHOB[m] = rhob + rng.normal(0, rhs, k)
        NPHI[m] = np.clip(nphi + rng.normal(0, nps, k), 0, 0.55)
        RT[m]   = np.clip(np.exp(lnrt + rng.normal(0, rts, k)), 0.3, 600)
    return pd.DataFrame({"GR": np.round(GR, 1), "RHOB": np.round(RHOB, 3),
                         "NPHI": np.round(NPHI, 3), "RT": np.round(RT, 2),
                         "facies": facies})


# ── Verified split + the two models the chapter builds (do not edit) ─────
rock = make_facies_dataset()                       # n=1500, seed=11
X = rock[["GR", "RHOB", "NPHI", "RT"]].values
y = rock["facies"].values
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.25, random_state=0, stratify=y)

baseline = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)).fit(X_train, y_train)
forest   = RandomForestClassifier(n_estimators=200, random_state=0).fit(X_train, y_train)
# ─────────────────────────────────────────────────────────────────────────


def most_confused_pair(model, X_test, y_test):
    """Name the single most-often-confused (true, predicted) facies pair.

    Build the 4x4 confusion matrix, zero the diagonal (correct calls), find the
    largest off-diagonal cell, and return (FACIES[true_idx], FACIES[pred_idx]).
    """
    cm = confusion_matrix(y_test, model.predict(X_test))
    off = cm - np.diag(np.diag(cm))                 # drop the correct calls
    true_idx, pred_idx = np.unravel_index(np.argmax(off), off.shape)
    return (FACIES[true_idx], FACIES[pred_idx])


baseline_pair = most_confused_pair(baseline, X_test, y_test)
forest_pair   = most_confused_pair(forest,   X_test, y_test)
same_pair     = (baseline_pair == forest_pair)

print("baseline most-confused pair:", baseline_pair)
print("forest   most-confused pair:", forest_pair)
print("same pair?", same_pair)

lockCopying code is a Full Access feature.