API Reference

Gaussian Mixture

class pyest.gm.gm.GaussianMixture(w, m, cov, cov_type='full', Seig=None)

Bases: object

__init__(w, m, cov, cov_type='full', Seig=None)

Parameters:

• w (array_like) – (nC,) array of mixand weights

• m (array_like) – (nC,nx) array of mixand means

• cov (array_like) – (nC,nx,nx) array of covariances. Covariances can be specified in different ways, as specified by the optional cov_type argument

• Optional

• --------

• cov_type (str) – type of covariance matrix. Defaults to ‘full’.

• Seig (3d array) – eigenvalues and eigenvectors of the covariance matrix. Defaults to None.

get_size()

Returns the number of components in the Gaussian mixture model.

get_w()

Returns the weights (1d array) of the Gaussian mixture model.

set_w(w)

Sets the weights (1d array) of the Gaussian mixture model.

get_m()

Returns the means (2d array (nC,nx)) of the Gaussian mixture model.

set_m(m)

Sets the means (2d array (nC,nx)) of the Gaussian mixture model.

get_P(ind=None)

Returns the covariance matrices (3d array (nC,nx,nx)) of the Gaussian mixture model.

set_P(P)

Sets the covariance matrices (3d array (nC,nx,nx)) of the Gaussian mixture model.

set_Seig(S)

Sets the covariance matrix eigenvalues and eigenvectors.

Schol

Seig

log_pdet

prec_sqrt

get_msize()

get vector length of mean

set_msize(msize)

set vector length of mean

get_comp(ind)

eval(x)

evaluate gm at n_eval points stored in numpy array x

Required

x: ndarray

(n_eval, nx) array of points at which to evaluate the GM

returns:

(n_eval,) array of GM evaluations at points x

rtype:

ndarray

mean()

return the mean of the distribution

Returns:

(nx,) distribution mean

Return type:

ndarray

cov()

Compute the covariance matrix of the distribution.

Returns:

(nx, nx) full covariance matrix of the distribution

Return type:

ndarray

cdf(x, allow_singular=False)

evaluate GM CDF at points x

comp_bounds(sigma_mult=3)

marginal_2d(dimensions=(0, 1))

return 2D marginalized GM

marginal_nd(dimensions)

return nD marginalization of GM over specified dimensions

pdf_2d(dimensions=(0, 1), res=100, xbnd=None, ybnd=None)

evaluate GM pdf in two dimensions

Optional

dimensions: tuple

dimensions to evaluate pdf over. Defaults to (0, 1)

res: int

resolution of the grid. Defaults to 100

xbnd: tuple

x-axis bounds. Determined automatically if not provided

ybnd: tuple

y-axis bounds. Determined automatically if not provided

returns:

• ndarray – (res, res) pdf values

• ndarray – (res, res) x values

• ndarray – (res, res) y values

pop(idx)

remove and return component by index

Required

idx: int

index of component to remove

returns:

(weight, mean, covariance) of removed component

rtype:

tuple

rvs(size=None, random_state=None)

Gaussian Mixture generated random variates

Parameters:

• size (int, optional) – Defining number of random variates (Default is 1)

• random_state (int, numpy.random.Generator, numpy.random.RandomState, optional) – Random state for reproducibility. Can be: - An integer seed - A numpy.random.Generator instance - A numpy.random.RandomState instance - None (default) for non-deterministic results

Returns:

Random variates

Return type:

ndarray

Notes

If size is not specified, a single (possibly-vector) random variate is returned. If size is specified, an ndarray of random variates is returned, even for size=1.

property size

Returns the number of components in the Gaussian mixture model.

property w

Returns the weights (1d array) of the Gaussian mixture model.

property m

Returns the means (2d array (nC,nx)) of the Gaussian mixture model.

property P

Returns the covariance matrices (3d array (nC,nx,nx)) of the Gaussian mixture model.

property dim

get vector length of mean

class pyest.gm.gm.GaussSplitOptions(L=3, lam=0.001, recurse_depth=inf, min_weight=0.001, state_idxs=(0, 1), spectral_direction_only=False, variance_preserving=True)

Bases: object

__init__(L=3, lam=0.001, recurse_depth=inf, min_weight=0.001, state_idxs=(0, 1), spectral_direction_only=False, variance_preserving=True)

Options for FoV splitting

Parameters:

• L (int, optional) – number of components to spit Gaussian into. Default is 3

• lam (float, optional) – Gauss split optimization regularization parameter. Lower values result in less component spread and larger component covariances. This value should be decreased as L increases. Default is 1e-3

• recurse_depth (int, optional) – maximum recursion depth. Default is inf

• min_weight (float, optional) – smallest component weight that is considered for splitting. Default is 1e-3

• state_idxs (tuple) – optional (2,) indices of state vector corresponding to the FoV space. Default is (0,1)

• spectral_direction_only (bool, optional) – if True, only use spectral direction for splitting. Default is False

• variance_preserving (bool, optional) – if True, preserve variance in splitting. Default is True

class pyest.gm.gm.Covariance

Bases: object

Representation of a covariance matrix

Calculations involving covariance matrices (e.g. data whitening, multivariate normal function evaluation) are often performed more efficiently using a decomposition of the covariance matrix instead of the covariance matrix itself. This class allows the user to construct an object representing a covariance matrix using any of several decompositions and perform calculations using a common interface.

Note

The Covariance class cannot be instantiated directly. Instead, use one of the factory methods (e.g. Covariance.from_diagonal).

Examples

The Covariance class is used by calling one of its factory methods to create a Covariance object, then pass that representation of the Covariance matrix as a shape parameter of a multivariate distribution.

For instance, the multivariate normal distribution can accept an array representing a covariance matrix:

>>> from scipy import stats
>>> import numpy as np
>>> d = [1, 2, 3]
>>> A = np.diag(d)  # a diagonal covariance matrix
>>> x = [4, -2, 5]  # a point of interest
>>> dist = stats.multivariate_normal(mean=[0, 0, 0], cov=A)
>>> dist.pdf(x)
4.9595685102808205e-08

but the calculations are performed in a very generic way that does not take advantage of any special properties of the covariance matrix. Because our covariance matrix is diagonal, we can use Covariance.from_diagonal to create an object representing the covariance matrix, and multivariate_normal can use this to compute the probability density function more efficiently.

