Calibrated CSSR-family tests, automatic Lmax, and reconstruction diagnostics

docs/generators/epsilon_inference.rst

@@ -41,7 +41,8 @@ CSSR :cite:`Shalizi2002` starts from an IID model and grows causal states in thr
 1. **Initialize** — one state for the empty history.
 2. **Homogenize** — extend each suffix one symbol into the past, up to ``Lmax``; a
    child suffix whose next-symbol distribution differs significantly from its
-   state's (G-test, :math:`\chi^2`, or total-variation threshold) moves to the best
+   state's (G-test, :math:`\chi^2`, Monte Carlo exact G-test, or total-variation
+   threshold) moves to the best
    matching state, or starts a new one. States keep suffixes of every length.
 3. **Determinize** — drop transient states, then split states until each state and
    symbol lead to a single successor, then keep the most-visited recurrent class.
@@ -59,8 +60,45 @@ split states by chance; lowering ``alpha`` counters this. A process that is not
 exactly synchronizable has no finite-``Lmax`` reconstruction, and CSSR returns
 extra states.
 
+Choosing ``Lmax`` and calibrating the tests
+-------------------------------------------
+
+``Lmax="auto"`` sets ``Lmax`` to :func:`suggest_lmax`, the Markov order estimated
+by :func:`dit.inference.select_markov_order`. Its default method tests order
+:math:`n` against :math:`n + 1` with surrogates that preserve the observed
+:math:`(n + 1)`-gram counts exactly, so the test holds its nominal size at any
+sample length, where the asymptotic chi-squared test is badly anti-conservative
+:cite:`Pethel2014`. For a Markov source this is its order, which is the
+synchronization length CSSR needs. A strictly sofic source such as the even
+process has infinite Markov order, so the suggestion grows with the sample: read
+it as the longest history the data supports, not as a synchronization length.
+
+The morph tests also rely on the chi-squared limit, which fails for the sparse
+counts of long suffixes. ``test="exact"`` compares the G statistic with tables
+drawn uniformly given the observed margins whenever an expected count is below 5
+(seeded from the table, so reconstruction stays deterministic).
+``correction="bonferroni"`` divides ``alpha`` by the number of suffixes eligible
+for testing, bounding the chance of any spurious split. Because CSSR chooses each
+test in light of earlier outcomes, false-discovery-rate step-up procedures do not
+apply directly.
+
+.. code-block:: python
+
+   inferred = EpsilonMachine.from_sequence(
+       observations, method="cssr", Lmax="auto", test="exact", correction="bonferroni"
+   )
+
 .. autofunction:: cssr
 
+.. autofunction:: suggest_lmax
+
+After reconstruction, check the result with
+:func:`~sofic.inference.diagnostics.goodness_of_fit` and
+:func:`~sofic.inference.diagnostics.structure_stability`
+(see :doc:`../inference/diagnostics`).
+
+.. autofunction:: morphs_differ
+
 Subtree merging
 ===============
 
@@ -70,6 +108,10 @@ to a unifilar presentation.  With ``delta=0``, two morphs are equivalent unless
 G-test at significance 0.01 tells them apart, a tolerance that scales with the
 sample. Transitions follow the same successor rule as CSSR.
 
+``subtree_merge`` accepts ``alpha``, ``test`` (including ``"exact"``) and
+``correction="bonferroni"``, which divides ``alpha`` over the history pairs
+compared, as well as ``L="auto"``.
+
 .. autofunction:: subtree_merge
 
 Spectral reconstruction

docs/generators/epsilon_transducer_inference.rst

@@ -14,6 +14,11 @@ Causal states are equivalence classes of joint ``(input, output)`` pasts that
 induce the same conditional next-output law ``P(y | history, x)`` for every input
 symbol ``x``. Rare histories inherit their parent's state (controlled by
 ``min_count``); the split decision uses a G-test at significance ``alpha``.
+As in :func:`~sofic.generators.epsilon_inference.cssr`, ``test="exact"`` uses a
+Monte Carlo exact G-test for tables with small expected counts, and
+``correction="bonferroni"`` divides ``alpha`` by the number of
+(history, input symbol) tests. ``Lmax="auto"`` sets the depth from the Markov
+order of the joint ``(input, output)`` sequence :cite:`Pethel2014`.
 
 .. ipython::
 

docs/generators/hmm_inference.rst

@@ -50,6 +50,14 @@ transition-graph topology fixed:
 
    In [13]: fitted, loglik_trace = baum_welch(golden_mean(0.6), data)
 
