Multivariate Normal MLE¶

Maximum likelihood estimation for multivariate normal distributions with missing data. Direct BFGS and EM algorithms. Little’s MCAR test. Missing data pattern analysis.

Multivariate normal maximum likelihood estimation with missing data.

Public API:: mlest(data, …) -> MVNSolution

pystatistics.mvnmle.mlest(data_or_design, *, method='direct', backend=None, solver=None, tol=None, max_iter=None, regularize=True, force=False, collinearity_tol=None, verbose=False)[source]¶

Maximum likelihood estimation for multivariate normal with missing data.

Accepts EITHER:

An MVNDesign object
Raw data array or DataFrame (convenience)

Parameters:

data_or_design (array-like or MVNDesign) – Data matrix with NaN for missing values, or MVNDesign object.
method (str) –
Estimation method: - ‘direct’ (default): gradient-based optimization on the

log-likelihood. The parameterization depends on the backend (see backend): the default CPU path uses a forward-Cholesky factorization; backend='cpu-reference' uses the R-exact inverse Cholesky parameterization.
- ’em’: Expectation-Maximization algorithm. Typically slower to converge but guaranteed monotone likelihood increase.
- ’monotone’: Closed-form MLE for monotone missingness patterns (Anderson 1957). Raises ValidationError if the data are not monotone — users should check with pystatistics.mvnmle.is_monotone() first, or use EM/direct for general patterns. When applicable, this is orders of magnitude faster than iterative algorithms.
backend (str or None) –
Backend selection. Default None → ‘cpu’. - ‘cpu’ (and the default): fast PyTorch forward-Cholesky FP64 path

when PyTorch is installed; otherwise falls back (with a warning) to the numpy inverse-Cholesky reference. Both match R; the PyTorch path is substantially faster.
- ’cpu-reference’: force the numpy inverse-Cholesky reference. This is the R-exact validation anchor and the only direct path that needs no PyTorch. Valid only with method='direct'.
- ’gpu’: require a GPU (CUDA or MPS), float32; raises if none available.
- ’gpu_fp64’: require CUDA, float64 (raises on MPS, which has no float64) — CPU-matching precision on the GPU.
- ’auto’: prefer CUDA (float32) when present, else the fast CPU path.
solver (str or None) – Numerical optimizer for the direct method. If None, auto-selected by backend. Ignored for EM.
tol (float or None) – Convergence tolerance. If None, uses algorithm-appropriate default: direct = 1e-5 (gradient tolerance), em = 1e-4 (parameter change).
max_iter (int or None) – Maximum iterations. If None, uses algorithm-appropriate default: direct = 100, em = 1000.
force (bool) – When False (default), a rank-deficient fit — caused by (near-)collinear variables, for which no interior maximum-likelihood estimate exists — raises SingularMatrixError instead of returning a meaningless result. When True, the degenerate result is returned anyway with converged=False and a warning attached.
collinearity_tol (float or None) – Full-rank detection threshold on the fitted correlation matrix’s minimum eigenvalue. If None, uses the calibrated default (1e-5). Smaller values make the collinearity check more permissive.
verbose (bool) – Print progress information.
regularize (bool)

Return type:

MVNSolution

Examples

>>> from pystatistics.mvnmle import mlest, datasets
>>> result = mlest(datasets.apple)
>>> result_em = mlest(datasets.apple, method='em')
>>> print(result.muhat)
>>> print(result.loglik)

class pystatistics.mvnmle.MVNDesign(_data, _n, _p)[source]¶

Bases: object

Design for multivariate normal MLE with missing data.

Wraps a data matrix (n observations x p variables) that may contain NaN values representing missing data. Immutable after construction.

Construction:: MVNDesign.from_array(data) MVNDesign.from_datasource(ds, columns=[‘a’, ‘b’, ‘c’])

Parameters:

_data (ndarray[tuple[Any, ...], dtype[floating[Any]]])
_n (int)
_p (int)

classmethod from_array(data)[source]¶

Build MVNDesign from array-like data.

Parameters:: data (array-like) – 2D data matrix. Can be numpy array, pandas DataFrame, or any array-like with .values attribute.
Return type:: MVNDesign

classmethod from_datasource(source, *, columns=None)[source]¶

Build MVNDesign from a DataSource.

Parameters:

source (DataSource) – Data source providing columns
columns (list of str, optional) – Column names to include. If None, uses all columns.

Return type:

MVNDesign

property data: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶: Data matrix (n x p), may contain NaN.

property n: int¶: Number of observations.

property p: int¶: Number of variables.

property n_missing: int¶: Total number of missing values.

property missing_rate: float¶: Overall missing rate (0.0 to 1.0).

property has_missing: bool¶: Whether data has any missing values.

class pystatistics.mvnmle.MVNSolution(_result, _design)[source]¶

Bases: SolutionReprMixin

User-facing MVN MLE results.

Wraps the backend Result and provides convenient accessors for all MVN estimation outputs.