>>> cov = stats.Covariance.from_diagonal(d)
>>> dist = stats.multivariate_normal(mean=[0, 0, 0], cov=cov)
>>> dist.pdf(x)
4.9595685102808205e-08

__init__()

static from_diagonal(diagonal)

Return a representation of a covariance matrix from its diagonal.

Parameters:

diagonal (array_like) – The diagonal elements of a diagonal matrix.

Notes

Let the diagonal elements of a diagonal covariance matrix \(D\) be stored in the vector \(d\).

When all elements of \(d\) are strictly positive, whitening of a data point \(x\) is performed by computing \(x \cdot d^{-1/2}\), where the inverse square root can be taken element-wise. \(\log\det{D}\) is calculated as \(-2 \sum(\log{d})\), where the \(\log\) operation is performed element-wise.

This Covariance class supports singular covariance matrices. When computing _log_pdet, non-positive elements of \(d\) are ignored. Whitening is not well defined when the point to be whitened does not lie in the span of the columns of the covariance matrix. The convention taken here is to treat the inverse square root of non-positive elements of \(d\) as zeros.

Examples

Prepare a symmetric positive definite covariance matrix A and a data point x.

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> n = 5
>>> A = np.diag(rng.random(n))
>>> x = rng.random(size=n)

Extract the diagonal from A and create the Covariance object.

>>> d = np.diag(A)
>>> cov = stats.Covariance.from_diagonal(d)

Compare the functionality of the Covariance object against a reference implementations.

>>> res = cov.whiten(x)
>>> ref = np.diag(d**-0.5) @ x
>>> np.allclose(res, ref)
True
>>> res = cov.log_pdet
>>> ref = np.linalg.slogdet(A)[-1]
>>> np.allclose(res, ref)
True

static from_precision(precision, covariance=None)

Return a representation of a covariance from its precision matrix.

Parameters:

• precision (array_like) – The precision matrix; that is, the inverse of a square, symmetric, positive definite covariance matrix.

• covariance (array_like, optional) – The square, symmetric, positive definite covariance matrix. If not provided, this may need to be calculated (e.g. to evaluate the cumulative distribution function of scipy.stats.multivariate_normal) by inverting precision.

Notes

Let the covariance matrix be \(A\), its precision matrix be \(P = A^{-1}\), and \(L\) be the lower Cholesky factor such that \(L L^T = P\). Whitening of a data point \(x\) is performed by computing \(x^T L\). \(\log\det{A}\) is calculated as \(-2tr(\log{L})\), where the \(\log\) operation is performed element-wise.

This Covariance class does not support singular covariance matrices because the precision matrix does not exist for a singular covariance matrix.

Examples

Prepare a symmetric positive definite precision matrix P and a data point x. (If the precision matrix is not already available, consider the other factory methods of the Covariance class.)

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> n = 5
>>> P = rng.random(size=(n, n))
>>> P = P @ P.T  # a precision matrix must be positive definite
>>> x = rng.random(size=n)

Create the Covariance object.

>>> cov = stats.Covariance.from_precision(P)

Compare the functionality of the Covariance object against reference implementations.

>>> res = cov.whiten(x)
>>> ref = x @ np.linalg.cholesky(P)
>>> np.allclose(res, ref)
True
>>> res = cov.log_pdet
>>> ref = -np.linalg.slogdet(P)[-1]
>>> np.allclose(res, ref)
True

static from_cholesky(cholesky)

Representation of a covariance provided via the (lower) Cholesky factor

Parameters:

cholesky (array_like) – The lower triangular Cholesky factor of the covariance matrix.

Notes

Let the covariance matrix be \(A\) and \(L\) be the lower Cholesky factor such that \(L L^T = A\). Whitening of a data point \(x\) is performed by computing \(L^{-1} x\). \(\log\det{A}\) is calculated as \(2tr(\log{L})\), where the \(\log\) operation is performed element-wise.

This Covariance class does not support singular covariance matrices because the Cholesky decomposition does not exist for a singular covariance matrix.

Examples

Prepare a symmetric positive definite covariance matrix A and a data point x.

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> n = 5
>>> A = rng.random(size=(n, n))
>>> A = A @ A.T  # make the covariance symmetric positive definite
>>> x = rng.random(size=n)

Perform the Cholesky decomposition of A and create the Covariance object.

>>> L = np.linalg.cholesky(A)
>>> cov = stats.Covariance.from_cholesky(L)

Compare the functionality of the Covariance object against reference implementation.

>>> from scipy.linalg import solve_triangular
>>> res = cov.whiten(x)
>>> ref = solve_triangular(L, x, lower=True)
>>> np.allclose(res, ref)
True
>>> res = cov.log_pdet
>>> ref = np.linalg.slogdet(A)[-1]
>>> np.allclose(res, ref)
True

static from_eigendecomposition(eigendecomposition)

Representation of a covariance provided via eigendecomposition

Parameters:

eigendecomposition (sequence) – A sequence (nominally a tuple) containing the eigenvalue and eigenvector arrays as computed by scipy.linalg.eigh or numpy.linalg.eigh.

Notes

Let the covariance matrix be \(A\), let \(V\) be matrix of eigenvectors, and let \(W\) be the diagonal matrix of eigenvalues such that V W V^T = A.

When all of the eigenvalues are strictly positive, whitening of a data point \(x\) is performed by computing \(x^T (V W^{-1/2})\), where the inverse square root can be taken element-wise. \(\log\det{A}\) is calculated as \(tr(\log{W})\), where the \(\log\) operation is performed element-wise.

This Covariance class supports singular covariance matrices. When computing _log_pdet, non-positive eigenvalues are ignored. Whitening is not well defined when the point to be whitened does not lie in the span of the columns of the covariance matrix. The convention taken here is to treat the inverse square root of non-positive eigenvalues as zeros.

Examples

Prepare a symmetric positive definite covariance matrix A and a data point x.

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> n = 5
>>> A = rng.random(size=(n, n))
>>> A = A @ A.T  # make the covariance symmetric positive definite
>>> x = rng.random(size=n)

Perform the eigendecomposition of A and create the Covariance object.

