Skip to content

Online Gaussian mixture

Some streams are not one distribution. A joint that switches between loaded and unloaded, a sensor that alternates between two modes, a return series that mixes a calm and a volatile regime: a single mean and a single variance average the two behaviours into a value that describes neither. A mixture keeps one Gaussian per behaviour and learns which one a sample belongs to.

What this is

A Gaussian mixture represents a distribution as a weighted sum of k Gaussian components. Each component has a weight omega_j, a mean mu_j and a variance sigma_j^2. The model is trained online: one sample at a time, an EM step updates the weights and the components. transform_one returns, for the sample, the log-density and the responsibilities (the posterior probability of belonging to each component).

The number of components k is fixed at construction. Full covariance is not implemented: the covariance is diagonal, positive by construction.

1. The first estimator

from dense_armor.utility.signals.stats import OnlineGaussianMixture

mx = OnlineGaussianMixture(n_components=2, learning_rate=0.1)
for v in [-2.0, -1.9, -2.1, 2.0, 2.1, 1.9] * 200:
    mx.learn_one({"x": v})
print(sorted(round(float(m), 1) for m in mx.means_.ravel()))

The output is the two recovered means; the class docstring example gives them. Two components on a bimodal stream end up on the two modes.

2. What a responsibility is

For each sample, the model computes the probability that the sample was generated by each component. This is the responsibility r_j(x), a number in [0, 1], with the responsibilities over all components summing to one. The EM update uses the responsibilities as weights: a sample that is clearly on the left mode gets r_left = 1.0, r_right = 0.0, and it moves only the left mean; a sample between the modes gets r_left = r_right = 0.5, and it moves both means a little.

from dense_armor.utility.signals.stats import OnlineGaussianMixture

mx = OnlineGaussianMixture(
    n_components=2, learning_rate=0.1, init_means=[[-2.0], [2.0]]
)
for v in [-2.0, 2.0] * 50:
    mx.learn_one({"x": v})
print([round(r, 3) for r in mx.transform_one({"x": -2.0})["responsibilities"]])
print([round(r, 3) for r in mx.transform_one({"x": 0.0})["responsibilities"]])

A sample on a mode gets a responsibility near 1 for that mode; a sample exactly between the modes gets 0.5 and 0.5. The numbers are in the table at the bottom of the page.

3. The formula

The Gaussian mixture density is

f(x) = sum_j omega_j N(x; mu_j, sigma_j^2)

with N(x; mu, sigma^2) the Gaussian density with mean mu and variance sigma^2. The online update follows Eq. (15) of Cappé & Moulines (2009), page 4, with the responsibilities from Eq. (18), page 5. The paper works the Poisson mixture explicitly; the Gaussian form is derived in the module docstring by substituting the Gaussian density for the Poisson density.

The M-step uses the weighted sufficient statistics, not a plain (1 - gamma) * mu + gamma * (r/w) * x. The two brackets are the zero-th, first and second moment of each component; the exact formulas are in the module docstring.

4. A stream from a robot

A joint position in radians, a slow sinusoid with rare spikes and a step in the mean. The stream is unimodal: there is no second mode to separate.

import numpy as np
from dense_armor.utility.signals.stats import OnlineGaussianMixture

