"""Parametric bootstrap goodness-of-fit for the BLIM.
The chi-squared reference for the likelihood-ratio statistic ``G2`` is
asymptotic and can be unreliable on sparse response tables (Koehler &
Larntz, 1980); sparsity alone is not a universal failure criterion. The
parametric bootstrap calibrates the test on the estimated model instead:
simulate ``n_replicates`` data sets of the observed size from the fitted
parameters, refit each with the same estimator settings, and locate the
observed ``G2`` within the bootstrap distribution of replicated ``G2``
values (Efron & Tibshirani, 1993). The Monte Carlo p-value uses the
add-one convention ``(1 + #{G2* >= G2_obs}) / (n_replicates + 1)``
(Davison & Hinkley, 1997, Eq. 4.12), so it is never exactly zero.
The computation is sequential by design — one process, one core; on the
data sets of the accompanying article a ``n_replicates=2000`` run takes
minutes. Refits that stop at the iteration cap are counted in
``n_capped`` and reported with a :class:`ConvergenceWarning`: their
log-likelihood may be below the attainable maximum, potentially inflating
replicated ``G2``. Convergence also does not guarantee a global maximum.
Parameter replicates retain their item/state alignment and convergence
status. Their dispersion describes the fitted resampling procedure; it
does not resolve non-identifiability or guarantee confidence coverage.
References:
Davison, A. C., & Hinkley, D. V. (1997). Bootstrap Methods and
their Application. Cambridge University Press.
Efron, B., & Tibshirani, R. J. (1993). An Introduction to the
Bootstrap. Chapman & Hall.
Koehler, K. J., & Larntz, K. (1980). An empirical investigation of
goodness-of-fit statistics for sparse multinomials. JASA, 75.
"""
from __future__ import annotations
import warnings
from dataclasses import dataclass, field
from operator import index
from typing import Literal
import numpy as np
from knowledgespaces._patterns import _positive_integer
from knowledgespaces.estimation._starts import _validate_start_options
from knowledgespaces.estimation.blim_em import (
BLIMEstimate,
ConvergenceWarning,
ResponseMatrix,
_log_likelihood_at_params,
estimate_blim,
estimate_blim_restarts,
)
from knowledgespaces.estimation.constraints import BLIMConstraints, _compile_constraints
from knowledgespaces.estimation.discrepancy import MDOptions
from knowledgespaces.estimation.simulate import simulate_blim
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
[docs]
@dataclass(frozen=True)
class BootstrapParameterSummary:
"""Means and sample standard deviations of parameter replicates.
Arrays follow ``items`` and ``states`` and are read-only. ``n_used``
records the selected replicates, with ``n_capped`` unconverged refits
among them. Standard deviations use ddof=1 and are NaN with fewer than
two selected replicates; means are NaN with none. These are empirical
bootstrap dispersions, not automatic confidence intervals or evidence
of identifiable parameters. Selection of converged fits can itself
change the resampling distribution.
"""
items: tuple[str, ...]
states: tuple[frozenset[str], ...]
n_used: int
n_capped: int
beta_mean: np.ndarray
beta_sd: np.ndarray
eta_mean: np.ndarray
eta_sd: np.ndarray
pi_mean: np.ndarray
pi_sd: np.ndarray
[docs]
@dataclass(frozen=True)
class BootstrapGOF:
"""Result of a parametric bootstrap of G2 or complete-table Pearson X2.
Attributes
----------
p_value : float
Monte Carlo p-value ``(1 + n_extreme) / (n_replicates + 1)`` when
observed/replicate fits converge and all selected statistics are
finite; NaN otherwise. descriptive_tail_fraction retains the finite
counting result for inspection even with capped fits.
g2_observed : float
``G2`` of the fit to the observed data.
g2_replicates : np.ndarray
The ``n_replicates`` bootstrap ``G2`` values.
n_replicates : int
Number of bootstrap replicates.
n_extreme : int
Replicates at least as extreme as the selected observed statistic.
n_capped : int
Replicate refits stopped by the iteration cap (not converged);
nonzero counts are also reported via :class:`ConvergenceWarning`.
method : str
Estimation method used for the observed fit and every refit.
seed : int | None
Seed of the replicate generator (``None`` = nondeterministic).
estimate : BLIMEstimate
The fit to the observed data that generated the replicates.
beta_replicates, eta_replicates, pi_replicates : np.ndarray
Rows correspond to ``g2_replicates``; columns match ``estimate.items``
and ``estimate.states``. No unconverged refit is silently removed.
converged_replicates, iterations_replicates : np.ndarray
Convergence status and iteration counts of each refit. All replicate
arrays returned by :func:`bootstrap_gof` are read-only snapshots.
statistic : {"G2", "X2"}
The selected calibration statistic. statistic_observed and
statistic_replicates expose it consistently. G2 fields remain G2
even when Pearson X2 determines the p-value. X2 fields are only
populated when requested. n_failed/replicate_errors retain failed
refits with NaN parameter/statistic rows, distinct from capped runs.
"""
p_value: float
g2_observed: float
g2_replicates: np.ndarray
n_replicates: int
n_extreme: int
n_capped: int
method: str
seed: int | None
estimate: BLIMEstimate
beta_replicates: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
eta_replicates: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
pi_replicates: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
converged_replicates: np.ndarray = field(default_factory=lambda: np.empty(0, dtype=bool))
iterations_replicates: np.ndarray = field(default_factory=lambda: np.empty(0, dtype=int))
statistic: Literal["G2", "X2"] = field(default="G2", kw_only=True)
x2_observed: float | None = field(default=None, kw_only=True)
x2_replicates: np.ndarray = field(default_factory=lambda: np.empty(0), kw_only=True)
n_failed: int = field(default=0, kw_only=True)
replicate_errors: tuple[str | None, ...] = field(default=(), kw_only=True)
n_restarts: int | None = field(default=None, kw_only=True)
init_strategy: Literal["uniform", "pks"] = field(default="uniform", kw_only=True)
init_range: tuple[float, float] = field(default=(0.01, 0.4), kw_only=True)
@property
def statistic_observed(self) -> float:
"""Observed value of the statistic selected for calibration."""
if self.statistic == "G2":
return self.g2_observed
if self.x2_observed is None:
raise ValueError("This result does not contain a Pearson statistic.")
return self.x2_observed
@property
def statistic_replicates(self) -> np.ndarray:
"""All replicate values of the selected statistic, including NaN failures."""
return self.g2_replicates if self.statistic == "G2" else self.x2_replicates
@property
def descriptive_tail_fraction(self) -> float:
"""Add-one tail fraction for all finite statistics, even capped fits.
This is descriptive when p_value is NaN; it is not calibrated
inference. Missing/nonfinite statistics give NaN, never deletion.
"""
if (
not np.isfinite(self.statistic_observed)
or not np.isfinite(self.statistic_replicates).all()
):
return float("nan")
return (1 + self.n_extreme) / (1 + self.n_replicates)
[docs]
def parameter_summary(self, *, converged_only: bool = False) -> BootstrapParameterSummary:
"""Summarize all refits, or explicitly select only converged refits."""
if len(self.converged_replicates) != self.n_replicates:
raise ValueError("This result does not contain parameter replicates.")
mask = (
self.converged_replicates if converged_only else np.ones(self.n_replicates, dtype=bool)
)
n = int(mask.sum())
summaries: list[np.ndarray] = []
for values in (self.beta_replicates, self.eta_replicates, self.pi_replicates):
selected = values[mask]
mean = selected.mean(axis=0) if n else np.full(values.shape[1], np.nan)
sd = (
(selected - selected[0]).std(axis=0, ddof=1)
if n > 1
else np.full(values.shape[1], np.nan)
)
mean.flags.writeable = sd.flags.writeable = False
summaries.extend((mean, sd))
return BootstrapParameterSummary(
tuple(self.estimate.items),
tuple(self.estimate.states),
n,
int(((~self.converged_replicates & (self.iterations_replicates > 0))[mask]).sum()),
*summaries,
)
[docs]
def bootstrap_gof(
structure: KnowledgeStructure,
data: ResponseMatrix,
*,
n_replicates: int = 1000,
method: Literal["ML", "MD", "MDML"] = "ML",
max_iter: int = 500,
tol: float = 1e-6,
seed: int | None = None,
estimate: BLIMEstimate | None = None,
constraints: BLIMConstraints | None = None,
discrepancy: MDOptions | None = None,
statistic: Literal["G2", "X2"] = "G2",
n_restarts: int | None = None,
init_strategy: Literal["uniform", "pks"] = "uniform",
init_range: tuple[float, float] = (0.01, 0.4),
max_patterns: int = 1_048_576,
max_memory_bytes: int = 512_000_000,
) -> BootstrapGOF:
"""Parametric bootstrap of G2 or Pearson X2 and BLIM parameter estimates.
Parameters
----------
structure : KnowledgeStructure
The knowledge structure under test.
data : ResponseMatrix
Observed response patterns with integer frequencies. Fractional
analysis weights cannot be resampled as counts of respondents.
n_replicates : int
Number of bootstrap replicates. Default 1000; the p-value
resolution is ``1 / (n_replicates + 1)``.
method : {"ML", "MD", "MDML"}
Estimator for the observed fit and every replicate refit.
max_iter, tol :
Estimator settings, applied identically to the observed fit
(unless ``estimate`` is supplied) and to every refit.
seed : int | None
Seed for the replicate generator. Fixing it makes the entire
bootstrap reproducible.
estimate : BLIMEstimate | None
A previously computed fit of ``structure`` to ``data``. When
given, the observed fit is not recomputed; its ``method`` must
match ``method``, since the refits must use the same estimator
as the fit that produced the observed ``G2``. State and item
alignment, probabilities, and likelihood/G2 on the supplied
observations are checked before any bootstrap simulation.
SLM estimates require :func:`bootstrap_slm`, which preserves their
restricted state distribution during refitting.
constraints : BLIMConstraints or None
Restrictions for the observed fit and every refit. If an estimate
is supplied, None inherits its constraints; an explicit value must
match. Simulating from a constrained fit and refitting an unrestricted
model would calibrate a different procedure.
discrepancy : MDOptions or None
Radius or inclusion weights for the observed fit and every refit.
None inherits a supplied estimate's options; explicit options must
match. Applies to MD/MDML, with hyperbolic weights limited to MD.
statistic : {"G2", "X2"}
Default G2. X2 sums squared Pearson residuals over the complete
response table, including unobserved cells. This does not assert an
asymptotic chi-squared reference distribution.
n_restarts, init_strategy, init_range :
None retains deterministic standard starts. A positive integer uses
estimate_blim_restarts with the same settings for observed/refit
samples, from one seeded stream. Only ML/MDML support this search.
Omit a supplied estimate to ensure an identical observed policy.
max_patterns, max_memory_bytes :
Bounds for Pearson full-table diagnostics and retained replicates;
max_memory_bytes is also forwarded to estimation. Estimates of array
allocations, not a hard process-memory cap.
Returns
-------
BootstrapGOF
Warns
-----
ConvergenceWarning
When calibration is unresolved: capped/failed fits or nonfinite
selected statistics. All rows remain present and p_value is NaN.
"""
if isinstance(n_replicates, bool):
raise ValueError("n_replicates must be a positive integer.")
try:
n_replicates = index(n_replicates)
except TypeError as error:
raise ValueError("n_replicates must be a positive integer.") from error
if n_replicates < 1:
raise ValueError(f"n_replicates must be >= 1, got {n_replicates}.")
if statistic not in ("G2", "X2"):
raise ValueError("statistic must be 'G2' or 'X2'.")
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
if 8 * n_replicates * (2 * data.n_items + len(structure) + 6) > max_memory_bytes:
raise MemoryError("Bootstrap replicate arrays exceed max_memory_bytes.")
if n_restarts is not None:
n_restarts = _positive_integer(n_restarts, "n_restarts")
_validate_start_options(init_range, init_strategy)
if method not in ("ML", "MDML"):
raise ValueError("Random restarts require ML or MDML.")
if estimate is not None:
raise ValueError(
"Omit estimate with n_restarts so the observed search uses the same policy."
)
elif init_strategy != "uniform" or init_range != (0.01, 0.4):
raise ValueError("Initialization options require n_restarts.")
rng = np.random.default_rng(seed)
def fit_sample(sample: ResponseMatrix) -> BLIMEstimate:
if n_restarts is not None:
assert method in ("ML", "MDML")
return estimate_blim_restarts(
structure,
sample,
method=method,
n_restarts=n_restarts,
seed=rng,
init_strategy=init_strategy,
init_range=init_range,
max_iter=max_iter,
tol=tol,
constraints=constraints,
discrepancy=discrepancy,
max_memory_bytes=max_memory_bytes,
)
return estimate_blim(
structure,
sample,
method=method,
max_iter=max_iter,
tol=tol,
constraints=constraints,
discrepancy=discrepancy,
max_memory_bytes=max_memory_bytes,
)
counts = data.effective_counts
if not np.equal(counts, np.floor(counts)).all():
raise ValueError("Parametric bootstrap requires integer response frequencies, not weights.")
if estimate is None:
estimate = fit_sample(data)
else:
if type(estimate) is not BLIMEstimate:
raise TypeError("bootstrap_gof refits BLIM estimates; use bootstrap_slm for an SLM.")
if discrepancy is None:
discrepancy = estimate.discrepancy
elif discrepancy != estimate.discrepancy:
raise ValueError("discrepancy must match the supplied estimate.")
if constraints is None:
constraints = estimate.constraints
elif constraints != estimate.constraints:
raise ValueError("constraints must match those used for the supplied estimate.")
if estimate.method != method:
raise ValueError(
f"estimate was fitted with method={estimate.method!r} but the "
f"bootstrap refits would use method={method!r}; pass matching "
f"values."
)
if list(estimate.items) != list(data.items):
raise ValueError(
"estimate.items does not match data.items; the estimate must "
"be a fit of this structure to this data."
)
states = list(structure)
if estimate.states != states:
raise ValueError("estimate.states does not match the structure in canonical order.")
if set(data.items) != structure.domain:
raise ValueError("data.items must match the structure's domain.")
pi = np.asarray(estimate.pi)
if (
pi.shape != (len(states),)
or not np.isfinite(pi).all()
or np.any(pi < 0)
or not np.isclose(pi.sum(), 1, rtol=0, atol=1e-8)
):
raise ValueError("estimate.pi must be a probability distribution over the states.")
for name in ("beta", "eta"):
values = np.asarray(getattr(estimate, name))
if (
values.shape != (len(data.items),)
or not np.isfinite(values).all()
or np.any(values < 0)
or np.any(values >= 1)
):
raise ValueError(
f"estimate.{name} must contain one finite value in [0, 1) per item."
)
if constraints is not None:
compiled = _compile_constraints(constraints, list(data.items), states)
projected_beta = estimate.beta.copy()
projected_eta = estimate.eta.copy()
compiled.beta.project(projected_beta, clip_free=False)
compiled.eta.project(projected_eta, clip_free=False)
if (
not np.allclose(projected_beta, estimate.beta, rtol=0, atol=1e-12)
or not np.allclose(projected_eta, estimate.eta, rtol=0, atol=1e-12)
or (
compiled.pi is not None and not np.allclose(compiled.pi, pi, rtol=0, atol=1e-12)
)
or estimate.gof.npar != compiled.n_parameters(len(states))
):
raise ValueError("estimate does not satisfy its declared constraints.")
# Matching labels alone can admit a fit from another dataset.
ll = _log_likelihood_at_params(structure, data, estimate.beta, estimate.eta, pi)
_, inverse = np.unique(data.patterns, axis=0, return_inverse=True)
frequencies = np.bincount(inverse, weights=counts)
frequencies = frequencies[frequencies > 0]
saturated_ll = float(np.sum(frequencies * np.log(frequencies / frequencies.sum())))
g2_check = 2 * (saturated_ll - ll)
if not np.isclose(ll, estimate.log_likelihood, rtol=0, atol=1e-7) or not np.isclose(
g2_check, estimate.gof.G2, rtol=0, atol=1e-7
):
raise ValueError("estimate does not match this data: log-likelihood or G2 differs.")
if not estimate.converged:
warnings.warn(
ConvergenceWarning(
"bootstrap_gof: the supplied observed estimate has not converged; "
"refit the observed data before interpreting the bootstrap p-value."
),
stacklevel=2,
)
n_total = round(float(np.sum(data.effective_counts)))
beta_map = {q: float(b) for q, b in zip(estimate.items, estimate.beta, strict=True)}
eta_map = {q: float(e) for q, e in zip(estimate.items, estimate.eta, strict=True)}
pi_raw = np.array(estimate.pi, dtype=float, copy=True)
pi_raw /= pi_raw.sum()
pi_map = {state: float(p) for state, p in zip(estimate.states, pi_raw, strict=True)}
x2_observed = (
estimate.residuals(data, max_patterns=max_patterns, max_memory_bytes=max_memory_bytes).X2
if statistic == "X2"
else None
)
g2, x2 = np.full(n_replicates, np.nan), np.full(n_replicates, np.nan)
beta_replicates = np.full((n_replicates, len(estimate.items)), np.nan)
eta_replicates = np.full_like(beta_replicates, np.nan)
pi_replicates = np.full((n_replicates, len(estimate.states)), np.nan)
converged_replicates = np.zeros(n_replicates, dtype=bool)
iterations_replicates = np.zeros(n_replicates, dtype=int)
errors: list[str | None] = [None] * n_replicates
n_capped = 0
with warnings.catch_warnings():
# The observed fit already carried any sparse-table or convergence
# warnings; repeating them for every replicate is noise.
warnings.simplefilter("ignore")
for b in range(n_replicates):
sim = simulate_blim(structure, n_total, beta=beta_map, eta=eta_map, pi=pi_map, seed=rng)
assert isinstance(sim, ResponseMatrix) # return_states=False
# Refit in observed item order so random starts attach to the same
# labels on observed and simulated data.
order = [sim.items.index(item) for item in estimate.items]
sim = ResponseMatrix(list(estimate.items), sim.patterns[:, order], sim.counts)
try:
refit = fit_sample(sim)
except ValueError as exc:
if str(exc) not in (
"Observed responses are impossible under the fixed parameters.",
"No possible minimum-discrepancy state for an observed response.",
):
raise
errors[b] = str(exc)
continue
except FloatingPointError as exc:
errors[b] = str(exc)
continue
if not np.isfinite(refit.log_likelihood):
errors[b] = "Non-finite fitted log likelihood."
continue
g2[b] = refit.gof.G2
if statistic == "X2":
x2[b] = refit.residuals(
sim, max_patterns=max_patterns, max_memory_bytes=max_memory_bytes
).X2
# Simulation orders columns canonically; the observed data need not.
order = [refit.items.index(item) for item in estimate.items]
beta_replicates[b] = refit.beta[order]
eta_replicates[b] = refit.eta[order]
pi_replicates[b] = refit.pi
converged_replicates[b] = refit.converged
iterations_replicates[b] = refit.n_iterations
n_capped += not refit.converged
g2_observed = float(estimate.gof.G2)
observed_stat = g2_observed if statistic == "G2" else x2_observed
assert observed_stat is not None
replicates = g2 if statistic == "G2" else x2
n_failed = sum(error is not None for error in errors)
n_extreme = int(np.sum(replicates >= observed_stat - 1e-12))
valid = (
estimate.converged
and not n_failed
and not n_capped
and np.isfinite(observed_stat)
and np.isfinite(replicates).all()
)
p_value = (1 + n_extreme) / (n_replicates + 1) if valid else float("nan")
if not valid:
warnings.warn(
ConvergenceWarning(
f"bootstrap_gof: {n_capped}/{n_replicates} replicate refits "
f"stopped at the iteration cap, {n_failed} failed, or the observed "
f"fit/statistic is unresolved; p-value is NaN. "
f"Inspect convergence and consider raising max_iter={max_iter}."
),
stacklevel=2,
)
for array in (
g2,
x2,
beta_replicates,
eta_replicates,
pi_replicates,
converged_replicates,
iterations_replicates,
):
array.flags.writeable = False
return BootstrapGOF(
p_value=p_value,
g2_observed=g2_observed,
g2_replicates=g2,
n_replicates=n_replicates,
n_extreme=n_extreme,
n_capped=n_capped,
method=method,
seed=seed,
estimate=estimate,
beta_replicates=beta_replicates,
eta_replicates=eta_replicates,
pi_replicates=pi_replicates,
converged_replicates=converged_replicates,
iterations_replicates=iterations_replicates,
statistic=statistic,
x2_observed=x2_observed,
x2_replicates=x2,
n_failed=n_failed,
replicate_errors=tuple(errors),
n_restarts=n_restarts,
init_strategy=init_strategy,
init_range=init_range,
)