>>> w, v = np.linalg.eigh(A)
>>> cov = stats.Covariance.from_eigendecomposition((w, v))

Compare the functionality of the Covariance object against reference implementations.

>>> res = cov.whiten(x)
>>> ref = x @ (v @ np.diag(w**-0.5))
>>> np.allclose(res, ref)
True
>>> res = cov.log_pdet
>>> ref = np.linalg.slogdet(A)[-1]
>>> np.allclose(res, ref)
True

whiten(x)

Perform a whitening transformation on data.

“Whitening” (“white” as in “white noise”, in which each frequency has equal magnitude) transforms a set of random variables into a new set of random variables with unit-diagonal covariance. When a whitening transform is applied to a sample of points distributed according to a multivariate normal distribution with zero mean, the covariance of the transformed sample is approximately the identity matrix.

Parameters:

x (array_like) – An array of points. The last dimension must correspond with the dimensionality of the space, i.e., the number of columns in the covariance matrix.

Returns:

x_ – The transformed array of points.

Return type:

array_like

References

[1]

“Whitening Transformation”. Wikipedia. https://en.wikipedia.org/wiki/Whitening_transformation

[2]

Novak, Lukas, and Miroslav Vorechovsky. “Generalization of coloring linear transformation”. Transactions of VSB 18.2 (2018): 31-35. :doi:`10.31490/tces-2018-0013`

Examples

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> n = 3
>>> A = rng.random(size=(n, n))
>>> cov_array = A @ A.T  # make matrix symmetric positive definite
>>> precision = np.linalg.inv(cov_array)
>>> cov_object = stats.Covariance.from_precision(precision)
>>> x = rng.multivariate_normal(np.zeros(n), cov_array, size=(10000))
>>> x_ = cov_object.whiten(x)
>>> np.cov(x_, rowvar=False)  # near-identity covariance
array([[0.97862122, 0.00893147, 0.02430451],
       [0.00893147, 0.96719062, 0.02201312],
       [0.02430451, 0.02201312, 0.99206881]])

colorize(x)

Perform a colorizing transformation on data.

“Colorizing” (“color” as in “colored noise”, in which different frequencies may have different magnitudes) transforms a set of uncorrelated random variables into a new set of random variables with the desired covariance. When a coloring transform is applied to a sample of points distributed according to a multivariate normal distribution with identity covariance and zero mean, the covariance of the transformed sample is approximately the covariance matrix used in the coloring transform.

Parameters:

x (array_like) – An array of points. The last dimension must correspond with the dimensionality of the space, i.e., the number of columns in the covariance matrix.

Returns:

x_ – The transformed array of points.

Return type:

array_like

References

[1]

“Whitening Transformation”. Wikipedia. https://en.wikipedia.org/wiki/Whitening_transformation

[2]

Novak, Lukas, and Miroslav Vorechovsky. “Generalization of coloring linear transformation”. Transactions of VSB 18.2 (2018): 31-35. :doi:`10.31490/tces-2018-0013`

Examples

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng(1638083107694713882823079058616272161)
>>> n = 3
>>> A = rng.random(size=(n, n))
>>> cov_array = A @ A.T  # make matrix symmetric positive definite
>>> cholesky = np.linalg.cholesky(cov_array)
>>> cov_object = stats.Covariance.from_cholesky(cholesky)
>>> x = rng.multivariate_normal(np.zeros(n), np.eye(n), size=(10000))
>>> x_ = cov_object.colorize(x)
>>> cov_data = np.cov(x_, rowvar=False)
>>> np.allclose(cov_data, cov_array, rtol=3e-2)
True

property log_pdet

Log of the pseudo-determinant of the covariance matrix

property rank

Rank of the covariance matrix

property covariance

Explicit representation of the covariance matrix

property shape

Shape of the covariance array

pyest.gm.gm.fit_1d_uniform(lb, ub, L=100)

approximates 1D uniform distribution with GM

Parameters:

• lb (float) – lower bound of uniform distribution

• ub (float) – upper bound of uniform distribution

• L (int, optional) – number of components to use in approximation. Default is 100

Return type:

GaussianMixture

pyest.gm.gm.gm_fit_1d(xvals, pvals, support=None, stdmax=None, L=None)

pyest.gm.gm.equal_weights(L)

generate L equal weights that sum to 1

Parameters:

L (int) – number of weights

Return type:

ndarray

pyest.gm.gm.eval_gmpdf(x, w, m, P)

evaluates Gaussian mixture at points x

Parameters:

• x (array_like) – (n_samp, vect_length) points at which to evaluate the GM

• w (array_like) – (nC,) Gaussian mixture weights

• m (array_like) – (nC, vect_length) Gaussian mixture means

• P (array_like) – (nC, vect_length, vect_length) Gaussian mixture covariances

Returns:

(n_samp,) Gaussian mixture evaluated at points x

Return type:

ndarray

pyest.gm.gm.eval_gmpdfchol(x, w, m, S)

evaluates Gaussian mixture at points x

Parameters:

• x (array_like) – (n_samp, vect_length) points at which to evaluate the GM

• w (array_like) – (nC,) Gaussian mixture weights

• m (array_like) – (nC, vect_length) Gaussian mixture means

• S (array_like) – (nC, vect_length, vect_length) Gaussian mixture lower-triangular Cholesky factors

Returns:

(n_samp,) Gaussian mixture evaluated at points x

Return type:

ndarray

pyest.gm.gm.eval_mvnpdf(x, m, P, allow_singular=False)

Multivariate normal probability density function.

Parameters:

• x (array_like) – Quantiles, with the last axis of x denoting the components.

• %(_mvn_doc_default_callparams)s

Returns:

pdf – Probability density function evaluated at x

Return type:

ndarray or scalar

Notes

%(_mvn_doc_callparams_note)s

pyest.gm.gm.eval_mvnpdfchol(x, m, S, allow_singular=False)

Multivariate normal probability density function.

Parameters:

• x (array_like) – Quantiles, with the last axis of x denoting the components.

• m (ndarray) – (nx,) mean vector

• S (ndarray) – (nx, nx) lower-triangular Cholesky factor of covariance matrix

Returns:

pdf – Probability density function evaluated at x