Parameters:

_result (Result[MVNParams])
_design (MVNDesign)

property muhat: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶: Estimated mean vector.

property sigmahat: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶: Estimated covariance matrix.

property loglik: float¶: Log-likelihood at the estimated parameters.

property converged: bool¶: Whether the optimization converged.

property n_iter: int¶: Number of iterations.

property gradient_norm: float | None¶: Final gradient norm, if available.

property correlation_matrix: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶: Correlation matrix derived from estimated covariance.

property standard_deviations: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶: Standard deviations from estimated covariance.

property aic: float¶: Akaike Information Criterion.

property bic: float¶: Bayesian Information Criterion.

property info: dict[str, Any]¶: Backend metadata.

property timing: dict[str, float] | None¶: Execution timing breakdown.

property backend_name: str¶: Name of the backend that produced this result.

property warnings: tuple[str, ...]¶: Non-fatal warnings from computation.

summary()[source]¶

Generate summary output.

Return type:: str

to_dict()[source]¶

Convert to dictionary for serialization.

Return type:: dict[str, Any]

class pystatistics.mvnmle.MVNParams(muhat, sigmahat, loglik, n_iter, converged, gradient_norm=None)[source]¶

Bases: object

Parameter payload for MVN MLE.

Immutable data computed by backends.

Parameters:

muhat (ndarray[tuple[Any, ...], dtype[floating[Any]]])
sigmahat (ndarray[tuple[Any, ...], dtype[floating[Any]]])
loglik (float)
n_iter (int)
converged (bool)
gradient_norm (float | None)

muhat: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶

sigmahat: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶

loglik: float¶

n_iter: int¶

converged: bool¶

gradient_norm: float | None = None¶

pystatistics.mvnmle.analyze_patterns(data)[source]¶

Analyze missingness patterns in the data.

Parameters:: data (array-like) – Data matrix with missing values as np.nan. Can be NumPy array or pandas DataFrame.
Returns:: Patterns sorted by frequency (most common first).
Return type:: List[PatternInfo]

pystatistics.mvnmle.pattern_summary(patterns, data_shape=None)[source]¶

Generate summary statistics for missingness patterns.

Parameters:

patterns (List[PatternInfo]) – Pattern information from analyze_patterns()
data_shape (Optional[Tuple[int, int]]) – Original data shape (n_obs, n_vars)

Return type:

PatternSummary

class pystatistics.mvnmle.PatternInfo(pattern_id, observed_indices, missing_indices, n_cases, data, pattern_vector)[source]¶

Bases: object

Information about a single missingness pattern.

Parameters:

pattern_id (int)
observed_indices (ndarray)
missing_indices (ndarray)
n_cases (int)
data (ndarray)
pattern_vector (ndarray)

pattern_id: int¶

observed_indices: ndarray¶

missing_indices: ndarray¶

n_cases: int¶

data: ndarray¶

pattern_vector: ndarray¶

property n_observed: int¶

property n_missing: int¶

property percent_cases: float¶

class pystatistics.mvnmle.PatternSummary(n_patterns, total_cases, overall_missing_rate, most_common_pattern, complete_cases, complete_cases_percent, variable_missing_rates)[source]¶

Bases: object

Summary statistics for all missingness patterns in a dataset.

Parameters:

n_patterns (int)
total_cases (int)
overall_missing_rate (float)
most_common_pattern (PatternInfo)
complete_cases (int)
complete_cases_percent (float)
variable_missing_rates (Dict[int, float])

n_patterns: int¶

total_cases: int¶

overall_missing_rate: float¶

most_common_pattern: PatternInfo¶

complete_cases: int¶

complete_cases_percent: float¶

variable_missing_rates: Dict[int, float]¶

pystatistics.mvnmle.little_mcar_test(data, alpha=0.05, backend=None, method='em', regularize=True, condition_threshold=1000000000000.0, drop_all_missing_rows=True, verbose=False)[source]¶

Little’s test for Missing Completely at Random (MCAR).

Parameters:

data (array-like, shape (n_observations, n_variables)) – Data matrix with missing values as np.nan.
alpha (float, default=0.05) – Significance level.
backend (str or None, default None) – Backend for the ML estimation step (the dominant cost). Default None → ‘cpu’ (R-reference path). Explicit: ‘cpu’, ‘gpu’, or ‘auto’ to prefer GPU when available. The per-pattern test-statistic accumulation runs on CPU regardless — it’s O(P × v³) for tiny v and is not the bottleneck.
method (str, default 'em') – ML method forwarded to mlest: ‘em’ (EM, default) or ‘direct’ (BFGS). EM is the default because BFGS convergence scales poorly with the number of missingness patterns — on realistic tabular data (e.g. 13 vars × 107 patterns) BFGS can take 400+ seconds while EM finishes in under a second. Little’s statistic depends only on the final ML estimates, not on which method produced them.
regularize (bool, default True) – When an observed-variable covariance submatrix is ill-conditioned (common on real tabular data — iris, wine, breast_cancer all trip the default threshold), fall back to the Moore-Penrose pseudo-inverse with a warning rather than raising. Matches R’s BaylorEdPsych::LittleMCAR and misty::na.test defaults. Pass False for strict fail-fast behaviour.
condition_threshold (float, default 1e12) – Maximum condition number before a per-pattern covariance submatrix is considered ill-conditioned.
drop_all_missing_rows (bool, default True) – Drop rows with no observed values before fitting. Such rows contribute nothing to either the MLE or the chi-square statistic; at realistic missingness rates on low-dimensional data they show up routinely and shouldn’t block the test. A warning reports how many rows were dropped.
verbose (bool, default False) – Print detailed progress.

