From 50b9d0cd5906555c74ee0c1af6c4607c4f2eed5f Mon Sep 17 00:00:00 2001 From: Remis Lape Date: Tue, 8 Sep 2026 16:11:34 +0100 Subject: [PATCH] Add ekdist.mixtures: two populations in one set of estimates Fit a mechanism to many simulated records and the estimates do not always form one cloud. Some records take the optimiser to a second maximum, typically one where most events are missed, and the results come back as two groups an order of magnitude apart. The mean of a set like that describes neither group. mixture2 fits two Gaussian components in log space, for groups that overlap or nearly do, and reports the proportion in each and their geometric means. gap_split finds the largest multiplicative step between neighbouring sorted values, for a few runaways cleanly separated from the body. The point of both is that the split comes from the data. Choosing a threshold by eye and counting either side makes the answer a property of the threshold; these let an exclusion be stated -- "16 values above a 55-fold gap" -- and checked by someone else. Log space because rate constants are positive with multiplicative errors, so a group of them is log-normal rather than normal; that also makes Mixture.means geometric means. Written for the reproduction of Colquhoun, Hatton & Hawkes (2003), whose Fig. 2 is this situation exactly, but nothing in the module is specific to that paper or to ion channels. numpy only. 23 tests, against data generated from a known mixture rather than against values the function produced itself, because a mixture fit cannot look wrong: it returns two plausible components whatever it is handed, and raises nothing. Proportions are recovered at 10, 27 and 50 per cent so the tolerances are not tuned to one case. Two tests keep the suite honest rather than checking arithmetic: one records that a unimodal set is split without complaint, since that is the mistake a user will make and the separation of the means is the only thing that tells them; the other shows the recovery tolerances would reject the most likely way for the fit to be broken while still looking reasonable. Suite goes from 75 to 98. The clean-install job needs no change -- it enumerates modules since #8, and picks this one up as the eighth. Co-Authored-By: Claude Opus 5 --- ekdist/mixtures.py | 251 +++++++++++++++++++++++++++++++++++++++++ tests/test_mixtures.py | 222 ++++++++++++++++++++++++++++++++++++ 2 files changed, 473 insertions(+) create mode 100644 ekdist/mixtures.py create mode 100644 tests/test_mixtures.py diff --git a/ekdist/mixtures.py b/ekdist/mixtures.py new file mode 100644 index 0000000..5404ec3 --- /dev/null +++ b/ekdist/mixtures.py @@ -0,0 +1,251 @@ +"""Two populations in one set of estimates, and how to separate them. + +Fit a mechanism to many simulated records and the estimates do not always form +one cloud. Some records take the optimiser to a second maximum -- typically one +where most events are missed -- and the results come back as two groups an +order of magnitude apart. Any summary that ignores this is meaningless: the +mean of a set like that describes neither group. + +Two ways of finding the split, for two shapes of problem: + +:func:`mixture2` + fits two Gaussian components in log space, for two groups that overlap or + nearly do. Reports the proportion in each and their geometric means. +:func:`gap_split` + finds the largest multiplicative step between neighbouring sorted values, + for a handful of runaways cleanly separated from the body. + +**Why not just pick a threshold.** Choosing a cut by eye and counting either +side makes the answer a property of the cut. Both functions here determine the +split from the data, so an exclusion can be *stated* -- "16 values above a +55-fold gap" -- rather than assumed, and someone else can check it. + +**Why log space.** Rate constants are positive and their errors are +multiplicative, so a group of them is roughly log-normal, not normal. Fitting +in logs keeps every fitted mean positive and makes each group closer to the +Gaussian the mixture assumes. It also means :attr:`Mixture.means` are geometric +means. + +These were written for the reproduction of Colquhoun, Hatton & Hawkes (2003) +*J Physiol* **547**:699-728, whose Fig. 2 is exactly this situation: 1000 fits +of the same mechanism falling into a group near the true rates and a second +group where alpha2 comes out near 17 000 instead of 2000. Their examples appear +below because a concrete case is easier to read than an abstract one, but +nothing here is specific to that paper or to ion channels. + +Usage:: + + from ekdist import mixtures + + m = mixtures.mixture2(np.column_stack([alpha2, beta2])) + print(m.means, m.proportion(1)) # per cent, with its binomial SE + good = m.cluster(alpha2, 0) + + keep, ratio, cut = mixtures.gap_split(beta1b) + if cut is not None: + print(f"excluded {(~keep).sum()} values above a {ratio:.0f}-fold gap") +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np + + +def proportion(k, n): + """A proportion and its binomial standard error, **both per cent**.""" + p = k / n + return 100.0 * p, 100.0 * float(np.sqrt(p * (1.0 - p) / n)) + + +@dataclass(frozen=True) +class Mixture: + """A two-component Gaussian mixture, fitted in log space by :func:`mixture2`. + + Attributes + ---------- + weights : ndarray, shape (2,) + Mixing proportions, summing to 1, ordered so that component 0 is the + one with the smaller first coordinate. + means : ndarray, shape (2, d) + Component means in the **original** units, obtained by exponentiating + the log-space means, so they are geometric means. + labels : ndarray of int, shape (n,) + Each point assigned to its most probable component. + responsibilities : ndarray, shape (n, 2) + The posterior probability of each component for each point. Use these + rather than ``labels`` when the two groups overlap enough that a hard + assignment throws away what the fit knows. + loglik : float + Mean log-likelihood per point at the last EM step. + n_iter : int + converged : bool + Whether the change in ``loglik`` fell below the tolerance before the + iteration cap. A fit that did not converge is not necessarily wrong, + but it should not be quoted without looking at it. + """ + + weights: np.ndarray + means: np.ndarray + labels: np.ndarray + responsibilities: np.ndarray + loglik: float + n_iter: int + converged: bool + + def cluster(self, x, which): + """The members of ``x`` assigned to component ``which``. + + ``x`` is indexed by the same order as the data that was fitted, so it + can be any per-point quantity, not only a fitted coordinate. + """ + return np.asarray(x)[self.labels == which] + + def proportion(self, which=1): + """Proportion in one component, per cent, with its binomial SE. + + Taken from the hard assignment rather than from the fitted weight, so + that it is a count of points and its error is the error of a count. + The two agree well within that error when the groups are separated; + when they disagree, the groups overlap and neither number means much on + its own. + """ + return proportion(int((self.labels == which).sum()), self.labels.size) + + +def mixture2(X, n_iter=500, tol=1e-8): + """Fit two Gaussian components to positive data, in log space. + + Parameters + ---------- + X : array_like, shape (n,) or (n, d) + Strictly positive. Several coordinates at once are better than one: + two groups that overlap in every single parameter can still separate + cleanly in two, and the correlations within each group help. + n_iter, tol : int, float + EM iteration cap, and the change in mean log-likelihood per point that + counts as converged. + + Returns + ------- + Mixture + + Notes + ----- + Started from a split at the median of the first coordinate, which for a + genuinely bimodal set lies inside the gap, so no random restarts are + needed and the result is deterministic. + + **A set that is not bimodal will still return two components.** Nothing + here tests whether two is the right number. The separation of the means and + the proportion are what say whether the answer means anything: two + components fitted to one population typically land within a factor of two + of each other, where a real split is an order of magnitude. + """ + X = np.asarray(X, dtype=float) + if X.ndim == 1: + X = X[:, None] + if X.shape[0] < 4: + raise ValueError("need at least four points to fit two components") + if not np.all(X > 0): + raise ValueError("mixture2 fits in log space; all values must be > 0") + + Y = np.log(X) + n, d = Y.shape + # a split at the median of the first coordinate, which lies in the gap + hard = (Y[:, 0] > np.median(Y[:, 0])).astype(int) + resp = np.zeros((n, 2)) + resp[np.arange(n), hard] = 1.0 + + ridge = 1e-9 * np.eye(d) # keeps a near-degenerate component safe + prev = -np.inf + ll = -np.inf + converged = False + it = 0 + w = mu = None + for it in range(1, int(n_iter) + 1): + # M step + nk = resp.sum(axis=0) + 1e-300 + w = nk / n + mu = (resp.T @ Y) / nk[:, None] + cov = np.empty((2, d, d)) + for k in range(2): + D = Y - mu[k] + cov[k] = (D * resp[:, k, None]).T @ D / nk[k] + ridge + + # E step + logp = np.empty((n, 2)) + for k in range(2): + sign, logdet = np.linalg.slogdet(cov[k]) + if sign <= 0: + raise np.linalg.LinAlgError("degenerate mixture component") + D = Y - mu[k] + m = np.einsum("ij,jk,ik->i", D, np.linalg.inv(cov[k]), D) + logp[:, k] = np.log(w[k]) - 0.5 * (d * np.log(2 * np.pi) + + logdet + m) + top = logp.max(axis=1, keepdims=True) + lse = top[:, 0] + np.log(np.exp(logp - top).sum(axis=1)) + resp = np.exp(logp - lse[:, None]) + + ll = float(lse.mean()) + if abs(ll - prev) < tol: + converged = True + break + prev = ll + + order = np.argsort(mu[:, 0]) # component 0 = the smaller first coord + w, mu, resp = w[order], mu[order], resp[:, order] + return Mixture(weights=w, means=np.exp(mu), + labels=resp.argmax(axis=1), responsibilities=resp, + loglik=ll, n_iter=it, converged=converged) + + +def gap_split(x, min_ratio=5.0): + """Find a high outlier population separated by a multiplicative gap. + + For the case a mixture is the wrong tool for: a few values far above the + rest with nothing in between. In one of the CHH 2003 scenarios, 234 + estimates of beta1b lie below 671 and 16 above 37 180, with nothing + between -- a 55-fold step. Excluding them at a round number would be an + arbitrary cut; this finds the **largest ratio between neighbouring sorted + values** and, if it exceeds ``min_ratio``, splits there. + + Parameters + ---------- + x : array_like + Strictly positive. + min_ratio : float + The smallest gap that counts as a separation. 5 is well above the + largest step inside a unimodal set of a few hundred values, and well + below the 55-fold step above. Raise it if the body of the distribution + is itself broad. + + Returns + ------- + body : ndarray of bool + True for the members below the gap. **All True when no gap is found**, + so that the caller can apply the mask unconditionally. + ratio : float + The largest neighbouring ratio, whether or not it was taken -- so that + "no gap" can be reported with a number rather than as a bare negative. + cut : float or None + The last value below the gap, or None when no gap was found. + + Notes + ----- + Only the largest gap is considered, and only a high tail is separated. A + set with two gaps, or with outliers at both ends, needs looking at rather + than a function. + """ + x = np.asarray(x, dtype=float).ravel() + if not np.all(x > 0): + raise ValueError("gap_split needs positive values") + s = np.sort(x) + ratios = s[1:] / s[:-1] + i = int(np.argmax(ratios)) + ratio = float(ratios[i]) + if ratio < min_ratio: + return np.ones(x.shape, dtype=bool), ratio, None + cut = float(s[i]) + return x <= cut, ratio, cut diff --git a/tests/test_mixtures.py b/tests/test_mixtures.py new file mode 100644 index 0000000..2694199 --- /dev/null +++ b/tests/test_mixtures.py @@ -0,0 +1,222 @@ +"""Tests for ekdist.mixtures. + +A mixture fit is one of those functions that cannot look wrong: it returns two +components whatever it is handed, with plausible means and a plausible +proportion, and there is no exception to notice. So it is tested against data +**generated from a known mixture**, where the proportion and both means are +known before the fit runs, rather than against values it produced itself. + +Two tests exist to keep the suite honest rather than to check the arithmetic: +:func:`test_a_unimodal_set_is_split_without_a_gap` records that the function +will happily split one population, because that is the mistake a user will +make; and :func:`test_the_recovery_tests_could_fail` shows the tolerances are +tight enough to catch a wrong answer. +""" + +import numpy as np +import pytest + +from ekdist import mixtures + + +# --------------------------------------------------------------- fixtures + +def two_groups(n=1000, frac=0.27, seed=7): + """Two log-normal groups in two dimensions, an order of magnitude apart. + + Shaped like Fig. 2 of Colquhoun, Hatton & Hawkes (2003): a majority near + (2000, 52000) and a minority near (17000, 300000). The majority is + correlated within itself, as a real one is. + """ + rng = np.random.default_rng(seed) + k = int(round(frac * n)) + major = np.exp(rng.multivariate_normal( + [np.log(2000.0), np.log(52000.0)], + [[0.08 ** 2, 0.9 * 0.08 * 0.07], + [0.9 * 0.08 * 0.07, 0.07 ** 2]], n - k)) + minor = np.exp(rng.multivariate_normal( + [np.log(17000.0), np.log(3.0e5)], + [[0.10 ** 2, 0.0], [0.0, 0.09 ** 2]], k)) + return np.vstack([major, minor]) + + +# --------------------------------------------------------------- mixture2 + +def test_recovers_a_known_proportion(): + """27 % was put in; 27 % must come out, within the binomial error.""" + m = mixtures.mixture2(two_groups(n=1000, frac=0.27)) + assert m.converged + p, se = m.proportion(1) + assert abs(p - 27.0) < 3 * se, f"{p:.1f} % +- {se:.1f} against 27 %" + + +@pytest.mark.parametrize("frac", [0.10, 0.27, 0.50]) +def test_recovers_the_proportion_at_several_splits(frac): + """Not tuned to one case: an even split and a small minority both work.""" + m = mixtures.mixture2(two_groups(n=1500, frac=frac, seed=3)) + p, se = m.proportion(1) + assert abs(p - 100 * frac) < 3 * se + + +def test_recovers_both_group_means(): + m = mixtures.mixture2(two_groups(n=2000, frac=0.27)) + assert m.means[0][0] == pytest.approx(2000.0, rel=0.05) + assert m.means[1][0] == pytest.approx(17000.0, rel=0.05) + assert m.means[0][1] == pytest.approx(52000.0, rel=0.05) + assert m.means[1][1] == pytest.approx(3.0e5, rel=0.05) + + +def test_component_zero_is_always_the_lower_one(): + """The ordering is a promise, so that means[0] can be indexed blind.""" + for seed in (11, 17, 23): + m = mixtures.mixture2(two_groups(seed=seed)) + assert m.means[0][0] < m.means[1][0] + assert m.weights.sum() == pytest.approx(1.0) + + +def test_assignment_and_weight_agree_when_the_groups_are_separated(): + m = mixtures.mixture2(two_groups(n=1000, frac=0.27, seed=13)) + p, se = m.proportion(1) + assert abs(p - 100 * m.weights[1]) < 3 * se + + +def test_works_on_one_dimension(): + m = mixtures.mixture2(two_groups(n=800)[:, :1]) + p, _ = mixtures.proportion(int((m.labels == 1).sum()), m.labels.size) + assert 22.0 < p < 32.0 + + +def test_responsibilities_are_a_probability_distribution(): + m = mixtures.mixture2(two_groups(n=400, seed=29)) + assert m.responsibilities.shape == (400, 2) + np.testing.assert_allclose(m.responsibilities.sum(axis=1), 1.0) + assert np.all(m.responsibilities >= 0.0) + + +def test_cluster_selects_by_label_and_takes_any_per_point_quantity(): + X = two_groups(n=500, frac=0.20, seed=31) + m = mixtures.mixture2(X) + tag = np.arange(X.shape[0]) + assert m.cluster(tag, 0).size + m.cluster(tag, 1).size == X.shape[0] + np.testing.assert_array_equal(m.cluster(X[:, 0], 1), + X[m.labels == 1, 0]) + + +def test_is_deterministic(): + """No random restarts, so two runs on the same data agree exactly.""" + X = two_groups(n=600, seed=37) + a, b = mixtures.mixture2(X), mixtures.mixture2(X) + np.testing.assert_array_equal(a.labels, b.labels) + np.testing.assert_allclose(a.means, b.means, rtol=0, atol=0) + + +def test_rejects_non_positive_values(): + with pytest.raises(ValueError): + mixtures.mixture2(np.array([[1.0], [2.0], [-1.0], [3.0]])) + with pytest.raises(ValueError): + mixtures.mixture2(np.array([[1.0], [2.0], [0.0], [3.0]])) + + +def test_needs_four_points(): + with pytest.raises(ValueError): + mixtures.mixture2([1.0, 2.0, 3.0]) + + +def test_a_unimodal_set_is_split_without_a_gap(): + """The mistake a user will make, recorded rather than prevented. + + Nothing here tests whether two components is the right number. One + population comes back as two whose means are within a factor of two of each + other -- nothing like the eightfold separation of a real split -- and it is + that separation, not the absence of an error, that tells the caller. + """ + rng = np.random.default_rng(5) + X = np.exp(rng.normal(np.log(2000.0), 0.08, size=(500, 1))) + m = mixtures.mixture2(X) + assert m.means[1][0] / m.means[0][0] < 2.0 + + +def test_the_recovery_tests_could_fail(): + """The tolerances above are tight enough to catch a wrong answer. + + If mixture2 returned the overall mean for both components -- the most + likely way for it to be broken while still looking reasonable -- the + checks on the group means would reject it. + """ + X = two_groups(n=2000, frac=0.27) + overall = float(np.exp(np.log(X[:, 0]).mean())) + assert not (abs(overall - 2000.0) / 2000.0 < 0.05) + assert not (abs(overall - 17000.0) / 17000.0 < 0.05) + + +# -------------------------------------------------------------- gap_split + +def test_finds_a_separated_high_tail(): + """234 values below 670, 16 above 37 000: a 55-fold step.""" + body = np.linspace(320.0, 670.0, 234) + runaway = np.linspace(37000.0, 52000.0, 16) + x = np.concatenate([body, runaway]) + + keep, ratio, cut = mixtures.gap_split(x) + assert keep.sum() == 234 + assert ratio > 50.0 + assert cut == pytest.approx(670.0) + assert x[keep].max() < x[~keep].min() + + +def test_leaves_a_clean_set_alone(): + rng = np.random.default_rng(9) + x = np.exp(rng.normal(np.log(500.0), 0.2, size=300)) + keep, ratio, cut = mixtures.gap_split(x) + assert keep.all() + assert cut is None + assert ratio < 5.0 + + +def test_the_mask_is_usable_unconditionally(): + """All True when nothing is found, so a caller needs no special case.""" + rng = np.random.default_rng(41) + x = np.exp(rng.normal(np.log(500.0), 0.2, size=100)) + keep, _, _ = mixtures.gap_split(x) + assert x[keep].size == x.size + + +def test_the_threshold_is_honoured(): + x = np.array([1.0, 2.0, 3.0, 12.0]) # a 4-fold gap + assert mixtures.gap_split(x, min_ratio=5.0)[0].all() + assert mixtures.gap_split(x, min_ratio=3.0)[0].sum() == 3 + + +def test_the_ratio_is_reported_even_when_no_split_is_taken(): + """So that "no gap" can be stated with a number.""" + x = np.array([1.0, 2.0, 3.0, 12.0]) + keep, ratio, cut = mixtures.gap_split(x, min_ratio=5.0) + assert keep.all() and cut is None + assert ratio == pytest.approx(4.0) + + +def test_the_mask_follows_the_input_order_not_the_sorted_order(): + x = np.array([500.0, 40000.0, 300.0, 45000.0, 400.0]) + keep, _, cut = mixtures.gap_split(x) + np.testing.assert_array_equal(keep, [True, False, True, False, True]) + assert cut == pytest.approx(500.0) + + +def test_rejects_non_positive_values(): + with pytest.raises(ValueError): + mixtures.gap_split([1.0, 0.0, 2.0]) + with pytest.raises(ValueError): + mixtures.gap_split([1.0, -2.0, 3.0]) + + +# ------------------------------------------------------------- proportion + +def test_proportion_and_its_binomial_error(): + p, se = mixtures.proportion(27, 100) + assert p == pytest.approx(27.0) + assert se == pytest.approx(100 * np.sqrt(0.27 * 0.73 / 100)) + + +def test_proportion_is_exact_at_the_ends(): + assert mixtures.proportion(0, 50) == (0.0, 0.0) + assert mixtures.proportion(50, 50) == (100.0, 0.0)