pmcprg.pmc¶

Pairwise Markov Chain models (HMC-IN, HMC-IN2, HMC-DN, PMC-IN, PMC) loaded from TOML, with simulation, MPM classification, exact handling of missing observations, and unsupervised estimation (ICE, SEM, GICE). See the README for the TOML model format, the CLI and the GUI.

Model¶

PMCModel ¶

PMCModel(path: str | Path)

Container for a PMC/HMC model loaded from a TOML configuration file.

Parameters:

Name	Type	Description	Default
`path`	`str or Path`	Path to the `.toml` configuration file.	required

Attributes (read-only)

variant : Variant K : int Number of hidden states. N_default : int Default sequence length for simulation. name : str Human-readable label from [model]. margin_structure : str "state" (K margins f_i) or "pair" (K² margins f_ij) — see the module docstring. missingness : None (ignorable, the default) or the non-ignorable mechanism of the [missingness] table (:mod:pmcprg.pmc.missingness).

prior_p `property` ¶

prior_p: ndarray

K×K joint distribution p[i,j] = P(X_n=i, X_{n+1}=j).

transition_A `property` ¶

transition_A: ndarray

K×K row-stochastic matrix A[i,j] = P(X_{n+1}=j | X_n=i).

stationary_pi `property` ¶

stationary_pi: ndarray

Marginal distribution π[i] = P(X_n=i).

path `property` ¶

path: Path

Path the model was loaded from (or 'memory' if built via from_dict).

raw `property` ¶

raw: dict

Deep copy of the raw TOML dict — safe to mutate.

missingness `property` ¶

missingness

Missingness mechanism: None (ignorable) or a mechanism object.

None — the default, no [missingness] table or mechanism = "ignorable": missing rows are integrated out of the observed-data likelihood p(y_obs). Otherwise a :class:~pmcprg.pmc.missingness.StateMissingness ("state") or :class:~pmcprg.pmc.missingness.StateMarkovMissingness ("state-markov"), frozen: the mask m is then evidence on the states and every inference function works with p(y_obs, m) (:mod:pmcprg.pmc.missingness, :mod:pmcprg.pmc.gaps). Use :meth:with_missingness for a model with another mechanism.

margin_structure `property` ¶

margin_structure: str

"state" (K margins f_i) or "pair" (K² margins f_ij).

State margins are exactly the SR-PMCs whose hidden chain X is Markov (Proposition of DerrodePieczynski_CSDA2013 §2.1); pair margins give the general PMC of DerrodePieczynski_CSDA2013 Eqs. 12–14 (PMC and PMC-IN variants only). See the module docstring.

from_dict `classmethod` ¶

from_dict(raw: dict, path: str | Path = 'memory') -> PMCModel

Build a PMCModel directly from a dict (useful for unit tests).

The dict is deep-copied so subsequent mutations of either raw or model._raw do not affect the other.

with_missingness ¶

with_missingness(missingness) -> PMCModel

A copy of this model with another missingness mechanism.

missingness is None (ignorable), a mechanism object or a [missingness] table (dict), validated against K. The copy is built like the estimators' iterates — from_dict on a deep copy of raw with the [missingness] table replaced (or removed) — so every other parameter is unchanged and self is not modified.

ice_config ¶

ice_config() -> dict

Return the [ice] TOML section as a dict (empty if absent).

sem_config ¶

sem_config() -> dict

Return the [sem] TOML section as a dict (empty if absent).

SEM (:func:pmcprg.pmc.sem.sem) merges this on top of [ice] so a model file can override shared estimator keys (or set SEM-specific ones like sem_seed / max_iter) in a dedicated [sem] block while still inheriting the common [ice] settings.

margin_blocks ¶

margin_blocks() -> list[dict]

Return a deep copy of the [[margins]] blocks, in a stable order.

margin_structure == "state" — K blocks keyed by i ({"i", "dist", "params"[, "candidates"]}), in K-format even when the source listed K² blocks collapsed by [model].margin_structure = "state";
margin_structure == "pair" — K² blocks keyed by i and j ({"i", "j", "dist", "params"[, "candidates"]}), sorted by (i, j).

copula_blocks ¶

copula_blocks() -> list[dict]

Return a deep copy of the [[copulas]] blocks.

margin ¶

margin(i: int, j: int | None = None) -> _MarginDist

Return the margin f_ij of the observation y_n given (x_n, x_{n+1}) = (i, j).

State margins (margin_structure == "state"): the density is f_i whatever the other state, so j is ignored and may be omitted — margin(i, j) is margin(i).
Pair margins (margin_structure == "pair"): j is required. margin(i, j) is f_ij, the left margin of the pair (i, j) in DerrodePieczynski_CSDA2013 Eq. 12; the right margin of that pair (the density of y_{n+1}) is margin(j, i) = f_ji. margin(i) raises ValueError: no density is attached to a state alone.

pdf ¶

pdf(i: int, j: int | None, y: float) -> float

f_ij(y) — PDF of margin(i, j) at y (f_i(y) for state margins).

cdf ¶

cdf(i: int, j: int | None, y: float) -> float

F_ij(y) — CDF of margin(i, j) at y (F_i(y) for state margins).

ppf ¶

ppf(i: int, j: int | None, q: float) -> float

F_ij^{-1}(q) — quantile of margin(i, j) (F_i^{-1} for state margins).

copula ¶

copula(i: int, j: int) -> CopulaVirt

Return the copula C_{ij} for the pair (i, j).

weight ¶

weight(i: int, j: int, y_n: float, y_n1: float) -> float

Un-normalised weight of the pair (x_n, x_{n+1}) = (i, j) given (y_n, y_{n+1}).

The margin of the left observation y_n is margin(i, j) = f_ij and that of the right observation y_{n+1} is margin(j, i) = f_ji — the index inversion of DerrodePieczynski_CSDA2013 Eq. 12. With state margins (margin_structure == "state") both reduce to the per-state densities f_i and f_j.

Formulas by variant

HMC-IN : A[i,j] · f_j(y_{n+1}) HMC-IN2 : A[i,j] · f_j(y_{n+1}) HMC-DN : A[i,j] · f_j(y_{n+1}) · c_ij(F_i(y_n), F_j(y_{n+1})) PMC-IN : p[i,j] · f_ij(y_n) · f_ji(y_{n+1}) PMC : p[i,j] · f_ij(y_n) · f_ji(y_{n+1}) · c_ij(F_ij(y_n), F_ji(y_{n+1}))

(HMC-* variants have state margins, so f_ji = f_j there.)

Relation to forward–backward

For HMC- the weight is the transition p(z_{n+1} | z_n) itself: precompute_weights(model, Y)[0][n, i, j] = weight(i, j, Y[n], Y[n+1]). For PMC- it is the summand of the pair density of DerrodePieczynski_CSDA2013 Eq. 12, p(x_n = i, x_{n+1} = j, y_n, y_{n+1}); the transition used by forward–backward divides it by p(x_n = i, y_n) = Σ_k p[i,k] f_ik(y_n) (DerrodePieczynski_CSDA2013 Eqs. 13–14):

W[n, i, j] = weight(i, j, y_n, y_{n+1}) / Σ_k p[i,k] · f_ik(y_n).

(precompute_weights clips the CDF values to [EPS, 1 − EPS] before the copula density; weight does not.)

Why the margins may depend on the pair

Reversibility (DerrodePieczynski_CSDA2013 §2.2) only gives p(y_1 | x_1, x_2) = p(y_2 | x_2, x_1), whence f_ji for the right observation. It does not make f_ij independent of j: by the Proposition of DerrodePieczynski_CSDA2013 §2.1, f_ij = f_i holds exactly when the hidden chain X is Markov, in which case the PMC is an SR HMC-DN (the transition p(x_{n+1} = j | x_n = i, y_n) ∝ p_ij f_i(y_n) no longer depends on y_n).

save ¶

save(path: str | Path | None = None) -> Path

Write the model back to a TOML file. Returns the path written.

Variant ¶

Bases: Enum

The five PMC/HMC model variants.

Properties

uses_copula : True if Y_n depends on Y_{n-1} via a copula (HMC-DN, PMC). has_markov_prior : True if the prior is a transition matrix A (HMC-). False if the prior is a joint p[i,j] (PMC-). allows_pair_margins : True if the variant admits pair-indexed margins f_ij (PMC, PMC-IN). per_class_margin : True if the variant forces state-indexed margins f_i (HMC-*); the negation of allows_pair_margins.

allows_pair_margins `property` ¶

allows_pair_margins: bool

True for PMC and PMC-IN, whose margins may be pair-indexed (f_ij).

For an HMC- variant X is a Markov chain, and by the Proposition of DerrodePieczynski_CSDA2013 §2.1 an SR-PMC has a Markov X iff p(y_2 | x_1, x_2) = p(y_2 | x_2): its margins are state-indexed. Whether a given model uses* pair margins is :attr:PMCModel.margin_structure.

per_class_margin `property` ¶

per_class_margin: bool

True if the variant forces state-indexed margins f_i (HMC-*).

For PMC and PMC-IN the margins are state-indexed or pair-indexed depending on the model (:attr:PMCModel.margin_structure), so this is False. (Versions 0.5.0–0.8.x returned True for every variant, wrongly attributing f_ij = f_i to reversibility; the Proposition of DerrodePieczynski_CSDA2013 §2.1 ties it to X being Markov.)

Simulation¶

simulate ¶

simulate.py — Sequence generator for all 5 PMC/HMC variants.

Public API

simulate(model, N=None, seed=None) -> (X, Y) Generate a synthetic sequence (X_{1:N}, Y_{1:N}) from a PMCModel.

Algorithm

State margins (model.margin_structure == "state", f_ij = f_i) — the hidden chain X is Markov (Proposition of DerrodePieczynski_CSDA2013 §2.1), so for every variant: 1. Sample the latent chain X_0, X_1, ..., X_N (X_0 is a virtual "warm-up" state, not returned). The chain is driven by the row-stochastic matrix A[i,j]. 2. Sample observations Y_1, ..., Y_N according to the variant-specific conditional:

 HMC-IN  : Y_n | X_n = j                    ~  f_j
 HMC-IN2 : Y_n | X_{n-1}=i, X_n=j           ~  f_j
 HMC-DN  : Y_1 | X_0=i, X_1=j               ~  f_j  (marginal, no copula)
           Y_n | X_{n-1}=i, X_n=j, Y_{n-1}  ~  f_j(·) · c_{ij}(F_i(Y_{n-1}), F_j(·))
 PMC-IN  : identical to HMC-IN2 (different prior, same conditional)
 PMC     : identical to HMC-DN  (different prior, same conditional)