+EM stops at a local maximum of the likelihood. ``n_restarts`` reruns it from
+random edge laws on the same topology and keeps the best fit;
+``return_restarts=True`` also returns each run's final log-likelihood:
+
+.. code-block:: python
+
+   fitted, trace, finals = baum_welch(start, data, n_restarts=10, rng=0, return_restarts=True)
+
 The score and observed information quantify the log-likelihood gradient and
 parameter uncertainty at the current parameters:
 

docs/generators/stack_inference.rst

@@ -48,6 +48,12 @@ can follow is decided by the stack top through matched call-return pairs, not by
 the finite control. Return edges are matched only to calls observed to close
 them.
 
+``stack_cssr`` accepts the same calibration options as
+:func:`~sofic.generators.epsilon_inference.cssr`: ``test="exact"``,
+``correction="bonferroni"`` (over eligible configurations), and ``Lmax="auto"``.
+Stack processes generally have infinite Markov order, so the automatic depth is
+a lower bound on the suffix length the data support.
+
 API
 ===
 

docs/inference/diagnostics.rst

@@ -0,0 +1,72 @@
+.. diagnostics.rst
+.. py:module:: sofic.inference.diagnostics
+
+*****************************
+Reconstruction diagnostics
+*****************************
+
+A reconstruction algorithm always returns *some* machine. These tools ask
+whether that machine fits the data, and whether its structure is supported by the
+data or produced by one particular sample and one setting of the tuning
+parameters.
+
+Goodness of fit
+===============
+
+:func:`goodness_of_fit` is a parametric bootstrap :cite:`Efron1993`. It
+simulates sequences as long as the data from the fitted machine and compares a
+length-``L`` word statistic of the data with its distribution over the
+simulations. Two statistics are available:
+
+* ``"g"`` — the G statistic of the observed word counts against the machine's
+  stationary word probabilities;
+* ``"entropy_rate"`` — the gap between the plug-in conditional entropy and the
+  machine's.
+
+Because the null distribution is simulated, overlapping windows need no
+correction. A small p-value means the machine misses structure. For CSSR that
+usually means ``Lmax`` is shorter than the source's synchronization length.
+
+.. code-block:: python
+
+   from sofic.generators.epsilon_inference import cssr
+   from sofic.inference.diagnostics import goodness_of_fit
+
+   machine = cssr(data, Lmax=1)
+   goodness_of_fit(machine, data, L=6).pvalue   # small for the even process
+   machine = cssr(data, Lmax=4)
+   goodness_of_fit(machine, data, L=6).pvalue   # large
+
+Observed words the machine forbids are listed in ``forbidden_words``.
+
+Structural stability
+====================
+
+:func:`structure_stability` reconstructs from many resamples and counts how
+often each topology (compared up to isomorphism by :func:`topology_key`)
+reappears. The default ``resample="subsample"`` uses random contiguous segments
+:cite:`Politis1999`, which contain no artificial junctions. The stationary
+bootstrap (``resample="block"``, :cite:`Politis1994`) joins blocks, which creates
+words the source never emits and can add spurious states. For example, on
+even-process data it returns 6–12-state machines where subsampling returns the
+true 2 states.
+
+:func:`reconstruction_sweep` reconstructs over a grid of ``alpha`` and ``Lmax``.
+A structure that persists over a range of settings is better supported than one
+that appears at a single setting.
+
+API
+===
+
+.. autofunction:: goodness_of_fit
+
+.. autoclass:: GoodnessOfFit
+
+.. autofunction:: structure_stability
+
+.. autoclass:: StructureStability
+   :members: reference_fraction, modal_topology, n_resamples
+
+.. autofunction:: reconstruction_sweep
+
+.. autofunction:: topology_key

docs/inference/inference.rst

@@ -44,6 +44,7 @@ The historical names ``InferMC`` and ``InferEM`` are retained as aliases for
    epsilon
    spectral
    model_selection
+   diagnostics
    hdp_hmm
    stack_hmm
    pymc

docs/inference/model_selection.rst

@@ -46,9 +46,18 @@ Cross-validation and WAIC
 
    cross_validated_log_likelihood(fit, data, folds=5)   # held-out log score (higher is better)
 