Return type:

MCARTestSolution

class pystatistics.mvnmle.MCARTestSolution(statistic, df, p_value, rejected, alpha, patterns, n_patterns, n_patterns_used, ml_mean, ml_cov, convergence_warnings, method='Little (MLE plug-in)')[source]¶

Bases: SolutionReprMixin

Result of an MCAR test.

Canonical result type for little_mcar_test (Little 1988, MLE plug-in). The dataclass is also reused by downstream packages that ship their own MCAR-test variants (e.g. Lacuna’s MoM plug-in at lacuna.analysis.mcar.mom_mcar_test) — the method field identifies which estimator produced the result so a caller can disambiguate without tracking the calling function.

Notes

ml_mean / ml_cov carry the plug-in estimates used in the chi-square statistic.

Parameters:

statistic (float)
df (int)
p_value (float)
rejected (bool)
alpha (float)
patterns (List[PatternInfo])
n_patterns (int)
n_patterns_used (int)
ml_mean (ndarray)
ml_cov (ndarray)
convergence_warnings (List[str])
method (str)

statistic: float¶

df: int¶

p_value: float¶

rejected: bool¶

alpha: float¶

patterns: List[PatternInfo]¶

n_patterns: int¶

n_patterns_used: int¶

ml_mean: ndarray¶

ml_cov: ndarray¶

convergence_warnings: List[str]¶

method: str = 'Little (MLE plug-in)'¶

summary()[source]¶

Generate human-readable summary of test results.

Return type:: str

pystatistics.mvnmle.is_monotone(data)[source]¶

Return True iff the missingness pattern is monotone.

Parameters:: data (array-like, shape (n, v)) – Data matrix with NaN for missing values.
Return type:: bool

pystatistics.mvnmle.monotone_permutation(data)[source]¶

Find a variable ordering that makes the missingness pattern monotone, or return None if no such ordering exists.

Algorithm: order columns by number of missing values ascending. The pattern is monotone under that order iff each column’s set of missing-row indices is a subset of the next column’s set. Equivalent to “at each prefix, the set of rows with at least one missing among the prefix columns is the set of rows where the last prefix column is missing.”

Parameters:: data (array-like, shape (n, v))
Returns:: order – Permutation that yields monotone patterns, or None.
Return type:: ndarray of int64, shape (v,) or None

pystatistics.mvnmle.mlest_monotone_closed_form(data)[source]¶

Closed-form MVN MLE for monotone missingness.

Works by ordering variables as per monotone_permutation so that missingness is a nested sequence, then performing a chain of OLS regressions. Specifically, with variables ordered $X_1, X_2, ldots, X_v$ such that every observation with $X_k$ observed also has $X_1, ldots, X_{k-1}$ observed:

$hatmu_1, hatsigma_{11}$ are the sample mean and variance of $X_1$ over all observations (it is never missing under the monotone order).

For $k = 2, ldots, v$: using observations where $X_k$ is observed, fit OLS regression $X_k = beta_{k,0} + sum_{j<k} beta_{k,j} X_j + epsilon_k$ with $mathrm{Var}(epsilon_k) = rho_k^2$. Back- substitute into the full $(mu, Sigma)$ via

$hatmu_k = beta_{k,0} + beta_{k,1:}^top hatmu_{1:k-1}$ $hatSigma_{k, 1:k-1} = beta_{k,1:}^top hatSigma_{1:k-1, 1:k-1}$ $hatSigma_{k, k} = rho_k^2 + beta_{k,1:}^top hatSigma_{1:k-1, 1:k-1} beta_{k,1:}$

Undo the variable permutation to return parameters in the user’s original variable order.

Matches Anderson (1957). Runs in $O(v^3 n)$ deterministic work; no iteration.

Parameters:

data (array-like, shape (n, v))

Returns:

mu (ndarray (v,))
sigma (ndarray (v, v))
n_effective (int) – Number of rows that participated in at least one regression (all rows that aren’t all-missing).

Raises:

ValidationError – If the data are not monotone, or if a regression step has fewer observations than predictors (identifiability failure).

Return type:

tuple[ndarray[tuple[Any, …], dtype[_ScalarT]], ndarray[tuple[Any, …], dtype[_ScalarT]], int]