awesomePMC

Copula-based pairwise and hidden Markov chains for the unsupervised classification of non-Gaussian time series and images.

awesomePMC implements pairwise Markov chains (PMC) — up to the general PMC of Derrode & Pieczynski (2013), whose observation margins f_ij are indexed by the pair of consecutive hidden states — and the hidden Markov chains they contain, with the dependence between consecutive observations modelled by bivariate copulas. A model is a plain TOML file: simulate it, classify a signal or an image with it, or estimate it from the data alone (ICE or SEM, with automatic selection of the copula and margin families). Missing observations are integrated out exactly. Python API, pmc command line and PyQt6 GUI.

Reference implementation of:

S. Derrode, W. Pieczynski. Unsupervised data classification using pairwise Markov chains with automatic copulas selection. Computational Statistics & Data Analysis 63 (2013), 81–98. doi:10.1016/j.csda.2013.01.027

S. Derrode, W. Pieczynski. Unsupervised classification using hidden Markov chain with unknown noise copulas and margins. Signal Processing 128 (2016), 8–17. doi:10.1016/j.sigpro.2016.03.008

Scientific lineage

awesomePMC builds on the pairwise Markov chain model introduced by Wojciech Pieczynski — a strict generalisation of the hidden Markov chain that drops the assumption that the hidden process alone is Markov, and lets the observation depend on a pair of consecutive hidden states:

W. Pieczynski. Pairwise Markov chains. IEEE Transactions on Pattern Analysis and Machine Intelligence 25(5) (2003), 634–639. doi:10.1109/TPAMI.2003.1195998

The two papers this library implements (above) are by Stéphane Derrode and Wojciech Pieczynski, and stand on further work led by Pieczynski and collaborators:

• Copulas in Markov chains, first proposed for the noise of a pairwise Markov chain by N. Brunel and W. Pieczynski: doi:10.1109/SSP.2003.1289350 (2003) and, in full, doi:10.1016/j.sigpro.2005.01.018 (Signal Processing 85 (2005), 2304–2315).
• Iterative Conditional Estimation (ICE), the unsupervised estimator this package's ice() implements: W. Pieczynski, Statistical image segmentation, Machine Graphics and Vision 1(1/2) (1992), 261–268 (no DOI); Y. Delignon, A. Marzouki & W. Pieczynski, doi:10.1109/83.624951 (1997).
• Segmentation with pairwise Markov chains: S. Derrode & W. Pieczynski, doi:10.1109/TSP.2004.832015 (IEEE Transactions on Signal Processing 52(9) (2004), 2477–2489).

sem() implements the complementary stochastic estimator, SEM, from Celeux & Diebolt (1985) — see REFERENCES.md for the full list, and Citation and references below for how to cite this work.

awesomePMC is built on Wojciech Pieczynski's introduction of pairwise Markov chains and of Iterative Conditional Estimation, and on the two papers above that he co-authored with Stéphane Derrode — the model and the estimator this package packages for reuse.

Why awesomePMC?