+   # Drop 20 symbols next to each held-out block, and keep folds finite when the
+   # fitted model forbids a held-out transition:
+   cross_validated_log_likelihood(fit, data, folds=5, gap=20, smoothing=1e-3)
+
    posterior = EpsilonMachinePosterior(golden_mean(0.3), data)
    waic_epsilon_machine(posterior, [data], n_samples=200)
 
+Contiguous blocks of one sequence are dependent, so without a ``gap`` the
+held-out score is optimistic :cite:`Burman1994`. ``smoothing`` mixes each held-out
+prediction with the uniform distribution, so that one forbidden transition no
+longer makes a whole fold ``-inf``.
+
 Ranking candidate topologies
 ============================
 

docs/references.bib

@@ -1060,3 +1060,53 @@ @misc{HardDropGameBoy
   url          = {https://harddrop.com/wiki/Tetris_(Game_Boy)},
   note         = {Reverse-engineered Game Boy bitwise-OR randomizer},
 }
+
+@article{Pethel2014,
+  author  = {Pethel, Shawn D. and Hahs, Daniel W.},
+  title   = {Exact significance test for {Markov} order},
+  journal = {Physica D: Nonlinear Phenomena},
+  volume  = {269},
+  pages   = {42--47},
+  year    = {2014},
+  doi     = {10.1016/j.physd.2013.11.014},
+}
+
+@book{Efron1993,
+  author    = {Efron, Bradley and Tibshirani, Robert J.},
+  title     = {An Introduction to the Bootstrap},
+  publisher = {Chapman \& Hall},
+  address   = {New York},
+  year      = {1993},
+  doi       = {10.1007/978-1-4899-4541-9},
+}
+
+@book{Politis1999,
+  author    = {Politis, Dimitris N. and Romano, Joseph P. and Wolf, Michael},
+  title     = {Subsampling},
+  publisher = {Springer},
+  address   = {New York},
+  year      = {1999},
+  doi       = {10.1007/978-1-4612-1554-7},
+}
+
+@article{Politis1994,
+  author  = {Politis, Dimitris N. and Romano, Joseph P.},
+  title   = {The stationary bootstrap},
+  journal = {Journal of the American Statistical Association},
+  volume  = {89},
+  number  = {428},
+  pages   = {1303--1313},
+  year    = {1994},
+  doi     = {10.1080/01621459.1994.10476870},
+}
+
+@article{Burman1994,
+  author  = {Burman, Prabir and Chow, Edmond and Nolan, Deborah},
+  title   = {A cross-validatory method for dependent data},
+  journal = {Biometrika},
+  volume  = {81},
+  number  = {2},
+  pages   = {351--358},
+  year    = {1994},
+  doi     = {10.2307/2336965},
+}

sofic/generators/__init__.py

@@ -16,7 +16,7 @@
     synergistic_information_flow,
     transfer_entropy,
 )
-from sofic.generators.epsilon_inference import cssr, spectral, subtree_merge
+from sofic.generators.epsilon_inference import cssr, spectral, subtree_merge, suggest_lmax
 from sofic.generators.epsilon_machine import EpsilonMachine
 from sofic.generators.epsilon_transducer import EpsilonTransducer
 from sofic.generators.lumping import LumpabilityError, is_lumpable, lump, normalize_partition
@@ -82,6 +82,7 @@
     "normalize_partition",
     "spectral",
     "subtree_merge",
+    "suggest_lmax",
     "fit_stack_hmm_mle",
     "learn_stack_hmm_papni",
     "stack_cssr",

sofic/generators/base.py

@@ -147,8 +147,13 @@ def baum_welch(
         max_iter: int = 100,
         tol: float = 1e-6,
         estimate_initial: bool = True,
+        n_restarts: int = 1,
+        rng: np.random.Generator | int | None = None,
     ) -> tuple[MealyHMM, list[float]]:
-        """Fit parameters by Baum-Welch EM, returning ``(fitted_model, loglik_trace)``."""
+        """Fit parameters by Baum-Welch EM, returning ``(fitted_model, loglik_trace)``.
+
+        ``n_restarts`` and ``rng`` are as in :func:`~sofic.generators.hmm_inference.baum_welch`.
+        """
         from sofic.generators.hmm_inference import baum_welch
 
         return baum_welch(
@@ -157,6 +162,8 @@ def baum_welch(
             max_iter=max_iter,
             tol=tol,
             estimate_initial=estimate_initial,
+            n_restarts=n_restarts,
+            rng=rng,
         )
 
     def score(self, observations: Sequence[Any]) -> dict[tuple[Hashable, Any, Hashable], float]: