nongaussian-mixtures

Mixture models whose components are not Gaussian, with a scikit-learn API.

Version 0.1 ships four estimators:

• DirichletMixture, a mixture of Dirichlet distributions, for compositional data (rows summing to one);
• BetaMixture, a mixture of products of Beta distributions, for data confined to the unit interval (one Beta per feature, features independent within a component). With a single feature it is the univariate Beta mixture;
• BayesianDirichletMixture, the same Dirichlet mixture fitted by variational inference, where n_components is an upper bound and the fit decides how many components the data supports;
• BinnedGaussianMixture, fitted to a histogram rather than to samples. Its components are Gaussian, which makes it the exception here; what is not Gaussian is the observation model, a multinomial over the cells of a grid.

Careful with the vocabulary: in scikit-learn, "Dirichlet" names a prior on the mixture weights (BayesianGaussianMixture with a Dirichlet-process prior over Gaussian components). Here the Dirichlet is the component density itself. There is no such estimator in scikit-learn.

Compositional data is anything whose samples are vectors of non-negative parts summing to one: normalised power spectra, topic proportions, relative abundances, word frequency profiles. A Gaussian mixture on such data ignores both the positivity and the sum-to-one constraint; a Dirichlet mixture is the natural model.

Install

pip install nongaussian-mixtures

Usage

import numpy as np
from scipy.stats import dirichlet
from nongaussian_mixtures import DirichletMixture

X = np.vstack(
    [
        dirichlet.rvs([10.0, 1.0, 1.0], size=500, random_state=0),
        dirichlet.rvs([1.0, 1.0, 10.0], size=500, random_state=1),
    ]
)

model = DirichletMixture(n_components=2, random_state=0).fit(X)

model.alphas_  # (2, 3) concentration parameters
model.weights_  # (2,) mixture proportions
model.score_samples(X)  # log-likelihood per sample
model.predict(X)  # most likely component
model.predict_proba(X)  # posterior over components

Rows are projected onto the simplex before fitting: zeros are floored to eps (the Dirichlet density is undefined on the boundary) and rows are renormalised, so unnormalised counts or energies can be passed directly. Negative values are rejected.

For data in the unit interval rather than on the simplex:

import numpy as np
from nongaussian_mixtures import BetaMixture

rng = np.random.default_rng(0)
X = np.vstack(
    [
        rng.beta([2.0, 8.0], [8.0, 2.0], size=(500, 2)),
        rng.beta([9.0, 2.0], [2.0, 9.0], size=(500, 2)),
    ]
)

model = BetaMixture(n_components=2, random_state=0).fit(X)

model.alphas_, model.betas_  # (2, 2) each: one Beta pair per component and feature

To let the fit choose the number of components instead of fixing it:

from nongaussian_mixtures import BayesianDirichletMixture

model = BayesianDirichletMixture(n_components=10, random_state=0).fit(X)

model.n_components_  # how many survived pruning
model.lower_bound_  # evidence lower bound per sample, for model comparison

When the samples are gone and only a histogram is left:

import numpy as np
from nongaussian_mixtures import BinnedGaussianMixture

samples = np.random.default_rng(0).normal(3.0, 2.0, size=(20_000, 1))
counts, edges = np.histogramdd(samples, bins=(np.arange(-8.0, 14.1, 1.0),))
centers = ((edges[0][:-1] + edges[0][1:]) / 2)[:, None]

model = BinnedGaussianMixture(bin_width=1.0).fit(centers, counts=counts.ravel())

model.means_, model.variances_  # variance 3.97, against 4.06 for the raw centres

X holds the cell centres, one row per cell, and counts how many samples fell in each. Feeding those centres to GaussianMixture instead inflates every variance by bin_width ** 2 / 12, the spread of the samples inside a cell, and distorts the split between components along the way.

Every estimator passes sklearn.utils.estimator_checks.check_estimator and works inside pipelines and GridSearchCV.

Algorithm

Expectation-maximisation. The E-step is computed in log-space; the M-step is a weighted maximum-likelihood Dirichlet fit per component, solving the digamma system for the Dirichlet

$$\psi(\alpha_k) - \psi!\Big(\sum_j \alpha_j\Big) = \overline{\log x_k}$$

by damped Newton iterations, following Minka's reference implementation (fastfit):