• The general pairwise Markov chain. Five variants, from the classical hidden Markov chain (HMC-IN) to the PMC of DerrodePieczynski_CSDA2013 Eqs. 12–14 with pair margins f_ij — the law of yₙ given (xₙ, xₙ₊₁) = (i, j), under which X is not a Markov chain — or state margins fᵢ. Margins are any scipy.stats family, or multivariate Gaussian (d > 1) for the copula-free variants.
• 41 copula families — elliptical (Gaussian, Student-t), Archimedean (Clayton, Gumbel–Hougaard, Frank, Joe, A12, A14, AMH, BB1, BB6, BB7, BB8) and their survival forms, extreme-value (Galambos, Hüsler–Reiss, Tawn types 1 and 2, the three-parameter Tawn, t-EV), explicit (Plackett, FGM, cubic section, independence), and the 90°/270° rotations that carry the one-sided families to τ < 0 — all parameterised by Kendall's τ, with numerically robust vectorised log-densities: log-space Archimedean kernels (Hofert, Mächler & McNeil 2012) and native elliptical log-densities, never floored, so strong dependence and the tails keep their likelihood.
• Supervised classification: normalised forward–backward, MPM decision, posterior marginals and forward-filtering backward-sampling draws.
• Unsupervised estimation: ICE and SEM (stochastic EM) with, for every pair of states, automatic copula-family selection among candidates (MLE, AIC, BIC, Huard's Bayesian evidence of DerrodePieczynski_CSDA2013 Eq. 20, Cramér–von Mises, cross-validated CIC xvcic); GICE margin-family selection of DerrodePieczynski_SP2016 (MLE, Kolmogorov, AIC, BIC); k-means warm start; multistart with parameter jitter or over copula families.
• Is a family choice significant? Analytic standard errors of copula parameters, weighted Vuong and Clarke tests with a HAC variance, confidence sets of families not significantly worse than the best — also for the pairs of states of an ICE run.
• Missing observations: exact inference and ICE/SEM estimation with NaN in the observations, impute (posterior mean, sd, quantiles, draws) and forecast (state probabilities and predictive law h steps ahead), missingness patterns with ImputeGAP's geometry and imputation / classification metrics (pmcprg.missing).
• Images: grey-level or multichannel images are linearised along a generalised Hilbert curve, then classified or estimated (classify_image, ice_image, sem_image, pmc classify-image, pmc estimate-image).
• Goodness of fit: multivariate Kolmogorov–Smirnov test (Naaman 2021), parametric-bootstrap calibration, Neyman smooth tests on margins.
• Three front-ends: Python API, pmc command line, PyQt6 GUI (model editor, simulation, classification, estimation and goodness-of-fit in background threads).
• Reproducible: report/ reproduces the experiments of CSDA 2013 (§3.2, §3.3, §4.3) with committed results and a LaTeX report; report/missing_benchmark/ compares exact marginalisation with naive fill-ins on 2500 simulated sequences. The test suite runs in CI on Python 3.11–3.13.

Table of Contents

• Installation
• Example notebooks
• Copula toolkit
  • Available families
  • Quick start
  • Fitting and model selection
  • Bivariate joint laws
• PMC models
  • Model variants
  • TOML model files
  • Python API
  • Missing observations, imputation and forecasting
  • Command-line interface
  • Graphical interface (PyQt6)
  • Logging
• Diagnostics — multivariate Kolmogorov–Smirnov test
• Reproducibility reports
• Folder structure
• Citation and references
• License

Installation

From PyPI (recommended)

pip install awesomepmc

Optional extras:

pip install "awesomepmc[gui]"
pip install "awesomepmc[gui,ml,image]"

From source

git clone https://github.com/SDerrode/awesomePMC.git
cd awesomePMC
pip install .

Development install

pip install -e ".[dev,gui]"
pytest -n auto -m "not slow"   # the fast suite, in parallel: a few minutes
pytest -n auto                 # everything, incl. the Monte-Carlo studies (~15 min)

Requirements

The code snippets below run from the root of a clone of the repository (they load the example models of pmcprg/pmc/models/). The models are also shipped in the installed package:

from importlib.resources import files
from pmcprg.pmc import PMCModel

mdl = PMCModel(files("pmcprg.pmc") / "models" / "pmc_gauss_k2.toml")

Example notebooks

examples/quickstart.ipynb walks through the full workflow — copula density, sampling and fitting, Sklar bivariate law, then PMC simulate → classify → ICE — in less than a minute:

pip install jupyter
jupyter notebook examples/quickstart.ipynb

examples/uci_har_smartphone.ipynb runs ICE and SEM on a real 3-D accelerometer recording of the Human Activity Recognition Using Smartphones dataset (UCI ID 240). The dataset is downloaded on first run and cached in data/uci_har/ (not versioned); if the download fails or is declined, the notebook falls back to a synthetic 3-D signal simulated from the bundled hmc_in_mvn_k2_d3.toml model and every later cell runs unchanged. It exercises the k-means warm start (init = "kmeans"), the SEM estimator (from pmcprg.pmc import sem) and the multivariate KS test (from pmcprg.diagnostics import mks_2samp).

Both notebooks are re-run by the test suite (tests marked slow).

Copula toolkit

Available families

All copulas are τ-parameterised: Kendall's τ is the first parameter of every family, so two families are always comparable at matched dependence. Twelve of them carry one or two extra shape parameters (df, delta, delta6, psi, nu, psi_u/psi_v), fitted jointly with τ by bounded MLE.

The registry CopulaEnum is the single source of truth — the tables below are a reading of it, not a second list:

from pmcprg.copulas import CopulaEnum

len(CopulaEnum.available())                           # 40
[e.value.SHORT_NAME for e in CopulaEnum.available()]  # 'Prod', 'Gauss', 'Student', …

41 families: 2 elliptical, 11 Archimedean, 4 survival (180° rotations), 14 rotations at 90°/270°, 6 extreme-value and 4 explicit. ID and SHORT_NAME are the registry's own; SHORT_NAME is what the name = key of a TOML model expects.

Elliptical

Archimedean

² τ is jointly constrained with θ: BB7 needs τ > τ_Joe(θ), i.e. θ < θ_Joe(τ), the Joe copula being its δ → 0 limit. BB1 and BB6 carry the same kind of constraint on their own second parameter. BB7 is the only family here that keeps θ rather than δ as its free parameter, and for a reason specific to it: τ is strictly increasing in δ at fixed θ, but not monotone in θ at fixed δ — it dips below the Clayton value once δ ≳ 3.44 — so (τ, δ) would be two-to-one. BB7 nests both Clayton (θ = 1) and Joe (δ → 0), with λ_U = 2 − 2^{1/θ} and λ_L = 2^{−1/δ}: one tail coefficient per parameter.

³ BB8 is Joe at delta8 = 1 and the independence copula at θ = 1 — not Frank, despite the usual name: Frank appears only as the joint limit θ → ∞, delta8 → 0 at fixed θ·delta8. Its upper-tail coefficient is 2 − 2^(1/θ) at delta8 = 1 and 0 for every delta8 < 1, a genuine discontinuity of the family. Unlike BB1, BB6 and the Tawn models, τ and the second parameter are not jointly constrained: the whole τ range is reached at every delta8. Its τ is the package's only Archimedean τ computed by quadrature (no elementary closed form; a ₃F₂ one exists and is unusable).

Survival (180° rotations)

A 180° rotation reflects both arguments: it leaves τ unchanged and swaps the two tails.

Extreme-value

Max-stable families, each given by a Pickands dependence function A. All are upper-tail dependent, and the three Tawn models are the only unrotated families here that model asymmetric dependence: C(u, v) ≠ C(v, u). Every elliptical, Archimedean, survival and explicit family is exchangeable.

¹ τ is jointly constrained with the weights — τ < ψ for types 1 and 2, τ < 1/(1/ψ_u + 1/ψ_v − 1) for the three-parameter model. Tawn3 is the package's first three-parameter family.

Explicit

Closed-form CDF, no generator or Pickands function.

Rotations at 90° and 270°

A 90° or 270° rotation reflects a single argument, turning a positively-dependent family into a negatively-dependent one. This gives negative dependence to families that cannot otherwise represent τ < 0 — Clayton, Gumbel-Hougaard, Joe, A12, A14 and BB1 are all one-sided. Both tail-dependence coefficients vanish for every rotation (the mass sits in the (0, 1) and (1, 0) corners), so the column is dropped here.

With the rotations, 20 of the 41 families reach τ < 0. BB7 (ID 40) and BB8 (ID 41) are positively dependent only: FR-8 closed on a fixed list that does not include the BB6/BB7/BB8 pairs, so no 90°/270° rotation is registered for them.

Nonparametric comparison tool (not a family)

EmpiricalBetaCopula (Segers, Sibuya & Tsukahara 2017) is not a 40th entry of CopulaEnum — it has no τ, no fit(), and ICE's family selection never sees it. Given a bivariate sample it is the nonparametric copula that sample suggests: an assumption-free benchmark to check whether a fitted parametric family's cdf/pdf track what the data alone say.

from pmcprg.copulas import EmpiricalBetaCopula

data = CopulaClayton(tau_k=0.6).sample(n=300, seed=0)
nonparam = EmpiricalBetaCopula(data)       # or .from_pseudo_obs(pseudo_obs)

nonparam.cdf(0.3, 0.7)                     # C_n^beta(u, v)
nonparam.pdf(0.3, 0.7)                     # density, exp(logpdf) internally
nonparam.h1(0.3, 0.7), nonparam.h2(0.3, 0.7)  # conditional cdfs (h-functions)
nonparam.sample(1000, seed=1)              # exact two-stage mixture sampler

Built from ranks (average-rank tie handling, matching CopulaVirt.fit's own rankdata(...)/(n+1) pseudo-observations); its margins are exactly uniform when there are no ties (a deterministic identity of the incomplete beta function, not just "asymptotic"), and only a small, ties-only, shrinking-with-n discrepancy otherwise. See pmcprg/copulas/_nonparametric.py for the full derivation and pmcprg/tests/test_empirical_beta_copula.py for the numerical checks.

Quick start

from pmcprg.copulas import CopulaGaussian, CopulaStudent, CopulaClayton

# Instantiate by Kendall's τ
gauss  = CopulaGaussian(tau_k=0.5)
stud   = CopulaStudent(tau_k=0.5, df=4.0)    # df defaults to 4.0
clay   = CopulaClayton(tau_k=0.6)

# Evaluate
u, v = 0.3, 0.7
print(gauss.pdf([u, v]))          # copula density c(u,v)
print(gauss.cdf([u, v]))          # copula CDF C(u,v)
print(gauss.conditional_cdf(v, u))# h-function h(v|u) = ∂C/∂u
print(gauss.tail_dependence())    # (λ_L, λ_U)

# Sampling  (Rosenblatt inversion)
samples = gauss.sample(n=1000, seed=42)   # shape (1000, 2) on [0,1]²

Fitting and model selection

Data are rank-transformed to pseudo-observations before fitting, so the margins need not be uniform.

from pmcprg.copulas import CopulaClayton, CopulaFrank, CopulaGaussian, CopulaGH
from pmcprg.copulas._base import CopulaVirt

data = CopulaClayton(tau_k=0.5).sample(n=500, seed=0)

# Single-family fit
result = CopulaGaussian.fit(data, method='tau')   # or method='mle'
print(result.tau_k, result.aic, result.bic)

# 5-fold cross-validation log-likelihood
print(result.cv_loglik(K=5))

# Parametric bootstrap GoF  (Cramér-von Mises)
gof = result.gof_test(B=200)
print(gof.p_value)

# Bootstrap CI on τ
print(result.bootstrap_ci(B=500))

# Analytic standard errors on τ and the native parameters (sandwich for 'mle')
se = CopulaClayton.fit(data, method='mle').standard_errors()
print(se, se.ci())

# Automatic family selection (every family if `families` is omitted)
ranked = CopulaVirt.fit_best(
    data, families=[CopulaGaussian, CopulaClayton, CopulaGH, CopulaFrank], method='mle')
best = ranked[0]                        # sorted by AIC ascending
print(best.copula, best.aic)

# Is the ranking significant?  Vuong tests (HAC variance) and the ties
print(ranked.compare())                 # or ranked.compare(test="clarke")
print(ranked.confidence_set().members)  # families not significantly worse than the best

Bivariate joint laws

from pmcprg.copulas import CopulaGaussian
from pmcprg.copulas.bivariate import BivariateLaw
import scipy.stats as ss

law = BivariateLaw(
    copula=CopulaGaussian(tau_k=0.6),
    left_margin=(ss.norm, -1, 1),    # (scipy_dist, *params)
    right_margin=(ss.norm, +1, 1),
)
print(law.pdf([-0.5, 0.8]))   # joint density f(x,y)
print(law.cdf([ 0.0, 1.0]))   # joint CDF F(x,y)
law.set_seed(0)               # reproducible sampling
samples = law.sample(1000)    # shape (1000, 2)

# IFM two-step fit
fit = BivariateLaw.fit(samples, copula_class=CopulaGaussian,
                       left_family=ss.norm, right_family=ss.norm)
print(fit.copula_fit.tau_k, fit.aic)

PMC models

Model variants

Five variants of the jointly Markovian pair process (X, Y) are supported, from the simplest to the most general. Every model has one of two margin structures ([model].margin_structure):

• "state" — K densities fᵢ, one per state (f_ij = fᵢ). By the Proposition of DerrodePieczynski_CSDA2013 §2.1, this is exactly the case where X is a Markov chain: a PMC with state margins is a stationary reversible HMC-DN, and a PMC-IN with state margins an HMC-IN. Every HMC-* variant has state margins.
• "pair" — K² densities f_ij, the law of yₙ given (xₙ, xₙ₊₁) = (i, j): the general PMC of DerrodePieczynski_CSDA2013 Eqs. 12–14, where X is not Markov. Allowed for PMC and PMC-IN only.

Stationarity and reversibility only make the right margin of the pair (i, j) the left margin of (j, i) — the index inversion f_ji of DerrodePieczynski_CSDA2013 Eq. 12; they do not make f_ij independent of j. Versions 0.5.0–0.8.x collapsed K² margins to fᵢ on that mistaken ground.

TOML model files

Models are stored as TOML files. State margins use the K-format (one block per state, indexed by i only); pair margins use the K²-format (one block per pair, keys i and j, see pmc_pair_gauss_k2.toml):

[model]
name      = "PMC Gaussien K=2"
variant   = "PMC"        # HMC-IN | HMC-IN2 | HMC-DN | PMC-IN | PMC
K         = 2
N_default = 5000
# margin_structure = "state"   # optional: "state" | "pair"; default inferred
#                                # (K blocks → state, K² blocks → pair)

[prior]
# HMC-* variants: key "A" — K×K row-stochastic transition matrix
# PMC-* variants: key "p" — K×K symmetric joint distribution (sums to 1)
p = [[0.45, 0.05],
     [0.05, 0.45]]

# K margin blocks (one per state) — state margins f_i.
[[margins]]
i      = 0
dist   = "norm"          # any scipy.stats distribution name
params = {loc = -1.0, scale = 1.0}
# candidates = ["norm", "gamma", "invgamma", "betaprime"]   # GICE: families to
#                        # select from when fit_margins = true

[[margins]]
i      = 1
dist   = "norm"
params = {loc = 1.0, scale = 1.0}

[[copulas]]              # only for HMC-DN and PMC — K² blocks (i, j)
i    = 0
j    = 0
name = "Gauss"           # SHORT_NAME — see "Available families" above
tau  = 0.6
# df = 4.0               # extra shape param of the 12 families that take one
                         # (df, delta, delta6, psi, nu, psi_u/psi_v)
# … blocks (0, 1), (1, 0) and (1, 1)

[ice]                    # optional — ICE estimator defaults (also consumed
                         # by SEM for the keys they share)
fit_margins = false
max_iter    = 50
tol         = 1e-4
candidates  = ["Gauss", "Clayton", "GH", "Frank", "Joe"]
selection_criterion = "mle"      # "mle" (default) | "aic" | "bic" | "huard" | "cvm" | "xvcic"
                                 # (+ "huard_common", "huard_global")
margin_selection_rule = "mle"    # GICE: "mle" (default) | "kolmogorov" | "aic" | "bic"
init        = "model"    # "model" (default) | "kmeans"
# kmeans_seed = 0        # RNG seed for sklearn.cluster.KMeans (init="kmeans")
# n_starts = 1           # multistart: best final log-likelihood among n_starts runs
# multistart_families = "none"   # "none" (default) | "random" | "sweep"
# return_best_iterate = false    # true: return the highest-log-likelihood iterate, not the last

# [sem]                  # optional — SEM-specific overrides
# max_iter = 30          # SEM does not converge; defaults to fewer iterations
# sem_seed = 0           # RNG seed for the per-iteration FFBS draw

ICE vs SEM. Both estimators share the same M-step and config keys. ICE (algorithm = "ice", default) uses soft posteriors from forward-backward and converges deterministically. SEM (algorithm = "sem") draws a single realisation X̃ ~ P(X | Y) at every iteration via Forward-Filter Backward-Sample, then runs the same M-step on the hard labels — the log-likelihood fluctuates around its stationary regime instead of converging. CLI: pmc estimate --algorithm sem --sem-seed 0 …. GUI: the Estimator combobox in the ICE-config tab. Note that the papers (DerrodePieczynski_CSDA2013 §4.2, Eqs. 22–24; DerrodePieczynski_SP2016 §3) estimate copulas and margins on one posterior draw (L = 1) and reserve the expectation for the prior — closer to SEM than to this package's ICE; see the docstring of pmcprg/pmc/ice.py.

Selection criteria. At every M-step, each pair (i, j) gets the copula family of candidates that is best under selection_criterion; its τ is always the maximum-likelihood estimate. A margin block that declares candidates (GICE, DerrodePieczynski_SP2016 §3) gets its family chosen by margin_selection_rule when fit_margins = true — see sp2016_gice_k2.toml.

K-means warm-start. Setting init = "kmeans" clusters Y with k-means++ and derives a hard-labelled warm-start model via a single supervised-style M-step (prior + copula τ always re-estimated; margins re-estimated only if fit_margins = true) — useful when the declared initial parameters are far from the data. Requires scikit-learn (pip install "awesomepmc[ml]"). Works for both ICE and SEM.

Multistart over copula families. With n_starts > 1 ICE/SEM keep the run with the highest final log-likelihood; the extra starts jitter the parameters (multistart_jitter) but keep the copula families of the initial model. Those families can decide the fixed point: on the CSDA-2013 Exp. 3 design, ICE started at independence picks Gumbel instead of the true Gaussian copula on pair (0, 0) and no amount of jitter changes it. multistart_families = "random" also redraws the family of every pair (i, j) among candidates at each extra start (τ uniform on the central 60 % of the family's range); "sweep" runs every combination of candidate families on the diagonal pairs (i, i) at τ = 0.5, after the model itself — set n_starts = 1 + |C|^K (10 for 3 candidates and K = 2; fewer drops combinations with a warning, more adds random starts; beyond 256 combinations use "random"). The winning start is trace.run_tag (family-sweep:Gauss/Clayton, family-random-4). Variants without copulas ignore the key. GUI: Start families in the ICE-config tab.

Best iterate. ICE is not monotone, and after patience consecutive regressions it stops on its last iterate, possibly several nats below its best one. return_best_iterate = true returns the iterate θ^q with the highest trace.log_liks[q] (computed with θ^q, before the M-step); a run that reaches max_iter evaluates the model of its last M-step once more, so its trace has max_iter + 1 entries. Multistart then ranks the starts by the log-likelihood of the model each one returns. For SEM it is a heuristic (the best point of a noisy chain, not its average). trace.best_iter is the argmax of log_liks and trace.returned_iter the iterate returned. Default false (unchanged results). CLI: pmc estimate --best-iterate. GUI: Return best iterate in the ICE-config tab.

Degenerate fits. The likelihood is unbounded: a variance collapsing on an atom of discretised data, or a copula driven to τ = ±1, gains nats without limit. ICE and SEM check the returned model and log one WARNING (also printed by pmc estimate, which adds a Degenerate summary line) when a state has a stationary weight below 0.5 %, a margin (state or pair) a standard deviation below 1 % of that of the data, or a copula τ within 1e-3 of ±1. The findings are in trace.degenerate ([] when none); the fit itself is unchanged. Thresholds and their justification: pmcprg.pmc._estim_common.degenerate_states.

K²-format and older files. Blocks indexed by (i, j) are kept as pair margins f_ij on PMC and PMC-IN (an INFO line says so once per process), and ICE/SEM estimate each f_ij separately. Older files that relied on the collapse to K state margins must add margin_structure = "state" under [model] (the (i, 0) block is then kept as fᵢ). HMC-* variants refuse pair margins.

Estimating pair margins. ICE and SEM fit f_ij by weighted maximum likelihood on both views of the pair density: yₙ with weight ξₙ(i, j) and yₙ₊₁ with weight ξₙ(j, i) (each halved, so the total weight stays N − 1). It is an ICE-style estimator, not an exact EM M-step. Copula pseudo-observations for c_ij are (F_ij(yₙ), F_ji(yₙ₊₁)) with weight ξₙ(i, j).

Nine example models are provided in pmcprg/pmc/models/:

Python API

from pmcprg.pmc import PMCModel, simulate, classify, ice

# ── Load model ────────────────────────────────────────────────────
mdl = PMCModel("pmcprg/pmc/models/pmc_gauss_k2.toml")
print(mdl)                    # PMCModel(name='PMC Gaussien K=2', …)
print(mdl.stationary_pi)      # [0.5, 0.5]
print(mdl.transition_A)       # [[0.9, 0.1], [0.1, 0.9]]
print(mdl.prior_p)            # [[0.45, 0.05], …]

# Margin/copula access
print(mdl.pdf(0, 1, -0.5))    # f_{01}(-0.5)
print(mdl.cdf(1, 0, 0.3))     # F_{10}(0.3)
cop = mdl.copula(0, 0)        # CopulaGaussian(τ=0.6)

# ── Simulate ──────────────────────────────────────────────────────
X, Y = simulate(mdl, N=5000, seed=42)
# X: (5000,) int — latent states 0/1
# Y: (5000,) float — observations

# ── Supervised classification (MPM) ───────────────────────────────
X_hat, gamma, log_lik = classify(mdl, Y)
# gamma: (5000, K) — posterior P(X_n=k | Y)
from pmcprg.pmc import error_rate
print(f"Error rate: {error_rate(X, X_hat):.3f}")

# ── Unsupervised estimation (ICE) ─────────────────────────────────
raw_init = mdl.raw            # deep copy of the TOML dict — safe to mutate
# … perturb parameters …
init_mdl = PMCModel.from_dict(raw_init)
fitted, trace = ice(init_mdl, Y,
    ice_cfg={"max_iter": 30, "candidates": ["Gauss", "Clayton", "GH"]})
print(trace.log_liks[-1])     # final log-likelihood
print(trace.best_iter, trace.returned_iter, trace.degenerate)
fitted.save("fitted_pmc.toml")

# ── Edit and save ─────────────────────────────────────────────────
mdl.save("my_model.toml")

sem has the same signature as ice (sem_cfg= instead of ice_cfg=), and classify_image, ice_image and sem_image take a 2-D array (H, W) or (H, W, d) instead of Y.

Missing observations, imputation and forecasting

Observations may contain NaN (empty, NaN or NA cells in CSV files). With a known model, supervised inference integrates the missing values out exactly (up to quadrature), for every variant and both margin structures:

import numpy as np
from pmcprg.pmc import PMCModel, simulate, classify, impute, forecast
from pmcprg.missing import patterns, metrics

mdl = PMCModel("pmcprg/pmc/models/pmc_pair_gauss_k2.toml")
X, Y = simulate(mdl, N=2000, seed=1)
Y_gap, mask = patterns.mcar(Y, 0.2, block_size=10, seed=2)   # ImputeGAP-like pattern

X_hat, gamma, log_lik = classify(mdl, Y_gap)     # states at every n, missing ones included
imp = impute(mdl, Y_gap, quantiles=(0.05, 0.5, 0.95), n_samples=100, rng=np.random.default_rng(0))
# imp.index, imp.mean, imp.sd, imp.quantile_values, imp.gamma, imp.x_samples, imp.y_samples
fc = forecast(mdl, Y_gap, h=10)                  # fc.state_probs, fc.mean, fc.sd, fc.quantile_values
print(metrics.error_rate_split(X, X_hat, mask))  # error on missing / observed positions

• Exact shortcut — HMC-IN, HMC-IN2 and PMC-IN with state margins: the transition of X does not depend on a missing value, whose density is set to 1.

• Augmented grid — HMC-DN, PMC and PMC-IN with pair margins: inside a run of missing values the forward–backward runs on (state, value) over gap_nodes (default 64) Gauss–Legendre nodes in the quantiles of the stationary law of y, with an endpoint transform for copula tails. Checked against exact references (brute force over state paths, Gaussian AR(1) bridges and forecasts, distinct Gaussian regimes, quad for Clayton/Gumbel pair models): log-likelihood and posteriors to 1e-7 or better at G = 64 except in extreme tails (test_gaps_references.py).

• A forecast is a trailing gap; forecast returns the state probabilities and the predictive law of the observation k = 1, …, h steps ahead.

• Estimation: ice, sem (and pmc estimate, which prints the number of missing values) estimate from the observed data; trace.log_liks is log p(y_obs). SEM draws the states and the missing values jointly at each iteration. ICE has two strategies (config key missing_strategy): "available" (default; exact posteriors given the observed data, margins fitted on observed values, copulas on pairs with both ends observed; deterministic) and "impute" (missing_draws completed series per iteration, estimates averaged; better when GICE selects margin families). Complete data give bit-identical results to the previous code.

  fitted, trace = ice(init_mdl, Y_gap, ice_cfg={"missing_strategy": "available"})
  

• pmcprg.missing.patterns generates missingness patterns with the geometry of ImputeGAP's GenGap (mcar, aligned, scattered, blackout, disjoint, overlap, gaussian, distribution); pmcprg.missing.metrics scores imputations (RMSE, MAE, MI, Pearson, CRPS, interval coverage) and classification on missing vs observed positions.

State-dependent missingness

Everything above treats the mask as ignorable (MCAR/MAR): it carries no information on the states. Use a state-dependent mechanism instead when dropouts track the hidden regime — e.g. sensor windows go missing far more often during high-intensity activity than at rest, as observed on PAMAP2 (report/real_series/README.md, "Missingness is state-dependent"). A [missingness] table on the model (pmcprg.pmc.missingness) then makes the mask evidence on the states, m ⟂ y | x:

[missingness]
mechanism = "state"          # per-state rates, independent rows
rates = [0.01, 0.3]          # π_i = P(row missing | x_n = i)

# or, for bursts of gaps whose length depends on the state:
# [missingness]
# mechanism = "state-markov"
# onset = [0.002, 0.02]        # a_i: P(gap starts | x_n = i)
# persistence = [0.8, 0.95]    # b_i: P(gap continues | x_n = i)

from pmcprg.pmc import PMCModel, StateMissingness

mdl = PMCModel("pmcprg/pmc/models/hmc_in_gauss_k2.toml")
mdl_mnar = mdl.with_missingness(StateMissingness(rates=(0.02, 0.3)))

ICE and SEM estimate the mechanism with the missingness config key ("model", default, carries it fixed; "ignorable" drops it; "state" / "state-markov" estimate it — also pmc estimate --missingness):

from pmcprg.pmc import ice

fitted, trace = ice(mdl, Y_gap, ice_cfg={"missingness": "state"})
print(fitted.missingness)              # StateMissingness(rates=(...))
print(trace.missingness_history[-1])   # its [missingness] table

missingness_lr_test is a likelihood-ratio test of a state-dependent mechanism (H1) against a state-independent one of the same kind (H0), fitted by ICE, with an asymptotic χ² p-value and an optional parametric bootstrap (also pmc missingness-lr-test):

from pmcprg.pmc import missingness_lr_test

test = missingness_lr_test(mdl, Y_gap, alternative="state", n_bootstrap=200, seed=0)
print(test.summary())

Documented limit: when the hidden states are close to i.i.d. (a transition matrix with near-equal rows) the per-state rates are not identified from the observed mixture alone, and the ICE estimates wander instead of converging (pmcprg/pmc/missingness.py, section "Estimation"). On real data (PAMAP2, report/missing_state/pamap2/README.md) the mask helped classification only where each state's dropout rate was homogeneous, and estimating the mechanism from the ignorable fit (init = "model") avoided the worse basins a fresh start could reach: compare segmentation and calibration with the ignorable model before relying on it.

Command-line interface

pmc COMMAND [options]            # or: python -m pmcprg.pmc COMMAND [options]

pmc COMMAND --help lists the options of each command.

Examples:

# Simulate 10 000 samples and save to CSV
pmc simulate \
    --model pmcprg/pmc/models/pmc_gauss_k2.toml \
    --N 10000 --seed 42 --out sim.csv

# Classify (reference labels in column X → prints error rate)
pmc classify \
    --model pmcprg/pmc/models/pmc_gauss_k2.toml \
    --data sim.csv --out cls.csv --ref X

# Unsupervised ICE estimation (--algorithm sem --sem-seed 0 for SEM)
pmc estimate \
    --model pmcprg/pmc/models/pmc_gauss_k2.toml \
    --data sim.csv \
    --candidates "Gauss,Clayton,GH" \
    --max-iter 50 \
    --out fitted.toml

# Image: estimate from the image, then segment it (needs the [image] extra)
pmc estimate-image --model init.toml --image photo.png --fit-margins --out fitted.toml
pmc classify-image --model fitted.toml --image photo.png --out segmentation.png

# State-dependent missingness: estimate the mechanism, then test it
pmc estimate --model pmcprg/pmc/models/hmc_in_gauss_k2.toml --data gaps.csv \
    --missingness state --out fitted.toml
pmc missingness-lr-test --model pmcprg/pmc/models/hmc_in_gauss_k2.toml \
    --data gaps.csv --alternative state --bootstrap 200

# Launch GUI (with optional startup model)
pmc gui pmcprg/pmc/models/pmc_gauss_k2.toml

# Verbose mode (DEBUG to console)
pmc --verbose simulate --model pmcprg/pmc/models/pmc_gauss_k2.toml --N 500

Graphical interface (PyQt6)

pmc gui [MODEL.toml]

The window is divided into a model editor (left) and a result viewer (right):

┌─────────────────────────────────────┬────────────────────────────────┐
│  Model info  (name, variant, K, N)  │  Matplotlib result panel       │
│  ─────────────────────────────────  │  • Simulate: scatter, histo-   │
│  Tabs                               │    gram, successive-pair plot  │
│    Prior  — editable K×K matrix     │  • Classify: MPM labels,       │
│    Margins — state / pair structure │    posterior curves, error map │
│              (dbl-click a cell)     │  • Estimate: ICE log-lik curve │
│    Copulas — K×K grid (dbl-click)   │                                │
│    ICE Config — estimator, criteria,│                                │
│              init, multistart       │  Log panel (INFO+ records)     │
│  ─────────────────────────────────  │                                │
│  [▶ Simulate]   [◈ Classify]        │                                │
│  [⟳ Estimate]   [✓ GoF test]        │                                │
└─────────────────────────────────────┴────────────────────────────────┘

• File menu: open / save TOML models; load data (CSV) or an image; save data, segmentation, results and plots
• Analysis menu: τ confidence intervals, margin adequacy (KS), missingness LR test
• Double-click any Margins or Copulas cell to edit the distribution or copula parameters
• Long-running operations (simulate, classify, estimate) run in a background thread — the GUI stays responsive

Logging

All pmcprg.* loggers route to three sinks:

The file handler captures every logger.exception(...) call with its full traceback, the logger.debug(...) records of the forward–backward weights, ICE iterations and copula fallbacks, and the logger.warning(...) records of numerical near-failures (zero normalisation constants, incompatible observations).

To configure programmatically:

import logging
from pmcprg.pmc.logging_setup import configure, add_widget_handler

# Point the file to a custom path; set console to DEBUG
configure(level=logging.DEBUG, log_file="run.log")

# In a PyQt6 application, once the QApplication exists, mirror the records
# in a QTextEdit:
# add_widget_handler(text_edit, level=logging.INFO)

Diagnostics — multivariate Kolmogorov–Smirnov test

The pmcprg.diagnostics sub-package collects standalone goodness-of-fit utilities that complement the estimation routines of pmcprg.pmc: the multivariate KS test below, parametric_bootstrap (calibration of any post-fit statistic by resampling whole series from the fitted model), Neyman smooth-test components of a fitted margin (neyman_components, neyman_test) and the Vuong / Clarke family comparisons (vuong_test, clarke_test, confidence_set, ice_pair_comparisons).

import numpy as np
from scipy import stats
from pmcprg.diagnostics import mks_1samp, mks_2samp

rng = np.random.default_rng(0)

# 1-sample: is `x` drawn from a 2-D standard normal?
x   = rng.standard_normal(size=(300, 2))
cdf = lambda t: float(stats.norm.cdf(t[0]) * stats.norm.cdf(t[1]))
res = mks_1samp(x, cdf, alpha=0.05)
print(res.statistic, res.critical_value, res.reject)

# 2-sample: do `a` and `b` share the same distribution?
a = rng.standard_normal(size=(300, 2))
b = rng.standard_normal(size=(300, 2)) + 2.0     # shifted mean → reject H0
res = mks_2samp(a, b, alpha=0.05)
assert res.reject

The test extends the classical 1-D Kolmogorov–Smirnov statistic to d > 1 dimensions via Naaman's construction (Statistics & Probability Letters 173, 2021). The default critical value uses the finite-sample union bound (safe but conservative); pass asymptotic=True for the tighter large-N approximation.

Reproducibility reports

• report/ reproduces the experiments of DerrodePieczynski_CSDA2013 §3.2 (supervised PMC, impact of the copula shape), §3.3 (i.i.d. PMM baseline) and §4.3 (unsupervised ICE-based copula selection), with the committed CSV results, LaTeX tables and the report csda2013_reproduction.pdf. Its README lists what the reproduction established, the commands and the provenance of every result.
• report/missing_benchmark/ compares classification and imputation with missing observations under known models: exact marginalisation against plug-in, linear, LOCF and mean fill-ins, on 5 models × 5 patterns × 5 missing rates × 20 replicates.

DerrodePieczynski_SP2016 (GICE, §5.1 setting) is exercised by the test suite on the bundled sp2016_gice_k2.toml model.

Folder structure

awesomePMC/
├── pmcprg/                 import package
│   ├── copulas/            41 copula families (archimedean/, elliptical/,
│   │                       explicit/, extreme_value/), fitting, standard
│   │                       errors, bivariate joint laws
│   ├── diagnostics/        multivariate KS, parametric bootstrap, smooth tests,
│   │                       Vuong/Clarke tests, pseudo-observations
│   ├── missing/            missingness patterns, metrics, CSV missing cells
│   ├── pmc/                model, simulate, inference, gaps (missing data),
│   │   │                   ice, sem, peano (images), cli
│   │   ├── gui/            PyQt6 interface
│   │   └── models/         example TOML models
│   └── tests/              pytest suite
├── examples/               quickstart.ipynb, uci_har_smartphone.ipynb
├── report/                 CSDA 2013 reproduction: scripts, results/, tables/,
│   │                       figures/, LaTeX report and PDF
│   └── missing_benchmark/  missing-data benchmark
├── scripts/                maintenance scripts
├── data/                   local data folders (UCI HAR cache, not versioned)
├── CHANGELOG.md
├── CITATION.cff
├── LICENSE
├── README.md
├── REFERENCES.md
├── pyproject.toml
└── requirements.txt

Citation and references

Please cite both DerrodePieczynski_CSDA2013 and DerrodePieczynski_SP2016 if you use this package in published work. Machine-readable citation metadata is in CITATION.cff (GitHub's Cite this repository button reads it). Every reference the code and the reports rely on is listed with its DOI and the module that uses it in REFERENCES.md.

The pmcprg.pmc sub-package implements the unsupervised classification methods of two papers by S. Derrode and W. Pieczynski. The CSDA paper (DerrodePieczynski_CSDA2013) underpins the PMC model family and ICE-based copula selection; the Signal Processing paper (DerrodePieczynski_SP2016) adds GICE — automatic margin family selection — on top of ICE. pmcprg/pmc/README.md maps every feature to its paper.

• DerrodePieczynski_CSDA2013 — Derrode S., Pieczynski W. Unsupervised data classification using pairwise Markov chains with automatic copulas selection. Computational Statistics & Data Analysis 63 (2013), pp. 81–98. doi:10.1016/j.csda.2013.01.027
• DerrodePieczynski_SP2016 — Derrode S., Pieczynski W. Unsupervised classification using hidden Markov chain with unknown noise copulas and margins. Signal Processing 128 (2016), pp. 8–17. doi:10.1016/j.sigpro.2016.03.008

@ARTICLE{DerrodePieczynski_CSDA2013,
  author    = {S. Derrode and W. Pieczynski},
  title     = {Unsupervised data classification using pairwise {M}arkov chains with automatic copulas selection},
  journal   = {Comput. Stat. Data Anal.},
  volume    = {63},
  pages     = {81--98},
  year      = {2013},
  doi       = {10.1016/j.csda.2013.01.027},
}

@ARTICLE{DerrodePieczynski_SP2016,
  author    = {S. Derrode and W. Pieczynski},
  title     = {Unsupervised classification using hidden {M}arkov chain with unknown noise copulas and margins},
  journal   = {Signal Process.},
  volume    = {128},
  pages     = {8--17},
  year      = {2016},
  doi       = {10.1016/j.sigpro.2016.03.008},
}

The general PMC is defined by DerrodePieczynski_CSDA2013 Eqs. 12–14 (pair margins f_ij):

• DerrodePieczynski_CSDA2013, Eq. 12: f(yₙ, yₙ₊₁ | Xₙ = i, Xₙ₊₁ = j) = f_ij(yₙ) · f_ji(yₙ₊₁) · c_ij(F_ij(yₙ), F_ji(yₙ₊₁))
• DerrodePieczynski_CSDA2013, Eq. 13: p(Xₙ₊₁ = j | Xₙ = i, Yₙ = y) ∝ p(i, j) · f_ij(y)
• DerrodePieczynski_CSDA2013, Eq. 14: p(Yₙ₊₁ | Xₙ = i, Xₙ₊₁ = j, Yₙ = y) = f_ji(yₙ₊₁) · c_ij(F_ij(y), F_ji(yₙ₊₁))

With state margins (f_ij = fᵢ) Eq. 13 no longer depends on y, X is a Markov chain (DerrodePieczynski_CSDA2013 §2.1, Proposition) and the model is a stationary reversible HMC-DN — see Model variants above.

Methodological building blocks used by the package:

• Huard D., Évin G., Favre A.-C. Bayesian copula selection. Computational Statistics & Data Analysis 51(2) (2006), pp. 809–822 — the huard selection criterion (DerrodePieczynski_CSDA2013, Eq. 20). doi:10.1016/j.csda.2005.08.010
• Naaman M. On the tight constant in the multivariate Dvoretzky–Kiefer–Wolfowitz inequality. Statistics & Probability Letters 173 (2021), 109088 — the multivariate KS test in pmcprg.diagnostics. doi:10.1016/j.spl.2021.109088
• Anguita D. et al. A Public Domain Dataset for Human Activity Recognition Using Smartphones. ESANN 2013, pp. 437–442 (the paper has no DOI) — the UCI HAR benchmark used in examples/. Dataset DOI (UCI, Reyes-Ortiz & Anguita 2013): doi:10.24432/C54S4K

Related work (a selection; the references the code and the reports rely on are listed in REFERENCES.md):

• Gorynin I., Gangloff H., Monfrini E., Pieczynski W. Assessing the segmentation performance of pairwise and triplet Markov models. Signal Processing 145 (2018), pp. 183–192 — quantifies the PMM/TMM gain over classical HMMs. doi:10.1016/j.sigpro.2017.12.006
• Gangloff H., Morales K., Petetin Y. Deep parameterizations of pairwise and triplet Markov models for unsupervised classification of sequential data. Computational Statistics & Data Analysis 180 (2023), 107663 — the deep-learning branch of the PMC/TMC lineage. doi:10.1016/j.csda.2022.107663
• Zimmerman R., Craiu R.V., Leos-Barajas V. Copula Modelling of Serially Correlated Multivariate Data with Hidden Structures. JASA 119(548) (2024), pp. 2598–2609 — the closest copula-HMM methodology outside the pairwise family. doi:10.1080/01621459.2023.2263202
• Nasri B.R., Rémillard B.N., Thioub M.Y. Goodness-of-fit for regime-switching copula models. Canadian Journal of Statistics 48(1) (2020), pp. 79–96 — regime-switching copulas and the R package HMMcopula. doi:10.1002/cjs.11534
• Grønneberg S., Hjort N.L. The Copula Information Criteria. Scandinavian Journal of Statistics 41(2) (2014), pp. 436–459 — why AIC on rank pseudo-observations needs correction (xv-CIC). doi:10.1111/sjos.12042

License

MIT — see LICENSE.