Multi-dimensional gaussian fit

MultiMin: Multivariate Gaussian fitting

This package provides tools for fitting Mixture of Gaussians (MoG) and other statistical utilities.

Main Features

Multivariate fitting (MoG)
Visualization tools (Density plots)
Statistical utilities

Usage

>>> import multimin as mn

For more information, visit: https://github.com/seap-udea/multimin

class multimin.FitFunctionMoG(data, ngauss=1)[source]

Bases: MultiMinBase

Fit a MoG to a function evaluated on a mesh (instead of to data points).

The objective is to approximate the function by a composed multivariate normal distribution (MoG) by minimizing a loss between the function values and the MoG PDF on the same grid.

Initialize with data=(Xmesh, Fmesh) and ngauss: - Univariate: Xmesh and Fmesh are 1D arrays (same length). - Multivariate: Xmesh is a tuple of ndarrays (e.g. from np.meshgrid),

Fmesh is an array of the same shape as the grid.

X

Grid points, shape (N, nvars).

Type:: numpy.ndarray

F

Function values at the grid points, shape (N,).

Type:: numpy.ndarray

ngauss, nvars, mog, minparams, scales, uparams: Same as in FitMoG (inherited).

normalization

Scale factor fitted so that normalization * mog.pdf(X) approximates F. Set by fit_data(); None before fitting.

Type:: float

Theory

The fitting procedure is based on the Maximum Likelihood Estimation (MLE) method. Given a dataset of $S$ points $\tilde U_k$ (e.g., unbound orbital elements), the likelihood $\mathcal{L}$ of the MoG parameters is given by:

\[\mathcal{L}(\{w_k\}_M, \{\mu_k\}_M, \{\Sigma_k\}_M | \{\tilde U_k\}_S) = \prod_{i=1}^{S} \mathcal{C}_M(\tilde U_i)\]

To find the optimal parameters, we minimize the negative normalized log-likelihood:

\[-\frac{\log \mathcal{L}}{S} = -\frac{1}{S} \sum_{i=1}^{S} \log \mathcal{C}_M(\tilde U_i)\]

This quantity represents the average negative log-probability density of the data under the model. Minimizing this value is equivalent to maximizing the likelihood that the observed data was generated by the MoG model.

ngauss

Number of fitting MND.

Type:: int

nvars

Number of variables in each MND.

Type:: int

mog

Fitting object of the class MoG. This object will have the result of the fitting procedure.

Type:: MixtureOfGaussians

solution

Once the fitting is completed the solution object is returned. Attributes include:

fun: value of the function in the minimum.
x: value of the minimization parameters in the minimum.
nit: number of iterations of the minimization algorithm.
nfev: number of evaluations of the function.
success: if True it implies that the minimization fullfills all the conditions.

Type:: scipy.optimize.OptimizeResult

Ndim

Number of mus.

Type:: int

Ncorr

Number of correlations.

Type:: int

minparams

Array of minimization parameters at any stage in minimization process.

Type:: numpy.ndarray

scales

Array of scales to convert minparams to uparams (unbound) and viceversa.

Type:: numpy.ndarray

uparams

Array of unbound minimization parameters.

Type:: numpy.ndarray

Notes

Objects of this class must be always initialized with the number of gaussians and the number of random variables.

Examples

Example: Fitting a synthetic dataset.

>>> # 1. Generate synthetic data
>>> np.random.seed(1)
>>> weights = [[1.0]]
>>> mus = [[1.0, 0.5, -0.5]]
>>> sigmas = [[1, 1.2, 2.3]]
>>> angles = [[10*Angle.Deg, 30*Angle.Deg, 20*Angle.Deg]]
>>> Sigmas = Stats.calc_covariance_from_rotation(sigmas, angles)
>>> MoG = MixtureOfGaussians(mus=mus, weights=weights, Sigmas=Sigmas)
>>> data = MoG.rvs(10000)

>>> # 2. Initialize the fitter
>>> F = FitMoG(ngauss=MoG.ngauss, nvars=MoG.nvars)

>>> # 3. Set bounds (optional but recommended)
>>> bounds = F.set_bounds(bounds=(0.1, 0.9*F._sigmax))

>>> # 4. Run the fit
>>> F.fit_data(data, verbose=0, tol=1e-3, options=dict(maxiter=100, disp=True), method=None, bounds=bounds)

>>> # 5. Inspect results
>>> print(F.mog)
>>> G = F.plot_fit(figsize=3, hargs=dict(bins=30, cmap='YlGn'), sargs=dict(s=0.5, edgecolor='None', color='r'))

>>> # 6. Save/Load results
>>> F.save_fit("/tmp/fit.pkl", useprefix=False)
>>> F._load_fit("/tmp/fit.pkl")

fit_data(**kwargs)

log_l(data=None)[source]

Value of the -log(Likelihood). Ex. F.log_l(data)

Parameters:: data (numpy.ndarray, optional) – Array with data (Nsam x nvars). If None, uses self.data.
Returns:: log_l – Value of the -log(Likelihood).
Return type:: float

norm_log_l(data=None)[source]

Calculate the normalized log-likelihood of the data.

Parameters:: data (numpy.ndarray, optional) – Data for which to calculate the log-likelihood. If None, uses self.data.
Returns:: Normalized log-likelihood per point (-log L/N).
Return type:: float

plot_fit(N=10000, figsize=2, properties=None, ranges=None, hargs=None, sargs=None, pargs=True, cargs=None, marginals=False)[source]

Plot the result of the fitting procedure. Ex. F.plot_fit(figsize=3, hargs=dict(bins=30), sargs=dict(s=0.5, color=’r’))

Parameters:

N (int, optional) – number of points used to build a representation of the marginal distributions (default 10000).
figsize (int, optional) – Size of each axis (default 2).
properties (list or dict, optional) – Property names or MultiPlot-style properties. List: each element as axis label, range=None. Dict: same as MultiPlot (keys with ‘label’ and optional ‘range’).
ranges (list, optional) – Ranges per variable; used only when properties is a list or None.
hargs (dict, optional) – Dictionary with options for the hist2d function (or 1D hist when nvars=1). If provided, the sample is shown as a histogram.
sargs (dict, optional) – Dictionary with options for the scatter plot. If provided, the sample is shown as a scatter plot.
pargs (dict, optional) – Dictionary with options for the PDF plot. If None (default), defaults are used (color=’k’, lw=2, label=’PDF’). To suppress PDF, pass an empty dict is not sufficient, but pargs is default.
cargs (dict, optional) – Dictionary with options for the contour plot (mog_contour). If None (default), contours are not plotted.
marginals (bool, optional) – Include marginal distributions on diagonal panels (default False).

