diff --git a/docs/generators/epsilon_inference.rst b/docs/generators/epsilon_inference.rst
index b263141..a9494fb 100644
--- a/docs/generators/epsilon_inference.rst
+++ b/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
diff --git a/docs/generators/epsilon_transducer_inference.rst b/docs/generators/epsilon_transducer_inference.rst
index c58097a..707be51 100644
--- a/docs/generators/epsilon_transducer_inference.rst
+++ b/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::
 
diff --git a/docs/generators/hmm_inference.rst b/docs/generators/hmm_inference.rst
index e91ddae..c430f29 100644
--- a/docs/generators/hmm_inference.rst
+++ b/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:
 
diff --git a/docs/generators/stack_inference.rst b/docs/generators/stack_inference.rst
index bc9395d..af85695 100644
--- a/docs/generators/stack_inference.rst
+++ b/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
 ===
 
diff --git a/docs/inference/diagnostics.rst b/docs/inference/diagnostics.rst
new file mode 100644
index 0000000..a674dc9
--- /dev/null
+++ b/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
diff --git a/docs/inference/inference.rst b/docs/inference/inference.rst
index b228b4f..e91692c 100644
--- a/docs/inference/inference.rst
+++ b/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
diff --git a/docs/inference/model_selection.rst b/docs/inference/model_selection.rst
index 2d9f8ee..6e66dd0 100644
--- a/docs/inference/model_selection.rst
+++ b/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
 ============================
 
diff --git a/docs/references.bib b/docs/references.bib
index 2d6a9ed..afbdd2f 100644
--- a/docs/references.bib
+++ b/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},
+}
diff --git a/sofic/generators/__init__.py b/sofic/generators/__init__.py
index 35a7388..d0866e3 100644
--- a/sofic/generators/__init__.py
+++ b/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",
diff --git a/sofic/generators/base.py b/sofic/generators/base.py
index 123108a..276ce6b 100644
--- a/sofic/generators/base.py
+++ b/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]:
diff --git a/sofic/generators/epsilon_inference.py b/sofic/generators/epsilon_inference.py
index 37adf36..488be70 100644
--- a/sofic/generators/epsilon_inference.py
+++ b/sofic/generators/epsilon_inference.py
@@ -8,6 +8,7 @@
 
 from __future__ import annotations
 
+import zlib
 from collections import Counter, defaultdict
 from collections.abc import Callable, Iterable, Sequence
 from dataclasses import dataclass, field
@@ -23,6 +24,9 @@
 
 History = tuple[Any, ...]
 
+#: Morph-equality tests: G-test, chi-squared, total-variation threshold, or Monte Carlo exact G-test.
+MorphTest = Literal["g", "chi2", "tv", "exact"]
+
 
 @dataclass
 class SuffixCounts:
@@ -175,16 +179,86 @@ def _chi2_critical(alpha: float, dof: int) -> float:
     return float(stats.chi2.ppf(1.0 - alpha, dof))
 
 
+#: Monte Carlo tables drawn per ``"exact"`` morph test.
+_EXACT_DRAWS = 999
+
+#: Smallest expected count at which the ``"exact"`` test trusts the chi-squared limit.
+_EXACT_MIN_EXPECTED = 5.0
+
+
+def _exact_g_pvalue(table: np.ndarray) -> float:
+    """Monte Carlo p-value of the G statistic among tables with the same margins.
+
+    Tables are drawn uniformly given both margins (``scipy.stats.random_table``),
+    the exact null of equal morphs. The generator is seeded from the table itself,
+    so reconstruction stays deterministic.
+    """
+    counts = table.astype(np.int64)
+    rows, cols = counts.sum(axis=1), counts.sum(axis=0)
+    expected = np.outer(rows, cols) / counts.sum()
+
+    def g(observed: np.ndarray) -> np.ndarray:
+        with np.errstate(divide="ignore", invalid="ignore"):
+            terms = np.where(observed > 0, observed * np.log(observed / expected), 0.0)
+        return 2.0 * terms.sum(axis=(-2, -1))
+
+    rng = np.random.default_rng(zlib.crc32(counts.tobytes()))
+    draws = stats.random_table(rows, cols, seed=rng).rvs(size=_EXACT_DRAWS)
+    observed = g(counts.astype(float))
+    extreme = np.sum(g(draws.astype(float)) >= observed - 1e-9 * max(1.0, observed))
+    return float((1 + extreme) / (1 + _EXACT_DRAWS))
+
+
+def _exact_significant(table: np.ndarray, alpha: float) -> bool:
+    """The ``"exact"`` decision for a contingency table of next-symbol counts.
+
+    Uses the Monte Carlo exact G-test when an expected count is below
+    ``_EXACT_MIN_EXPECTED`` and the asymptotic G-test otherwise.
+    """
+    if np.any(table.sum(axis=1) == 0):
+        return False
+    expected = np.outer(table.sum(axis=1), table.sum(axis=0)) / table.sum()
+    if expected.min() < _EXACT_MIN_EXPECTED:
+        return _exact_g_pvalue(table) < alpha
+    statistic = _g_statistic(table)
+    if statistic is None or not np.isfinite(statistic):
+        return False
+    return statistic > _chi2_critical(alpha, max(1, table.shape[1] - 1))
+
+
+def _bonferroni_alpha(
+    counts: SuffixCounts,
+    alpha: float,
+    *,
+    max_length: int,
+    min_count: int,
+    suffix_length: Callable[[History], int] = len,
+) -> float:
+    """``alpha`` divided by the number of suffixes eligible for a split test."""
+    eligible = sum(
+        1
+        for history, following in counts.next_counts.items()
+        if 0 < suffix_length(history) <= max_length and sum(following.values()) >= max(1, min_count)
+    )
+    return alpha / max(1, eligible)
+
+
 def morphs_differ(
     counts: SuffixCounts,
     left_histories: set[History],
     right_histories: set[History],
     *,
     alpha: float = 0.05,
-    test: Literal["g", "chi2", "tv"] = "g",
+    test: MorphTest = "g",
     delta: float = 0.0,
 ) -> bool:
-    """Return whether two history sets have significantly different morphs."""
+    """Return whether two history sets have significantly different morphs.
+
+    ``"g"`` and ``"chi2"`` use the chi-squared limit, which is unreliable when
+    expected counts are small. ``"exact"`` instead compares the G statistic with
+    tables drawn uniformly given the observed margins whenever an expected count is
+    below 5, and uses the G-test otherwise.
+    """
     if test == "tv":
         left = counts.state_morph(left_histories)
         right = counts.state_morph(right_histories)