• the Hessian is diagonal plus rank-one, H = -diag(ψ'(α)) + ψ'(Σα) 11ᵀ, so H⁻¹g is obtained in O(D) by Sherman-Morrison instead of O(D³) by a dense solve. At D = 513 (an audio spectrum) that is what makes the fit usable;
• a Levenberg-Marquardt damping is applied to the diagonal, and a step is accepted only if it keeps every α_k > 0 and increases the weighted log-likelihood. Ronning (1989) showed the undamped Newton step can leave the admissible region;
• initialisation is by the method of moments, each dimension contributing an independent estimate of the precision.

Because every accepted M-step increases the weighted log-likelihood, the EM lower bound is monotone by construction (there is a test for it).

BetaMixture needs no separate solver: Beta(a, b) on x is Dirichlet(a, b) on (x, 1 - x), so each (alpha, beta) pair is fitted by the same damped Newton iteration on the sufficient statistic (mean log x, mean log(1 - x)).

BayesianDirichletMixture replaces maximum likelihood by mean-field variational inference, after Ma & Leijon (2014). Each concentration parameter gets a Gamma prior and the weights a symmetric Dirichlet one, so components the data does not support see their weight collapse and are pruned. The obstacle is that

$$\mathbb{E}_q\Big[\ln\Gamma\Big(\sum_d \alpha_d\Big) - \sum_d \ln\Gamma(\alpha_d)\Big]$$

has no closed form; it is expanded to second order in ln α around the posterior mean, which keeps the Gamma posteriors conjugate. That expansion is checked against a Monte-Carlo estimate in the test suite, and the resulting bound is verified to be monotone.

Being a local optimum, the bound can settle on redundant components. Fit a few random_state values and keep the largest lower_bound_, exactly as with BayesianGaussianMixture.

BinnedGaussianMixture maximises the multinomial likelihood of the histogram, Σ_j n_j ln P_j with P_j the probability the mixture assigns to cell j (McLachlan & Peel, chapter 9). EM treats the samples as the missing data, so the E-step needs the probability of each cell and the first two moments of each component truncated to it. Covariances are diagonal and cells are boxes, so everything factorises over features and those moments are closed-form, in terms of Φ and φ alone. The reference implementation of the thesis carries a full covariance in two dimensions and computes the same moments by numerical quadrature, one call per cell, component and matrix entry; diagonal covariances buy the closed form, and with it an arbitrary number of features.

The cell probabilities are computed in log space, which is not decoration: a cell thirty sigma out from a component underflows to exactly zero in float64, and the moment ratios are then 0 / 0. The reference implementation stops there. In log space the ratios stay finite and say the sensible thing, the conditional mean sitting on the near edge of the cell.

Development

uv run --extra dev ruff check .
uv run --extra dev ruff format --check .
uv run --extra dev mypy
uv run --extra dev pytest

CI runs the same four commands, the tests on Python 3.10 to 3.13. Type checking is pinned to 3.13: what mypy sees depends on the resolved scipy-stubs, which differs between the oldest and the newest supported dependency set.

Roadmap

• full covariances for BinnedGaussianMixture. A box then has no factorised probability, so both the cell probability and its moments go back to numerical integration;
• open-ended cells, for histograms whose extreme bins collect everything beyond the grid (censored tails).

References

• T. P. Minka, Estimating a Dirichlet distribution, 2000 (rev. 2012).
• Z. Ma, A. Leijon, Bayesian estimation of Dirichlet mixture model with variational inference, Pattern Recognition 47(9), 2014, 3143-3157.
• C. M. Bishop, Pattern Recognition and Machine Learning, 2006, chapter 10.
• G. Ronning, Maximum likelihood estimation of Dirichlet distributions, Journal of Statistical Computation and Simulation 32(4), 1989, 215-221.
• N. Wicker, J. Muller, R. K. R. Kalathur, O. Poch, A maximum likelihood approximation method for Dirichlet's parameter estimation, Computational Statistics & Data Analysis 52(3), 2008, 1315-1322.
• G. McLachlan, D. Peel, Finite Mixture Models, Wiley, 2000, chapter 9, and G. McLachlan, P. Jones, Fitting mixture models to grouped and truncated data via the EM algorithm, Biometrics 44(2), 1988, 571-578.
• M. Baelde, Modèles génératifs pour la classification et la séparation de sources sonores en temps-réel, PhD thesis, Université de Lille, 2019, appendix B.2, and mvbetapdf.m for the product-of-Betas component; appendix B.1 and gmm2d_binned.m for the binned mixture. The fitter here is the standalone version of the one used in generative-audio-source-models.

License

BSD 3-Clause.