rng = np.random.default_rng(20261026)
n = 2000
t = np.linspace(0.0, 4.0 * np.pi, n)
q = 0.1 * np.sin(t)
spikes = rng.random(n) < 0.02
q[spikes] += rng.normal(0.0, 0.5, int(spikes.sum()))
q[n // 2:] += 0.2

mx = OnlineGaussianMixture(
    n_components=2, learning_rate=0.05, init_means=[[0.0], [0.2]]
)
for v in q:
    mx.learn_one({"x": float(v)})
print([round(float(w), 3) for w in mx.weights_])

The two weights end up at [0.0, 1.0]: one component takes the whole stream, the other one has weight zero. This is the correct answer on a unimodal domain. A mixture that returned [0.5, 0.5] on a stream with one mode would be inventing a second mode that the data does not have. The benchmark records the result as a declared behaviour, not a failure.

5. A stream from a drone

An IMU acceleration in m/s², a random walk with Student-t noise and a variance change at the midpoint.

import numpy as np
from dense_armor.utility.signals.stats import OnlineGaussianMixture

rng = np.random.default_rng(20261026)
n = 2000
noise = rng.standard_t(df=4.0, size=n) * 0.5
walk = np.cumsum(rng.normal(0.0, 0.05, n))
walk[n // 2:] *= 3.0
acc = walk + noise

mx = OnlineGaussianMixture(
    n_components=2, learning_rate=0.05, init_means=[[-1.0], [1.0]]
)
for v in acc:
    mx.learn_one({"x": float(v)})
print([round(float(m), 3) for m in mx.means_.ravel()])
print([round(float(w), 3) for w in mx.weights_])

On this stream the two components end up close to each other: a random walk is not a mixture of two modes. The benchmark records the collapse, and the same answer as the robot case: a mixture on a unimodal stream reduces to a single component, which is the correct behaviour.

6. A bimodal stream

Two balanced Gaussian modes at -2 and +2, noise 0.3, shuffled. This is the domain the mixture exists for.

import numpy as np
from dense_armor.utility.signals.stats import OnlineGaussianMixture

rng = np.random.default_rng(20261026)
n = 2000
half = n // 2
x = np.concatenate([
    rng.normal(-2.0, 0.3, half),
    rng.normal(+2.0, 0.3, half),
])
rng.shuffle(x)

mx = OnlineGaussianMixture(
    n_components=2, learning_rate=0.05, init_means=[[-2.0], [2.0]]
)
for v in x:
    mx.learn_one({"x": float(v)})
print(sorted(round(float(m), 2) for m in mx.means_.ravel()))
print([round(float(w), 3) for w in mx.weights_])

The two means are recovered to -2.05 and +1.96, the two weights to 0.548 and 0.452: the model separates the two modes and assigns a plausible proportion to each. This is the result the mixture is designed for, and the benchmark records it as the positive case.

Details

References

  • Cappé, O., Moulines, E. (2009). On-line expectation-maximization algorithm for latent data models. Journal of the Royal Statistical Society: Series B 71(3), 593-613. Eq. (15) page 4; Eq. (18) page 5; Eq. (19) page 6. The paper works the Poisson mixture explicitly; the Gaussian form is derived in the module docstring.
  • Neal, R. M., Hinton, G. E. (1998). A view of the EM algorithm that justifies incremental, sparse, and other variants. NIPS 1998. Not in the folder, bibliographic only.

Choices stated explicitly

  • Covariance is diagonal, not full. Full covariance would make the M-step O(k d^3) per sample. Declared in the class docstring.
  • The default initialisation places the k component means at small deterministic offsets along the first axis. For arbitrary data, pass init_means; the default can leave a component with weight near zero.
  • On a sequentially ordered stream the components collapse on the first mode. On an interleaved stream (the realistic case) they separate given init_means. Two tests document the two behaviours with the same seed.
  • n_components is fixed. Adding or removing components on the fly is not implemented.
  • merge is a weighted average of the two states by n_; the online update is not additive, so the merge of two partial states is not identical to the sequential pass over the union.

Measured cost

On the host of this session, median per-sample latency of learn_one is 106 µs, 90th percentile 140 µs, maximum 260 µs (measured by benchmarks/statistics_online.py --smoke --robust). The cost is O(k d) per sample. The state is weights_ (k floats), means_ (k d floats), variances_ (k d floats), n_, n_missing_, and the initialisation flag. Memory class: O(k d) with k and d fixed at construction.

Numeric results on this page

Every number quoted above was produced by running the code on this page, in the session of benchmarks/statistics_online.py --smoke --robust. The JSON file is benchmarks/statistics_online.json.

Quantity Value
Robot, mixture weights [0.0, 1.0]
Robot, means [-0.027, 0.188]
Drone, mixture weights [0.045, 0.955]
Drone, means [-1.103, -0.710]
Bimodal, mixture means [-2.048, 1.956]
Bimodal, weights [0.548, 0.452]
Bimodal, true modes -2.0, +2.0
Median per-sample latency of learn_one 106 µs

The unit tests on test/stats/test_mixture.py (24 tests) check the density against scipy.stats.norm.logpdf for one component within 1e-10, the positive definiteness of every variance, the weights summing to one, the ValueError on invalid constructor arguments, merge, state_dict roundtrip, fit vs learn_one parity, the two declared behaviours on shuffled and sequential streams, and check_estimator with the skipped list empty.