Nonparametric Latent Class Analysis

Model

The nonparametric project of the Latent Class Analysis arc, and the direct answer to the one before it. Choosing the number of classes spent an entire notebook selecting $C$ with six criteria and a bootstrap test. This one removes the choice: a Dirichlet-process latent class model places a prior over a countably infinite set of classes, of which any finite sample occupies only finitely many — so the number of occupied classes $K$ becomes an ordinary posterior quantity, inferred rather than selected.

\text{partition}\sim\text{CRP}(\alpha),\qquad \delta_{kj}\sim\text{Beta}(a,b),\qquad x_{ij}\mid\text{class }k\sim\text{Bernoulli}(\delta_{kj})

From a finite mixture to an infinite one

The construction is a limit. LCA is a finite mixture with a $\text{Dirichlet}(\alpha/C,\dots,\alpha/C)$ prior on the class weights; let $C\to\infty$ and that prior converges to the stick-breaking weights of a Dirichlet process, with class profiles drawn from a base measure. The Chinese Restaurant Process then governs how subjects cluster — a subject joins an occupied class with probability proportional to its size, or opens a new one with probability proportional to the concentration $\alpha$ . The notebook opens by simulating the CRP directly to expose its two signatures: the logarithmic growth of $K$ , matched against the analytic $\mathbb{E}[K]=\alpha(\psi(\alpha+n)-\psi(\alpha))$ , and the rich-get-richer law where a few tables hold most customers and a long tail hold one or two.

The collapsed sampler

Because the base measure is conjugate, the class profiles are integrated out and the partition is sampled directly — Neal's Algorithm 3, using the Beta-Bernoulli posterior predictive. Class labels being arbitrary, every summary is label-invariant: the posterior of $K$ , and the co-clustering matrix giving $\Pr(i,j\ \text{share a class})$ . Four independent implementations agree on the carcinoma result — the from-scratch Python collapsed Gibbs, PyMC's truncated stick-breaking, a from-scratch base-R sampler, and R's dirichletprocess package with a custom Beta-Bernoulli component.

\Pr(\text{join }k)\propto n_k\prod_j\text{Pred}\big(x_{ij}\mid a+s_{kj},\,b+n_k-s_{kj}\big),\qquad \Pr(\text{new})\propto\alpha\prod_j\text{Pred}\big(x_{ij}\mid a,b\big)

Results

On the carcinoma ratings the posterior puts $P(K=3)\approx0.84$ — three latent slide types, reached in a single fit with no information criteria at all, and exactly the answer the previous project assembled from six criteria plus a bootstrap test. The concentration is learned rather than fixed: the data pull $\alpha$ down to $\approx0.48$ , an operating point whose prior expectation is about three classes. Parsimony emerges from the fit instead of being imposed on it.

Where it breaks — the Miller–Harrison inconsistency

The project is unusually honest about where this breaks, and that is its most valuable section. On simulated data with three true classes, the DP recovers the partition well — adjusted Rand index $\approx0.83$ — but its $K$ -posterior sits above the truth, with the mode at four. That is not sampling noise: holding the truth at three classes and growing $N$ through 200, 800 and 3200, the posterior mass on $K>3$ increases. This is the Miller–Harrison inconsistency — the Dirichlet process is not a consistent estimator of the number of components, however much data you give it.

Question	Dirichlet process	Finite $C$ with selection
How is $K$ obtained?	a posterior quantity, from one fit	chosen, by criteria across many fits
Consistent for the number of types?	no — over-extracts as $N$ grows	yes, under the usual conditions
Best used for	flexible clustering, density estimation	a small, interpretable count of types

Which to use

So the two views are complementary rather than competing. The DP is the right tool for flexible clustering, where a growing, data-adaptive number of components is a feature. When $K$ is meant to be a small, interpretable count of types, the finite- $C$ selection of the previous project — or a mixture-of-finite-mixtures prior, which restores consistency while keeping the conjugate machinery — is the principled route. The nonparametric model proposes the structure; the finite, selected model commits to a count.

Notebooks

Python DP-LCA from scratch — the CRP and its two signatures, a collapsed Gibbs sampler (Neal Algorithm 3), co-clustering matrices, the concentration α as a granularity dial, the Miller–Harrison inconsistency, and PyMC truncated stick-breaking (Python) Open → R Independent R re-implementation — a base-R collapsed CRP Gibbs sampler, cross-checked against the dirichletprocess package with a custom Beta-Bernoulli component (R) Open →

Downloads

dplca.py Collapsed CRP Gibbs sampler with concentration updating, co-clustering, representative partitions, cluster profiles and the adjusted Rand index (NumPy only) carcinoma.csv 118 slides rated by 7 pathologists — the same data the model-selection project analysed

DP-LCA Module — Source Code