Returns:

G – Graphic handle.

Return type:

matplotlib.figure.Figure or MultiPlot

Examples

>>> F = FitMoG(1, 3)
>>> F.fit_data(data, verbose=0, tol=1e-3, options=dict(maxiter=100, disp=True))
>>> G = F.plot_fit(figsize=3, hargs=dict(bins=30, cmap='YlGn'), sargs=dict(s=0.5, edgecolor=None, color='r'))

pmap(minparams)[source]

Mapping routine used in sample_mog_likelihood. Mapping may change depending on the complexity of the parameters to be minimized. Ex. stdcorr = F.pmap(minparams)

Parameters:: minparams (numpy.ndarray) – Minimization parameters.
Returns:: stdcorr – Flatten parameters with correlations.
Return type:: numpy.ndarray

quality_of_fit(data=None, n_sim=None, resample=False, plot_qq=False, figsize=5, ax=None, plot_line=True, title='Q-Q Plot', verbose=False)[source]

Calculate goodness-of-fit statistics using the Kolmogorov-Smirnov distance and the Coefficient of Determination ($R^2$) relative to the identity line.

Optionally generates a Q-Q plot to visually assess the quality of the fit.

The Kolmogorov-Smirnov Statistic ($D_{KS}$) measures the maximum absolute distance between the Cumulative Distribution Function (CDF) of the observed data and that of the model. It represents the “worst point” of local deviation.

\[D_{KS} = \sup_{x} | F_{obs}(x) - F_{model}(x) |\]

The range of the statistic is $D_{KS} \in [0, 1]$, where an acceptable value is $D_{KS} \approx 0$ (ideally < 0.05).

The Q-Q plot is a non-parametric visual diagnostic tool used to assess whether a set of observed data follows a specific theoretical distribution. It compares the distribution of the Negative Log-Likelihood ($S = -2 \ln \mathcal{L}$) of the real data against the model’s prediction.

If the fit is ideal, the points should align on the identity line $y=x$. Deviations indicate:

Curved ends: tail deviations (“Heavy” or “Light” tails).
Shift or slope change: systematic errors in mean or variance.

The Coefficient of Determination on Identity ($R^2$) measures the variance explained by the model with respect to the ideal line $y=x$. It penalizes both dispersion and bias.

\[R^2 = 1 - \frac{\sum (y_{obs} - x_{model})^2}{\sum (y_{obs} - \bar{y}_{obs})^2}\]

The range is $R^2 \in (-\infty, 1]$, where an acceptable value is $R^2 \approx 1$ (ideally > 0.9).

Parameters:

data (numpy.ndarray, optional) – Observed data. If None, uses self.data.
n_sim (int, optional) – Number of synthetic data points to generate. Default is 10 * len(data).
resample (bool, optional) – If True, regenerate the synthetic sample even if one exists.
plot_qq (bool, optional) – If True, generate the Q-Q plot. Default is False.
figsize (int or tuple, optional) – Size of the figure (if creating a new one).
ax (matplotlib.axes.Axes, optional) – Axes object to plot on. If None and plot_qq=True, creates a new figure.
plot_line (bool, optional) – If True, plot the identity line y=x.
title (str, optional) – Title of the plot.
verbose (bool, optional) – If True, print the results.

Returns:

stats – Dictionary containing: - ‘ks_dist’: Kolmogorov-Smirnov distance between S_obs and S_sim (lower is better) - ‘ks_pval’: p-value from KS test (note: not strictly valid as parameters are estimated) - ‘r2_identity’: R-squared coefficient regarding identity line (closer to 1 is better) - ‘ax’: The axes object if plot_qq=True, else None.

Return type:

dict

save_fit(objfile=None, useprefix=True, myprefix=None)[source]

Pickle the result of a fit. Ex. F.save_fit(objfile=’fit.pkl’)

Parameters:

objfile (str, optional) – Name of the file where the fit will be stored. If None, the name is set by the routine as FitMoG.pkl.
useprefix (bool, optional) – Use a prefix in the filename of the pickle file (default True). The prefix is normally {ngauss}mog_{hash}.
myprefix (str, optional) – Custom prefix.

Examples

If objfile=”fit.pkl”, the final filename will be fit-1mnd_asa33.pkl

set_bounds(boundw=None, bounds=None, boundr=None, boundsm=None)[source]

Set the minimization parameters. Ex. F.set_bounds(boundsm=((-2, 1), (-3, 0)))

Parameters:

boundw (tuple, optional) – Bound of weights (default (-np.inf, np.inf)).
bounds (tuple or list, optional) – Bounds of weights, mus, sigmas and rhos of each variable.
boundr (tuple, optional) – Bound of rhos (default (-np.inf, np.inf)).
boundsm (tuple of tuples, optional) – Bounds of averages (default (-np.inf, np.inf)). Normally the bounds on averages must be expressed in this way: boundsm = ((-min_1, max_1), (-min_2, max_2), …) Example for nvars = 2: boundsm = ((-2, 1), (-3, 0))

Returns:

bounds – Formatted bounds for minimization.

Return type:

tuple

set_initial_params(mus=None, sigmas=None, rhos=None, from_mog=None)[source]

Set initial values for the minimization parameters (means, standard deviations, correlations). Only the arguments provided are updated; the rest keep their current values. Ex. F.set_initial_params(mus=[[0.2, 0.3], [0.8, 0.7]], sigmas=[[0.1, 0.2], [0.1, 0.2]])

Parameters:

mus (array-like, optional) – Initial means. Shape (ngauss, nvars) or (nvars,) to use the same means for all components. Example: [0.2, 0.3] or [[0.2, 0.3], [0.8, 0.7]].
sigmas (array-like, optional) – Initial standard deviations. Shape (ngauss, nvars) or (nvars,) to use the same for all components.
rhos (array-like, optional) – Initial correlation coefficients (upper triangle). Shape (ngauss, Ncorr) or (Ncorr,) to use the same for all. Ncorr = nvars*(nvars-1)/2.
from_mog (MixtureOfGaussians, optional) – If provided, initialize parameters from this MoG object. Must have same ngauss and nvars as check. Any other provided arguments (mus, sigmas, rhos) will override values from from_mog.

Return type:

None

Examples

>>> F = FitMoG(ngauss=2, nvars=2)
>>> F.set_initial_params(mus=[0.2, 0.3], sigmas=[0.1, 0.1])  # same for both components
>>> F.set_initial_params(mus=[[0.2, 0.3], [0.8, 0.7]], sigmas=[[0.1, 0.2], [0.1, 0.2]])
>>> F.fit_data(data)

set_params(mu=0.5, sigma=1.0, rho=0.5)[source]

Set the value of the basic params (minparams, scales, etc.). It updates minparams, scales and uprams. Ex. F.set_params(mu=0.5, sigma=1.0, rho=0.5)

Parameters:

mu (float, optional) – Value of all initial mus (default 0.5).
sigma (float, optional) – Value of all initial sigmas (default 1.0).
rho (float, optional) – Value of all initial rhos (default 0.5).

Return type:

None

update_params(weights=None, mus=None, sigmas=None, rhos=None)[source]

Update the parameters of the MoG.

Parameters:

weights (array_like, optional) – Weights of the gaussians.
mus (array_like, optional) – Means of the gaussians.
sigmas (array_like, optional) – Standard deviations of the gaussians.
rhos (array_like, optional) – Correlation coefficients of the gaussians.

class multimin.MixtureOfGaussians(ngauss=0, nvars=0, params=None, stdcorr=None, weights=None, mus=None, Sigmas=None, domain=None, normalize_weights=True)[source]

Bases: MultiMinBase

The Mixture of Gaussians (MoG).

We conjecture that any multivariate distribution function $p(\tilde U):\Re^{N}\rightarrow\Re$, where $\tilde U:(u_1,u_2,u_3,\ldots,u_N)$ and $u_i$ are random variables, can be approximated with an arbitrary precision by a normalized linear combination of $M$ Multivariate Normal Distributions (MND):

\[p(\tilde U) \approx \mathcal{C}_{M,k}(\tilde U; \{w_k\}_M, \{\mu_k\}_M, \{\Sigma_k\}_M) \equiv \sum_{i=1}^{M} w_i\;\mathcal{N}(\tilde U; \tilde \mu_i, \Sigma_i)\]

where the multivariate normal $\mathcal{N}(\tilde U; \tilde \mu, \Sigma)$ with mean vector $\tilde \mu$ and covariance matrix $\Sigma$ is given by:

\[\mathcal{N}(\tilde U; \tilde \mu, \Sigma) = \frac{1}{\sqrt{(2\pi)^{k} \det \Sigma}} \exp\left[-\frac{1}{2}(\tilde U - \tilde \mu)^{\rm T} \Sigma^{-1} (\tilde U - \tilde \mu)\right]\]

The covariance matrix $\Sigma$ elements are defined as $\Sigma_{ij} = \rho_{ij}\sigma_{i}\sigma_{j}$, where $\sigma_i$ is the standard deviation of $u_i$ and $\rho_{ij}$ is the correlation coefficient among variable $u_i$ and $u_j$ ($-1<\rho_{ij}<1$, $\rho_{ii}=1$).

The normalization condition on $p(\tilde U)$ implies that the set of weights $\{w_k\}_M$ are also normalized, i.e., $\sum_i w_i=1$.

Truncated multivariate distributions

In real problems the domain of the variables is often bounded. Starting from the unbounded multivariate normal $\mathcal{N}_k(\tilde U; \tilde \mu, \Sigma)$, let $T \subset \{1,\ldots,k\}$ be the set of indices of truncated variables and $a_i < b_i$ the truncation bounds for $i \in T$. The truncation region is

\[A_T = \left\{ \tilde U \in \mathbb{R}^k \;:\; a_i \le \tilde U_i \le b_i \;\; \forall i \in T \right\},\]

with the remaining coordinates $i \notin T$ unbounded. The partially truncated multivariate normal is

\[\mathcal{TN}_T(\tilde U; \tilde\mu, \Sigma, \mathbf{a}_T, \mathbf{b}_T) = \frac{\mathcal{N}(\tilde U; \tilde\mu, \Sigma) \, \mathbf{1}_{A_T}(\tilde U)} {Z_T(\tilde\mu, \Sigma, \mathbf{a}_T, \mathbf{b}_T)},\]

where $\mathbf{1}_{A_T}$ is the indicator of $A_T$ and the normalization constant is

\[Z_T(\tilde\mu, \Sigma, \mathbf{a}_T, \mathbf{b}_T) = \int_{A_T} \mathcal{N}(\tilde T; \tilde\mu, \Sigma) \, d\tilde T = \mathbb{P}_{\tilde T \sim \mathcal{N}(\tilde\mu,\Sigma)}\!\left( \tilde T \in A_T \right).\]

When a finite domain is set (e.g. domain=[[0, 1], None]), the MoG uses truncated normals for the bounded variables so that the PDF is zero outside the domain and the mixture remains normalized.

ngauss

Number of composed MND.

Type:: int

nvars

Number of random variables.

Type:: int

mus

Array with average (mu) of random variables (ngauss x nvars).

Type:: numpy.ndarray

weights

Array with weights of each MND (ngauss). NOTE: These weights are normalized at the end.

Type:: numpy.ndarray

sigmas

Standard deviation of each variable(ngauss x nvars).

Type:: numpy.ndarray

rhos

Elements of the upper triangle of the correlation matrix (ngauss x nvars x (nvars-1)/2).

Type:: numpy.ndarray

Sigmas

Array with covariance matrices for each MND (ngauss x nvars x nvars).

Type:: numpy.ndarray

params

Parameters of the distribution in flatten form including symmetric elements of the covariance matrix (ngauss*(1+nvars+nvars*(nvars+1)/2)).