Return type:

ndarray or scalar

pyest.gm.gm.v_eval_mvnpdf(x, m, P)

pyest.gm.gm.v_eval_mvnpdfchol(x, m, S)

pyest.gm.gm.integral_gauss_product(m1, P1, m2, P2, allow_singular=False)

compute integral of product of two Gaussians

Parameters:

• m1 (np.ndarray) – mean of Gaussian 1

• P1 (np.ndarray) – covariance of Gaussian 1

• m2 (np.ndarray) – mean of Gaussian 2

• P2 (np.ndarray) – covariance of Gaussian 2

Returns:

integral – integral of the product of the two Gaussians

Return type:

float

pyest.gm.gm.integral_gauss_product_chol(m1, S1, m2, S2)

compute integral of product of two Gaussians

Parameters:

• m1 (np.ndarray) – mean of Gaussian 1

• S1 (np.ndarray) – Cholesky factor of covariance of Gaussian 1

• m2 (np.ndarray) – mean of Gaussian 2

• S2 (np.ndarray) – Cholesky factor of covariance of Gaussian 2

Returns:

integral – integral of the product of the two Gaussians

Return type:

float

pyest.gm.gm.integral_squared_gm(p)

compute integral of squared Gaussian mixture

Parameters:

p (GaussianMixture) – Gaussian mixture

Returns:

integral – integral of the squared Gaussian mixture

Return type:

float

pyest.gm.gm.marginal_2d(m, P, dimensions=[0, 1])

compute 2D marginal of GM

pyest.gm.gm.marginal_nd(m, P, dimensions)

returns nd-mariginal GM over specified dimensions

pyest.gm.gm.comp_bounds(m, P, sigma_mult=3)

find lower and upper sigma bounds for each component

pyest.gm.gm.bounds(m, P, sigma_mult=3)

find bounds of GM

pyest.gm.gm.gm_pdf_2d(w, m, P, dimensions=(0, 1), res=100, xbnd=None, ybnd=None)

evaluate GM pdf in two dimensions

Parameters:

• w (ndarray) – (nC,) weights

• m (ndarray) – (nC, nx) means

• P (ndarray) – (nC, nx, nx) covariances

• Optional

• --------

• dimensions (tuple) – dimensions to evaluate pdf over. Defaults to (0, 1)

• res (int) – resolution of the grid. Defaults to 100

• xbnd (tuple) – x-axis bounds. Determined automatically if not provided

• ybnd (tuple) – y-axis bounds. Determined automatically if not provided

Returns:

• ndarray – (res, res) pdf values

• ndarray – (res, res) x values

• ndarray – (res, res) y values

pyest.gm.gm.sigma_contour(m, P, sig_mul=1, num_pts=100)

compute sigma contours of a 2D Gaussian

Parameters:

• m (ndarray) – (2,) mean vector

• P (ndarray) – (2, 2) covariance matrix

• sig_mul (float, optional) – sigma multiplier factor. For example, sig_mul=2 will return points corresponding to the 2-sigma contour. Defaults to sig_mul=1

• num_pts (int, optional) – number of points to compute. Defaults to num_pts=100

Returns:

(num_pts, 2) array of points along the contour

Return type:

ndarray

pyest.gm.gm.distribute_mean_centers(lb, ub, L)

distribute L mean centers between lb and ub

pyest.gm.gm.optimal_homoscedastic_std(L, width=1)

generates optimal std for an L component GM under homoscedasticity constraint

Split Methods

pyest.gm.split._reflect_weights(w, L)

reflect weights so that the L-element set is symmetric

Parameters:

• w (arraylike) – (ceil(L/2),) array of weights

• L (int) – length of symmetric weight set

Returns:

L-element symmetric weight array

Return type:

ndarray

pyest.gm.split.obj_l2_gauss_split(eps, sig, w, L, lam)

objective function for optimal Gausian split

Parameters:

• L (int) – number of components in split

• lam (float) – optimization parameter that specifies the importance of standard deviation size and overall L2 distance. lam=0 will place zero importance on standard deviation, making objective equivalent to L2 distance

• eps (float) – spacing between means

• sig (float) – standard deviation of components

• w (ndarray) – (ceil(L/2),) weights of candidate Gaussian mixture. To enforce symmetry, only the left-half weights are used

Returns:

value of objective function J

Return type:

float

pyest.gm.split.equally_spaced_centered_means(L, eps)

generate equally spaced means centered at zero

Parameters:

• L (int) – number of components

• eps (float) – spacing between means

Returns:

(L,) array of means

Return type:

ndarray

pyest.gm.split.optimize_gauss_split_variance_preserving(L, lam)

optimize L-wise split of standard univariate Gaussian, preserving variance

Parameters:

• L (int) – number of components to split into

• lam (float) – optimization parameter that specifies the importance of standard deviation size and overall L2 distance. lam=0 will place zero importance on standard deviation, making objective equivalent to L2 distance

Returns:

L-component GM approximation of standard normal

Return type:

GaussianMixture

Notes

This optimization process uses the following constraints to restrict the dimensionality of the search space: 1. Means are evenly spaced 2. The components are homoscedastic 3. The distribution variance is preserved

pyest.gm.split.optimize_gauss_split(L, lam)

optimize L-wise split of standard univariate Gaussian

Parameters:

• L (int) – number of components to split into

• lam (float) – optimization parameter that specifies the importance of standard deviation size and overall L2 distance. lam=0 will place zero importance on standard deviation, making objective equivalent to L2 distance

Returns:

L-component GM approximation of standard normal

Return type:

GaussianMixture

Notes

This optimization process uses the following constraints to restrict the dimensionality of the search space: 1. Means are evenly spaced 2. The components are homoscedastic

pyest.gm.split.split_1d_standard_gaussian(L, lam, variance_preserving=True)

split 1D standard univariate Gaussian into L-component GM

Parameters:

• L (int) – number of components to split into

• lam (float) – optimization parameter that specifies the importance of standard deviation size and overall L2 distance. lam=0 will place zero importance on standard deviation, making objective equivalent to L2 distance

• variance_preserving (bool) – if True, the L2 distance will be minimized while preserving the variance of the original Gaussian. If False, the L2 distance will be minimized without regard to the variance of the original Gaussian

