Metadata-Version: 2.4
Name: nongaussian-mixtures
Version: 0.4.0
Summary: Mixture models with non-Gaussian components, scikit-learn compatible: mixtures of Dirichlet distributions for compositional data, of Beta distributions for data in the unit interval, and Gaussian mixtures fitted to binned data.
Project-URL: Homepage, https://github.com/mbaelde/nongaussian-mixtures
Project-URL: Issues, https://github.com/mbaelde/nongaussian-mixtures/issues
Author: Maxime Baelde
License-Expression: BSD-3-Clause
License-File: LICENSE
Keywords: beta-distribution,binned-data,compositional-data,density-estimation,dirichlet,mixture-model,scikit-learn
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Python: >=3.10
Requires-Dist: numpy>=1.26
Requires-Dist: scikit-learn>=1.6
Requires-Dist: scipy>=1.12
Provides-Extra: dev
Requires-Dist: mypy>=1.10; extra == 'dev'
Requires-Dist: pytest>=8; extra == 'dev'
Requires-Dist: ruff>=0.6; extra == 'dev'
Requires-Dist: scipy-stubs; extra == 'dev'
Description-Content-Type: text/markdown

# 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

```bash
pip install nongaussian-mixtures
```

## Usage

```python
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:

```python
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:

```python
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:

```python
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

```bash
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](https://github.com/mbaelde/generative-audio-source-models).

## License

BSD 3-Clause.