Type:: numpy.ndarray

stdcorr

Parameters of the distribution in flatten form, including upper triangle of the correlation matrix (ngauss*(1+nvars+nvars*(nvars+1)/2)).

Type:: numpy.ndarray

Notes

There are several ways of initialize a MoG:

Providing: ngauss and nvars In this case the class is instantiated with zero means, unitary dispersion and covariance matrix equal to Ngasus identity matrices nvars x nvars.
Providing: params, nvars In this case you have a flatted version of the parametes (weights, mus, Sigmas) and want to instantiate the system. All parameters are set and no other action is required.
Providing: stdcorr, nvars In this case you have a flatted version of the parametes (weights, mus, sigmas, rhos) and want to instantiate the system. All parameters are set and no other action is required.
Providing: weights, mus, Sigmas (optional), domain (optional), normalize_weights (optional) In this case the basic properties of the MoG are set. For univariate (one variable), mus may be a 1-D array of means, e.g. [0, 2], and Sigmas a 1-D array of variances, e.g. [1.0, 0.25]. domain: list of length nvars; each element is None (unbounded) or [low, high] (finite support). Example: [None, [0, 1], None] for variable 1 in [0, 1]. normalize_weights: if True (default), weights are scaled so sum(weights)=1 (proper PDF). If False, weights are used as-is (sum may != 1); useful for function fitting with an overall scale (e.g. FitFunctionMoG).

Examples

Example 1: Initialization using explicit arrays for means and weights.

>>> # Define means for 2 Gaussian components in 2D
>>> mus = [[0, 0], [1, 1]]
>>> # Define weights (normalization is handled automatically)
>>> weights = [0.1, 0]
>>> # Create the MoG object
>>> MND1 = MixtureOfGaussians(mus=mus, weights=weights)
>>> # Set covariance matrices explicitly
>>> MND1.set_sigmas([[[1, 0.2], [0, 1]], [[1, 0], [0, 1]]])
>>> print(MND1)

Example 2: Initialization using a flattened parameter array.

>>> # Flattened parameters: [weights..., mus..., flattened_covariances...]
>>> params = [0.1, 0.9, 0, 0, 1, 1, 1, 0.2, 0.2, 1, 1, 0, 0, 1]
>>> # Create MoG object specifying number of variables
>>> MND2 = MixtureOfGaussians(params=params, nvars=2)
>>> print(MND2)

drop(index)[source]

Drop one of the gaussians. Ex. mog.drop(0)

Parameters:: index (int) – Index of the gaussian to be dropped.

get_function(print_code=True, decimals=6, type='python', properties=None, cmog=False)[source]

Return the source code of mog(X) and an executable function (type=’python’), or LaTeX code with parameters in begin{array} (type=’latex’). Ex. code, mog = MoG.get_function()

Parameters:

print_code (bool, optional) – If True (default), print the string to screen so it can be copied.
decimals (int, optional) – Number of decimal places for numeric literals (default 6).
type (str, optional) –
- 'python' (default): return Python source and callable.
- 'latex': return LaTeX formula with explicit parameters and matrices in begin{array}.
LaTeX output is single-line per formula so it can be pasted into README or other RST-derived long descriptions without triggering PyPI/twine markup errors.
properties (dict or sequence, optional) – If None (default), variable subscripts in parameters are numeric (mu_1, sigma_1, …). If a dict (e.g. from MultiPlot), its keys are used as variable names (mu_x, sigma_x, …). If a sequence, its elements are used in order. Length must match the number of variables.
cmog (bool, optional) – If True and type=’python’, generate code using the C-optimized routines cmog.nmd_c and cmog.tnmd_c. Default False.

Returns:

First element: source code (Python or LaTeX). Second: callable mog if type=’python’, else None.

Return type:

tuple (str, callable or None)

Examples

>>> code, mog = MoG.get_function()
>>> latex_str, _ = MoG.get_function(type='latex')

pdf(X)[source]

Compute the PDF. Ex. mog.pdf([1.0, 0.5])

Parameters:: X (float or numpy.ndarray) – Point in the nvars-dimensional space (scalar when nvars=1, or length nvars).
Returns:: p – PDF value at X.
Return type:: float

plot_pdf(properties=None, ranges=None, figsize=3, grid_size=200, cmap='Spectral_r', colorbar=False, marginals=False)[source]

Plot only the PDF of the MoG.

For univariate distributions (nvars==1) this produces a single curve. For multivariate distributions (nvars>=2) this produces density panels using MultiPlot. For nvars>2 the density shown in each panel is a 2D slice of the full PDF, evaluated at the weighted mean for the remaining variables.

Parameters:

properties (list or dict, optional) – Axis specification, same format accepted by plot_sample().
ranges (list, optional) – Ranges per variable; used when properties is a list or None.
figsize (int, optional) – Size of each axis (default 3).
grid_size (int, optional) – Number of grid points per axis for density evaluation (default 200).
cmap (str, optional) – Matplotlib colormap name (default 'Spectral_r').
colorbar (bool, optional) – Include a colorbar for the first density panel (default False).
marginals (bool, optional) – Include marginal distributions on diagonal panels (default False).

Returns:

G – Handle to the density plot grid.

Return type:

MultiPlot

plot_sample(data=None, N=10000, properties=None, ranges=None, figsize=2, sargs=None, hargs=None, marginals=False)[source]

Plot a sample of the MoG. Ex. plot_sample(N=1000, sargs=dict(s=0.5))

Parameters:

data (numpy.ndarray, optional) – Data to plot. If None it generate a sample.
N (int, optional) – Number of points to generate the sample (default 10000).
properties (list or dict, optional) – Property names or MultiPlot-style properties. List (e.g. properties=[“x”,”y”]): each element is used as axis label, range=None. Dict: same as MultiPlot, e.g. properties=dict(x=dict(label=r”$x$”, range=None), y=dict(label=r”$y$”, range=[-1,1])).
ranges (list, optional) – Ranges per variable; used only when properties is a list or None. Ex. ranges=[[-3,3],[-5,5]].
figsize (int, optional) – Size of each axis (default 2).
sargs (dict, optional) – Dictionary with options for the scatter plot. Default: dict(s=0.5, edgecolor=None, color=’b’).
hargs (dict, optional) – Dictionary with options for the hist2d function. Ex. hargs=dict(bins=50).
marginals (bool, optional) – Include marginal distributions on diagonal panels (default False).