Returns:

L-component GM approximation of standard normal

Return type:

GaussianMixture

pyest.gm.split.split_gaussian(w, m, cov, split_options, cov_type='full', direction=None)

split multivariate Gaussian into smaller Gaussians

Parameters:

• w (float) – component weight

• m (ndarray) – (nx,) component mean

• cov (ndarray) – (nx, nx) component covariance

• split_options (GaussSplitOptions) – splitting options

• cov_type (str, optional) – type of covariance provided. Default is ‘full’

• direction (ndarray, optional) – (nx,) desired direction along which to split. The covariance eigenvector that closest matches this direction is used. If direction is None, the direction of largest variance will be used.

pyest.gm.split.split_spectral_direction(w, m, S, V, D, eig_idx, wt, mt, Pt)

split Gaussian along spectral direction

Parameters:

• w (float) – component weight

• m (ndarray) – (nx,) component mean

• S (ndarray) – (nx, nx) component covariance lower Cholesky factor

• V (ndarray) – (nx, nx) spectral decomposition eigenvectors

• D (ndarray) – (nx,,nx) spectral decomposition eigenvalues diagonal matrix

• eig_idx (int) – index of eigenvector to split along

• wt (ndarray) – (L,) weights of the standard split library

• mt (ndarray) – (L, nx) means of the standard split library

• Pt (ndarray)

Returns:

L-component GM approximation of the input weighted Gaussian

Return type:

GaussianMixture

pyest.gm.split.identify_split_components(p, fovs, split_opts)

Identify what components of a GM should be split for a FoV

Parameters:

• p (GaussianMixture) – Gaussian mixture to split

• fovs (PolygonalFieldOfView or list) – FoV, the geometry of which determines which components to split

• split_opts (GaussSplitOptions) – splitting options

Returns:

• split_mask (ndarray) – boolean array, where true elements denote a component that should be split

• split_dir (ndarray) – (nC, nX) array where each row indicates the best direction to split along. Non-positional state elements are left as zero

pyest.gm.split.identify_split_components_3d_fov(p, fovs, split_opts)

Identify what components of a GM should be split for a set of 3D FoVs

Parameters:

• p (GaussianMixture) – Gaussian mixture to split

• fovs (FieldOfView or list) – FoV, the geometry of which determines which components to split

• split_opts (GaussSplitOptions) – splitting options

Returns:

• split_mask (ndarray) – boolean array, where true elements denote a component that should be split

• split_dir (ndarray) – (nC, nX) array where each row indicates the best direction to split along. Non-positional state elements are left as zero

pyest.gm.split.group_preserving_split_dir(intact_rows, intact_cols, eigvals)

choose direction for FoV splitting to minimize number of downstream components

Parameters:

• intact_cols (ndarray) – (L,) logical, where intact_cols[i]=True indicates that all grid points in the ith column are either entirely included in the FoV or entirely excluded

• intact_rows (ndarray) – (L,) logical, where intact_rows[i]=True indicates that all grid points in the ith row are either entirely included in the FoV or entirely excluded

• eigvals (ndarray) – (2,) positional covariance eigenvalues. In the case that the direction cannot be determined by point inclusion, the direction will be chosen according to the largest positional covariance eigenvalue.

Returns:

index of direction to split along. 0 corresponds to a split along the horizontal axis, and 1 corresponds to a split along the vertical axis

Return type:

int

pyest.gm.split.find_intact_slices(in_mask_tensor, num_pts_in_slice)

find slices that are totally included/excluded

Parameters:

• in_mask_tensor (ndarray) – (L,L,L) boolean mask indicating point-wise inclusion by the FoV, where L is the number of equally spaced points per dimension

• num_pts_in_slice (int) – (L,) number of valid collocation points in each slice

Returns:

• intact_xy_plane (ndarray) – (L,) logical, where intact_xy[i]=True indicates that all grid points in the ith horizontal slice [i,:,:] are either entirely included in the FoV or entirely excluded

• intact_xz_plane (ndarray) – (L,) logical, where intact_xz_plane[j]=True indicates that all grid points in the jth lateral slice [:,j,:] are either entirely included in the FoV or entirely excluded

• intact_yz_plane (ndarray) – (L,) logical, where intact_xz_plane[k]=True indicates that all grid points in the kth frontal slice [:,:,k] are either entirely included in the FoV or entirely excluded

pyest.gm.split.find_intact_rows_cols(in_mask_mat, num_pts_in_slice)

find rows and columns that are totally included/excluded

Parameters:

• in_mask_mat (ndarray) – (L,L) boolean mask indicating point-wise inclusion by the FoV, where L is the number of equally spaced points per dimension

• num_pts_in_slice (int) – (L,1) number of collocation points in each slice

Returns:

• intact_cols (ndarray) – (L,) logical, where intact_cols[i]=True indicates that all grid points in the ith column are either entirely included in the FoV or entirely excluded

• intact_rows (ndarray) – (L,) logical, where intact_rows[i]=True indicates that all grid points in the ith row are either entirely included in the FoV or entirely excluded

pyest.gm.split.split_for_fov(p, fovs, split_opts)

SPLIT_FOR_FOV Splits Gaussians that are close to FoV edges

Assumes that first two dimensions correspond to the same coordinates that the FOV is expressed in

p: GaussianMixture

mixture to be split

fovs: list

FieldOfView, where the boundaries are used to determine where to split the distribution

split_opts: GaussSplitOptions

splitting options

p_split: GaussianMixture

split Gaussian mixture

Written by Keith LeGrand

pyest.gm.split.recursive_split(p, split_opts, identify_split_components, *args)

RECURSIVE_SPLIT Splits Gaussians recursively based on a splitting criterion

p: GaussianMixture

mixture to be split

split_opts: MixtureSplittingOptions

splitting options

indentify_split_components: callable

function to identify components to split, of the form split_mask, split_dir = identify_split_components(p, *args)

p_split: GaussianMixture

split Gaussian mixture

pyest.gm.split.id_variance(p, tol)

identify components and split directions based on variance

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• tol (float) – tolerance for splitting, given as a standard deviation. If the weighted standard deviation of the considered component is greater than tol in any spectral direction, the component is marked for splitting

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{P_x}\) be the initial covariance matrix.