Pair margins (model.margin_structure == "pair", general PMC, PMC and PMC-IN only) — X is not Markov, and the pairs z_n = (x_n, y_n) are drawn alternately, as in DerrodePieczynski_CSDA2013 §3.1 (Eqs. 13–14):

 x_1 ~ π,  π_i = Σ_j p[i,j];   y_1 | x_1 = i  ~  Σ_j A[i,j] f_ij
 x_{n+1} | x_n = i, y_n               ∝  p[i, x_{n+1}] · f_{i,x_{n+1}}(y_n)      (Eq. 13)
 y_{n+1} | x_n = i, x_{n+1} = j, y_n  ~  f_ji(·) · c_ij(F_ij(y_n), F_ji(·))     (Eq. 14)

PMC-IN drops the copula factor (y_{n+1} ~ f_ji). The mixture for y_1 is drawn through an auxiliary state j ~ A[x_1, ·], and y_{n+1} by conditional inversion of c_ij at u = F_ij(y_n), then F_ji^{-1}.

In both cases the right margin of the pair (i, j) is model.margin(j, i) = f_ji (the index inversion of DerrodePieczynski_CSDA2013 Eq. 12), which is f_j for state margins.

Returns:

Name	Type	Description
`X`	`np.ndarray, shape (N,), dtype int`	Hidden state sequence X_{1:N} with values in {0, ..., K-1}.
`Y`	`np.ndarray, shape (N,), dtype float`	Observed sequence Y_{1:N}.

simulate ¶

simulate(model: PMCModel, N: int | None = None, seed: int | None = None) -> tuple[np.ndarray, np.ndarray]

Generate a synthetic sequence from a PMCModel.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`	The model to simulate from.	required
`N`	`int`	Sequence length. Defaults to `model.N_default`.	`None`
`seed`	`int`	Seed for the random number generator (reproducibility).	`None`

Returns:

Name	Type	Description
`X`	`np.ndarray, shape (N,), dtype int`	Hidden state sequence X_{1:N} ∈ {0, …, K-1}.
`Y`	`ndarray`	Observed sequence Y_{1:N}. Shape `(N,)` for scalar models (`model.d == 1`); `(N, d)` when observations are vectors.

Inference¶

Forward-backward, MPM classification, and Forward-Filter Backward-Sample.

classify ¶

classify(model: PMCModel, Y: ndarray, *, gap_nodes: int | None = None) -> tuple[np.ndarray, np.ndarray, float]

Supervised MPM classification of the observation sequence Y.

Runs the full forward-backward-smooth-MPM pipeline.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel — fully specified model (parameters known).`		required
`Y`	`ndarray`	`rows are missing observations, integrated out (module docstring).`	required
`gap_nodes`	`quadrature nodes for the missing rows of the grid variants`	(default 64, see :mod:`pmcprg.pmc.gaps`); ignored otherwise.	`None`

Returns:

Name	Type	Description
`X_hat`	`np.ndarray, shape (N,) int — MPM class labels (every n,`	missing rows included).
`gamma`	`np.ndarray, shape (N, K) — posterior marginals P(x_n \| y_obs).`
`log_lik`	`float — log p(y_obs \| model), the`	observed-data log-likelihood (log p(y_{1:N}) without gaps); log p(y_obs, m \| model) when `model.missingness` is not None, γ being then P(x_n \| y_obs, m).

classify_image ¶

classify_image(model: PMCModel, img: ndarray) -> tuple[np.ndarray, np.ndarray, float]

Supervised MPM classification of a 2D image (mono- or multi-channel).

The image is linearised along the Generalized Hilbert ("gilbert") path, classified with :func:classify, then re-folded into 2D maps.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel — fully specified model (parameters known).`		required
`img`	`ndarray`	Shape `(H, W)` for grayscale (when `model.d == 1`) or `(H, W, d)` for multi-channel (when `model.d > 1`).	required

Returns:

Name	Type	Description
`X_hat_2d`	`np.ndarray, shape (H, W), int — MPM class-label map.`
`gamma_2d`	`np.ndarray, shape (H, W, K) — per-pixel posterior marginals.`
`log_lik`	`float — log p(image \| model).`

Notes

The class-label map is always 2D (segmentation is per-pixel, scalar) regardless of the input dimensionality. Posterior marginals are over states, so their last axis is always K — independent of model.d.

forward ¶

forward(model: PMCModel, Y: ndarray, W: ndarray | None = None, f_pdf: ndarray | None = None, *, gap_nodes: int | None = None) -> tuple[np.ndarray, float]

Normalized forward pass (Devijver 1985 / Baum-Welch).

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`Y`	`ndarray`		required
`W`	`ndarray \| None`		`None`
`f_pdf`	`pre-computed margin-PDF tensor (N, K, K). Computed with W if None.`		`None`
`gap_nodes`	`quadrature nodes for the missing rows of the grid variants`	(default 64, see :mod:`pmcprg.pmc.gaps`); ignored otherwise.	`None`

Returns:

Name	Type	Description
`alpha_hat`	`np.ndarray, shape (N, K) — normalized forward variables`
`log_lik`	`float — log-likelihood log p(y_{1:N})`
`Underflow`	`when a step's normaliser C is not finite or below`
	``MIN_POSITIVE`` — every weight reachable from α̂_n underflowed — the
	`pass is recomputed in log space from the model (see the module`
	docstring); ``W`` must therefore be ``precompute_weights(model, Y)``.
	`A step impossible even in log space raises`
	class:`IncompatibleObservationError`.
	`Missing rows (module docstring): the exact-shortcut variants run this`
	recursion on the marginalised ``precompute_weights``; the grid variants
	run the augmented chain of :mod:`pmcprg.pmc.gaps` (``W`` must then be None),
	``alpha_hat[n, i]`` being P(x_n = i \| observations up to n) at every n and
	``log_lik`` the observed-data log-likelihood.
	Non-ignorable missingness (``model.missingness`` not None): ``alpha_hat``
	and ``log_lik`` are those of p(y_obs, m) — the mask up to n is part of
	the conditioning; the factors live in ``W`` and ``f_pdf``
	(:func:`precompute_weights`), which must then come from the same model.

backward ¶

backward(model: PMCModel, Y: ndarray, W: ndarray | None = None, *, gap_nodes: int | None = None) -> np.ndarray

Normalized backward pass.

Each β̂n is rescaled by its own sum (Rabiner 1989 §V.A-type scaling), not by the forward constants C{n+1} as in Devijver (1985) — harmless because :func:smooth and :func:joint_posteriors renormalise (audit K-13).

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`Y`	`ndarray`		required
`W`	`ndarray \| None`		`None`
`gap_nodes`	`quadrature nodes for the missing rows of the grid variants`	(default 64, see :mod:`pmcprg.pmc.gaps`); ignored otherwise.	`None`

Returns:

Name	Type	Description
`beta_hat`	`(ndarray, shape(N, K))`
`Underflow`	a step whose sum D is not finite or below ``MIN_POSITIVE``
	`makes the pass restart in log space from the model, as in`
	func:`forward` (the former floor zeroed β̂ for every earlier step).
	Missing rows: exact-shortcut variants as in :func:`forward`. For the grid
	variants (``W`` must be None) β̂ at a missing n depends on the forward
	messages — ``β̂_n(i) ∝ Σ_g α̃_n(i, g) β̃_n(i, g) / Σ_g α̃_n(i, g)`` — so the
	forward pass is run too; ``smooth(forward(...)[0], backward(...))`` is
	then exact at every n (:mod:`pmcprg.pmc.gaps`).
	`Non-ignorable missingness: β̂ carries the factors of the later masks, as`
	``W`` does (:func:`precompute_weights`).

smooth ¶

smooth(alpha_hat: ndarray, beta_hat: ndarray) -> np.ndarray

Compute the posterior marginals γ_n(j) = P(X_n=j | y_{1:N}).

Parameters:

Name	Type	Description	Default
`alpha_hat`	`(N, K)`		required
`beta_hat`	`ndarray`		required

Returns:

Name	Type	Description
`gamma`	`(N, K) — each row sums to 1.`

joint_posteriors ¶

joint_posteriors(alpha_hat: ndarray, W: ndarray, beta_hat: ndarray) -> np.ndarray

Pairwise posteriors ξ_n(i, j) = P(X_n=i, X_{n+1}=j | Y), n = 0 … N-2.

The pairwise companion of :func:smooth (which returns the marginal γ_n). Fully vectorised::

ξ_n(i, j) ∝ α̂_n(i) · W[n, i, j] · β̂_{n+1}(j)

with a per-step normalisation. A step whose un-normalised mass is zero — its transition weights underflowed to 0.0 — falls back to the product γ_n(i) · γ_{n+1}(j) of the smoothed marginals (:func:smooth), the coupling that keeps both marginal-consistency invariants γ_n(i) = Σ_j ξ_n(i, j) and γ_{n+1}(j) = Σ_i ξ_n(i, j) (audit N-11). The former uniform 1/K² kept them only while γ itself was uniform, which was the case when an underflowed step zeroed α̂ and β̂; since :func:forward and :func:backward recover from underflow in log space, γ around such a step is informative. With uniform γ the two coincide.

Parameters:

Name	Type	Default
`alpha_hat`	(N, K) — normalised forward variables (:func:`forward`).	required
`W`	`ndarray`	required
`beta_hat`	`ndarray`	required

Returns:

Name	Type	Description
`xi`	`(N - 1, K, K)`

mpm ¶

mpm(gamma: ndarray) -> np.ndarray

Marginal Posterior Mode (MPM) classification.

Parameters:

Name	Type	Description	Default
`gamma`	`(N, K) — posterior marginals (rows sum to 1).`		required

Returns:

Name	Type	Description
`X_hat`	`(N,) int — state sequence maximising the marginal posterior.`

sample_posterior ¶

sample_posterior(model: PMCModel, Y: ndarray, rng: Generator, *, W: ndarray | None = None, f_pdf: ndarray | None = None, alpha_hat: ndarray | None = None, gap_nodes: int | None = None, return_y: bool = False) -> np.ndarray | tuple[np.ndarray, np.ndarray]

Draw one realisation of the state sequence X̃ ~ P(X | Y).

Implements Forward-Filter Backward-Sample (FFBS) — the standard posterior sampler for state-space models (Carter & Kohn 1994; Frühwirth-Schnatter 1994; and, for the discrete-state case used here, Chib 1996). DerrodePieczynski_CSDA2013 Eq. 23 samples forward instead, along p(x_{n+1} | x_n, y_{1:N}) — the same posterior law, drawn in the other direction.

Filter: compute (or accept pre-computed) α̂_n(j) = P(X_n=j | Y_{1:n}).
Sample X̃_{N-1} ~ α̂_{N-1} (the last filtered posterior is the smoothed posterior because there is no future observation).
Backward-sample, for n = N-2, …, 0::

P(X_n=i | X̃{n+1}, Y{1:N}) ∝ α̂n(i) · W[n, i, X̃{n+1}]