Returns:

G – Graphic handle. If nvars = 2, it is a figure object, otherwise is a MultiPlot instance.

Return type:

matplotlib.figure.Figure or MultiPlot

Examples

Example 1: Plotting a generated sample.

>>> # Generate and plot 10000 points
>>> G = MoG.plot_sample(N=10000, sargs=dict(s=1, c='r'))
>>> # Plot with histogram bins
>>> G = MoG.plot_sample(N=1000, sargs=dict(s=1, c='r'), hargs=dict(bins=20))

Example 2: Plotting for a 2D distribution.

>>> MoG = MixtureOfGaussians(ngauss=1, nvars=2)
>>> fig = MoG.plot_sample(N=1000, hargs=dict(bins=20), sargs=dict(s=1, c='r'))

Example 3: Defining and plotting a complex MoG.

>>> # Define parameters for 2 Gaussians in 2D
>>> params = [0.1, 0.9, 0.0, 0.0, 1.0, 1.0, 1.0, 0.2, 1.0, 1.0, 0.0, 1.0]
>>> MND2 = MixtureOfGaussians(params=params, nvars=2)
>>> # Print the object details
>>> print(MND2)
>>> # Calculate PDF at a point
>>> print(MND2.pdf([1, 1]))

rvs(**kwargs)

sample_mog_likelihood(uparams, data=None, pmap=None, tset='stdcorr', scales=[], verbose=0)[source]

Compute the negative value of the logarithm of the likelihood of a sample. Ex. mog.sample_mog_likelihood(uparams, data=data, pmap=pmap)

Parameters:

uparams (numpy.ndarray) – Minimization parameters (unbound).
data (numpy.ndarray, optional) – Data for which log_l is computed.
pmap (function, optional) – Routine to map from minparams to params or stdcorr. Example: >>> def pmap(minparams): … stdcorr = np.array([1] + list(minparams)) … stdcorr[-1:] -= 1 … return stdcorr
tset (str, optional) – Type of minimization parameters. Values “params”, “stdcorr” (default “stdcorr”).
scales (list, optional) – List of scales for transforming uparams (unbound) in minparams (natural scale).
verbose (int, optional) – Verbosity level (0, none, 1: input parameters, 2: full definition of the MoG) (default 0).

Returns:

log_l – Negative log-likelihood.

Return type:

float

set_domain(domain)[source]

Set the support domain for each variable (finite or infinite). Ex. set_domain([None, [0, 1], None])

Parameters:: domain (list of length nvars) – Each element is None (unbounded) or [low, high] (finite interval). Example: [None, [0, 1], None] for variable 1 bounded to [0, 1].

set_params(params, nvars)[source]

Set the properties of the MoG from flatten params. After setting it generate flattend stdcorr and normalize weights. Ex. set_params(params, nvars=2)

Parameters:

params (list or numpy.ndarray) – Flattened parameters.
nvars (int) – Number of variables.

set_sigmas(Sigmas)[source]

Set the value of list of covariance matrices. After setting Sigmas it update params and stdcorr. Ex. set_sigmas([[[1, 0.2], [0.2, 1]]])

Parameters:: Sigmas (list or numpy.ndarray) – Array of covariance matrices. For univariate (nvars=1), may be a 1-D array of variances, e.g. [1.0, 0.25] for two components.

set_stdcorr(stdcorr, nvars)[source]

Set the properties of the MoG from flatten stdcorr. After setting it generate flattened params and normalize weights. Ex. set_stdcorr(stdcorr, nvars=2)

Parameters:

stdcorr (list or numpy.ndarray) – Flattened standard deviations and correlations.
nvars (int) – Number of variables.

tabulate(sort_by='weight', return_df=True, type='df', properties=None)[source]

Build a table of MoG parameters: weights, means (mu_i), standard deviations (sigma_i), and correlations (rho_ij, i < j). Diagonal rho_ii are 1 by definition and omitted. Ex. df = MoG.tabulate(sort_by=’weight’)

Parameters:

sort_by (str, optional) –
How to order the rows (one row per Gaussian component): - 'weight': by weight descending (heaviest first). - 'distance': by Euclidean distance of the mean vector to the origin ascending

(closest to origin first).
return_df (bool, optional) – If True (default), return a pandas DataFrame. If False, print the table and return None. Ignored when type=’latex’.
type (str, optional) –
- 'df' (default): return pandas DataFrame (or print and return None).
- 'latex': return a LaTeX tabular string suitable for papers.
properties (dict or sequence, optional) – If None (default), column headers use numeric subscripts (mu_1, sigma_1, …). If a dict (e.g. from MultiPlot), its keys are used (mu_x, sigma_x, …). If a sequence, its elements are used in order.

Returns:

DataFrame when type=’df’ and return_df=True; LaTeX string when type=’latex’; None when type=’df’ and return_df=False.

Return type:

pandas.DataFrame, str, or None

update_params(weights=None, mus=None, sigmas=None, rhos=None)[source]

Update MoG parameters in-place using FitMoG-like syntax.

This method updates the internal parameters of the composed distribution (means, standard deviations, and correlations) using the same broadcasting rules as FitMoG.set_initial_params().

Only arguments provided are updated; other parameters keep their current values.

Parameters:

weights (array-like, optional) – Mixture weights. Must have length ngauss (or length 1 when ngauss==1). If normalize_weights is True (default), weights are scaled so that sum(weights)=1; in all cases weights are kept non-negative.
mus (array-like, optional) – Means. Shape (ngauss, nvars) or (nvars,) (same for all components).
sigmas (array-like, optional) – Standard deviations. Shape (ngauss, nvars) or (nvars,).
rhos (array-like, optional) – Correlation coefficients (upper triangle). Shape (ngauss, Ncorr) or (Ncorr,) where Ncorr = nvars*(nvars-1)/2.

Return type:

None

class multimin.MultiPlot(properties, figsize=3, fontsize=10, direction='out', marginals=False)[source]

Bases: MultiMinBase

Create a grid of plots showing the projection of a N-dimensional data.

Parameters:

properties (dict) –
List of properties to be shown, dictionary of dictionaries (N entries). Keys are label of attribute, ex. “q”. Dictionary values:
- label: label used in axis, string
- range: range for property, tuple (2)
figsize (int, optional) – Base size for panels (the size of figure will be M x figsize), default 3.
fontsize (int, optional) – Base fontsize, default 10.
direction (str, optional) – Direction of ticks in panels, default ‘out’.

N

Number of properties.

Type:: int

M

Size of grid matrix (M=N-1).

Type:: int

fw

Figsize.

Type:: int

fs

Fontsize.

Type:: int

fig

Figure handle.

Type:: matplotlib.figure.Figure

axs

Matrix with subplots, axes handles (MxM).

Type:: numpy.ndarray

axp

Matrix with subplots, dictionary of dictionaries.

Type:: dict

properties

List of properties labels, list of strings (N).

Type:: list

tight_layout()[source]: Tight layout if no constrained_layout was used.

set_labels(**args)[source]: Set labels parameters.

set_ranges()[source]: Set ranges in panels according to ranges defined in dparameters.

set_tick_params(**args)[source]: Set tick parameters.

sample_hist(data, colorbar=False, \*\*args)[source]: Create a 2d-histograms of data on all panels of the MultiPlot.

sample_scatter(data, \*\*args)[source]: Scatter plot on all panels of the MultiPlot.

mog_pdf(mog, \*\*args)[source]: Plot the PDF of a MoG on all panels of the MultiPlot.

mog_contour(mog, \*\*args)[source]: Plot the contours of a MoG on all panels of the MultiPlot.

mog_contour(mog, grid_size=200, **args)[source]

Plot the contours of a MoG on all panels of the MultiPlot. Ex. G.mog_contour(mog, levels=5, cmap=’Reds’, margs=dict(color=’blue’))

Parameters:

mog (MixtureOfGaussians) – MoG object to plot.
grid_size (int, optional) – Number of points for the grid (default 200).
**args (dict) – Arguments for contour function. Can include ‘margs’ dict with arguments for marginal plots. If margs=None, marginals are not drawn.

mog_pdf(mog, grid_size=200, **args)[source]

Plot the PDF of a MoG on all panels of the MultiPlot. Ex. G.mog_pdf(mog, color=’k’, lw=2, margs=dict(color=’blue’))

Parameters:

mog (MixtureOfGaussians) – MoG object to plot.
grid_size (int, optional) – Number of points for the grid (default 200).
**args (dict) – Arguments for the plot function (e.g. color, linewidth). Can include ‘margs’ dict with arguments for marginal plots. If margs=None, marginals are not drawn.

reset_ranges()[source]: Reset ranges to match the data limits.

sample_hist(data, colorbar=False, **args)[source]

Create a 2d-histograms of data on all panels of the MultiPlot. Ex. G.sample_hist(data, bins=100, cmap=’viridis’, margs=dict(color=’blue’))

Parameters:

data (numpy.ndarray) – Data to be histogramed (n=len(data)), numpy array (nxN).
colorbar (bool, optional) – Include a colorbar? (default False).
**args (dict) – All arguments of hist2d method. Can include ‘margs’ dict with arguments for marginal plots. If margs=None, marginals are not drawn.

Returns:

hist – List of histogram instances.

Return type:

list

Examples

>>> properties = {
...     'Q': {'label': r"$Q$", 'range': None},
...     'E': {'label': r"$C$", 'range': None},
...     'I': {'label': r"$I$", 'range': None},
... }
>>> G = mm.MultiPlot(properties, figsize=3)
>>> hargs = dict(bins=100, cmap='viridis', margs=dict(color='blue'))
>>> hist = G.sample_hist(udata, **hargs)

sample_scatter(data, nbins=20, **args)[source]

Scatter plot on all panels of the MultiPlot. Ex. G.sample_scatter(data, s=0.2, color=’r’, margs=dict(color=’b’))

Parameters:

data (numpy.ndarray) – Data to be histogramed (n=len(data)), numpy array (nxN).
nbins (int, optional) – Number of bins for marginal histograms (default 20).
**args (dict) – All arguments of scatter method. Can include ‘margs’ dict with arguments for marginal plots. If margs=None, marginals are not drawn.

Returns:

scatter – List of scatter instances.

Return type:

list

Examples

>>> # With marginals (blue histogram)
>>> sargs = dict(s=0.2, edgecolor='None', color='r', margs=dict(color='b'))
>>> hist = G.sample_scatter(udata, **sargs)
>>> # Without marginals
>>> sargs = dict(s=0.2, edgecolor='None', color='r', margs=None)
>>> hist = G.sample_scatter(udata, **sargs)

set_labels(**args)[source]

Set labels parameters. Ex. set_labels(fontsize=12)

Parameters:: **args (dict) – Common arguments of set_xlabel, set_ylabel and text.

set_ranges()[source]: Set ranges in panels according to ranges defined in dparameters.

set_tick_params(**args)[source]

Set tick parameters. Ex. set_tick_params(labelsize=10)

Parameters:: **args (dict) – Same arguments as tick_params method.

tight_layout()[source]: Tight layout if no constrained_layout was used.

class multimin.Stats[source]

Bases: MultiMinBase

Abstract class with useful routines

calc_correlations_from_covariances()[source]

Compute the standard deviations and corresponding correlation coefficients given a set of covariance matrices. Ex. sigmas, rhos = Stats.calc_correlations_from_covariances(Sigmas)

Parameters:

Sigmas (numpy.ndarray) – Array of covariance matrices (ngauss x nvars x nvars).

Returns:

sigmas (numpy.ndarray) – Array of standard deviations (ngauss x nvars).
rhos (numpy.ndarray) – Array of correlation coefficients (ngauss x nvars * (nvars-1) / 2).

Examples

>>> Sigmas = [
...     [[1. , 0.2, 0.6],
...      [0.2, 4. , 1.8],
...      [0.6, 1.8, 9. ]]
... ]
>>> sigmas, rhos = Stats.calc_correlations_from_covariances(Sigmas)
>>> print(sigmas)
[1. 2. 3.]
>>> print(rhos)
[[0.1 0.2 0.3]]