variance measure is given by

\[\max_{|\mathbf{x}\|_{P_x^{-1}} = 1} \|\mathbf{x}\|\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

pyest.gm.split.id_fos(p, jacobian_func, tol)

identify components and split directions based on first-order stretching

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• jacobian_func (callable) – function that returns the Jacobian of the nonlinear function

• tol (float) – tolerance for splitting. If the weighted norm exceeds tol, the component is marked for splitting

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

The first order stretching (FOS) measure is given by

\[\max_{|\mathbf{x}\|_{2} = 1} \|\mathbf{G}\mathbf{x}\|_{2}\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

See also

id_usfos

uncertainty scaled first-order stretching

id_safos

spherical average first-order stretching

pyest.gm.split.id_safos(p, jacobian_func, tol)

identify components and split directions based on spherical average first-order stretching

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• jacobian_func (callable) – function that returns the Jacobian of the measurement model

• tol (float) – tolerance for splitting.

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{\\varphi(u)}\) be the (n-1)-dimensional surface measure on the ellipsoid \(\{u : u^TP_x^{-1}u = 1\}\)

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

\(\mathbf{P_x}\) be the initial covariance matrix.

\(\mathbf{u}\) be .

The spherical-average first order stretching (SAFOS) measure is given by

\[\begin{split}\max_{\|\mathbf{x}\|_{2} = 1}\int\limits_{\mathbf{u}^T\mathbf{P}_x^{-1} \mathbf{u}}(\mathbf{u}^{T}\mathbf{x})^{2}\Vert \mathbf{G}\mathbf{u}\Vert_{2}^{2} \mathrm{d}\\varphi(\mathbf{u})\end{split}\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

pyest.gm.split.id_usfos(p, jacobian_func, tol)

identify components and split directions based on scaled first-order stretching

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• jacobian_func (callable) – function that returns the Jacobian of the measurement model

• tol (float) – tolerance for splitting. If the weighted norm exceeds tol, the component is marked for splitting

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

\(\mathbf{P_x}\) be the initial covariance matrix.

The uncertainty-scaled first order stretching (USFOS) measure is given by

\[\max_{\|\mathbf{x}\|_{P_x^{-1}} = 1} \Vert\mathbf{Gx}\Vert_2\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

pyest.gm.split.id_sos(p, pdt_func, jacobian_func, tol, single_fn=False)

identify components and split directions based on nonlinear stretching

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• pdt_func (callable) – nonlinear function partial derivative tensor of the form pdt_func(x1, x2, …, xn)

• jacobian_func (callable) – nonlinear function Jacobian of the form jacobian_func(x1, x2, …, xn)

• tol (float) – tolerance for splitting, given as a Mahalanobis distance. If e^T Pf^-1 e > tol, the component is marked for splitting, where Pf is the linearly-mapped covariance of the considered component

• single_fn (bool) – interpret pdt_func argument as returning the tuple (x_f, jac, hess) and disregard jacobian function completely

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

The second order stretching (SOS) measure is given by

\[\max_{|\mathbf{x}\|_{2} = 1} \|\mathbf{G}^{(2)}\mathbf{x}^2\|_{2}\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

pyest.gm.split.id_wussos(p, pdt_func, jacobian_func, tol, single_fn=False)

identify components and split directions based on scaled nonlinear stretching

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• pdt_func (callable) – nonlinear function second-order partial derivative tensor of the form pdt_func(x1, x2, …, xn)

• jacobian_func (callable) – nonlinear function Jacobian of the form jacobian_func(x1, x2, …, xn)

• tol (float) –

  tolerance for splitting, given as a Mahalanobis distance. If

  e^T Pf^-1 e > tol, the component is marked for splitting, where Pf is the linearly-mapped covariance of the considered component

• single_fn (bool) – interpret pdt_func argument as returning the tuple (x_f, jac, hess) and disregard jacobian function completely

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{\|x\|_{P_z^{-1}}}\) be the norm of x induced by the final precision matrix.

\(\mathbf{\|x\|_{P_x^{-1}}}\) be the norm of x induced by the initial precision matrix.

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

The whitened uncertainty-scaled second order stretching (WUSSOS) measure is given by

\[\max_{\|\mathbf{x}\|_{P_x^{-1}} = 1}\Vert\mathbf{G}^{(2)}\mathbf{x}^2\Vert_{\mathbf{P}_z^{-1}}\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

pyest.gm.split.id_solc(p, pdt_func, tol)

identify components and split directions based on second-order linearization change

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• pdt_func (callable) – nonlinear function partial derivative tensor of the form pdt_func(x1, x2, …, xn)

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{\|x\|_F}\) be the Frobenius norm of x.

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

The second order linearization change (SOLC) measure is given by

\[\max_{\|\mathbf{x}\|_{2} = 1} \Vert\mathbf{G}^{(2)}\mathbf{x}\Vert_F\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

[2]

K. Tuggle and R. Zanetti, “Automated Splitting Gaussian mixture Nonlinear Measurement Update,” Journal of Guidance, Control, and Dynamics, vol. 41, no. 3, pp. 725–734, 2018.

[3]

K. Tuggle, “Model Selection for Gaussian Mixture Model Filtering and Sensor Scheduling,” Ph.D. dissertation, 2020.

See also

id_ussolc

uncertainty scaled second-order linearization change

id_wussolc

whitened uncertainty scaled second-order linearization change

pyest.gm.split.id_ussolc(p, pdt_func, tol)

identify components and split directions based on uncertainty scaled second-order linearization change

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• pdt_func (callable) – function that returns the second-order partial derivative tensor of the nonlinear function

• tol (float) – tolerance for splitting. If the weighted norm exceeds tol, the component is marked for splitting

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{\|x\|_F}\) be the Frobenius norm of x.

\(\mathbf{\|x\|_{P_x^{-1}}}\) be the norm of x induced by the initial precision matrix.

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

The uncertainty-scaled second order linearization change (USSOLC) measure is given by

\[\max_{\|\mathbf{x}\|_{P_x^{-1}} = 1} \Vert\mathbf{G^{(2)}x}\Vert_F\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

[2]

K. Tuggle and R. Zanetti, “Automated Splitting Gaussian mixture Nonlinear Measurement Update,” Journal of Guidance, Control, and Dynamics, vol. 41, no. 3, pp. 725–734, 2018.

[3]

K. Tuggle, “Model Selection for Gaussian Mixture Model Filtering and Sensor Scheduling,” Ph.D. dissertation, 2020.

See also

id_solc

second-order linearization change

id_wussolc

whitened uncertainty scaled second-order linearization change

pyest.gm.split.id_wussolc(p, pdt_func, jacobian_func, tol, single_fn=False)

identify components and split directions based on output-whitened uncertainty scaled second-order linearization change

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• pdt_func (callable) – function that returns the nonlinear function’s second-order partial derivative tensor of the nonlinear function of the form pdt_func(x1, x2, …, xn)

• jacobian_func (callable) – function that returns the Jacobian of the nonlinear function

• tol (float) – tolerance for splitting. If the weighted norm exceeds tol, the component is marked for splitting

• single_fn (bool) – interpret pdt_func argument as returning the tuple (x_f, jac, pdt) and disregard jacobian function completely

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{\|x\|_F}\) be the Frobenius norm of x.

\(\mathbf{\|x\|_{P_x^{-1}}}\) be the norm of x induced by the initial precision matrix.

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

\(\mathbf{P_x}\) be the initial covariance matrix.

\(\mathbf{P_z}\) be the linear prediction of the covariance matrix.

The whitened uncertainty-scaled second-order linearization change (WUSSOLC) measure is given by

\[\max_{\|\mathbf{x}\|_{P_x^{-1}} = 1}\Vert \mathbf{P}_z^{-1/2}(\mathbf{G}^{(2)}\mathbf{x})\mathbf{P}_x^{1/2}\Vert_{F}^2\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

See also

id_solc

second-order linearization change

id_ussolc

uncertainty scaled second-order linearization change

pyest.gm.split.id_sasos(p, pdt_func, tol)

identify components and split directions based on scaled nonlinear stretching

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• pdt_func (callable) – nonlinear function second-order partial derivative tensor of the form pdt_func(x1, x2, …, xn)

• tol (float) –

  tolerance for splitting, given as a Mahalanobis distance. If

  e^T Pf^-1 e > tol, the component is marked for splitting, where Pf is the linearly-mapped covariance of the considered component

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{u}\) be the .

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

\(\mathbf{P_x}\) be the initial covariance matrix.

\(\mathbf{\\varphi(u)}\) be the (n-1)-dimensional surface measure on the ellipsoid \(\{u : u^TP_x^{-1}u = 1\}\)

The spherical-average second order stretching (SASOS) measure is given by

\[\begin{split}\max_{\|\mathbf{x}\|_{2} = 1}\int\limits_{\mathbf{u}^T\mathbf{P}_x^{-1}\mathbf{u}}(\mathbf{u}^{T}\mathbf{x})^{2}\Vert\mathbf{G^{(2)}}\mathbf{u^{2}}\Vert_{2}^{2}\mathrm{d}\\varphi(\mathbf{u})\end{split}\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

pyest.gm.split.id_wsasos(p, pdt_func, jacobian_func, tol, single_fn=False)

identify components and split directions based on scaled nonlinear stretching

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• pdt_func (callable) – nonlinear function second-order partial derivative tensor of the form pdt_func(x1, x2, …, xn)

• jacobian_func (callable) – nonlinear function Jacobian of the form jacobian_func(x1, x2, …, xn)

• tol (float) –

  tolerance for splitting, given as a Mahalanobis distance. If

  e^T Pf^-1 e > tol, the component is marked for splitting, where Pf is the linearly-mapped covariance of the considered component

• single_fn (bool) – interpret pdt_func argument as returning the tuple (x_f, jac, pdt) and disregard jacobian function completely

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{\|x\|_{P_x^{-1}}}\) be the norm of x induced by the initial precision matrix.

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

\(\mathbf{P_x}\) be the initial covariance matrix.

\(\mathbf{\\varphi(u)}\) be the (n-1)-dimensional surface measure on the ellipsoid \(\{u : u^TP_x^{-1}u = 1\}\)

\(\mathbf{u}\) be the .

The whitened spherical-average second order stretching (WSASOS) measure is given by

\[\begin{split}\max_{\|\mathbf{x}\|_{2} = 1}\int\limits_{\mathbf{u}^T\mathbf{P}_x^{-1} \mathbf{u}}(\mathbf{u}^{T}\mathbf{x})^{2}\Vert \mathbf{G}^{(2)}\mathbf{u}^2\Vert_{\mathbf{P}_x^{-1}}^{2} \mathrm{d} \\varphi(\mathbf{u})\end{split}\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

pyest.gm.split.id_alodt(p, g, sigma_pt_opts, tol)

identify components and split directions based on sigma point curvature

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• g (callable) – nonlinear function through which to propagate the sigma points

• sigma_pt_opts (SigmaPointOptions) – options for sigma point generation

• tol (float) – tolerance for splitting based on the deviation of the sigma points from a linear fit

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

This method is referred to by the original authors as the “Adaptive Level of Detail” method [1].

References

[1] F. Faubel and D. Klakow, “Further improvement of the adaptive level of

detail transform: Splitting in direction of the nonlinearity,” in 2010 18th European Signal Processing Conference, Aug. 2010, pp. 850–854.

pyest.gm.split.id_sadl(p, jacobian_func, g, sigma_pt_opts, tol)

identify components and split directions based on the difference in deterministic and statistical linearization

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• jacobian_func (callable) – function that returns the Jacobian of the nonlinear function

• g (callable) – nonlinear function through which to propagate the sigma points

• sigma_pt_opts (SigmaPointOptions) – options for sigma point generation

• tol (float) – tolerance for splitting based on the difference between the deterministic and statistical linearizations

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

\(\mathbf{G}^{(SL)}\) be the statistical linearization of G.

The statistical and deterministic linearization (SADL) measure is given by

\[\max_{\|\mathbf{x}\|_{2} = 1} \Vert\mathbf({G}^{(SL)} - G)\mathbf{x}\Vert_2\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

pyest.gm.split.id_wussadl(p, jacobian_func, g, sigma_pt_opts, tol, deterministic_whitening=True)

identify components and split directions based on output-whitened uncertainty scaled statistical and deterministic linearization difference

Parameters:

• p (GaussianMixture) – input Gaussian mixture to be considered for splitting

• jacobian_func (callable) – function that returns the Jacobian of the nonlinear function

• g (callable) – nonlinear function through which to propagate the sigma points

• sigma_pt_opts (SigmaPointOptions) – options for sigma point generation

• tol (float) – tolerance for splitting. If the weighted norm exceeds tol, the component is marked for splitting

Returns:

• split_mask (np.ndarray) – (nC,) boolean array indicating which components are marked for splitting

• split_dir (np.ndarray) – (nC, nX) array of split directions for each component

Notes

Let

\(\mathbf{G}\) be the Jacobian of the nonlinear function.

\(\mathbf{G}^{(SL)}\) be the statistical linearization of G.

\(\mathbf{P_z}\) be the linear prediction of the covariance matrix.

\(\mathbf{\|x\|_{P_x^{-1}}}\) be the norm of x induced by the initial precision matrix.

The whitened uncertainty-scaled statistical and deterministic linearization (WUSSADL) measure is given by

\[\max_{\|\mathbf{x}\|_{P_x^{-1}} = 1}\Vert \mathbf{P}_z^{-1/2}(\mathbf{G}^{(SL)}-\mathbf{G})\mathbf{x}\Vert_{2}\]

References

[1]

Jackson Kulik and Keith A. LeGrand, “Nonlinearity and Uncertainty Informed Moment-Matching Gaussian Mixture Splitting,” https://arxiv.org/abs/2412.00343, 2024

Filters

Linear Algebra

pyest.linalg.is_pos_def(x)

check that x is positive definite

pyest.linalg.is_square_matrix(A)

returns true if A is scalar or 2D square matrix

pyest.linalg.is_valid_chol(x)

check the cholesky factor is valid

pyest.linalg.choldowndate(A, x)

downdate Cholesky factorization of R with x

Parameters:

• A (ndarray) – (n,n) upper-triangular Cholesky factorization of covariance matrix

• x (ndarray) – (n,) vector to downdate R

Returns:

(n,n) downdated upper-triangular Cholesky factorization of covariance matrix

Return type:

ndarray

Raises:

ValueError – if the downdated Cholesky factorization is not valid

pyest.linalg.make_chol_diag_positive(S)

make the diagonal of the Cholesky factorization positive

pyest.linalg.triangularize(A, upper=False)

triangularize a square root factor A

Required

Anumpy.ndarray

square root factor of a symmetric positive definite matrix

returns:

lower triangular matrix

rtype:

numpy.ndarray

pyest.linalg.cholesky_from_sqrt_precision(U)

Compute the Cholesky factor of the covariance matrix from the square root of the precision matrix.

Required

Unumpy.ndarray

square root of the precision matrix: UU.T = inv(P)

returns:

Cholesky factor of the covariance matrix: LL.T = P

rtype:

numpy.ndarray

Metrics

pyest.metrics.l2_dist(p1, p2)

compute L2 distance between GMs p1 and p2

Parameters:

• p1 (GaussianMixture) – first Gaussian mixture

• p2 (GaussianMixture) – second Gaussian mixture

Returns:

L2 distance between the two input GMs

Return type:

float

pyest.metrics.max_covariance_ratio(S, S_ref)

compute the maximum covariance ratio between two distributions

Required:

Snp.ndarray

covariance matrix lower-triangular Cholesky square-root factor of the test distribution

S_refnp.ndarray

covariance matrix lower-triangular Cholesky square-root factor of the reference distribution

Returns:

float

maximum covariance ratio

pyest.metrics.madem(m, S, m_ref)

Mahalanobis distance of the error of the mean

Required:

mnp.ndarray

mean of the test distribution

Snp.ndarray

covariance matrix lower-triangular Cholesky square-root factor of the test distribution

m_refnp.ndarray

mean of the reference distribution

Returns:

float

Mahalanobis distance of the error of the mean (MaDEM)

pyest.metrics.integral_squared_error_2d(p1, p2, a, b, c, d, epsabs=0.0149, epsrel=0.0149)

compute integral squared error between two 2D densities

Parameters:

• p1 (callable) – first density, p1([x,y])

• p2 (callable) – second density, p2([x,y])

• a (float) – The limits of integration in x: a<b

• b (float) – The limits of integration in x: a<b

• c (float) – The limits of integration in y: c<d

• d (float) – The limits of integration in y: c<d

• epsabs (float, optional) – absolute error tolerance for numerical integration

• epsrel (float, optional) – relative error tolerance for numerical integration

Returns:

• ise (foat) – kld

• int_error – numerical integration estimated error

See also

normalized_integral_squared_error_2d

compute normalized integral squared error between two 2D densities

l2_dist

compute L2 distance (ISE) between two Gaussian mixtures

Notes

This function is intended for use with generic callable densities and makes no assumptions about the form of the densities. If both p1 and p2 are Gaussian mixtures, use l2_dist instead, which is exact and more efficient.

pyest.metrics.normalized_integral_squared_error_2d(p1, p2, a, b, c, d, epsabs=0.0149, epsrel=0.0149)

compute normalized integral squared error between two 2D, densities

Parameters:

• p1 (callable) – first density

• p2 (callable) – second density

• a (float) – The limits of integration in x: a<b

• b (float) – The limits of integration in x: a<b

• c (float) – The limits of integration in y: c<d

• d (float) – The limits of integration in y: c<d

• epsabs (float, optional) – absolute error tolerance for numerical integration

• epsrel (float, optional) – relative error tolerance for numerical integration

Returns:

• nise (float) – normalized integral squared error

• ise – integral squared error

• err – numerical integration estimated error in ISE computation

See also

integral_squared_error_2d

compute integral squared error between two 2D densities

l2_dist

compute L2 distance (ISE) between two Gaussian mixtures

Notes

This function is intended for use with generic callable densities and makes no assumptions about the form of the densities. If both p1 and p2 are Gaussian mixtures, use metrics.l2_dist and gm.integral_squared_gm instead for the numerator and denominator terms separately, which is exact and more efficient.

Sensors