Z = (X, Y) being Markov (DerrodePieczynski_CSDA2013 Eq. 4), X_n and Y_{n+2:N} are independent given Z_{n+1} = (X_{n+1}, Y_{n+1}), which makes this sampler exact for pair margins (general PMC, X alone not Markov) as well as for state margins.

Used by SEM (Stochastic EM) to obtain the hard pseudo-labels on which the M-step is then run.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`Y`	`ndarray`		required
`rng`	`Generator`		required
`W`	`optional pre-computed tensors (skip re-computation).`		`None`
`f_pdf`	`optional pre-computed tensors (skip re-computation).`		`None`
`alpha_hat`	`optional pre-computed tensors (skip re-computation).`		`None`
`gap_nodes`	`quadrature nodes for the missing rows of the grid variants`	`(default 64); ignored otherwise.`	`None`
`return_y`	`bool`		`False`

Returns:

Name	Type	Description
`X_sample`	`np.ndarray, shape (N,), dtype int — one draw from P(X \| Y).`
`Y_sample`	(only with ``return_y``) copy of Y whose missing rows are drawn	from P(y_miss \| X̃, y_obs) — jointly with X̃.
	`Missing rows (non-finite rows of Y): the draw is from P(X \| y_obs). The`
	`exact-shortcut variants (HMC-IN, HMC-IN2, PMC-IN with state margins) run`
	`the FFBS below on the marginalised weights, then draw each missing`
	`y_n ~ f_{X̃_n}. The grid variants run FFBS on the augmented chain of`
	mod:`pmcprg.pmc.gaps` (``W``/``alpha_hat`` must be None); the missing y are
	`then drawn on its quadrature nodes.`
	Non-ignorable missingness (``model.missingness`` not None): the draw is
	from P(X \| y_obs, m) — the factors live in ``W`` (:func:`precompute_weights`)
	`and in the augmented chain; given X̃ the missing y are drawn as above`
	`(y_miss does not depend on m given the states).`

error_rate ¶

error_rate(X_true: ndarray, X_hat: ndarray) -> float

Classification error rate (fraction of misclassified labels), invariant under label permutation.

For unsupervised classifiers, predicted class indices are arbitrary. We therefore find the optimal one-to-one relabeling π̂ of predicted labels that minimises mean(X_true != π̂(X_hat)) and return that minimum.

Implementation: build the K×K confusion matrix C[k, ℓ] = #{n : X_true=k, X_hat=ℓ} then solve the linear assignment problem (Kuhn-Munkres / Hungarian) on -C to maximise total agreement.

Parameters:

Name	Type	Description	Default
`X_true`	`np.ndarray, shape (N,) — ground-truth labels in {0, …, K-1}.`		required
`X_hat`	`ndarray`		required

Returns:

Type	Description
`float — minimum error rate over all label permutations, in [0, 1].`

Missing observations: exact inference, imputation, forecasting¶

Missing values (NaN) are integrated out exactly rather than imputed before inference.

gap_posterior ¶

gap_posterior(model: PMCModel, Y: ndarray, *, gap_nodes: int | None = None, xi: bool = True) -> GapPosterior

Exact (up to quadrature) posterior quantities of a NaN-bearing Y.

Works for every variant and for a Y without missing values too (then identical to the functions of :mod:pmcprg.pmc.inference). The exact shortcut variants use precompute_weights / forward / backward / smooth / joint_posteriors on the marginalised weights; the grid variants run the augmented chain (module docstring). With a non-ignorable model.missingness every quantity is given (y_obs, m), the missingness factors included (module docstring, "Missingness mechanisms").

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`Y`	`ndarray`		required
`gap_nodes`	`number G of quadrature nodes (default 64; grid variants only).`	`Memory grows with G: about 0.3 kB per missing row and per node (K = 3; messages, grids, posteriors), plus 8·(K·G)² bytes per step of a gap that touches a local grid (295 kB at G = 64, 4.7 MB at G = 256), of which at most 2 GiB are kept and the rest rebuilt when needed (module docstring, "Memory").`	`None`
`xi`	`bool`		`True`

GapPosterior `dataclass` ¶

GapPosterior(miss: ndarray, log_lik: float, alpha_hat: ndarray, beta_hat: ndarray, gamma: ndarray, xi: ndarray | None, method: str, grid: QuadratureGrid | None = None, node_post: ndarray | None = None, nodes: ndarray | None = None, quad_error: float | None = None)

Posterior quantities of a sequence with missing observations.

Attributes:

Name	Type	Description
`miss`	`(N,) bool — missing rows.`
`log_lik`	`float — log p(y_obs), the observed-data log-likelihood`	(log p(y_obs, m) with a non-ignorable `model.missingness`, every posterior below being then given the mask m too).
`alpha_hat`	`(N, K) — P(x_n = i \| observations up to n).`
`beta_hat`	`(N, K) — normalised backward messages (module docstring for`	their definition at a missing n).
`gamma`	`(N, K) — P(x_n = i \| y_obs).`
`xi`	`(N-1, K, K) or None — P(x_n = i, x_{n+1} = j \| y_obs).`
`method`	``"exact"`` (K-state shortcut) or ``"grid"``.
`grid`	`QuadratureGrid or None (exact shortcut) — the reference grid.`
`node_post`	`(M, K, G) or None — P(x_n = i, y_n ∈ node g \| y_obs) at the`	M missing positions `np.nonzero(miss)[0]` (grid variants).
`nodes`	`(M, G) or None — the quadrature nodes y_g of each missing`	position: the reference nodes, or those of its local grid (module docstring, "Local grids").
`quad_error`	`float or None — the quadrature diagnostic of`	:func:`_quadrature_report` (grid variants): an estimate of Σ_runs \|error of the run's log-likelihood factor\| in nats, from the local errors of the steps into the missing rows (each discrete transition against the exact one-step law, on the backward message), weighted by the posterior; 0 for a trailing gap, whose factor is exactly 1. Inside a leading gap it also counts the errors of the posterior there (module docstring, "Convergence diagnostic"). Above `QUAD_WARN` (0.05 nats) a WARNING is logged.

index `property` ¶

index: ndarray

Positions of the missing rows.

impute ¶

impute(model: PMCModel, Y: ndarray, *, gap_nodes: int | None = DEFAULT_GAP_NODES, quantiles=DEFAULT_QUANTILES, n_samples: int = 0, rng=None) -> Imputation

Posterior law of every missing y_n given all the observed data.

Parameters:

Name	Type	Default
`model`	`PMCModel`	required
`Y`	`ndarray`	required
`gap_nodes`	`quadrature nodes G for the grid variants (default 64).`	`DEFAULT_GAP_NODES`
`quantiles`	`levels in (0, 1) of the reported posterior quantiles.`	`DEFAULT_QUANTILES`
`n_samples`	`number S of joint FFBS draws of (x_{1:N}, y_miss) (default 0).`	`0`
`rng`		`None`

Returns:

Type	Description
`Imputation — see its docstring. A Y without missing rows gives empty`
`per-position arrays.`

Imputation `dataclass` ¶

Imputation(index: ndarray, mean: ndarray, sd: ndarray, quantiles: tuple, quantile_values: ndarray, gamma: ndarray, log_lik: float, Y_mean: ndarray, method: str, nodes: ndarray | None = None, density: ndarray | None = None, x_samples: ndarray | None = None, y_samples: ndarray | None = None, grid_nodes: ndarray | None = None, grid_mass: ndarray | None = None)

Posterior law of the missing observations given all the observed ones.

Attributes:

Name	Type	Description
`index`	`(M,) positions of the missing rows.`
`mean, sd`	`(M,) — or (M, d) — posterior mean and standard deviation`	of y_n given y_obs.
`quantiles`	`the requested levels.`
`quantile_values`	`(M, Q) — or (M, Q, d), per component — posterior quantiles.`
`gamma`	`(M, K) — P(x_n = i \| y_obs) at the missing positions.`
`log_lik`	`float — log p(y_obs).`
`Y_mean`	`Y with every missing row replaced by its posterior mean.`
`method`	``"exact"`` or ``"grid"``.
`nodes`	`(G,) or None — the reference nodes y_g (d = 1).`
`density`	(M, G) or None — posterior density of y_n at ``nodes``	(the exact mixture density for the shortcut variants; mass / ω on the reference grid, the Nyström density on a local grid, see :func:`_grid_laws`).
`grid_nodes`	`(M, G) or None — the quadrature nodes of each missing`	position (grid variants: reference or local grid).
`grid_mass`	(M, G) or None — the posterior masses on ``grid_nodes``,	the discrete law whose moments are `mean` and `sd`.
`x_samples`	`(S, N) int or None — FFBS draws of the whole state path.`
`y_samples`	`(S, M) — or (S, M, d) — or None: the matching draws of the`	missing observations (exact continuous draws for the shortcut variants; drawn on the quadrature nodes — the discrete law whose moments are `mean` and `sd` — for the grid variants).

forecast ¶

forecast(model: PMCModel, Y: ndarray, h: int, *, gap_nodes: int | None = DEFAULT_GAP_NODES, quantiles=DEFAULT_QUANTILES, check_nodes: bool | int = False) -> Forecast

Predictive laws of the next h observations, a trailing gap of length h.

Runs the forward filter on Y followed by h missing rows: the normalised filter at N + k is P(x_{N+k}, y_{N+k} | y_obs) (no backward message is needed for a trailing gap). Y may contain missing rows.

With a non-ignorable model.missingness the laws are given (y_obs, m_{1:N}): the N rows of Y carry their missingness factors, the h appended rows none — their mask is unknown, and summing p(m_n | m_{n-1}, x_n) over it gives 1.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`Y`	`ndarray`		required
`h`	`int`		required
`gap_nodes`	`int \| None`	Memory, as for :func:`gap_posterior` over the missing rows of Y and the h appended ones: about 0.3 kB per missing row and per node (K = 3), plus at most 2 GiB of transitions between local grids, 8·(K·G)² bytes per step (module docstring, "Memory").	`DEFAULT_GAP_NODES`
`quantiles`			`DEFAULT_QUANTILES`
`check_nodes`	`False (default), True, or an integer G_ref > G. Checks that`	the predictive laws are converged in ``gap_nodes``, which ``GapPosterior.quad_error`` does not cover (a forecast horizon is a trailing gap; module docstring, "Convergence diagnostic"). The laws are computed a second time at G_ref nodes (True: 2G), and ``Forecast.node_check`` compares them (:class:`NodeCheck`). When they are not converged (:data:`FORECAST_CHECK_TOL`, :data:`FORECAST_CHECK_TOL_TAIL`), a WARNING is logged ("Forecast not stable in gap_nodes: … Increase gap_nodes."). The returned forecast is the one at G, the same bit for bit as without the check. The reference is a separate pass, run after it and without its own quadrature diagnostic (nor its WARNING). The exact variants are not recomputed: their laws do not depend on G. Cost: a whole forecast at 2G costs 1.6–3.6 times the one at G, mostly for its densities and quantiles. The check's pass computes only the moments and the node masses: 0.04– 0.64 times the forecast at G (measured: N = 500 with 25 missing rows, h = 24, G = 32–128, K = 2–3, τ = 0.6–0.99, the default quantiles). It runs once the forecast's chain is released; its transitions between local grids take 4 times the memory per step, 8·(K·2G)² bytes.	`False`