"""
dplca.py -- From-scratch NONPARAMETRIC (Dirichlet-process) LATENT CLASS ANALYSIS.

Backs the notebooks in  "Nonparametric Latent Class Analysis"  (the nonparametric project of
the Latent Class Analysis arc).

A DP-LCA replaces the finite Dirichlet prior on C class weights with a DIRICHLET
PROCESS, giving a countably infinite number of latent classes of which only finitely
many are occupied by any finite sample. The number of occupied classes K is then a
posterior quantity -- inferred, not chosen. Model:

    (partition of subjects)   ~ CRP(alpha)                         [Chinese Restaurant Process]
    delta_{k,j}               ~ Beta(a, b)          i.i.d.         [base measure G0]
    x_ij | class k            ~ Bernoulli(delta_{k,j})            [local independence]

We sample it with the COLLAPSED CRP GIBBS sampler (Neal 2000, Algorithm 3): the
class-specific delta are integrated out, so each subject is reassigned using the
Beta-Bernoulli posterior predictive,

    P(join class k) proportional to  n_k  * prod_j  Pred(x_ij ; a+s_kj, b+n_k-s_kj)
    P(open a new class) proportional to alpha * prod_j Pred(x_ij ; a, b),

with Pred(x=1 ; A, B) = A/(A+B). The concentration alpha can be given a Gamma prior
and updated by the Escobar-West (1995) auxiliary-variable step.

Because the class labels are arbitrary, results are summarised by LABEL-INVARIANT
quantities: the posterior of K, and the co-clustering (posterior similarity) matrix.
"""

import numpy as np


# --------------------------------------------------------------------------- #
#  Collapsed CRP Gibbs sampler for the DP mixture of product-Bernoulli          #
# --------------------------------------------------------------------------- #

def dp_lca_gibbs(X, rng, alpha=1.0, a=1.0, b=1.0, draws=3000, burn=1000,
                 sample_alpha=False, alpha_prior=(2.0, 4.0), init_k=1):
    """Collapsed Gibbs (Neal Alg. 3) for DP-LCA. Returns the trace of the number of
    occupied classes K, the concentration alpha, and the (canonically relabelled)
    subject assignments at each retained sweep."""
    X = np.asarray(X, float); N, J = X.shape
    # initialise: init_k clusters at random
    z = rng.integers(0, init_k, N)
    clusters = {}
    for k in range(init_k):
        m = z == k
        if m.any():
            clusters[k] = dict(n=int(m.sum()), s=X[m].sum(0))
    # relabel to compact ids
    z, clusters, next_id = _compact(z, clusters)
    K_tr, al_tr, Z_tr = [], [], []
    for it in range(draws + burn):
        for i in range(N):
            xi = X[i]; ci = z[i]
            clusters[ci]["n"] -= 1; clusters[ci]["s"] = clusters[ci]["s"] - xi
            if clusters[ci]["n"] == 0:
                del clusters[ci]
            keys = list(clusters.keys())
            logp = np.empty(len(keys) + 1)
            for idx, k in enumerate(keys):
                nk = clusters[k]["n"]; sk = clusters[k]["s"]
                p1 = (a + sk) / (a + b + nk)
                logp[idx] = np.log(nk) + np.sum(xi * np.log(p1) + (1 - xi) * np.log1p(-p1))
            p1n = a / (a + b)
            logp[-1] = np.log(alpha) + np.sum(xi * np.log(p1n) + (1 - xi) * np.log(1 - p1n))
            logp -= logp.max(); pr = np.exp(logp); pr /= pr.sum()
            ch = rng.choice(len(keys) + 1, p=pr)
            if ch == len(keys):
                clusters[next_id] = dict(n=1, s=xi.copy()); z[i] = next_id; next_id += 1
            else:
                k = keys[ch]; clusters[k]["n"] += 1; clusters[k]["s"] = clusters[k]["s"] + xi; z[i] = k
        K = len(clusters)
        if sample_alpha:
            alpha = _update_alpha(alpha, K, N, rng, *alpha_prior)
        if it >= burn:
            K_tr.append(K); al_tr.append(alpha); Z_tr.append(_canonical(z))
    return dict(K=np.array(K_tr), alpha=np.array(al_tr), Z=np.array(Z_tr),
                a=a, b=b, X=np.asarray(X, float))


def _compact(z, clusters):
    old = sorted(clusters.keys()); remap = {o: i for i, o in enumerate(old)}
    z = np.array([remap[v] for v in z]); cl = {remap[o]: clusters[o] for o in old}
    return z, cl, len(old)

def _canonical(z):
    """Relabel a partition by order of first appearance (a label-invariant canonical form)."""
    seen = {}; out = np.empty(len(z), int); nxt = 0
    for i, v in enumerate(z):
        if v not in seen: seen[v] = nxt; nxt += 1
        out[i] = seen[v]
    return out