calc_covariance_from_correlations(rhos)[source]

Compute covariance matrices from the standard deviations and correlations (rho). Ex. Stats.calc_covariance_from_correlations(sigmas, rhos)

Parameters:

sigmas (numpy.ndarray) – Array of values of standard deviation for variables (ngauss x nvars).
rhos (numpy.ndarray) – Array with correlations (ngauss x nvars x (nvars-1)/2).

Returns:

Sigmas – Array with covariance matrices corresponding to these sigmas and rhos (ngauss x nvars x nvars).

Return type:

numpy.ndarray

Examples

>>> import numpy as np
>>> sigmas = np.array([[1, 2, 3]])
>>> # rho_12, rho_13, rho_23
>>> rhos = np.array([[0.1, 0.2, 0.3]])
>>> S = Stats.calc_covariance_from_correlations(sigmas, rhos)
>>> print(S)
[[[1.  0.2 0.6]
  [0.2 4.  1.8]
  [0.6 1.8 9. ]]]

This is equivalent to:

>>> rho = rhos[0]
>>> sigma = sigmas[0]
>>> R = np.eye(3)
>>> Stats.set_matrix_off_diagonal(R, rho)
>>> M = np.zeros((3, 3))
>>> for i in range(3):
...     for j in range(3):
...         M[i,j] = R[i,j] * sigma[i] * sigma[j]
>>> print(M)
[[1.  0.2 0.6]
 [0.2 4.  1.8]
 [0.6 1.8 9. ]]

calc_covariance_from_rotation(angles)[source]

Compute covariance matrices from the stds and the angles. Ex. Stats.calc_covariance_from_rotation(sigmas, angles)

Parameters:

sigmas (numpy.ndarray) – Array of values of standard deviation for variables (ngauss x 3).
angles (numpy.ndarray) – Euler angles expressing the directions of the principal axes of the distribution (ngauss x 3).

Returns:

Sigmas – Array with covariance matrices corresponding to these sigmas and angles (ngauss x 3 x 3).

Return type:

numpy.ndarray

flatten_symmetric_matrix()[source]

Given a symmetric matrix the routine returns the flatten version of the Matrix. Ex. Stats.flatten_symmetric_matrix(M)

Parameters:: M (numpy.ndarray) – Matrix (n x n).
Returns:: F – Flatten array (nx(n+1)/2).
Return type:: numpy.ndarray

Examples

>>> M = np.array([[1, 0.2], [0.2, 3]])
>>> F = Stats.flatten_symmetric_matrix(M)
>>> print(F)
[1.  0.2 3. ]

gen_index()[source]

Given a set of (normalized) probabilities, randomly generate an index n following the probabilities. For instance if we have 3 events with probabilities 0.1, 0.7, 0.2, gen_index will generate a number in the set (0,1,2) having those probabilities, ie. 1 will have 70%. Ex. Stats.gen_index([0.1, 0.7, 0.2])

Parameters:: probs (numpy.ndarray) – Probabilities (N), adimensional. NOTE: It should be normalized, ie. sum(probs)=1
Returns:: n – Index in the set [0,1,2,… len(probs)-1].
Return type:: int

Examples

>>> n = Stats.gen_index([0.1, 0.7, 0.2])

phi = 1.618033988749895

set_matrix_off_diagonal(off)[source]

Set a matrix with the terms of the off diagonal. Ex. Stats.set_matrix_off_diagonal(M, [0.1, 0.2, 0.3])

Parameters:

M (numpy.ndarray) – Matrix (n x n).
off (list or numpy.ndarray) – Terms off diagonal (n x (n-1) / 2).

Returns:

Implicitly the matrix M has now the off diagonal terms.

Return type:

None

Examples

>>> M = np.eye(3)
>>> off = [0.1, 0.2, 0.3]
>>> Stats.set_matrix_off_diagonal(M, off)
>>> print(M)
[[1. , 0.1, 0.2],
 [0.1, 1. , 0.3],
 [0.2, 0.3, 1. ]]

unflatten_symmetric_matrix(M)[source]

Given a flatten version of a matrix, returns the symmetric matrix. Ex. Stats.unflatten_symmetric_matrix(F, M)

Parameters:

F (numpy.ndarray) – Flatten array (n x (n+1)/2).
M (numpy.ndarray) – Matrix where the result will be stored (n x n).

Returns:

It return the results in matrix M.

Return type:

None

Examples

>>> F = [1, 0.2, 3]
>>> M = np.zeros((2, 2))
>>> Stats.unflatten_symmetric_matrix(F, M)
>>> print(M)
[[1.  0.2]
 [0.2 3. ]]

class multimin.Util[source]

Bases: MultiMinBase

This abstract class contains useful methods for the package.

DTIME = -1

DUTIME = []

TIME = 1771179770.4275477

TIMESTART = 1771179770.4275472

cos = <ufunc 'cos'>

static docstring_summary(doc)[source]

Extract the first line or paragraph of a docstring as summary.

Stops at common section headers (Parameters, Returns, etc.) so that the description does not include parameter lists.

Parameters:: doc (str or None) – The docstring.
Returns:: One-line summary or empty string if no docstring.
Return type:: str

el_time(start=False)[source]

Compute the time elapsed since last call of this routine. The displayed time is preseneted in the more convenient unit, ns (nano seconds), us (micro seconds), ms (miliseconds), s (seconds), min (minutes), h (hours), d (days). Ex. el_time(), el_time(verbose=0), el_time(start=True)

Parameters:

verbose (int or bool, optional) – Show the time in screen (default 1).
start (int or bool, optional) – Compute time from program start (default 0).

Returns:

dt (float) – Elapsed time in seconds.
dtu_unit (list) – List containing [time in units, unit string].

Examples

>>> Util.el_time() # basic usage (show output)
>>> Util.el_time(verbose=0) # no output
>>> Util.el_time(start=True) # measure elapsed time since program start
>>> print(Util.DTIME, Util.DUTIME) # show values of elapsed time

error_msg(msg)[source]