Returns:

Type	Description
`Forecast — see its docstring.`

Forecast `dataclass` ¶

Forecast(h: int, state_probs: ndarray, mean: ndarray, sd: ndarray, quantiles: tuple, quantile_values: ndarray, log_lik: float, method: str, nodes: ndarray | None = None, density: ndarray | None = None, grid_nodes: ndarray | None = None, grid_mass: ndarray | None = None, node_check: NodeCheck | None = None)

h-step predictive law of (x_{N+k}, y_{N+k}) given the observed part of Y.

Attributes:

Name	Type	Description
`h`	`horizon.`
`state_probs`	`(h, K) — P(x_{N+k} = j \| y_obs), k = 1..h.`
`mean, sd`	`(h,) — or (h, d) — predictive mean and standard deviation.`
`quantiles`	`the requested levels.`
`quantile_values`	`(h, Q) — or (h, Q, d) — predictive quantiles.`
`log_lik`	`float — log p(y_obs) of the conditioning sequence.`
`method`	``"exact"`` or ``"grid"``.
`nodes`	`(G,) or None — the reference nodes y_g (d = 1).`
`density`	(h, G) or None — predictive density at ``nodes``.
`grid_nodes`	`(h, G) or None — the quadrature nodes of each horizon`	(grid variants: reference or local grid).
`grid_mass`	(h, G) or None — the predictive masses on ``grid_nodes``.
`node_check`	:class:`NodeCheck` or None — the convergence of these	laws in `gap_nodes` (`check_nodes` of :func:`forecast`); None when not requested.

NodeCheck `dataclass` ¶

NodeCheck(gap_nodes: int, gap_nodes_ref: int, mean_sd: float, tails: float, converged: bool)

Convergence of a forecast's predictive laws in gap_nodes.

The result of :func:forecast's check_nodes: the forecast at gap_nodes = G against a second one at gap_nodes_ref = G_ref > G (2G by default), at every horizon, with the definitions of the forecasting study (report/forecasting/fc_common.py, quad_diff).

Attributes:

Name	Type	Description
`gap_nodes`	`int — G, the nodes of the returned forecast.`
`gap_nodes_ref`	`int — G_ref, the nodes of the reference forecast.`
`mean_sd`	`float — the largest \|mean(G) − mean(G_ref)\| and`	\|sd(G) − sd(G_ref)\| over the horizons, in units of sd(G_ref) at the same horizon.
`tails`	`float — the largest change of the 2.5 % and 97.5 %`	quantiles of the node law over the horizons, in the same unit. The node law is the discrete law on `grid_nodes` / `grid_mass`, its CDF linear between the cell edges (the midpoints of the sorted nodes).
`converged`	bool — mean_sd ≤ :data:`FORECAST_CHECK_TOL` (0.01) and	tails ≤ :data:`FORECAST_CHECK_TOL_TAIL` (0.05).
`The node-law quantiles are first order in the node spacing, and the`
returned ``quantile_values`` (the Legendre interpolant of the Nyström
`density) converge much faster`	``tails`` can flag a G whose returned
`quantiles are converged (τ = 0.3, G = 16`	`tails 0.28, while the returned`
`2.5 % and 97.5 % quantiles are within 1.1e-3 sd of those at G = 256; at`
`the default G = 64, tails 1.5e-3). It is the study's criterion`	`the`
`study scores its interval forecasts on the node law.`
``mean_sd`` and ``tails`` are inf, and ``converged`` False, when a
`predictive sd collapses`	`below 1e-3 of the sd of the observed values of`
`Y, at some horizon, at G or at G_ref. When Y has fewer than two distinct`
`observed values, the sd of the reference law g_ref is used instead. A law`
`that is genuinely that narrow is flagged too. The exact variants`
(``Forecast.method == "exact"``) get 0.0, 0.0 and True	`their laws do`
not depend on ``gap_nodes``, and nothing is recomputed.

Non-ignorable missingness mechanisms¶

The [missingness] block of a model's TOML file: the missingness mask itself carries evidence on the hidden states.

StateMissingness `dataclass` ¶

StateMissingness(rates: tuple[float, ...])

Mechanism "state": P(m_n = 1 | x_n = i) = π_i, independently given X.

Attributes:

Name	Type	Description
`rates`	`tuple of K floats in [0, 1] — π_i.`

evidence ¶

evidence(miss) -> np.ndarray

(N, K) factors e_n(i) = π_i if miss[n] else 1 − π_i.

log_evidence ¶

log_evidence(miss) -> np.ndarray

(N, K) log e_n(i) (−∞ where e_n(i) = 0), log(1 − π) by log1p.

to_table ¶

to_table() -> dict

The [missingness] TOML table of this mechanism.

StateMarkovMissingness `dataclass` ¶

StateMarkovMissingness(onset: tuple[float, ...], persistence: tuple[float, ...])

Mechanism "state-markov": a Markov mask whose transition depends on x_n.

Attributes:

Name	Type	Description
`onset`	`tuple of K floats in [0, 1] — a_i = P(m_n = 1 \| m_{n-1} = 0, x_n = i).`
`persistence`	`tuple of K floats in [0, 1] — b_i = P(m_n = 1 \| m_{n-1} = 1, x_n = i).`
`The initial mask follows the stationary law of the mask chain under a`
`constant state`	`P(m_0 = 1 \| x_0 = i) = s_i = a_i / (1 − b_i + a_i)`
`(`	attr:`stationary`). a_i = 0 with b_i = 1 is refused (s_i = 0/0).

stationary `property` ¶

stationary: ndarray

(K,) s_i = a_i / (1 − b_i + a_i) — the law of m_0 given x_0 = i.

evidence ¶

evidence(miss) -> np.ndarray

(N, K) factors e_n(i) = p(m_n | m_{n-1}, x_n = i); e_0(i) = p(m_0 | x_0 = i).

log_evidence ¶

log_evidence(miss) -> np.ndarray

(N, K) log e_n(i) (−∞ where e_n(i) = 0), log(1 − ·) by log1p.

to_table ¶

to_table() -> dict

The [missingness] TOML table of this mechanism.

missingness_lr_test ¶

missingness_lr_test(model, Y, *, alternative: str = 'state', null: str = 'common', ice_cfg: dict | None = None, n_bootstrap: int = 0, seed: int = 0, progress_cb=None) -> MissingnessLRTest