@@ -194,6 +268,8 @@ def morphs_differ(
     table = _contingency_rows(counts, left_histories, right_histories)
     if table is None:
         return False
+    if test == "exact":
+        return _exact_significant(table, alpha)
     if test == "g":
         statistic = _g_statistic(table)
         if statistic is None or not np.isfinite(statistic):
@@ -213,7 +289,7 @@ def morph_test_score(
     left_histories: set[History],
     right_histories: set[History],
     *,
-    test: Literal["g", "chi2", "tv"] = "g",
+    test: MorphTest = "g",
 ) -> float:
     """Score for matching morphs (lower is more similar)."""
     if test == "tv":
@@ -224,7 +300,7 @@ def morph_test_score(
     if table is None:
         return 0.0
     try:
-        if test == "g":
+        if test in ("g", "exact"):
             statistic = _g_statistic(table)
             return statistic if statistic is not None and np.isfinite(statistic) else 0.0
         statistic, _p, _dof, _expected = stats.chi2_contingency(table)
@@ -249,7 +325,7 @@ def _cssr_homogenize(
     *,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2", "tv"],
+    test: MorphTest,
     successor_fn: Callable[[History, Any], History] = _grow_history,
 ) -> tuple[dict[int, set[History]], dict[History, int]]:
     """Return state id -> histories and history -> state id."""
@@ -361,7 +437,7 @@ def _merge_similar_states(
     counts: SuffixCounts,
     *,
     alpha: float,
-    test: Literal["g", "chi2", "tv"],
+    test: MorphTest,
 ) -> dict[int, set[History]]:
     """Merge inferred states whose pooled morphs are statistically indistinguishable.
 
@@ -465,7 +541,7 @@ def _suffix_homogenize(
     *,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2", "tv"],
+    test: MorphTest,
     min_count: int = 1,
 ) -> list[set[History]]:
     """CSSR homogenization: grow suffixes one symbol into the past, up to length ``Lmax``.
@@ -517,7 +593,7 @@ def _suffix_edges(
     *,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2", "tv"],
+    test: MorphTest,
     resolve: bool = True,
 ) -> dict[int, dict[Any, dict[int, set[History]]]]:
     """Successor states of each alive state, by symbol, with the suffixes that lead there.
@@ -587,7 +663,7 @@ def _suffix_determinize(
     *,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2", "tv"],
+    test: MorphTest,
     resolve: bool = True,
 ) -> tuple[list[set[History]], set[int]]:
     """Split alive states until each (state, symbol) pair has a single alive successor.
@@ -624,7 +700,7 @@ def _suffix_machine(
     *,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2", "tv"],
+    test: MorphTest,
     resolve: bool = True,
 ) -> EpsilonMachine:
     """Build the ε-machine on the most-visited recurrent class of the alive states."""
@@ -691,7 +767,7 @@ def _suffix_reconstruct(
     *,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2", "tv"],
+    test: MorphTest,
 ) -> EpsilonMachine:
     """Prune transient states, determinize, and build the machine from homogeneous ``states``."""
 
@@ -710,14 +786,69 @@ def reconstruct(resolve: bool) -> EpsilonMachine:
     return machine
 
 
+def suggest_lmax(
+    sequence: Sequence[Any],
+    *,
+    alpha: float = 0.01,
+    max_order: int | None = None,
+    method: Literal["exact", "chi2", "aic", "bic"] = "exact",
+    n_surrogates: int = 999,
+    seed: int = 0,
+) -> int:
+    """A data-driven ``Lmax`` for :func:`cssr`: the estimated Markov order, at least 1.
+
+    Orders ``0, 1, ...`` are tested against the next order with
+    :func:`dit.inference.select_markov_order`. The default ``"exact"`` method
+    compares the conditional block entropy against surrogates that preserve the
+    observed ``(order + 1)``-gram counts, which is valid at any sample size, unlike
+    the asymptotic chi-squared test :cite:`Pethel2014`.
+
+    Parameters
+    ----------
+    sequence
+        Observed symbols.
+    alpha
+        Significance level of each order test.
+    max_order
+        Largest order considered; by default the largest ``L`` whose
+        ``(L + 1)``-words are seen about 5 times each on average, at most 10.
+    method
+        ``"exact"`` or ``"chi2"`` (sequential tests), or ``"aic"`` / ``"bic"``.
+    n_surrogates
+        Surrogates per test for ``"exact"``.
+    seed
+        Seed for the surrogates, so the suggestion is reproducible.
+
+    Notes
+    -----
+    For a Markov source this recovers 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 keeps growing with the sample; treat it
+    as a lower bound on the history length the data can support, not as the source's
+    synchronization length.
+    """
+    import dit.inference
+
+    select_markov_order = getattr(dit.inference, "select_markov_order", None)
+    if select_markov_order is None:  # pragma: no cover - depends on the installed dit
+        raise ImportError("suggest_lmax requires a dit release with dit.inference.select_markov_order")
+    seq = [repr(symbol) for symbol in sequence]
+    if max_order is None:
+        k = max(2, len(set(seq)))
+        max_order = max(1, min(10, int(np.log(max(len(seq), 1) / 5) / np.log(k)) - 1))
+    order = select_markov_order(seq, max_order, method=method, alpha=alpha, n_surrogates=n_surrogates, prng=seed)
+    return max(1, int(order))
+
+
 def cssr(
     sequence: Sequence[Any],
     *,
     alphabet: Sequence[Any] | None = None,
-    Lmax: int | None = None,
+    Lmax: int | Literal["auto"] | None = None,
     alpha: float = 0.01,
-    test: Literal["g", "chi2", "tv"] = "g",
+    test: MorphTest = "g",
     min_count: int = 5,
+    correction: Literal["bonferroni"] | None = None,
 ) -> EpsilonMachine:
     """Reconstruct an ε-machine by Causal-State Splitting Reconstruction :cite:`Shalizi2004`.
 
@@ -736,15 +867,23 @@ def cssr(
         Longest suffix considered; by default a third of ``log_k len(sequence)``,
         between 1 and 10. It should be at least the synchronization length of the
         source (for a Markov source, its order). Larger values run many more
-        significance tests, and some of them split states by chance.
+        significance tests, and some of them split states by chance. ``"auto"``
+        uses :func:`suggest_lmax`, the Markov order estimated by exact tests.
     alpha
         Significance level of each morph-equality test. The worked example of
         :cite:`Shalizi2002` uses 0.01; smaller values guard against spurious states
         when ``Lmax`` is large.
     test
-        ``"g"`` (G-test), ``"chi2"``, or ``"tv"`` (total-variation threshold).
+        ``"g"`` (G-test), ``"chi2"``, ``"tv"`` (total-variation threshold), or
+        ``"exact"`` (Monte Carlo exact G-test when expected counts are small; see
+        :func:`morphs_differ`).
     min_count
         Suffixes seen fewer than this many times are not tested or placed in a state.
+    correction
+        ``"bonferroni"`` divides ``alpha`` by the number of suffixes eligible for
+        testing, bounding the chance of any spurious split. CSSR decides each test
+        in light of earlier ones, so step-up procedures that control the false
+        discovery rate (Benjamini–Hochberg) do not apply directly.
 
     Notes
     -----
@@ -757,10 +896,17 @@ def cssr(
     if len(seq) < 2:
         raise ValueError("sequence must contain at least two symbols")
     alphabet_size = len(set(seq)) if alphabet is None else len(tuple(alphabet))
-    max_length = Lmax if Lmax is not None else _cssr_default_lmax(len(seq), alphabet_size)
+    if Lmax == "auto":
+        max_length = suggest_lmax(seq, alpha=alpha)
+    else:
+        max_length = Lmax if Lmax is not None else _cssr_default_lmax(len(seq), alphabet_size)
     if max_length < 0:
         raise ValueError("Lmax must be non-negative")
     counts = SuffixCounts.from_sequence(seq, alphabet=alphabet, max_length=max_length + 1)
+    if correction == "bonferroni":
+        alpha = _bonferroni_alpha(counts, alpha, max_length=max_length, min_count=min_count)
+    elif correction is not None:
+        raise ValueError(f"unknown correction {correction!r}")
 
     homogeneous = _suffix_homogenize(counts, Lmax=max_length, alpha=alpha, test=test, min_count=min_count)
     return _suffix_reconstruct(homogeneous, counts, seq, Lmax=max_length, alpha=alpha, test=test)
@@ -788,10 +934,12 @@ def _morphs_equivalent(
     right: History,
     *,
     delta: float,
+    alpha: float = _SUBTREE_ALPHA,
+    test: MorphTest = "g",
 ) -> bool:
     if delta > 0.0:
         return _morph_distance(counts, left, right, delta=delta) <= delta
-    return not morphs_differ(counts, {left}, {right}, alpha=_SUBTREE_ALPHA, test="g")
+    return not morphs_differ(counts, {left}, {right}, alpha=alpha, test=test)
 
 
 def _cluster_histories_by_morph(
@@ -799,6 +947,8 @@ def _cluster_histories_by_morph(
     histories: set[History],
     *,
     delta: float,
+    alpha: float = _SUBTREE_ALPHA,
+    test: MorphTest = "g",
 ) -> dict[int, set[History]]:
     parent: dict[History, History] = {history: history for history in histories}
 
@@ -818,7 +968,7 @@ def union(left: History, right: History) -> None:
     history_list = sorted(histories)
     for index, left in enumerate(history_list):
         for right in history_list[index + 1 :]:
-            if _morphs_equivalent(counts, left, right, delta=delta):
+            if _morphs_equivalent(counts, left, right, delta=delta, alpha=alpha, test=test):
                 union(left, right)
 
     clusters: dict[History, set[History]] = defaultdict(set)
@@ -834,17 +984,27 @@ def union(left: History, right: History) -> None:
 def subtree_merge(
     sequence: Sequence[Any],
     *,
-    L: int,
+    L: int | Literal["auto"],
     delta: float = 0.0,
     alphabet: Sequence[Any] | None = None,
+    alpha: float = _SUBTREE_ALPHA,
+    test: MorphTest = "g",
+    correction: Literal["bonferroni"] | None = None,
 ) -> EpsilonMachine:
     """Reconstruct an ε-machine by merging depth-``L`` subtrees (Crutchfield--Young).
 
     Histories up to length ``L`` are clustered by next-symbol distribution: within
-    total-variation distance ``delta``, or, when ``delta = 0``, unless a G-test at
-    significance 0.01 tells them apart. The clusters are then determinized as in
-    :func:`cssr`.
+    total-variation distance ``delta``, or, when ``delta = 0``, unless a morph test
+    (``test``, at level ``alpha``) tells them apart. The clusters are then
+    determinized as in :func:`cssr`.
+
+    ``L="auto"`` uses :func:`suggest_lmax`. ``correction="bonferroni"`` divides
+    ``alpha`` by the number of history pairs compared, so that no pair is split
+    apart by chance; since a rejected test *separates* histories, this makes the
+    reconstruction more conservative (fewer states).
     """
+    if L == "auto":
+        L = suggest_lmax(sequence, alpha=alpha)
     if L < 0:
         raise ValueError("L must be non-negative")
     seq = tuple(sequence)
@@ -854,9 +1014,13 @@ def subtree_merge(
 
     histories = {history for history in counts.history_counts if len(history) <= L}
     histories.add(())
-    states = list(_cluster_histories_by_morph(counts, histories, delta=delta).values())
+    if correction == "bonferroni":
+        alpha /= max(1, len(histories) * (len(histories) - 1) // 2)
+    elif correction is not None:
+        raise ValueError(f"unknown correction {correction!r}")
+    states = list(_cluster_histories_by_morph(counts, histories, delta=delta, alpha=alpha, test=test).values())
 
-    return _suffix_reconstruct(states, counts, seq, Lmax=L, alpha=_SUBTREE_ALPHA, test="g")
+    return _suffix_reconstruct(states, counts, seq, Lmax=L, alpha=alpha, test=test)
 
 
 def spectral(
diff --git a/sofic/generators/epsilon_transducer_inference.py b/sofic/generators/epsilon_transducer_inference.py
index b782b8d..2466b1b 100644
--- a/sofic/generators/epsilon_transducer_inference.py
+++ b/sofic/generators/epsilon_transducer_inference.py
@@ -126,7 +126,7 @@ def aggregates_differ(
     input_alphabet: tuple[Any, ...],
     output_alphabet: tuple[Any, ...],
     alpha: float,
-    test: Literal["g", "chi2"] = "g",
+    test: Literal["g", "chi2", "exact"] = "g",
 ) -> bool:
     """Return whether two aggregated morphs differ on ``P(output | ., input)`` for some input."""
     for input_symbol in input_alphabet:
@@ -168,7 +168,11 @@ def _aggregate_score(
     return total
 
 
-def _table_significant(table: np.ndarray, *, alpha: float, test: Literal["g", "chi2"]) -> bool:
+def _table_significant(table: np.ndarray, *, alpha: float, test: Literal["g", "chi2", "exact"]) -> bool:
+    if test == "exact":
+        from sofic.generators.epsilon_inference import _exact_significant
+
+        return _exact_significant(table, alpha)
     try:
         if test == "g":
             with np.errstate(invalid="ignore", divide="ignore"):
@@ -201,7 +205,7 @@ def _homogenize(
     *,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2"],
+    test: Literal["g", "chi2", "exact"],
     min_count: int,
 ) -> list[set[JointHistory]]:
     """transCSSR homogenization: grow joint suffixes one ``(input, output)`` pair into the past.
@@ -268,7 +272,7 @@ def _edges(
     *,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2"],
+    test: Literal["g", "chi2", "exact"],
 ) -> dict[int, dict[tuple[Any, Any], dict[int, set[JointHistory]]]]:
     """Successor states by ``(input, output)`` pair, with the histories that lead there.
 
@@ -349,7 +353,7 @@ def _determinize(
     *,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2"],
+    test: Literal["g", "chi2", "exact"],
 ) -> tuple[list[set[JointHistory]], set[int]]:
     """Split alive states until each ``(input, output)`` pair has one alive successor."""
     states = [set(h) for h in states]
@@ -384,7 +388,7 @@ def _build_transducer(
     *,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2"],
+    test: Literal["g", "chi2", "exact"],
 ) -> EpsilonTransducer:
     edges = _edges(states, counts, alive, Lmax=Lmax, alpha=alpha, test=test)
     history_to_state = {h: index for index in alive for h in states[index]}
@@ -462,10 +466,11 @@ def transcssr(
     *,
     input_alphabet: Sequence[Any] | None = None,
     output_alphabet: Sequence[Any] | None = None,
-    Lmax: int | None = None,
+    Lmax: int | Literal["auto"] | None = None,
     alpha: float = 0.001,
-    test: Literal["g", "chi2"] = "g",
+    test: Literal["g", "chi2", "exact"] = "g",
     min_count: int = 5,
+    correction: Literal["bonferroni"] | None = None,
 ) -> EpsilonTransducer:
     """Reconstruct an ε-transducer from paired input/output sequences (transCSSR).
 
@@ -473,6 +478,13 @@ def transcssr(
     decision; the transCSSR/CSSR default of ``0.001`` favors fewer, more robust
     states. ``Lmax`` bounds the joint-history depth and ``min_count`` the minimum
     occurrences before a history is eligible to seed a new state.
+
+    ``Lmax="auto"`` applies :func:`~sofic.generators.epsilon_inference.suggest_lmax`
+    to the joint ``(input, output)`` sequence. ``test="exact"`` uses the Monte
+    Carlo exact G-test when expected counts are small (see
+    :func:`~sofic.generators.epsilon_inference.morphs_differ`), and
+    ``correction="bonferroni"`` divides ``alpha`` by the number of
+    (history, input symbol) tests that can split a state.
     """
     xs = tuple(inputs)
     ys = tuple(outputs)
@@ -485,7 +497,12 @@ def transcssr(
         if input_alphabet is None or output_alphabet is None
         else len(tuple(input_alphabet)) * len(tuple(output_alphabet))
     )
-    max_length = Lmax if Lmax is not None else _default_lmax(len(xs), joint_alphabet_size, min_count)
+    if Lmax == "auto":
+        from sofic.generators.epsilon_inference import suggest_lmax
+
+        max_length = suggest_lmax(list(zip(xs, ys, strict=True)), alpha=alpha)
+    else:
+        max_length = Lmax if Lmax is not None else _default_lmax(len(xs), joint_alphabet_size, min_count)
     counts = JointSuffixCounts.from_sequences(
         xs,
         ys,
@@ -493,6 +510,15 @@ def transcssr(
         output_alphabet=output_alphabet,
         max_length=max_length + 1,
     )
+    if correction == "bonferroni":
+        eligible = sum(
+            1
+            for history in counts.next_counts
+            if 0 < len(history) <= max_length and _observed(counts, history) >= max(1, min_count)
+        )
+        alpha /= max(1, eligible * len(counts.input_alphabet))
+    elif correction is not None:
+        raise ValueError(f"unknown correction {correction!r}")
 
     states = _homogenize(counts, Lmax=max_length, alpha=alpha, test=test, min_count=min_count)
     everything = set(range(len(states)))
diff --git a/sofic/generators/hmm_inference.py b/sofic/generators/hmm_inference.py
index eb86703..f83728a 100644
--- a/sofic/generators/hmm_inference.py
+++ b/sofic/generators/hmm_inference.py
@@ -370,7 +370,10 @@ def baum_welch(
     max_iter: int = 100,
     tol: float = 1e-6,
     estimate_initial: bool = True,
-) -> tuple[Any, list[float]]:
+    n_restarts: int = 1,
+    rng: np.random.Generator | int | None = None,
+    return_restarts: bool = False,
+) -> tuple[Any, list[float]] | tuple[Any, list[float], list[float]]:
     r"""Fit HMM parameters by Baum-Welch (EM) expectation-maximization.
 
     Re-estimates the Mealy joint edge law
@@ -390,6 +393,13 @@ def baum_welch(
     Returns ``(fitted_model, loglik_trace)`` where ``loglik_trace`` is the
     non-decreasing sequence of total natural-log likelihoods observed before each
     parameter update.
+
+    EM converges to a local maximum of the likelihood. With ``n_restarts > 1`` the
+    first run starts from ``hmm``'s parameters and each further run from edge laws
+    drawn uniformly (Dirichlet(1)) over each state's structurally allowed edges;
+    the fit with the highest final log-likelihood is returned. Pass
+    ``return_restarts=True`` to also get every run's final log-likelihood, which
+    shows whether near-equal optima exist.
     """
     from sofic.generators.mealy import MealyHMM
 
@@ -409,6 +419,63 @@ def baum_welch(
         if matrix[i, j] > 0.0
     }
 
+    if n_restarts < 1:
+        raise ValueError("n_restarts must be at least 1")
+    generator = rng if isinstance(rng, np.random.Generator) else np.random.default_rng(rng)
+    runs = []
+    for restart in range(n_restarts):
+        start_joint = joint if restart == 0 else _random_edge_law(joint, support, n_states, generator)
+        runs.append(
+            _baum_welch_run(pi, start_joint, seqs, max_iter=max_iter, tol=tol, estimate_initial=estimate_initial)
+        )
+    finals = [trace[-1] if trace else float("-inf") for _pi, _joint, trace in runs]
+    pi, joint, loglik_trace = runs[int(np.argmax(finals))]
+
+    fitted = MealyHMM(
+        initial_distribution={states[i]: float(pi[i]) for i in range(n_states) if pi[i] > 0.0},
+        observation_alphabet=alphabet,
+    )
+    for state in states:
+        fitted.graph.add_state(state)
+    for i, symbol, j in sorted(support, key=lambda edge: (edge[0], str(edge[1]), edge[2])):
+        prob = float(joint[symbol][i, j])
+        if prob > 0.0:
+            fitted.add_transition(states[i], states[j], symbol, prob)
+    fitted.validate()
+    if return_restarts:
+        return fitted, loglik_trace, finals
+    return fitted, loglik_trace
+
+
+def _random_edge_law(
+    joint: dict[Any, np.ndarray],
+    support: set[tuple[int, Any, int]],
+    n_states: int,
+    rng: np.random.Generator,
+) -> dict[Any, np.ndarray]:
+    """Edge laws drawn uniformly over each state's allowed ``(symbol, target)`` edges."""
+    new_joint = {symbol: np.zeros((n_states, n_states), dtype=float) for symbol in joint}
+    for i in range(n_states):
+        edges = sorted(((symbol, j) for (source, symbol, j) in support if source == i), key=lambda e: (str(e[0]), e[1]))
+        if not edges:
+            continue
+        weights = rng.dirichlet(np.ones(len(edges)))
+        for (symbol, j), weight in zip(edges, weights, strict=True):
+            new_joint[symbol][i, j] = weight
+    return new_joint
+
+
+def _baum_welch_run(
+    pi: np.ndarray,
+    joint: dict[Any, np.ndarray],
+    seqs: list[Any],
+    *,
+    max_iter: int,
+    tol: float,
+    estimate_initial: bool,
+) -> tuple[np.ndarray, dict[Any, np.ndarray], list[float]]:
+    """One EM run from ``(pi, joint)``; returns the final parameters and trace."""
+    n_states = len(pi)
     loglik_trace: list[float] = []
     prev_ll: float | None = None
     for _iteration in range(max_iter):
@@ -443,19 +510,7 @@ def baum_welch(
             mass = float(gamma0_sum.sum())
             if mass > 0.0:
                 pi = gamma0_sum / mass
-
-    fitted = MealyHMM(
-        initial_distribution={states[i]: float(pi[i]) for i in range(n_states) if pi[i] > 0.0},
-        observation_alphabet=alphabet,
-    )
-    for state in states:
-        fitted.graph.add_state(state)
-    for i, symbol, j in sorted(support, key=lambda edge: (edge[0], str(edge[1]), edge[2])):
-        prob = float(joint[symbol][i, j])
-        if prob > 0.0:
-            fitted.add_transition(states[i], states[j], symbol, prob)
-    fitted.validate()
-    return fitted, loglik_trace
+    return pi, joint, loglik_trace
 
 
 def score(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> dict[tuple[Hashable, Any, Hashable], float]:
diff --git a/sofic/generators/stack_inference.py b/sofic/generators/stack_inference.py
index 84bc97f..6028214 100644
--- a/sofic/generators/stack_inference.py
+++ b/sofic/generators/stack_inference.py
@@ -10,7 +10,9 @@
 from sofic.exceptions import StochasticValidationError
 from sofic.generators.epsilon_inference import (
     History,
+    MorphTest,
     SuffixCounts,
+    _bonferroni_alpha,
     _cluster_histories_by_morph,
     _cssr_default_lmax,
     _cssr_determinize,
@@ -18,6 +20,7 @@
     _merge_similar_states,
     morph_test_score,
     morphs_differ,
+    suggest_lmax,
 )
 from sofic.generators.stack_hmm import HiddenMarkovStackModel
 from sofic.graph import ATTR_SYMBOL
@@ -166,7 +169,7 @@ def _stack_homogenize(
     alphabet: DyckAlphabet,
     Lmax: int,
     alpha: float,
-    test: Literal["g", "chi2", "tv"],
+    test: MorphTest,
     max_stack_depth: int,
     min_count: int = 1,
 ) -> tuple[dict[int, set[ConfigurationHistory]], dict[ConfigurationHistory, int]]:
@@ -241,7 +244,7 @@ def _stack_merge(
     counts: StackSuffixCounts,
     *,
     alpha: float,
-    test: Literal["g", "chi2", "tv"],
+    test: MorphTest,
     alphabet: DyckAlphabet,
 ) -> dict[int, set[ConfigurationHistory]]:
     proxy = _control_counts(counts, alphabet).restricted_to(set(history_to_state))
@@ -378,23 +381,45 @@ def stack_cssr(
     sequence: Sequence[Any],
     *,
     alphabet: DyckAlphabet,
-    Lmax: int | None = None,
+    Lmax: int | Literal["auto"] | None = None,
     max_stack_depth: int = 8,
     alpha: float = 0.05,
-    test: Literal["g", "chi2", "tv"] = "g",
+    test: MorphTest = "g",
     min_count: int = 5,
+    correction: Literal["bonferroni"] | None = None,
 ) -> HiddenMarkovStackModel:
-    """Reconstruct a stack HMM via configuration-lifted CSSR."""
+    """Reconstruct a stack HMM via configuration-lifted CSSR.
+
+    ``Lmax="auto"`` uses :func:`~sofic.generators.epsilon_inference.suggest_lmax`
+    on the observed symbols. Stack processes generally have infinite Markov
+    order, so treat it as a lower bound on the suffix length the data support.
+    ``test="exact"`` and ``correction="bonferroni"`` are as in
+    :func:`~sofic.generators.epsilon_inference.cssr`; the correction counts
+    eligible (suffix, stack) configurations.
+    """
     seq = tuple(sequence)
     if len(seq) < 2:
         raise ValueError("sequence must contain at least two symbols")
-    max_length = Lmax if Lmax is not None else _cssr_default_lmax(len(seq), len(alphabet.symbol_alphabet))
+    if Lmax == "auto":
+        max_length = suggest_lmax(seq, alpha=alpha)
+    else:
+        max_length = Lmax if Lmax is not None else _cssr_default_lmax(len(seq), len(alphabet.symbol_alphabet))
     counts = StackSuffixCounts.from_sequence(
         seq,
         alphabet=alphabet,
         max_length=max_length + 1,
         max_stack_depth=max_stack_depth,
     )
+    if correction == "bonferroni":
+        alpha = _bonferroni_alpha(
+            counts,
+            alpha,
+            max_length=max_length,
+            min_count=min_count,
+            suffix_length=lambda history: len(history[0]),
+        )
+    elif correction is not None:
+        raise ValueError(f"unknown correction {correction!r}")
     states, history_to_state = _stack_homogenize(
         counts,
         alphabet=alphabet,
diff --git a/sofic/inference/__init__.py b/sofic/inference/__init__.py
index 05d7664..76fdf41 100644
--- a/sofic/inference/__init__.py
+++ b/sofic/inference/__init__.py
@@ -1,6 +1,14 @@
 """Inference algorithms for stochastic generators."""
 
 from sofic.inference import bayesian
+from sofic.inference.diagnostics import (
+    GoodnessOfFit,
+    StructureStability,
+    goodness_of_fit,
+    reconstruction_sweep,
+    structure_stability,
+    topology_key,
+)
 from sofic.inference.model_selection import (
     ModelScores,
     WAICResult,
@@ -25,6 +33,12 @@
 
 __all__ = [
     "bayesian",
+    "GoodnessOfFit",
+    "StructureStability",
+    "goodness_of_fit",
+    "reconstruction_sweep",
+    "structure_stability",
+    "topology_key",
     "ModelScores",
     "WAICResult",
     "compare_information_criteria",
diff --git a/sofic/inference/diagnostics.py b/sofic/inference/diagnostics.py
new file mode 100644
index 0000000..9b97732
--- /dev/null
+++ b/sofic/inference/diagnostics.py
@@ -0,0 +1,374 @@
+"""Diagnostics for reconstructed machines: goodness of fit and structural stability.
+
+:func:`goodness_of_fit` asks whether a fitted machine reproduces the observed
+word statistics, via a parametric bootstrap from the machine itself
+:cite:`Efron1993`. :func:`structure_stability` and :func:`reconstruction_sweep`
+ask whether the reconstructed *structure* survives resampling the data or
+changing the reconstruction's tuning parameters.
+"""
+
+from __future__ import annotations
+
+from collections import Counter
+from collections.abc import Callable, Hashable, Iterable, Sequence
+from dataclasses import dataclass, field
+from typing import Any, Literal
+
+import numpy as np
+
+from sofic.exceptions import StochasticValidationError
+from sofic.graph import ATTR_EMISSION
+
+__all__ = [
+    "GoodnessOfFit",
+    "StructureStability",
+    "goodness_of_fit",
+    "reconstruction_sweep",
+    "structure_stability",
+    "topology_key",
+]
+
+
+@dataclass(frozen=True)
+class GoodnessOfFit:
+    """The result of :func:`goodness_of_fit`.
+
+    Attributes
+    ----------
+    statistic
+        ``"g"`` or ``"entropy_rate"``.
+    L
+        The word length compared.
+    value
+        The statistic on the observed data.
+    pvalue
+        ``(1 + #{simulated >= observed}) / (1 + n_samples)``. Small values mean
+        the machine does not reproduce the data's length-``L`` statistics.
+    null
+        The statistic on each sequence simulated from the machine.
+    forbidden_words
+        Observed words that the machine assigns probability zero.
+    """
+
+    statistic: str
+    L: int
+    value: float
+    pvalue: float
+    null: np.ndarray
+    forbidden_words: tuple[tuple[Any, ...], ...] = ()
+
+
+def _default_word_length(n: int, alphabet_size: int) -> int:
+    """The longest ``L`` with about ten observations per possible word, from 1 to 6."""
+    k = max(2, alphabet_size)
+    return max(1, min(6, int(np.log(max(n, 1) / 10) / np.log(k))))
+
+
+def _word_counts(sequence: Sequence[Any], L: int) -> Counter[tuple[Any, ...]]:
+    seq = tuple(sequence)
+    return Counter(seq[i : i + L] for i in range(len(seq) - L + 1))
+
+
+def _conditional_entropy(counts: Counter[tuple[Any, ...]]) -> float:
+    """Plug-in ``H[X_{L-1} | X_{0:L-1}]`` in bits from length-``L`` word counts."""
+    total = sum(counts.values())
+    prefixes: Counter[tuple[Any, ...]] = Counter()
+    for word, count in counts.items():
+        prefixes[word[:-1]] += count
+    h = -sum(c * np.log2(c) for c in counts.values()) + sum(c * np.log2(c) for c in prefixes.values())
+    return float(h / total)
+
+
+def _model_conditional_entropy(probabilities: dict[tuple[Any, ...], float]) -> float:
+    prefixes: Counter[tuple[Any, ...]] = Counter()
+    for word, p in probabilities.items():
+        prefixes[word[:-1]] += p
+    joint = -sum(p * np.log2(p) for p in probabilities.values() if p > 0)
+    marginal = -sum(p * np.log2(p) for p in prefixes.values() if p > 0)
+    return float(joint - marginal)
+
+
+def goodness_of_fit(
+    machine: Any,
+    data: Sequence[Any],
+    *,
+    L: int | None = None,
+    statistic: Literal["g", "entropy_rate"] = "g",
+    n_samples: int = 199,
+    burn_in: int = 100,
+    rng: np.random.Generator | int | None = None,
+) -> GoodnessOfFit:
+    """Parametric-bootstrap test that ``machine`` generated ``data``.
+
+    Sequences as long as ``data`` are simulated from ``machine`` (after
+    ``burn_in`` steps, so they start near stationarity), and a length-``L`` word
+    statistic of the data is compared with its distribution over the simulations.
+    Because the null distribution is simulated, the overlap between successive
+    words is accounted for; no chi-squared approximation is used.
+
+    Parameters
+    ----------
+    machine
+        A fitted generator with ``sample``, ``word_probabilities`` and
+        ``stationary_distribution`` (e.g. an ε-machine from :func:`cssr`).
+    data
+        The observed sequence the machine was fitted to.
+    L
+        Word length; by default the longest with about ten observations per
+        possible word (between 1 and 6).
+    statistic
+        ``"g"`` is the G statistic of the observed length-``L`` word counts
+        against the machine's stationary word probabilities. ``"entropy_rate"`` is
+        ``|h_hat - h_L|``, the gap between the plug-in conditional entropy
+        ``H[X_{L-1} | X_{0:L-1}]`` and the machine's value.
+    n_samples
+        Number of simulated sequences.
+    burn_in
+        Steps discarded at the start of each simulation.
+    rng
+        Seed or generator.
+
+    Returns
+    -------
+    GoodnessOfFit
+
+    Notes
+    -----
+    Fitting and testing on the same data makes the test conservative, as in any
+    parametric bootstrap without refitting. A small p-value is still evidence
+    that the reconstruction misses structure. For CSSR that usually means
+    ``Lmax`` is shorter than the source's synchronization length, which happens
+    for strictly sofic sources. An observed word that the machine forbids gives
+    ``G = inf`` and the smallest possible p-value.
+    """
+    generator = rng if isinstance(rng, np.random.Generator) else np.random.default_rng(rng)
+    seq = tuple(data)
+    n = len(seq)
+    alphabet = set(seq) | set(machine.observation_alphabet)
+    if L is None:
+        L = _default_word_length(n, len(alphabet))
+    if L < 1 or n < L:
+        raise ValueError("L must be between 1 and len(data)")
+    pi = np.asarray(machine.stationary_distribution(), dtype=float)
+    probabilities = {tuple(w): float(p) for w, p in machine.word_probabilities(L, start=pi).items()}
+
+    if statistic == "g":
+
+        def compute(sample: Sequence[Any]) -> float:
+            counts = _word_counts(sample, L)
+            total = sum(counts.values())
+            g = 0.0
+            for word, count in counts.items():
+                p = probabilities.get(word, 0.0)
+                if p <= 0.0:
+                    return float("inf")
+                g += count * np.log(count / (total * p))
+            return 2.0 * g
+    elif statistic == "entropy_rate":
+        target = _model_conditional_entropy(probabilities)
+
+        def compute(sample: Sequence[Any]) -> float:
+            return abs(_conditional_entropy(_word_counts(sample, L)) - target)
+    else:
+        raise ValueError(f"unknown statistic {statistic!r}")
+
+    value = compute(seq)
+    null = np.empty(n_samples)
+    for i in range(n_samples):
+        simulated, _ = machine.sample(n + burn_in, generator)
+        null[i] = compute(simulated[burn_in:])
+    tol = 1e-12 * max(1.0, abs(value)) if np.isfinite(value) else 0.0
+    pvalue = float((1 + np.sum(null >= value - tol)) / (1 + n_samples))
+    forbidden = tuple(sorted((w for w in _word_counts(seq, L) if probabilities.get(w, 0.0) <= 0.0), key=repr))
+    return GoodnessOfFit(statistic, int(L), float(value), pvalue, null, forbidden)
+
+
+def topology_key(machine: Any) -> tuple[int, tuple[tuple[int, str, int], ...]]:
+    """An isomorphism-invariant signature of a unifilar machine's labeled topology.
+
+    States are relabeled in breadth-first order from each start state, following
+    symbols in ``repr`` order, and the lexicographically smallest edge list is
+    kept. Two unifilar machines have the same key exactly when their labeled
+    transition graphs are isomorphic; transition probabilities are ignored.
+    """
+    states = list(machine.states())
+    out = {
+        state: sorted(
+            ((repr(t.data.get(ATTR_EMISSION)), t.target) for t in machine.graph.out_transitions(state)),
+            key=lambda item: (item[0], repr(item[1])),
+        )
+        for state in states
+    }
+    best: tuple[tuple[int, str, int], ...] | None = None
+    for start in states:
+        labels: dict[Hashable, int] = {start: 0}
+        order = [start]
+        index = 0
+        while index < len(order):
+            for _symbol, target in out[order[index]]:
+                if target not in labels:
+                    labels[target] = len(labels)
+                    order.append(target)
+            index += 1
+        for state in sorted(states, key=repr):
+            labels.setdefault(state, len(labels))
+        key = tuple(sorted((labels[s], symbol, labels[t]) for s in states for symbol, t in out[s]))
+        if best is None or key < best:
+            best = key
+    return len(states), best or ()
+
+
+@dataclass
+class StructureStability:
+    """The result of :func:`structure_stability`.
+
+    Attributes
+    ----------
+    reference
+        :func:`topology_key` of the reconstruction from the full data.
+    topologies
+        How often each topology was reconstructed across resamples.
+    state_counts
+        How often each number of states was reconstructed.
+    failures
+        Resamples on which reconstruction raised an error.
+    """
+
+    reference: tuple[int, tuple[tuple[int, str, int], ...]]
+    topologies: Counter = field(default_factory=Counter)
+    state_counts: Counter = field(default_factory=Counter)
+    failures: int = 0
+
+    @property
+    def n_resamples(self) -> int:
+        return sum(self.topologies.values()) + self.failures
+
+    @property
+    def reference_fraction(self) -> float:
+        """Fraction of resamples reproducing the full-data topology."""
+        return self.topologies.get(self.reference, 0) / max(1, self.n_resamples)
+
+    @property
+    def modal_topology(self) -> tuple[int, tuple[tuple[int, str, int], ...]] | None:
+        return self.topologies.most_common(1)[0][0] if self.topologies else None
+
+
+Method = Literal["cssr", "subtree", "spectral"]
+
+
+def _reconstruct(sequence: Sequence[Any], method: Method | Callable[..., Any], kwargs: dict[str, Any]) -> Any:
+    if isinstance(method, str):
+        from sofic.generators.epsilon_machine import EpsilonMachine
+
+        return EpsilonMachine.from_sequence(sequence, method=method, **kwargs)
+    return method(sequence, **kwargs)
+
+
+def structure_stability(
+    sequence: Sequence[Any],
+    *,
+    method: Method | Callable[..., Any] = "cssr",
+    n_resamples: int = 50,
+    resample: Literal["subsample", "block"] = "subsample",
+    fraction: float = 0.5,
+    mean_block_length: float | None = None,
+    rng: np.random.Generator | int | None = None,
+    **kwargs: Any,
+) -> StructureStability:
+    """How often the reconstructed topology survives resampling the data.
+
+    Parameters
+    ----------
+    sequence
+        The observed sequence.
+    method
+        Reconstruction method passed to
+        :meth:`~sofic.generators.epsilon_machine.EpsilonMachine.from_sequence`,
+        or a callable ``method(sequence, **kwargs)`` returning a machine.
+    n_resamples
+        Number of resampled reconstructions.
+    resample
+        ``"subsample"`` reconstructs from random contiguous segments of
+        ``fraction * len(sequence)`` symbols, which contain no artificial
+        junctions :cite:`Politis1999`. ``"block"`` uses the stationary bootstrap
+        :cite:`Politis1994` (via :func:`dit.inference.stationary_bootstrap`).
+        Block junctions create words the source never emits, which can add
+        spurious states.
+    fraction
+        Segment length for ``"subsample"``, as a fraction of the data.
+    mean_block_length
+        Mean block length for ``"block"``.
+    rng
+        Seed or generator.
+    **kwargs
+        Forwarded to the reconstruction (e.g. ``Lmax``, ``alpha``).
+
+    Returns
+    -------
+    StructureStability
+
+    Notes
+    -----
+    Subsamples are shorter than the data, so they detect less structure; a
+    topology that appears in most half-length subsamples is well supported, while
+    one that rarely reappears reflects the particular sample. Statistical
+    confidence in an inferred *structure* is not otherwise quantified by CSSR.
+    """
+    generator = rng if isinstance(rng, np.random.Generator) else np.random.default_rng(rng)
+    seq = list(sequence)
+    n = len(seq)
+    result = StructureStability(reference=topology_key(_reconstruct(seq, method, kwargs)))
+    if resample == "subsample":
+        if not 0.0 < fraction <= 1.0:
+            raise ValueError("fraction must be in (0, 1]")
+        m = max(2, int(round(fraction * n)))
+        starts = generator.integers(0, n - m + 1, size=n_resamples)
+        replicates: Iterable[Sequence[Any]] = (seq[s : s + m] for s in starts)
+    elif resample == "block":
+        import dit.inference
+
+        stationary_bootstrap = getattr(dit.inference, "stationary_bootstrap", None)
+        if stationary_bootstrap is None:  # pragma: no cover - depends on the installed dit
+            raise ImportError('resample="block" requires a dit release with dit.inference.stationary_bootstrap')
+        array = np.empty(n, dtype=object)
+        array[:] = seq
+        replicates = (list(r) for r in stationary_bootstrap(array, n_resamples, mean_block_length, generator))
+    else:
+        raise ValueError(f"unknown resample {resample!r}")
+    for replicate in replicates:
+        try:
+            machine = _reconstruct(replicate, method, kwargs)
+        except (StochasticValidationError, ValueError):
+            result.failures += 1
+            continue
+        key = topology_key(machine)
+        result.topologies[key] += 1
+        result.state_counts[key[0]] += 1
+    return result
+
+
+def reconstruction_sweep(
+    sequence: Sequence[Any],
+    *,
+    alphas: Sequence[float] = (0.05, 0.01, 0.001),
+    lmaxes: Sequence[int] = (1, 2, 3, 4),
+    method: Literal["cssr"] | Callable[..., Any] = "cssr",
+    **kwargs: Any,
+) -> dict[tuple[float, int], tuple[int, tuple[tuple[int, str, int], ...]] | None]:
+    """Reconstruct over a grid of significance levels and history lengths.
+
+    Returns ``{(alpha, Lmax): topology_key}``, with ``None`` where reconstruction
+    failed. A structure that persists across a range of ``alpha`` and ``Lmax`` is
+    better supported than one that appears at a single setting. For a Markov
+    source, it should be stable for every ``Lmax`` at or above the source's order.
+    """
+    results: dict[tuple[float, int], tuple[int, tuple[tuple[int, str, int], ...]] | None] = {}
+    for alpha in alphas:
+        for lmax in lmaxes:
+            try:
+                machine = _reconstruct(sequence, method, {**kwargs, "alpha": alpha, "Lmax": lmax})
+            except (StochasticValidationError, ValueError):
+                results[alpha, lmax] = None
+                continue
+            results[alpha, lmax] = topology_key(machine)
+    return results
diff --git a/sofic/inference/model_selection.py b/sofic/inference/model_selection.py
index 970fe4f..c108140 100644
--- a/sofic/inference/model_selection.py
+++ b/sofic/inference/model_selection.py
@@ -116,6 +116,39 @@ def _total_log_likelihood(model: HiddenMarkovModel, sequences: Sequence[Sequence
     return total
 
 
+def _smoothed_log_likelihood(
+    model: HiddenMarkovModel,
+    sequence: Sequence[Any],
+    *,
+    smoothing: float,
+    alphabet_size: int,
+) -> float:
+    """Natural-log likelihood with each one-step prediction mixed with the uniform law.
+
+    ``P'(x_t | x_{0:t}) = (1 - smoothing) P(x_t | x_{0:t}) + smoothing / alphabet_size``,
+    computed by forward filtering. After a symbol the model forbids, the belief is
+    propagated without conditioning on it.
+    """
+    from sofic.generators.hmm_inference import _emission_transition_tensors_from_mealy
+
+    pi, joint = _emission_transition_tensors_from_mealy(model.to_mealy())
+    total_step = sum(joint.values())
+    belief = np.asarray(pi, dtype=float)
+    belief = belief / belief.sum()
+    total = 0.0
+    for symbol in sequence:
+        matrix = joint.get(symbol)
+        unnormalized = belief @ matrix if matrix is not None else np.zeros_like(belief)
+        predicted = float(unnormalized.sum())
+        total += float(np.log((1.0 - smoothing) * predicted + smoothing / alphabet_size))
+        if predicted > 0.0:
+            belief = unnormalized / predicted
+        else:
+            belief = belief @ total_step
+            belief = belief / belief.sum()
+    return total
+
+
 def score_model(
     model: HiddenMarkovModel,
     data: Iterable[Any],
@@ -193,6 +226,8 @@ def cross_validated_log_likelihood(
     *,
     folds: int = 5,
     rng: np.random.Generator | int | None = None,
+    gap: int = 0,
+    smoothing: float = 0.0,
 ) -> float:
     """Return the total held-out natural-log likelihood under ``folds``-fold CV.
 
@@ -201,8 +236,26 @@ def cross_validated_log_likelihood(
     sequences; a single long sequence is split into ``folds`` contiguous blocks.
     Each held-out block is scored under a model trained on the remaining data and
     the contributions are summed (higher is better). A fold whose held-out data
-    has zero probability contributes ``-inf``.
+    has zero probability contributes ``-inf`` unless ``smoothing > 0``.
+
+    Parameters
+    ----------
+    gap
+        For contiguous blocks of one long sequence, drop this many symbols from
+        each training block on the side adjacent to the held-out block. Neighboring
+        blocks of a dependent sequence are correlated, so without a gap the
+        held-out score is optimistic (buffered or "h-block" cross-validation
+        :cite:`Burman1994`). A gap of the order of the process's memory suffices.
+    smoothing
+        Mix each one-step held-out prediction with the uniform distribution over
+        the observed alphabet, with this weight. Then a single transition that the
+        fitted model forbids costs ``log(smoothing / |A|)`` instead of making the
+        whole fold ``-inf``, so models can still be compared.
     """
+    if gap < 0:
+        raise ValueError("gap must be non-negative")
+    if not 0.0 <= smoothing < 1.0:
+        raise ValueError("smoothing must be in [0, 1)")
     generator = rng if isinstance(rng, np.random.Generator) else np.random.default_rng(rng)
     sequences = _normalize_sequences(data)
     if folds < 2:
@@ -219,16 +272,31 @@ def cross_validated_log_likelihood(
         blocks = [list(chunk) for chunk in np.array_split(np.array(flat, dtype=object), folds)]
         blocks = [[list(block)] for block in blocks]
 
+    contiguous = len(sequences) < folds
+    alphabet_size = max(1, len({symbol for seq in sequences for symbol in seq}))
     total = 0.0
     for held_out_index in range(len(blocks)):
         train: list[Any] = []
         for index, block in enumerate(blocks):
             if index == held_out_index:
                 continue
+            if contiguous and gap:
+                (segment,) = block
+                if index == held_out_index - 1:
+                    segment = segment[: max(0, len(segment) - gap)]
+                elif index == held_out_index + 1:
+                    segment = segment[gap:]
+                block = [segment] if segment else []
             train.extend(block)
-        held_out = blocks[held_out_index]
+        held_out = _normalize_sequences(blocks[held_out_index])
         model = fit(train)
-        total += _total_log_likelihood(model, _normalize_sequences(held_out))
+        if smoothing > 0.0:
+            total += sum(
+                _smoothed_log_likelihood(model, seq, smoothing=smoothing, alphabet_size=alphabet_size)
+                for seq in held_out
+            )
+        else:
+            total += _total_log_likelihood(model, held_out)
     return float(total)
 
 
diff --git a/tests/test_epsilon_inference.py b/tests/test_epsilon_inference.py
index 47a2773..532b49d 100644
--- a/tests/test_epsilon_inference.py
+++ b/tests/test_epsilon_inference.py
@@ -273,3 +273,94 @@ def test_subtree_merge_default_delta_recovers_process(name: str, L: int, n_state
     inferred.validate()
     assert len(list(inferred.states())) == n_states
     assert inferred.entropy_rate() == pytest.approx(oracle.entropy_rate(), abs=0.02)
+
+
+def _has_markov_order_selection() -> bool:
+    import dit.inference
+
+    return hasattr(dit.inference, "select_markov_order")
+
+
+needs_markov_order = pytest.mark.skipif(
+    not _has_markov_order_selection(), reason="needs dit.inference.select_markov_order"
+)
+
+
+@needs_markov_order
+def test_suggest_lmax_markov_sources():
+    from sofic.examples import processes
+    from sofic.generators.epsilon_inference import suggest_lmax
+
+    observations, _ = sample(golden_mean(0.5), 4000, np.random.default_rng(1))
+    assert suggest_lmax(observations) == 1
+    observations, _ = sample(processes.RkGM(3, 2), 20000, np.random.default_rng(2))
+    assert suggest_lmax(observations, method="bic") == 3
+
+
+@needs_markov_order
+def test_suggest_lmax_grows_for_even_process():
+    """The even process has infinite Markov order, so the suggestion grows with data."""
+    from sofic.generators.epsilon_inference import suggest_lmax
+
+    short, _ = sample(even_process(0.5), 300, np.random.default_rng(3))
+    long, _ = sample(even_process(0.5), 30000, np.random.default_rng(3))
+    assert suggest_lmax(long, method="bic") > suggest_lmax(short, method="bic")
+
+
+@needs_markov_order
+def test_cssr_auto_lmax_golden_mean(rng: np.random.Generator):
+    observations, _ = sample(golden_mean(0.5), 8000, rng)
+    inferred = cssr(observations, Lmax="auto", alpha=0.001)
+    assert len(list(inferred.states())) == 2
+    assert _signatures_isomorphic(inferred, golden_mean(0.5), prob_tol=0.1)
+
+
+def test_exact_morph_test_small_counts():
+    """With tiny counts the exact test is calibrated where the chi-squared limit is not."""
+    from sofic.generators.epsilon_inference import SuffixCounts, morphs_differ
+
+    rng = np.random.default_rng(4)
+    rejections = {"g": 0, "exact": 0}
+    trials = 300
+    for _ in range(trials):
+        counts = SuffixCounts(alphabet=(0, 1, 2))
+        for history in ((0,), (1,)):
+            for symbol in rng.choice(3, size=6, p=[0.8, 0.1, 0.1]):
+                counts.next_counts[history][int(symbol)] += 1
+                counts.history_counts[history] += 1
+        for test in rejections:
+            rejections[test] += morphs_differ(counts, {(0,)}, {(1,)}, alpha=0.05, test=test)
+    assert rejections["exact"] / trials <= 0.08
+    assert rejections["exact"] <= rejections["g"]
+
+
+def test_exact_morph_test_is_deterministic(rng: np.random.Generator):
+    observations, _ = sample(golden_mean(0.5), 3000, rng)
+    first = cssr(observations, Lmax=3, alpha=0.01, test="exact")
+    second = cssr(observations, Lmax=3, alpha=0.01, test="exact")
+    assert _transition_signature(first) == _transition_signature(second)
+    assert len(list(first.states())) == 2
+
+
+def test_cssr_bonferroni_reduces_spurious_states():
+    """An i.i.d. source with a long Lmax: the corrected test keeps a single state."""
+    observations, _ = sample(bernoulli(0.3), 3000, np.random.default_rng(6))
+    inferred = cssr(observations, Lmax=6, alpha=0.05, correction="bonferroni")
+    assert len(list(inferred.states())) == 1
+    with pytest.raises(ValueError, match="unknown correction"):
+        cssr(observations, Lmax=2, correction="holm")
+
+
+@pytest.mark.parametrize("kwargs", [{"test": "exact"}, {"correction": "bonferroni"}, {"alpha": 0.001}])
+def test_subtree_merge_options_golden_mean(kwargs):
+    observations, _ = sample(golden_mean(0.5), 6000, np.random.default_rng(7))
+    inferred = subtree_merge(observations, L=2, **kwargs)
+    assert len(list(inferred.states())) == 2
+    with pytest.raises(ValueError, match="unknown correction"):
+        subtree_merge(observations, L=2, correction="holm")
+
+
+@needs_markov_order
+def test_subtree_merge_auto_depth():
+    observations, _ = sample(golden_mean(0.5), 6000, np.random.default_rng(8))
+    assert len(list(subtree_merge(observations, L="auto").states())) == 2
diff --git a/tests/test_epsilon_transducer_inference.py b/tests/test_epsilon_transducer_inference.py
index 66094c9..6e7137f 100644
--- a/tests/test_epsilon_transducer_inference.py
+++ b/tests/test_epsilon_transducer_inference.py
@@ -114,3 +114,35 @@ def test_recovers_two_step_delay():
     assert len(list(eps.states())) == 4
     test_xs, test_ys = _paired_samples(Delay(2), 3000, seed=1)
     assert _held_out_bits_per_symbol(eps, test_xs, test_ys) == pytest.approx(0.0, abs=1e-9)
+
+
+def _has_markov_order_selection() -> bool:
+    import dit.inference
+
+    return hasattr(dit.inference, "select_markov_order")
+
+
+@pytest.mark.parametrize("seed", [0, 1])
+def test_exact_and_bonferroni_recover_delay_memory(seed):
+    xs, ys = _paired_samples(Delay(1), 6000, seed=seed)
+    for kwargs in ({"test": "exact"}, {"correction": "bonferroni"}):
+        eps = transcssr(xs, ys, input_alphabet=("0", "1"), output_alphabet=("0", "1"), Lmax=2, **kwargs)
+        eps.validate()
+        assert len(list(eps.states())) == 2
+
+
+def test_bonferroni_keeps_memoryless_channel_single_state():
+    xs, ys = _paired_samples(BinaryChannel(0.1, 0.2), 4000, seed=3)
+    eps = transcssr(
+        xs, ys, input_alphabet=("0", "1"), output_alphabet=("0", "1"), Lmax=4, alpha=0.05, correction="bonferroni"
+    )
+    assert len(list(eps.states())) == 1
+    with pytest.raises(ValueError, match="unknown correction"):
+        transcssr(xs, ys, Lmax=1, correction="holm")
+
+
+@pytest.mark.skipif(not _has_markov_order_selection(), reason="needs dit.inference.select_markov_order")
+def test_auto_lmax_delay():
+    xs, ys = _paired_samples(Delay(1), 6000, seed=4)
+    eps = transcssr(xs, ys, input_alphabet=("0", "1"), output_alphabet=("0", "1"), Lmax="auto")
+    assert len(list(eps.states())) == 2
diff --git a/tests/test_hmm_inference.py b/tests/test_hmm_inference.py
index 2cf6756..db6668c 100644
--- a/tests/test_hmm_inference.py
+++ b/tests/test_hmm_inference.py
@@ -347,3 +347,36 @@ def test_seeded_sample_is_reproducible_across_hash_seeds():
         for seed in ("1", "2", "3")
     }
     assert len(outputs) == 1
+
+
+def _symmetric_two_state() -> MealyHMM:
+    """A fully connected 2-state binary HMM whose symmetric start is an EM fixed point."""
+    hmm = MealyHMM(observation_alphabet=frozenset({0, 1}), initial_distribution={"A": 0.5, "B": 0.5})
+    for state in ("A", "B"):
+        hmm.graph.add_state(state)
+    for source in ("A", "B"):
+        for target in ("A", "B"):
+            for symbol in (0, 1):
+                hmm.add_transition(source, target, symbol, 0.25)
+    hmm.validate()
+    return hmm
+
+
+def test_baum_welch_restarts_escape_symmetric_fixed_point():
+    data = sample(golden_mean(0.5), 2000, rng=np.random.default_rng(4))[0]
+    single, _ = baum_welch(_symmetric_two_state(), data, max_iter=300)
+    best, _trace, finals = baum_welch(
+        _symmetric_two_state(), data, max_iter=300, n_restarts=5, rng=0, return_restarts=True
+    )
+    assert len(finals) == 5
+    assert log_likelihood(best, data) > log_likelihood(single, data) + 10
+    assert log_likelihood(best, data) == pytest.approx(max(finals), abs=1.0)
+
+
+def test_baum_welch_restarts_reproducible_and_validated():
+    data = sample(golden_mean(0.4), 500, rng=np.random.default_rng(5))[0]
+    first, _ = baum_welch(_symmetric_two_state(), data, n_restarts=3, rng=1)
+    second, _ = baum_welch(_symmetric_two_state(), data, n_restarts=3, rng=1)
+    assert log_likelihood(first, data) == pytest.approx(log_likelihood(second, data))
+    with pytest.raises(ValueError):
+        baum_welch(_symmetric_two_state(), data, n_restarts=0)
diff --git a/tests/test_inference_diagnostics.py b/tests/test_inference_diagnostics.py
new file mode 100644
index 0000000..f87c31b
--- /dev/null
+++ b/tests/test_inference_diagnostics.py
@@ -0,0 +1,87 @@
+"""Tests for sofic.inference.diagnostics."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from sofic.examples.epsilon_machines import even_process, golden_mean
+from sofic.generators.epsilon_inference import cssr
+from sofic.generators.hmm_inference import sample
+from sofic.inference.diagnostics import (
+    goodness_of_fit,
+    reconstruction_sweep,
+    structure_stability,
+    topology_key,
+)
+
+
+@pytest.fixture(scope="module")
+def even_sample():
+    observations, _ = sample(even_process(0.5), 4000, np.random.default_rng(0))
+    return observations
+
+
+@pytest.mark.parametrize("statistic", ["g", "entropy_rate"])
+def test_goodness_of_fit_accepts_correct_machine(even_sample, statistic):
+    machine = cssr(even_sample, Lmax=4, alpha=0.001)
+    result = goodness_of_fit(machine, even_sample, L=5, statistic=statistic, n_samples=49, rng=1)
+    assert result.pvalue > 0.05
+    assert result.null.shape == (49,)
+    assert result.forbidden_words == ()
+
+
+@pytest.mark.parametrize("statistic", ["g", "entropy_rate"])
+def test_goodness_of_fit_rejects_short_lmax(even_sample, statistic):
+    """CSSR with Lmax=1 cannot capture the even process's parity."""
+    machine = cssr(even_sample, Lmax=1, alpha=0.001)
+    result = goodness_of_fit(machine, even_sample, L=6, statistic=statistic, n_samples=49, rng=1)
+    assert result.pvalue <= 0.05
+
+
+def test_goodness_of_fit_forbidden_word():
+    observations = [0, 1, 1, 0, 1, 0, 0, 1] * 20
+    result = goodness_of_fit(golden_mean(0.5), observations, L=2, n_samples=9, rng=0)
+    assert result.value == float("inf")
+    assert (1, 1) in result.forbidden_words
+    assert result.pvalue == pytest.approx(0.1)
+    with pytest.raises(ValueError):
+        goodness_of_fit(golden_mean(0.5), observations, L=2, statistic="nope")
+
+
+def test_topology_key_is_isomorphism_invariant(even_sample):
+    inferred = cssr(even_sample, Lmax=4, alpha=0.001)
+    assert topology_key(inferred) == topology_key(even_process(0.5))
+    assert topology_key(golden_mean(0.3)) == topology_key(golden_mean(0.7))
+    assert topology_key(golden_mean(0.5)) != topology_key(even_process(0.5))
+
+
+def test_structure_stability_subsample(even_sample):
+    result = structure_stability(even_sample, n_resamples=12, rng=0, Lmax=4, alpha=0.001)
+    assert result.reference == topology_key(even_process(0.5))
+    assert result.reference_fraction >= 0.6
+    assert result.n_resamples == 12
+    assert result.modal_topology == result.reference
+
+
+def _has_stationary_bootstrap() -> bool:
+    import dit.inference
+
+    return hasattr(dit.inference, "stationary_bootstrap")
+
+
+@pytest.mark.skipif(not _has_stationary_bootstrap(), reason="needs dit.inference.stationary_bootstrap")
+def test_structure_stability_block_runs(even_sample):
+    result = structure_stability(
+        even_sample, n_resamples=4, rng=0, resample="block", mean_block_length=500, Lmax=3, alpha=0.001
+    )
+    assert result.n_resamples == 4
+    with pytest.raises(ValueError):
+        structure_stability(even_sample, n_resamples=1, resample="nope", Lmax=2)
+
+
+def test_reconstruction_sweep_golden_mean():
+    observations, _ = sample(golden_mean(0.5), 4000, np.random.default_rng(2))
+    sweep = reconstruction_sweep(observations, alphas=(0.01, 0.001), lmaxes=(1, 2, 3))
+    target = topology_key(golden_mean(0.5))
+    assert all(key == target for key in sweep.values())
diff --git a/tests/test_model_selection.py b/tests/test_model_selection.py
index 22d1950..39d661b 100644
--- a/tests/test_model_selection.py
+++ b/tests/test_model_selection.py
@@ -187,3 +187,40 @@ def test_rank_topological_epsilon_machines_prefers_two_states():
     best = ranked[0]
     assert len(list(best.machine.states())) == 2
     assert best.criterion_value < ranked[-1].criterion_value
+
+
+def test_cross_validation_smoothing_keeps_forbidden_folds_finite():
+    """The golden mean forbids 11; a held-out '11' makes an unsmoothed fold -inf."""
+    rng = np.random.default_rng(6)
+    data, _ = golden_mean(0.5).sample(1000, rng=rng)
+    data = list(data)
+    data[500:502] = [1, 1]
+
+    def fit_golden(train):
+        model, _trace = golden_mean(0.6).baum_welch(train)
+        return model
+
+    assert cross_validated_log_likelihood(fit_golden, data, folds=4) == float("-inf")
+    smoothed = cross_validated_log_likelihood(fit_golden, data, folds=4, smoothing=0.01)
+    assert np.isfinite(smoothed)
+    clean, _ = golden_mean(0.5).sample(1000, rng=np.random.default_rng(7))
+    exact = cross_validated_log_likelihood(fit_golden, list(clean), folds=4)
+    nearly = cross_validated_log_likelihood(fit_golden, list(clean), folds=4, smoothing=1e-9)
+    assert nearly == pytest.approx(exact, rel=1e-6)
+
+
+def test_cross_validation_gap_drops_boundary_symbols():
+    seen: list[int] = []
+
+    def fit_counting(train):
+        seen.append(sum(len(seq) for seq in train))
+        return _iid_binary(0.5)
+
+    data = [0, 1] * 200
+    cross_validated_log_likelihood(fit_counting, data, folds=4)
+    cross_validated_log_likelihood(fit_counting, data, folds=4, gap=10)
+    plain, gapped = seen[:4], seen[4:]
+    # Interior folds lose two gaps, the end folds one.
+    assert [p - g for p, g in zip(plain, gapped, strict=True)] == [10, 20, 20, 10]
+    with pytest.raises(ValueError):
+        cross_validated_log_likelihood(fit_counting, data, folds=4, gap=-1)
diff --git a/tests/test_stack_inference.py b/tests/test_stack_inference.py
index e87af43..ce3d78e 100644
--- a/tests/test_stack_inference.py
+++ b/tests/test_stack_inference.py
@@ -244,3 +244,20 @@ def test_stack_cssr_matches_motzkin_likelihood(Lmax: int):
     assert _held_out_bits_per_symbol(inferred, held_out) == pytest.approx(
         _held_out_bits_per_symbol(oracle, held_out), abs=0.05
     )
+
+
+def test_stack_cssr_exact_and_bonferroni_run():
+    shift = motzkin_shift()
+    oracle = HiddenMarkovStackModel.from_sofic_dyck_shift(shift, _uniform_probabilities(shift))
+    observations, _ = oracle.sample(3000, rng=np.random.default_rng(1))
+    alphabet = DyckAlphabet(
+        call_alphabet=shift.call_alphabet,
+        return_alphabet=shift.return_alphabet,
+        internal_alphabet=shift.internal_alphabet,
+    )
+    for kwargs in ({"test": "exact"}, {"correction": "bonferroni"}):
+        inferred = stack_cssr(observations, alphabet=alphabet, Lmax=3, max_stack_depth=4, alpha=0.001, **kwargs)
+        inferred.validate()
+        assert oracle.word_probability(tuple(observations[:12])) > 0.0
+    with pytest.raises(ValueError, match="unknown correction"):
+        stack_cssr(observations, alphabet=alphabet, Lmax=2, correction="holm")