Add a custom message msg to an error handle. Ex. error_msg(error, ‘custom message’)

Parameters:

error (Exception) – Error handle (eg. except ValueError as error).
msg (str) – Message to add to error.

exp = <ufunc 'exp'>

static f2u(x, s)[source]

Convert from a finite interval [0,s] to an unbound one [-inf,inf]. Ex. f2u(0.5, 1.0)

Parameters:

x (float or array_like) – Value in the interval [0,s].
s (float) – Scale (upper limit of the interval).

Returns:

u – Unbound value.

Return type:

float or array_like

static get_data(filename)[source]

Get the full path of the filename which is one of the datafiles provided with the package. Ex. get_data(“nea_extended.json.gz”)

Parameters:: filename (str) – Name of the data file.
Returns:: path – Full path to package datafile.
Return type:: str

Examples

>>> from multimin.util import Util
>>> path=Util.get_data("nea_extended.json.gz")

log = <ufunc 'log'>

static logit_transform(data, scales, inverse=False)[source]

Apply logit (log-odds) transformation to data, or its inverse. This provides a vectorized alternative to utilizing t_if in a loop.

Parameters:

data (array_like) – Data to transform (N, nvars) or (nvars,).
scales (array_like) – Scales/bounds for each variable. If scale[i] > 0, the variable is transformed. If scale[i] <= 0, variable is unchanged.
inverse (bool, optional) – If True, apply inverse transformation (False is finite -> unbounded). Default False.

Returns:

out – Transformed data (same shape as input).

Return type:

ndarray

mantisa_exp()[source]

Calculate the mantisa and exponent of a number. Ex. m, e = Util.mantisa_exp(234.5)

Parameters:

x (float) – Number.

Returns:

man (float) – Mantisa.
exp (float) – Exponent.

Examples

>>> m, e = Util.mantisa_exp(234.5)
# returns m=2.345, e=2
>>> m, e = Util.mantisa_exp(-0.000023213)
# return m=-2.3213, e=-5

static nmd(X, mu, Sigma)[source]

PDF of a multivariate normal at x; mu=mean vector, Sigma=covariance matrix. Ex. Util.nmd(X, mu, Sigma)

Parameters:

X (array-like) – Point(s) at which to evaluate the PDF. Shape (n,) or (N, n).
mu (array-like) – Mean vector, shape (n,).
Sigma (array-like) – Covariance matrix, shape (n, n). For univariate (n=1), a 1x1 matrix or scalar variance.

Returns:

PDF value(s). Scalar if x is 1D, array if x is 2D.

Return type:

float or ndarray

static props_to_properties(properties, nvars, ranges=None)[source]

Convert properties (list or MultiPlot-style dict) to properties dict for MultiPlot. Ex. Util.props_to_properties(properties, nvars, ranges)

If properties is a dict: use as-is (MultiPlot-style); each value must have ‘label’ and optionally ‘range’. Keys define variable names; first nvars keys are used.
If properties is a list or sequence: use each element as variable name and as axis label, with range=None unless ranges is provided (backward compatible).
If properties is None: use ascii_letters and ranges from the ranges argument.

Parameters:

properties (dict, list or None) – Properties specification.
nvars (int) – Number of variables.
ranges (list, optional) – Ranges for each variable.

Returns:

Properties dictionary.

Return type:

dict

sin = <ufunc 'sin'>

static t_if(p, s, f)[source]

Transform a set of parameters using a transformation function f and scales s. This routine allows the conversion from a finite interval [0,s] to an unbound one [-inf,inf] (using f=Util.f2u) or vice versa (using f=Util.u2f). Ex. Util.t_if(minparams, scales, Util.f2u)

Parameters:

p (list or array_like) – Parameters to transform.
s (list or array_like) – Scales for each parameter. If s[i] > 0, the transformation is applied. If s[i] == 0, the parameter is unchanged.
f (function) – Transformation function (e.g. Util.f2u or Util.u2f).

Returns:

tp – Transformed parameters.

Return type:

list

Examples

>>> scales = [0, 0, 10, 10, 1]
>>> minparams = [0.0, 0.0, 1, 1, 0.7]
>>> uparams = Util.t_if(minparams, scales, Util.f2u)
>>> print(uparams)
[0.0, 0.0, -2.197224577336219, -2.197224577336219, 0.8472978603872034]

static timer(func)[source]

static tnmd(X, mu, Sigma, a, b, Z=None)[source]

PDF of a truncated (multivariate) normal at X; single evaluation like nmd. Uses scipy.stats.truncnorm for 1D (one call). For nD, returns normal PDF in the box [a, b] divided by the normalization constant Z. Ex. Util.tnmd(X, mu, Sigma, a, b)

Parameters:

X (array-like) – Point(s). Shape (n,) or (N, n).
mu (array-like) – Mean vector, shape (n,).
Sigma (array-like) – Covariance matrix (n, n). For n=1, scalar variance or 1x1.
a (array-like or float) – Truncation bounds (lower, upper). For 1D, scalars; for nD, vectors of length n.
b (array-like or float) – Truncation bounds (lower, upper). For 1D, scalars; for nD, vectors of length n.
Z (float, optional) – Normalization constant P(a < X < b). For nD, pass to avoid extra computation.

Returns:

PDF value(s). Zero outside [a, b].

Return type:

float or ndarray

static true_anomaly_to_mean_anomaly(e, f)[source]

Convert true anomaly to mean anomaly. Ex. true_anomaly_to_mean_anomaly(0.1, 0.5)

Parameters:

e (float) – Eccentricity.
f (float) – True anomaly.

Returns:

M – Mean anomaly.

Return type:

float

static u2f(t, s)[source]

Convert from an unbound interval [-inf,inf] to a finite one [0,s]. Ex. u2f(0.0, 1.0)

Parameters:

t (float or array_like) – Unbound value.
s (float) – Scale (upper limit of the interval).

Returns:

x – Value in the interval [0,s].

Return type:

float or array_like

multimin.multimin_watermark(ax, frac=0.25, alpha=1)[source]

Add a water mark to a 2d or 3d plot.

Parameters:

ax: Class axes:
Axe where the pryngles mark will be placed.