Likelihood-ratio test of state-independent against state-dependent missingness.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`	Starting point of the null fit (its own mechanism, if any, is only a start: the null fit replaces it). Its `[ice]` table and `ice_cfg` configure both ICE fits (`fit_margins`, `missing_strategy`, …).	required
`Y`	`array(N) or (N, d)`	Observations; the rows with a non-finite value are the mask m. Y must have at least one missing and one observed row.	required
`alternative`	``"state"`` or ``"state-markov"`` — H1.		`'state'`
`null`	``"common"`` (a state-independent mechanism of the same kind,	default) or `"state"` (only against `"state-markov"`) — H0.	`'common'`
`ice_cfg`	`dict`	ICE settings of the fits, over the model's `[ice]` table and the test's defaults (`tol = 1e-8`, `max_iter = 500`, `patience = 500`: ICE's fixed points). `missingness` is set by the test, and the alternative starts from the null fit (`init = "model"`), a `"state"` null from the ignorable fit.	`None`
`n_bootstrap`	`int`	Parametric-bootstrap replicates B (0: asymptotic p-value only). Each costs two ICE fits.	`0`
`seed`	`int`	Seed of the bootstrap (separate streams for the paths and the masks).	`0`
`progress_cb`	`callable`	`progress_cb(b, B, 0.0, "bootstrap")` before each replicate.	`None`

Returns:

Type	Description
class:`MissingnessLRTest` (see the module docstring for the statistic,
`the null log-likelihood and the two p-values).`

MissingnessLRTest `dataclass` ¶

MissingnessLRTest(alternative: str, null: str, statistic: float, df: int, p_value: float, p_value_bootstrap: float, n_bootstrap: int, n_bootstrap_valid: int, bootstrap_statistics: tuple, log_lik_null: float, log_lik_alt: float, null_model: object, alt_model: object, null_params: dict, alt_params: dict, n_obs: int, n_missing: int, sup_log_lik_null: float, sup_log_lik_alt: float, profile_log_liks: dict, statistic_fits: float)

Result of :func:missingness_lr_test.

Fields

alternative : "state" or "state-markov" (H1). null : "common" (state-independent) or "state" (H0). statistic : LR = 2 (sup_log_lik_alt − sup_log_lik_null) ≥ 0 (module docstring, "The statistic"). df : degrees of freedom of the asymptotic χ². p_value : asymptotic, χ²(df) survival function at statistic. p_value_bootstrap : parametric-bootstrap p-value, NaN without bootstrap. n_bootstrap : replicates requested (B). n_bootstrap_valid : replicates that produced a statistic. bootstrap_statistics : the replicate statistics LR_b (valid ones). log_lik_null, log_lik_alt : log p(y_obs, m) of the two ICE fits. null_model, alt_model : the fitted models (:class:~pmcprg.pmc.model.PMCModel), each carrying its mechanism (missingness). null_params, alt_params : the [missingness] tables of the two fitted mechanisms. n_obs, n_missing : N and M. sup_log_lik_null, sup_log_lik_alt : the approximations of the two suprema, max over the θ of both fits of the profile log-likelihoods profile_log_liks. profile_log_liks : {"null_at_null_theta": ℓ0(θ̂0), "null_at_alt_theta": ℓ0(θ̂1), "alt_at_null_theta": ℓ1(θ̂0), "alt_at_alt_theta": ℓ1(θ̂1)} — ℓ_h(θ) the log-likelihood maximised over hypothesis h's mechanism at θ. statistic_fits : 2 (log_lik_alt − log_lik_null), the difference of the two ICE fits alone — a diagnostic: it can be negative where ICE is not EM (module docstring).

summary ¶

summary() -> str

A few lines for a log or a console.

Erroneous observations: predictive PIT, flags, flag-and-mask estimation¶

One-step-ahead predictive PIT of every observed row (exact per variant, missing rows integrated out), outlier flags with innovation gating and multiplicity corrections, calibration checks, and a flag-and-mask ICE/SEM that trims the flagged rows through the missing-data machinery. The formulas, the multiplicity discussion and the design of the robust loop are in the docstring of pmcprg/pmc/outliers.py; the simulation study is in report/erroneous_data.

predictive_pit ¶

predictive_pit(model: PMCModel, Y, *, gap_nodes: int | None = None) -> PredictivePIT

One-step-ahead predictive PIT, p-values and normal scores of Y.

PIT_n = P(Y_n ≤ y_n | y_1:n−1) from the forward filter and the conditional CDFs of the transition (module docstring, per variant); missing rows before n are integrated out (exact shortcut, or the quadrature grid of :mod:pmcprg.pmc.gaps), missing rows get NaN. Every observed row enters the filter: to keep a flagged row out of the prediction of the next ones, use :func:flag_outliers.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`Y`			required
`gap_nodes`	`quadrature nodes G for the missing rows of the grid variants`	`(default 64).`	`None`

Returns:

Type	Description
`PredictivePIT`

PredictivePIT `dataclass` ¶

PredictivePIT(pit: ndarray, pvalue: ndarray, z: ndarray, log_pred: ndarray, log_lik: float, alpha_hat: ndarray, miss: ndarray, method: str)

One-step-ahead predictive PIT of a series (module docstring).

Attributes:

Name	Type	Description
`pit`	`(N,) PIT_n = P(Y_n ≤ y_n \| past observed rows); NaN at a`	missing row.
`pvalue`	`(N,) two-sided p-value 2 min(PIT_n, 1 − PIT_n), capped at 1;`	NaN at a missing row.
`z`	`(N,) normal score Φ⁻¹(PIT_n), from the smaller tail (finite:`	a tail probability is floored at the smallest normal float).
`log_pred`	`(N,) log p(y_n \| past), the log predictive density; NaN at a`	missing row.
`log_lik`	`float — Σ of the log normalisers of the pass: the`	observed-data log-likelihood of the rows the filter used.
`alpha_hat`	`(N, K) — P(x_n = i \| rows up to n used by the filter).`
`miss`	`(N,) bool — missing rows of Y (non-finite values).`
`method`	``"exact"`` (K-state filter) or ``"grid"`` (quadrature grid	for the missing rows).

flag_outliers ¶

flag_outliers(model: PMCModel, Y, *, alpha: float = DEFAULT_ALPHA, sequential: bool = True, correction: str | None = None, gap_nodes: int | None = None) -> OutlierFlags

Flag the observed rows whose predictive p-value is too small.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`Y`			required
`alpha`	`float`	``correction``, family-wise with ``"bonferroni"``, false discovery rate with ``"bh"``.	`DEFAULT_ALPHA`
`sequential`	`gate a flagged row out of the filter of the rows after it`	(default True; module docstring, "Sequential gating"). False: every row enters the filter, p-values are those of :func:`predictive_pit`.	`True`
`correction`	``None`` (default), ``"bonferroni"`` or ``"bh"`` — module	`docstring, "Flagging and multiplicity".`	`None`
`gap_nodes`	`int \| None`	Memory: the grids of the missing rows come from a pass of :func:`pmcprg.pmc.gaps.gap_posterior`'s chain — about 0.3 kB per missing row and per node (K = 3), plus at most 2 GiB of transitions between local grids, 8·(K·G)² bytes per step (``gaps`` module docstring, "Memory").	`None`

Returns:

Type	Description
OutlierFlags — ``flagged ⇔ pvalue < threshold`` at every observed row.

OutlierFlags `dataclass` ¶

OutlierFlags(flagged: ndarray, pvalue: ndarray, pit: ndarray, z: ndarray, log_pred: ndarray, threshold: float, alpha: float, correction: str | None, sequential: bool, n_tests: int, log_lik: float, alpha_hat: ndarray, miss: ndarray, method: str)

Rows flagged by :func:flag_outliers and the statistics used.

Attributes:

Name	Type	Description
`flagged`	`(N,) bool — flagged rows (never a missing row).`
`pvalue`	`(N,) two-sided p-values the decision used (with`	`sequential`, each from the filter gated by the flags before it).
`pit, z`	`(N,) the matching PITs and normal scores.`
`log_pred`	`(N,) log predictive densities of the rows that entered the`	filter; NaN at missing and flagged rows when `sequential`.
`threshold`	`float — the per-row threshold: flagged ⇔ p < threshold.`
`alpha`	`float — the requested level (per row, family-wise or FDR).`
`correction`	``None``, ``"bonferroni"`` or ``"bh"``.
`sequential`	`bool — whether flagged rows were gated out of the filter.`
`n_tests`	`int — number of tested (observed) rows.`
`log_lik`	`float — log-likelihood of the rows that entered the filter`	(flagged rows integrated out when `sequential`).
`alpha_hat`	`(N, K) — the filter P(x_n = i \| rows up to n it used).`
`miss`	`(N,) bool — missing rows of Y.`
`method`	``"exact"`` or ``"grid"``.

index `property` ¶

index: ndarray

Positions of the flagged rows.

pit_checks ¶

pit_checks(pit, *, lags: int = 10, exclude=None) -> PitChecks

Uniformity and independence checks of a PIT sequence.

Parameters:

Name	Type	Description	Default
`pit`		`array of PITs (NaN rows are skipped).`	required
`lags`	`int`		`10`
`exclude`	`optional (N,) bool mask of rows to leave out (e.g. the flagged`	`rows).`	`None`

Returns:

Type	Description
`PitChecks — KS test of U(0, 1) on the PITs; Ljung–Box on the normal`
`scores and on their squares, over the kept rows in order (the PITs of`
`the observed rows of a correct model are iid, whatever the missing rows`
`between them).`

PitChecks `dataclass` ¶

PitChecks(n: int, ks_stat: float, ks_pvalue: float, lb_stat: float, lb_pvalue: float, lb2_stat: float, lb2_pvalue: float, lags: int, z_mean: float, z_sd: float)

Calibration checks of a PIT sequence (:func:pit_checks).

Attributes:

Name	Type	Description
`n`	`number of PITs used (observed, not excluded).`
`ks_stat, ks_pvalue`	`Kolmogorov–Smirnov test of U(0, 1) on the PITs.`
`lb_stat, lb_pvalue`	Ljung–Box Q on the normal scores, ``lags`` lags.
`lb2_stat, lb2_pvalue`	`Ljung–Box Q on the squared normal scores`	(McLeod–Li).
`lags`	`number of lags.`
`z_mean, z_sd`	`mean and standard deviation of the normal scores (0 and 1`	under the model).

robust_estimate ¶

robust_estimate(model: PMCModel, Y, estim_cfg: dict | None = None, *, algorithm: str = 'ice', flag_cfg: dict | None = None, max_rounds: int = 10, initial_mask=None, progress_cb=None) -> RobustFit

Flag-and-mask estimation: fit, flag, set the flags to NaN, refit.

The design (masking as trimming, iteration against the masking effect, flags recomputed on the original Y, restarts from the initial model, breakdown, stopping rule) is in the module docstring, "Robust estimation".

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`Y`			required
`estim_cfg`	`dict \| None`	:func:`pmcprg.pmc.sem`, used unchanged by every fit.	`None`
`algorithm`	`str`		`'ice'`
`flag_cfg`	`dict \| None`	``sequential``, ``correction``, ``gap_nodes``); by default its defaults, with the ``gap_nodes`` of ``estim_cfg``.	`None`
`max_rounds`	`int`		`10`
`initial_mask`	`optional (N,) bool — rows set to NaN for the first fit`	(a pre-screen such as the Hampel filter of ``report/erroneous_data/run_study.py``, sensor-level flags); re-tested by the model from the first round on. Needed when the spikes are frequent enough to form a state of the first fit (module docstring, "Breakdown").	`None`
`progress_cb`			`None`

Returns:

Type	Description
`RobustFit`

RobustFit `dataclass` ¶

RobustFit(model: PMCModel, trace: object, mask: ndarray, flags: OutlierFlags, masks: list = list(), traces: list = list(), converged: bool = False)

Result of :func:robust_estimate.

Attributes:

Name	Type	Description
`model`	`PMCModel — the last fit.`
`trace`	`IceTrace or SemTrace of the last fit.`
`mask`	`(N,) bool — rows set to NaN for the last fit.`
`flags`	`OutlierFlags — the flags of the last fit on the original Y`	(`flags.flagged == mask` when `converged`).
`masks`	list of (N,) bool — the mask of every fit, ``masks[0]`` the	`initial_mask` (empty by default: the fit on Y as given).
`traces`	`list — the trace of every fit.`
`converged`	`bool — the last flags equal the last mask.`

ICE / SEM unsupervised estimation¶

Both estimators share the same M-step; ICE uses the soft forward-backward posteriors and converges deterministically, SEM draws one sample_posterior realisation per iteration. See "ICE vs SEM" in the README, and State labelling for how the hidden-state indices are numbered and identified across a fit.

ice ¶

ice.py — Iterative Conditional Estimation (ICE) for unsupervised PMC/HMC fitting.

Public API

ice(model, Y, ice_cfg=None) -> (PMCModel, IceTrace) Fit a PMCModel to the observation sequence Y using ICE. Returns the fitted model and an :class:IceTrace capturing per-iteration diagnostics (log-lik, τ, family, prior matrix, margin params; plus the losing runs when multistart is enabled). Use trace.log_liks for the plain log-likelihood history.

Algorithm (ICE for SR-PMC)

Given an initial model θ^(0) and observations Y = y_{1:N}:

Iterate until convergence:

E-step: Compute forward-backward quantities: α̂n(j), β̂_n(j) — normalized forward/backward variables γ_n(j) = P(X_n=j | Y) — marginal posteriors ξ_n(i,j) = P(X_n=i, X{n+1}=j | Y) — joint posteriors

M-step (Conditional estimation): 1. Prior p̂[i,j] = (1/N-1) Σ_n ξ_n(i,j)

2. For each pair (i,j) — copula selection + τ estimation:
   a. Collect pseudo-observations (u_n, v_n) = (F_{ij}(y_n), F_{ji}(y_{n+1}))
      weighted by ξ_n(i,j)  for n = 1, …, N-1  (F_{ij} = F_i for state
      margins). With ``copula_margins = "empirical"`` the F are the
      posterior-weighted empirical margins of :func:`_empirical_margin_cdfs`
      (a weighted rank pseudo-likelihood; the E-step keeps the model's F).
   b. Select best copula family from 'candidates' by weighted log-likelihood.
   c. Estimate τ_{ij} by weighted MLE on τ ∈ [τ_min, τ_max].

3. (Optional, ``fit_margins``) Margin re-estimation by weighted MLE, the
   declared family kept (or selected among ``candidates`` — GICE):
   * state margins f_i : y_n with weight γ_n(i);
   * pair margins f_ij (general PMC, DerrodePieczynski_CSDA2013
     Eqs. 12–14) — "dual view": y_n with weight ξ_n(i,j) and
     y_{n+1} with weight ξ_n(j,i), since the right margin of the
     pair (j,i) is f_ij. For a Gaussian margin,
     μ̂_{ij} = Σ w y / Σ w and σ̂²_{ij} = Σ w (y − μ̂)² / Σ w
     over that sample. An ICE-style estimator, not an exact EM
     M-step (see ``_pair_margin_sample``).

Relation to the papers' ICE (DerrodePieczynski_CSDA2013 §4.2, DerrodePieczynski_SP2016 §3)

This is a responsibility-weighted variant of ICE, not the scheme of Derrode & Pieczynski (2013, §4.2, Eqs. 21–24) and (2016, §3 steps (b)–(c), Remark 3.1). There, only the prior p_ij is updated through the conditional expectation of Eq. 22; the copula (and margin) parameters, whose conditional expectation is not computable, are estimated by the complete-data estimator applied to one posterior draw x^(q) ~ p(x | y, θ^q) — "replacing x_{1:N} by x^q_{1:N} (L = 1)" — i.e. on the hard sub-samples y^{ij}(x^(q)). Here every ξ_n(i,j) enters as a weight instead of a 0/1 membership, which makes the M-step deterministic and EM-like (for known margins it is EM's M-step for the copula parameters). :func:pmcprg.pmc.sem.sem is the hard-draw estimator and therefore the closest thing in the package to the papers' ICE with L = 1 — except that it also takes p_ij from the draw rather than from the expectation. Neither reproduces the papers' scheme to the letter; the reproduction report (report/csda2013_reproduction.tex, §Exp. 3) says which one it ran. (Audit K-4.)

ICE configuration (TOML [ice] section or dict)

fit_margins : bool (default False) — re-estimate margin parameters. max_iter : int (default 50) — maximum EM iterations. tol : float (default 1e-4) — relative log-likelihood convergence threshold. candidates : list[str] — SHORT_NAMEs of candidate copula families. default: all 1-parameter available families except Product. n_starts, multistart_seed, multistart_jitter, multistart_workers, multistart_families — multistart (see :func:_parse_ice_cfg for the full list). return_best_iterate : bool (default False) — return the iterate with the highest log-likelihood instead of the last one. missing_strategy, missing_draws, missing_seed, gap_nodes — missing observations (NaN rows of Y): see :func:ice, section "Missing observations". copula_margins : str (default "parametric") — "empirical" computes the copula step's pseudo-observations from posterior-weighted empirical margins instead of the model's F (AUDIT_COPULES FR-7 a; see :func:_parse_ice_cfg and :func:_empirical_margin_cdfs). missingness : str (default "model") — the missingness mechanism: the model's, held fixed; "ignorable"; or "state" / "state-markov", estimated (see :func:ice, section "Missingness mechanism").

References

ICE itself (the references DerrodePieczynski_SP2016 gives for the method):

Pieczynski, W. (1992). Statistical image segmentation. Machine Graphics and Vision 1(1/2), 261–268.
Delignon, Y., Marzouki, A. & Pieczynski, W. (1997). Estimation of generalized mixtures and its application in image segmentation. IEEE Trans. Image Processing 6(10), 1364–1375.
Giordana, N. & Pieczynski, W. (1997). Estimation of generalized multisensor hidden Markov chains and unsupervised image segmentation. IEEE Trans. PAMI 19(5), 465–475.

The two papers this package reproduces:

DerrodePieczynski_CSDA2013 — Derrode, S. & Pieczynski, W. (2013). Unsupervised data classification using pairwise Markov chains with automatic copulas selection. CSDA 63, 81–98.
DerrodePieczynski_SP2016 — Derrode, S. & Pieczynski, W. (2016). Unsupervised classification using hidden Markov chain with unknown noise copulas and margins. Signal Processing 128, 8–17.

IceTrace `dataclass` ¶

IceTrace(log_liks: list[float] = list(), tau_history: ndarray = (lambda: np.empty((0, 0, 0), dtype=float))(), family_history: list[list[list[str]]] = list(), p_history: ndarray = (lambda: np.empty((0, 0, 0), dtype=float))(), margin_history: list[list[dict]] = list(), multistart_runs: list[IceTrace] = list(), run_tag: str = '', candidates: list[str] = list(), best_iter: int = -1, returned_iter: int = -1, degenerate: list | None = None, missingness_history: list[dict | None] = list())

Per-iteration ICE history, captured for diagnostic visualisations.

Snapshot semantics

All time-indexed arrays/lists have length T = number of completed E-steps. The snapshot at index t reflects the model state as the log-likelihood log_liks[t] was computed — i.e. just before the M-step that produces the model used for log_liks[t + 1]::

   E-step ──→  M-step  ──→  E-step ──→ …
      ▲              ▲
      │              │
snapshot[t]     snapshot[t+1]
log_liks[t]     log_liks[t+1]

For T = 0 (no completed iteration), the array fields are empty arrays of shape (0, K, K) — never None.

Fields

log_liks : list[float] — log-likelihood trace. tau_history : np.ndarray shape (T, K, K), dtype float — τ_K of the copula at pair (i, j) for each iteration. NaN if the variant has no copula at that pair (e.g. HMC-IN) or the candidate selection failed. family_history : list[list[list[str]]] — same shape (T, K, K), each entry is the SHORT_NAME of the selected family (or "" if no copula). p_history : np.ndarray shape (T, K, K) — joint prior p[i,j] = P(X_n=i, X_{n+1}=j) (computed as π·A for HMC variants). margin_history : list[list[dict]] — at index t a list of margin blocks {"i": ..., "j": ..., "dist": ..., "params": {...}} in declaration order, mirroring model.margin_blocks(). multistart_runs: list[IceTrace] — when n_starts > 1, the traces from the non-best runs (the chosen run is the trace itself). Empty list when single-start. run_tag : str — label of this run ("unperturbed", "perturbed-3", …) when multistart was used. candidates : list[str] — copula candidate SHORT_NAMEs that ICE iterated over, captured for the family-ribbon legend. best_iter : int — index of the highest log_liks value (first one on ties, NaN ranked lowest); -1 when the trace is empty. Set whatever return_best_iterate. returned_iter : int — iterate q of the returned parameter set θ^q, so that log_liks[returned_iter] is its log-likelihood when returned_iter < len(log_liks). With return_best_iterate it equals best_iter. Without it, len(log_liks) - 1 when ICE stopped early (convergence, patience) and len(log_liks) when the run used all max_iter iterations: the model of the last M-step was returned without being evaluated (always the case for SEM). -1 when not recorded. degenerate : list[DegenerateFinding] | None — degenerate states of the returned model (:func:pmcprg.pmc._estim_common.degenerate_states), [] when none; None when not checked (the traces of multistart_runs, whose models are not kept). missingness_history : list[dict | None] — length T: at index t the [missingness] table of the iterate ({"mechanism": "state", "rates": [...]} or {"mechanism": "state-markov", "onset": [...], "persistence": [...]}, :mod:pmcprg.pmc.missingness), None when it is ignorable. Constant unless the missingness config key estimates the mechanism.

IceResult `dataclass` ¶

IceResult(initial_model: PMCModel, fitted_model: PMCModel, Y: ndarray, trace: IceTrace)

Bundle returned by GUI/diagnostics layers around an ICE run.

Carries everything the View-selector needs to draw the 12 ICE views (the trace alone is not sufficient — view J needs the initial model, view C and H need the observation sequence Y).

Fields

initial_model : :class:PMCModel — the model handed to ice(). fitted_model : :class:PMCModel — the best-run model returned by ICE. Y : np.ndarray — observations ICE was fitted to. trace : :class:IceTrace — per-iteration history.

ice ¶

ice(model: PMCModel, Y: ndarray, ice_cfg: dict | None = None, progress_cb=None) -> tuple[PMCModel, IceTrace]

Iterative Conditional Estimation (ICE) for unsupervised PMC/HMC fitting.

With n_starts > 1 (config), runs ICE several times from random perturbations of model and keeps the run with the highest final log-likelihood — this hardens the fit against multimodal log-likelihoods (a real concern for archimedean copula mixtures and HMC variants with many similar margins). The first start is always the unperturbed model. Parameter jitter keeps the copula families of model; with multistart_families = "random" or "sweep" the starts also change them, which is what leaves a wrong-family basin (see :func:_parse_ice_cfg).

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`Y`	`ndarray`	`non-finite rows are missing observations (below).`	required
`ice_cfg`	`dict \| None`	See :func:`_parse_ice_cfg` for the recognised keys.	`None`
`progress_cb`		``progress_cb(it, max_iter, log_lik, run_tag)`` after every E-step. Designed for the GUI's progress bar; safe to leave as ``None`` for headless / scripted use.	`None`

Returns:

Name	Type	Description
`fitted_model`	`PMCModel — model with updated parameters (the best run`	when multistart is enabled; its best iterate with `return_best_iterate`).
`trace`	`IceTrace — captured per-iteration diagnostics. Use`	`trace.log_liks` for the plain log-likelihood history. When multistart was active, `trace.multistart_runs` holds the traces of the non-best runs. `trace.best_iter` / `trace.returned_iter` locate the best and the returned iterates.

Degenerate states

The returned model is checked by :func:pmcprg.pmc._estim_common.degenerate_states (a state of stationary weight below 0.5 %, a margin standard deviation below 1 % of that of the observed data, a copula τ within 1e-3 of ±1): the findings are stored in trace.degenerate and, when there are any, logged as one WARNING line. On such a fit the likelihood may be unbounded — a variance collapsing on a data atom gains nats without limit — and should not be compared with others. Nothing in the estimate changes.

Missing observations

A row of Y with a non-finite value is missing (ignorable missingness, MCAR/MAR — :mod:pmcprg.pmc.gaps). A Y without missing rows runs the historical algorithm unchanged. Otherwise every E-step uses the missing-data inference of :mod:pmcprg.pmc.gaps (exact K-state shortcut for HMC-IN, HMC-IN2 and PMC-IN with state margins, quadrature grid of gap_nodes nodes otherwise), trace.log_liks is the observed-data log-likelihood log p(y_obs), and the M-step follows missing_strategy:

"available" (default; deterministic) — γ, ξ exact given y_obs at every n (gap_posterior); prior from ξ over all n; state margins on the observed y_n with weight γ_n(i); pair margins on the observed endpoints of their dual view (same ξ weights); copulas on the pairs whose two endpoints are observed, weight ξ_n(i, j). See :func:_m_step.
"impute" — multiple imputation: missing_draws completions (x, y_mis) ~ P(· | y_obs, θ^q) per iteration (seed missing_seed), the complete-data ICE estimator on each completed series, parameter estimates averaged (DerrodePieczynski_CSDA2013 Eq. 24 with L = missing_draws), prior from the exact ξ. See :func:_m_step_impute.

init = "kmeans" clusters the observed rows only (:func:_warmstart_from_kmeans); GICE margin selection runs on the samples above (observed values, or each completed series). Why "available" is the default (bias/RMSE against the missing rate) is documented at :data:DEFAULT_MISSING_STRATEGY.

Missingness mechanism

The config key missingness (:data:MISSINGNESS_MODES, also in the TOML [ice] table) says what ICE does with the law of the mask m (:mod:pmcprg.pmc.missingness, section "Estimation"):

"model" (default) — model.missingness is carried unchanged to every iterate and to the returned model, its parameters held fixed: ignorable when the model has none (the text above), otherwise every E-step is given (y_obs, m) — trace.log_liks is log p(y_obs, m), γ and ξ are P(· | y_obs, m), the completions of "impute" are drawn given (y_obs, m) and the complete-data E-step of each completed series keeps the observed mask. A Y without missing rows still has a mask (m = 0), which the E-step uses too. The M-step formulas are those above: the factors of p(m | x) do not involve the other parameters.
"ignorable" — any mechanism is dropped: every iterate and the returned model have missingness is None.
"state" / "state-markov" — the mechanism of that kind is estimated with the other parameters. E-steps as for "model"; the M-step adds :func:pmcprg.pmc.missingness.estimate_mechanism on the exact γ given (y_obs, m) — for "impute" too, as the prior takes the exact ξ. Start: the model's mechanism of that kind if it has one, otherwise the state-independent MLE of the mask, whose first E-step is the ignorable one. Each M-step adds one pseudo-observation at the pooled rate (the boundary guard: a rate at 0 or 1 would be absorbing). The returned model carries the estimate. A Y without missing rows identifies no rate (MLE 0: the ignorable model): WARNING and "ignorable".

trace.missingness_history records the [missingness] table of every iterate (None when ignorable), like the other parameters. Testing a state-dependent mechanism against a state-independent one: :func:pmcprg.pmc.missingness_lr.missingness_lr_test.

ice_image ¶

ice_image(model: PMCModel, img: ndarray, ice_cfg: dict | None = None, progress_cb=None) -> tuple[PMCModel, IceTrace]

ICE on a 2D image — linearises along the gilbert path then calls :func:ice.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`img`	`ndarray`	Shape `(H, W)` for grayscale (when `model.d == 1`) or `(H, W, d)` for multi-channel (when `model.d > 1`).	required
`ice_cfg`	`dict \| None`		`None`
`progress_cb`	callable, optional — same as :func:`ice`.		`None`

Returns:

Name	Type	Description
`fitted_model`	`PMCModel`
`trace`	`IceTrace`

ice_image ¶

ice_image(model: PMCModel, img: ndarray, ice_cfg: dict | None = None, progress_cb=None) -> tuple[PMCModel, IceTrace]

ICE on a 2D image — linearises along the gilbert path then calls :func:ice.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`img`	`ndarray`	Shape `(H, W)` for grayscale (when `model.d == 1`) or `(H, W, d)` for multi-channel (when `model.d > 1`).	required
`ice_cfg`	`dict \| None`		`None`
`progress_cb`	callable, optional — same as :func:`ice`.		`None`

Returns:

Name	Type	Description
`fitted_model`	`PMCModel`
`trace`	`IceTrace`

IceResult `dataclass` ¶

IceResult(initial_model: PMCModel, fitted_model: PMCModel, Y: ndarray, trace: IceTrace)

Bundle returned by GUI/diagnostics layers around an ICE run.

Carries everything the View-selector needs to draw the 12 ICE views (the trace alone is not sufficient — view J needs the initial model, view C and H need the observation sequence Y).

Fields

initial_model : :class:PMCModel — the model handed to ice(). fitted_model : :class:PMCModel — the best-run model returned by ICE. Y : np.ndarray — observations ICE was fitted to. trace : :class:IceTrace — per-iteration history.

IceTrace `dataclass` ¶

IceTrace(log_liks: list[float] = list(), tau_history: ndarray = (lambda: np.empty((0, 0, 0), dtype=float))(), family_history: list[list[list[str]]] = list(), p_history: ndarray = (lambda: np.empty((0, 0, 0), dtype=float))(), margin_history: list[list[dict]] = list(), multistart_runs: list[IceTrace] = list(), run_tag: str = '', candidates: list[str] = list(), best_iter: int = -1, returned_iter: int = -1, degenerate: list | None = None, missingness_history: list[dict | None] = list())

Per-iteration ICE history, captured for diagnostic visualisations.

Snapshot semantics

All time-indexed arrays/lists have length T = number of completed E-steps. The snapshot at index t reflects the model state as the log-likelihood log_liks[t] was computed — i.e. just before the M-step that produces the model used for log_liks[t + 1]::

   E-step ──→  M-step  ──→  E-step ──→ …
      ▲              ▲
      │              │
snapshot[t]     snapshot[t+1]
log_liks[t]     log_liks[t+1]

For T = 0 (no completed iteration), the array fields are empty arrays of shape (0, K, K) — never None.

Fields

log_liks : list[float] — log-likelihood trace. tau_history : np.ndarray shape (T, K, K), dtype float — τ_K of the copula at pair (i, j) for each iteration. NaN if the variant has no copula at that pair (e.g. HMC-IN) or the candidate selection failed. family_history : list[list[list[str]]] — same shape (T, K, K), each entry is the SHORT_NAME of the selected family (or "" if no copula). p_history : np.ndarray shape (T, K, K) — joint prior p[i,j] = P(X_n=i, X_{n+1}=j) (computed as π·A for HMC variants). margin_history : list[list[dict]] — at index t a list of margin blocks {"i": ..., "j": ..., "dist": ..., "params": {...}} in declaration order, mirroring model.margin_blocks(). multistart_runs: list[IceTrace] — when n_starts > 1, the traces from the non-best runs (the chosen run is the trace itself). Empty list when single-start. run_tag : str — label of this run ("unperturbed", "perturbed-3", …) when multistart was used. candidates : list[str] — copula candidate SHORT_NAMEs that ICE iterated over, captured for the family-ribbon legend. best_iter : int — index of the highest log_liks value (first one on ties, NaN ranked lowest); -1 when the trace is empty. Set whatever return_best_iterate. returned_iter : int — iterate q of the returned parameter set θ^q, so that log_liks[returned_iter] is its log-likelihood when returned_iter < len(log_liks). With return_best_iterate it equals best_iter. Without it, len(log_liks) - 1 when ICE stopped early (convergence, patience) and len(log_liks) when the run used all max_iter iterations: the model of the last M-step was returned without being evaluated (always the case for SEM). -1 when not recorded. degenerate : list[DegenerateFinding] | None — degenerate states of the returned model (:func:pmcprg.pmc._estim_common.degenerate_states), [] when none; None when not checked (the traces of multistart_runs, whose models are not kept). missingness_history : list[dict | None] — length T: at index t the [missingness] table of the iterate ({"mechanism": "state", "rates": [...]} or {"mechanism": "state-markov", "onset": [...], "persistence": [...]}, :mod:pmcprg.pmc.missingness), None when it is ignorable. Constant unless the missingness config key estimates the mechanism.

sem ¶

sem.py — Stochastic EM (SEM) for unsupervised PMC/HMC fitting.

Public API

sem(model, Y, sem_cfg=None) -> (PMCModel, SemTrace) Fit a :class:PMCModel to the observation sequence Y using SEM. Returns the fitted model and a :class:SemTrace capturing per-iteration diagnostics.

sem_image(model, img, sem_cfg=None) -> (PMCModel, SemTrace) Convenience wrapper that linearises a 2D image along the Generalized Hilbert ("gilbert") path and runs SEM on it.

Algorithm

Given an initial model θ⁽⁰⁾ and observations Y = y_{1:N}, SEM alternates:

E-step (filter) Compute α̂n(j) = P(X_n=j | Y{1:n}) via :func:forward.

S-step (stochastic completion) Draw one realisation X̃ ~ P(X | Y) via Forward-Filter Backward-Sample (see :func:sample_posterior). Build hard one-hot posteriors

    γ̃_n(i) = 1[X̃_n = i]
    ξ̃_n(i, j) = 1[X̃_n = i, X̃_{n+1} = j]

M-step (supervised-style) Update prior, margins (optional), copula τ from the hard posteriors — delegated to :func:pmcprg.pmc._estim_common.m_step so the implementation stays identical to ICE's M-step body. State margins f_i are fitted on {y_n : X̃n = i}; pair margins f_ij (general PMC, DerrodePieczynski_CSDA2013 Eqs. 12–14) on the dual-view sub-sample {y_n : (X̃_n, X̃{n+1}) = (i, j)} ∪ {y_{n+1} : (X̃n, X̃{n+1}) = (j, i)}, and the copula c_ij on the pseudo-observations (F_ij(y_n), F_ji(y_{n+1})) of the pairs drawn in (i, j) — the model's F, or with copula_margins = "empirical" the rescaled empirical CDFs of those drawn sub-samples. A pair that does not occur in the draw keeps its margin (and copula) unchanged.

Where SEM differs from ICE

Stochastic E-step. Each iteration draws a single X̃ ~ P(X|Y) instead of computing the soft posteriors ξ, γ. Successive iterations are random; the log-likelihood trace fluctuates rather than monotonically increasing.
Convergence detection. The log-lik does NOT converge in the deterministic sense. We therefore stop only at max_iter (no patience / tol early-stop). For inference, prefer averaging the post burn-in estimates externally (the full trace is exposed).
Multistart and K-means warm-start are supported with the same semantics as ICE — see :func:pmcprg.pmc._estim_common.warmstart_from_kmeans.

Missing observations

A row of Y with a non-finite value is missing (ignorable missingness, :mod:pmcprg.pmc.gaps); a Y without missing rows runs the algorithm above unchanged. Otherwise SEM is exact data augmentation (Tanner & Wong 1987): the S-step draws the states and the missing values jointly, (X̃, ỹ_mis) ~ P(x, y_mis | y_obs, θ^q), by forward-filter backward-sample on the missing-data chain of :mod:pmcprg.pmc.gaps (the exact K-state shortcut for HMC-IN, HMC-IN2 and PMC-IN with state margins — then ỹ_n ~ f_{X̃_n}; the quadrature grid of gap_nodes nodes otherwise — then ỹ_n is the drawn node). The M-step is the complete-data one on the completed series (ỹ_mis filled in), GICE margin selection included; log_liks is the observed-data log-likelihood log p(y_obs) of the forward pass. The k-means warm start clusters the observed rows only.

Missingness mechanism

The config key missingness is ICE's (:func:pmcprg.pmc.ice.ice, section "Missingness mechanism"; :mod:pmcprg.pmc.missingness, "Estimation"):

"model" (default) — model.missingness is carried unchanged to every iterate and to the returned model, its parameters held fixed. With a mechanism the S-step draws from P(x, y_mis | y_obs, m, θ^q) and log_liks is log p(y_obs, m) — the missingness factors are in the forward messages the draw comes from, for a complete Y (mask m = 0) too.
"ignorable" — any mechanism is dropped.
"state" / "state-markov" — estimated: the M-step adds ICE's formulas (:func:pmcprg.pmc.missingness.estimate_mechanism) on the one-hot drawn path γ̃, with the same start (the model's mechanism of that kind, or the state-independent MLE of the mask) and the same boundary guard. A Y without missing rows falls back to "ignorable" with a WARNING.

trace.missingness_history records the mechanism of every iterate.

Drawing ỹ_n on the quadrature nodes rather than from the continuous law is an approximation of the grid variants; measured against a continuous within-cell jitter and a 4× finer grid, it leaves the estimates unbiased where the jitter is not: see :data:pmcprg.pmc.ice._GAP_DRAW.

Relation to the papers' ICE

SEM is the closest estimator in the package to the ICE of the papers with L = 1 (DerrodePieczynski_CSDA2013 §4.2, DerrodePieczynski_SP2016 §3): both estimate the copula and margin parameters on one posterior draw x^(q) ~ p(x | y, θ^q). It differs in taking the prior p_ij from the draw too, where the papers use its conditional expectation. The :mod:pmcprg.pmc.ice module docstring ("Relation to the papers' ICE") sets the two side by side.

Reference

Celeux, G. & Diebolt, J. (1985). The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem. Computational Statistics Quarterly 2, 73–82.

Tanner, M. A. & Wong, W. H. (1987). The calculation of posterior distributions by data augmentation. Journal of the American Statistical Association 82(398), 528–540.

The awesomePMC implementation reuses the shared M-step / warm-start helpers in :mod:pmcprg.pmc._estim_common to avoid duplicating the update logic.

SemTrace `dataclass` ¶

SemTrace(log_liks: list[float] = list(), tau_history: ndarray = (lambda: np.empty((0, 0, 0), dtype=float))(), family_history: list[list[list[str]]] = list(), p_history: ndarray = (lambda: np.empty((0, 0, 0), dtype=float))(), margin_history: list[list[dict]] = list(), multistart_runs: list[IceTrace] = list(), run_tag: str = '', candidates: list[str] = list(), best_iter: int = -1, returned_iter: int = -1, degenerate: list | None = None, missingness_history: list[dict | None] = list(), sampled_X_history: ndarray = (lambda: np.empty((0, 0), dtype=int))())

Bases: IceTrace

Per-iteration SEM history — :class:pmcprg.pmc.ice.IceTrace plus the SEM-specific stochastic completions.

Inherits every :class:IceTrace field (log_liks, tau_history, family_history, p_history, margin_history, multistart_runs, run_tag, candidates, best_iter, returned_iter, degenerate, missingness_history) and the n_iters / __len__ helpers, so the two share all GUI/plotting code. Adds:

sampled_X_history : (T, N) int — the FFBS draw X̃ used at each M-step (no counterpart in :class:IceTrace).

SemResult `dataclass` ¶

SemResult(initial_model: PMCModel, fitted_model: PMCModel, Y: ndarray, trace: IceTrace)

Bases: IceResult

Bundle returned by GUI/diagnostics layers around an SEM run.

Identical structure to :class:pmcprg.pmc.ice.IceResult (initial_model, fitted_model, Y, trace); kept as a distinct class so callers can tell an SEM result from an ICE one by type.

sem ¶

sem(model: PMCModel, Y: ndarray, sem_cfg: dict | None = None, progress_cb=None) -> tuple[PMCModel, SemTrace]

Stochastic EM (SEM) for unsupervised PMC/HMC fitting.

With n_starts > 1 (config), runs SEM several times from random perturbations of model and keeps the run with the highest final log-likelihood — the same multistart strategy as :func:pmcprg.pmc.ice.ice, multistart_families included.

Parameters:

Name	Type	Description	Default
`model`	`PMCModel`		required
`Y`	`ndarray`	`non-finite rows are missing observations (module docstring, "Missing observations").`	required
`sem_cfg`	`dict \| None`	:func:`_parse_sem_cfg` for recognised keys.	`None`
`progress_cb`		``progress_cb(it, max_iter, log_lik, run_tag)`` after every iteration.	`None`

Returns:

Name	Type	Description
`fitted_model`	`PMCModel — the model at the last iteration (or the best`	run when multistart is enabled).
`trace`	`SemTrace — per-iteration history.`

Best iterate

With return_best_iterate = True the model returned is the iterate θ^q with the highest observed-data log-likelihood log_liks[q] (computed with θ^q before the S- and M-steps that give θ^(q+1); after the last M-step, θ^max_iter is evaluated once more, so the trace has max_iter + 1 log-likelihoods and max_iter draws). This is a heuristic, not an SEM estimator: SEM's iterates are a Markov chain whose stationary spread is the Monte-Carlo noise of one draw, and the likelihood-best iterate of that chain is not its mean — averaging the post-burn-in iterates is the usual estimate. It picks the point of the chain that fits best, and inherits its noise. trace.best_iter and trace.returned_iter locate it; multistart ranks each start by the log-likelihood of the model it returns.

Degenerate states

Checked on the returned model as for :func:pmcprg.pmc.ice.ice (trace.degenerate, one WARNING line).

sem_image ¶

sem_image(model: PMCModel, img: ndarray, sem_cfg: dict | None = None, progress_cb=None) -> tuple[PMCModel, SemTrace]

SEM on a 2D image — linearises along the gilbert path then calls :func:sem.

Same API as :func:pmcprg.pmc.ice.ice_image.

sem_image ¶

sem_image(model: PMCModel, img: ndarray, sem_cfg: dict | None = None, progress_cb=None) -> tuple[PMCModel, SemTrace]

SEM on a 2D image — linearises along the gilbert path then calls :func:sem.

Same API as :func:pmcprg.pmc.ice.ice_image.

SemResult `dataclass` ¶

SemResult(initial_model: PMCModel, fitted_model: PMCModel, Y: ndarray, trace: IceTrace)

Bases: IceResult

Bundle returned by GUI/diagnostics layers around an SEM run.

Identical structure to :class:pmcprg.pmc.ice.IceResult (initial_model, fitted_model, Y, trace); kept as a distinct class so callers can tell an SEM result from an ICE one by type.

SemTrace `dataclass` ¶

SemTrace(log_liks: list[float] = list(), tau_history: ndarray = (lambda: np.empty((0, 0, 0), dtype=float))(), family_history: list[list[list[str]]] = list(), p_history: ndarray = (lambda: np.empty((0, 0, 0), dtype=float))(), margin_history: list[list[dict]] = list(), multistart_runs: list[IceTrace] = list(), run_tag: str = '', candidates: list[str] = list(), best_iter: int = -1, returned_iter: int = -1, degenerate: list | None = None, missingness_history: list[dict | None] = list(), sampled_X_history: ndarray = (lambda: np.empty((0, 0), dtype=int))())

Bases: IceTrace

Per-iteration SEM history — :class:pmcprg.pmc.ice.IceTrace plus the SEM-specific stochastic completions.

Inherits every :class:IceTrace field (log_liks, tau_history, family_history, p_history, margin_history, multistart_runs, run_tag, candidates, best_iter, returned_iter, degenerate, missingness_history) and the n_iters / __len__ helpers, so the two share all GUI/plotting code. Adds:

sampled_X_history : (T, N) int — the FFBS draw X̃ used at each M-step (no counterpart in :class:IceTrace).

pmcprg.pmc¶

Model¶

PMCModel ¶

prior_p property ¶

transition_A property ¶

stationary_pi property ¶

path property ¶

raw property ¶

missingness property ¶

margin_structure property ¶

from_dict classmethod ¶

with_missingness ¶

ice_config ¶

sem_config ¶

margin_blocks ¶

copula_blocks ¶

margin ¶

pdf ¶

cdf ¶

ppf ¶

copula ¶

weight ¶

save ¶

Variant ¶

allows_pair_margins property ¶

per_class_margin property ¶

Simulation¶

simulate ¶

simulate ¶

Inference¶

classify ¶

classify_image ¶

forward ¶

backward ¶

smooth ¶

joint_posteriors ¶

mpm ¶

sample_posterior ¶

error_rate ¶

Missing observations: exact inference, imputation, forecasting¶

gap_posterior ¶

GapPosterior dataclass ¶

index property ¶

impute ¶

Imputation dataclass ¶

forecast ¶

Forecast dataclass ¶

NodeCheck dataclass ¶

Non-ignorable missingness mechanisms¶

StateMissingness dataclass ¶

evidence ¶

log_evidence ¶

to_table ¶

StateMarkovMissingness dataclass ¶

stationary property ¶

evidence ¶

log_evidence ¶

to_table ¶

missingness_lr_test ¶

MissingnessLRTest dataclass ¶

summary ¶

Erroneous observations: predictive PIT, flags, flag-and-mask estimation¶

predictive_pit ¶

PredictivePIT dataclass ¶

flag_outliers ¶

OutlierFlags dataclass ¶

index property ¶

pit_checks ¶

PitChecks dataclass ¶

robust_estimate ¶

RobustFit dataclass ¶

ICE / SEM unsupervised estimation¶

ice ¶

IceTrace dataclass ¶

IceResult dataclass ¶

ice ¶

ice_image ¶

ice_image ¶

IceResult dataclass ¶

IceTrace dataclass ¶

prior_p `property` ¶

transition_A `property` ¶

stationary_pi `property` ¶

path `property` ¶

raw `property` ¶

missingness `property` ¶

margin_structure `property` ¶

from_dict `classmethod` ¶

allows_pair_margins `property` ¶

per_class_margin `property` ¶

GapPosterior `dataclass` ¶

index `property` ¶

Imputation `dataclass` ¶

Forecast `dataclass` ¶

NodeCheck `dataclass` ¶

StateMissingness `dataclass` ¶

StateMarkovMissingness `dataclass` ¶

stationary `property` ¶

MissingnessLRTest `dataclass` ¶

PredictivePIT `dataclass` ¶

OutlierFlags `dataclass` ¶

index `property` ¶

PitChecks `dataclass` ¶

RobustFit `dataclass` ¶

IceTrace `dataclass` ¶

IceResult `dataclass` ¶

IceResult `dataclass` ¶

IceTrace `dataclass` ¶

SemTrace `dataclass` ¶

SemResult `dataclass` ¶

SemResult `dataclass` ¶

SemTrace `dataclass` ¶