def _update_alpha(alpha, K, N, rng, a0, b0):
    """Escobar-West (1995) update of the DP concentration under a Gamma(a0, b0) prior."""
    eta = rng.beta(alpha + 1, N)
    c = a0 + K - 1; d = N * (b0 - np.log(eta))
    pi = c / (c + d)
    shape = a0 + K if rng.random() < pi else a0 + K - 1
    return rng.gamma(shape, 1.0 / (b0 - np.log(eta)))


# --------------------------------------------------------------------------- #
#  Label-invariant summaries                                                    #
# --------------------------------------------------------------------------- #

def coclustering(Z):
    """Posterior co-clustering (similarity) matrix: P(subjects i,j share a class)."""
    Z = np.asarray(Z); S, N = Z.shape
    P = np.zeros((N, N))
    for s in range(S):
        same = (Z[s][:, None] == Z[s][None, :])
        P += same
    return P / S

def representative_partition(Z, target_K=None):
    """A representative clustering: the retained draw whose partition is closest (in
    co-clustering Hamming distance) to the posterior mean, optionally restricted to
    draws with K = target_K (e.g. the modal K)."""
    Z = np.asarray(Z); P = coclustering(Z); S, N = Z.shape
    Ks = np.array([len(np.unique(z)) for z in Z])
    idx = np.arange(S) if target_K is None else np.where(Ks == target_K)[0]
    best, bd = None, np.inf
    for s in idx:
        same = (Z[s][:, None] == Z[s][None, :]).astype(float)
        d = np.abs(same - P).sum()
        if d < bd: bd, best = d, s
    return _canonical(Z[best])

def cluster_profiles(X, z, a=1.0, b=1.0):
    """Posterior-mean item-response probabilities delta_k for each cluster in partition z,
    ordered by cluster size (descending)."""
    X = np.asarray(X, float); z = np.asarray(z)
    ks, counts = np.unique(z, return_counts=True)
    order = ks[np.argsort(counts)[::-1]]
    lam = np.array([np.mean(z == k) for k in order])
    delta = np.array([(a + X[z == k].sum(0)) / (a + b + (z == k).sum()) for k in order])
    return lam, delta


def adjusted_rand(labels_true, labels_pred):
    """Adjusted Rand Index between two clusterings (label-invariant agreement; 1 = identical)."""
    from math import comb
    a = np.asarray(labels_true); b = np.asarray(labels_pred)
    ca = {v: i for i, v in enumerate(np.unique(a))}; cb = {v: i for i, v in enumerate(np.unique(b))}
    A = np.array([ca[v] for v in a]); B = np.array([cb[v] for v in b])
    cont = np.zeros((len(ca), len(cb)), int)
    for x, y in zip(A, B): cont[x, y] += 1
    sum_c = sum(comb(v, 2) for v in cont.sum(1))
    sum_k = sum(comb(v, 2) for v in cont.sum(0))
    sum_ij = sum(comb(v, 2) for v in cont.flat)
    n2 = comb(len(a), 2); exp = sum_c * sum_k / n2; mx = 0.5 * (sum_c + sum_k)
    return (sum_ij - exp) / (mx - exp) if mx != exp else 1.0


def simulate_lca(N, lam, delta, rng):
    lam = np.asarray(lam, float); delta = np.asarray(delta, float); C, J = delta.shape
    T = rng.choice(C, size=N, p=lam); X = (rng.random((N, J)) < delta[T]).astype(int)
    return X, T

References

Neal, R. M. (2000). Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics 9(2), 249–265. — Algorithm 3, the collapsed conjugate sampler implemented here
Ferguson, T. S. (1973). A Bayesian analysis of some nonparametric problems. Annals of Statistics 1(2), 209–230. — the Dirichlet process as a prior over distributions
Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statistica Sinica 4(2), 639–650. — stick-breaking, the representation PyMC's truncated model uses
Miller, J. W. & Harrison, M. T. (2014). Inconsistency of Pitman–Yor process mixtures for the number of components. Journal of Machine Learning Research 15, 3333–3370. — the over-extraction demonstrated here by growing the sample size at fixed truth
Miller, J. W. & Harrison, M. T. (2018). Mixture models with a prior on the number of components. Journal of the American Statistical Association 113(521), 340–356. — the mixture-of-finite-mixtures alternative that restores consistency
Escobar, M. D. & West, M. (1995). Bayesian density estimation and inference using mixtures. Journal of the American Statistical Association 90(430), 577–588. — the Gamma-prior update for the concentration parameter used here
Hubert, L. & Arabie, P. (1985). Comparing partitions. Journal of Classification 2(1), 193–218. — the adjusted Rand index used to score the recovered partition