"""Simulation and bootstrap with an explicit observation-mask model."""
from __future__ import annotations
import warnings
from collections.abc import Callable
from dataclasses import dataclass
from typing import Literal
import numpy as np
from knowledgespaces._patterns import _positive_integer
from knowledgespaces.estimation.blim_em import ConvergenceWarning
from knowledgespaces.estimation.constraints import BLIMConstraints
from knowledgespaces.estimation.incomplete import (
IncompleteBLIMEstimate,
IncompleteResponseMatrix,
estimate_blim_incomplete,
)
from knowledgespaces.estimation.incomplete_gof import IncompleteGOF, incomplete_gof
from knowledgespaces.estimation.incomplete_restarts import (
IncompleteBLIMRestart,
IncompleteBLIMRestartError,
estimate_blim_incomplete_restarts,
)
from knowledgespaces.estimation.observation_masks import (
independent_mask_pearson,
sample_observation_masks,
)
from knowledgespaces.estimation.prediction import ItemParameter, StateParameter, _item_vector
from knowledgespaces.estimation.simulate import simulate_blim
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
MaskGenerator = Callable[[np.ndarray, tuple[str, ...], np.random.Generator], np.ndarray]
def _mask(values: np.ndarray, shape: tuple[int, ...]) -> np.ndarray:
result = np.asarray(values)
if result.shape != shape or result.dtype != np.bool_:
raise ValueError(
"observed_mask must be a boolean matrix matching respondent rows and item columns."
)
return result
[docs]
def simulate_blim_incomplete(
structure: KnowledgeStructure,
observed_mask: np.ndarray,
*,
items: list[str] | tuple[str, ...] | None = None,
beta: ItemParameter = 0.1,
eta: ItemParameter = 0.1,
pi: StateParameter = None,
seed: int | np.random.Generator | None = None,
max_memory_bytes: int = 512_000_000,
) -> IncompleteResponseMatrix:
"""Generate complete BLIM responses, then apply a fixed exogenous mask.
True indicates an observed cell. Rows remain individual, in mask order;
columns follow items (sorted by default). This is a fixed-design/MCAR
simulation, not a generic MAR missingness model. No values are imputed.
"""
labels = tuple(sorted(structure.domain)) if items is None else tuple(items)
if set(labels) != structure.domain or len(labels) != len(set(labels)):
raise ValueError("items must match the structure domain exactly once.")
mask = np.asarray(observed_mask)
if mask.ndim != 2 or not len(mask):
raise ValueError("observed_mask must be a nonempty 2D matrix.")
mask = _mask(mask, (len(mask), len(labels)))
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
if 8 * (8 * mask.size + 4 * len(structure) * len(labels) + 4 * len(mask)) > max_memory_bytes:
raise MemoryError("Incomplete simulation exceeds max_memory_bytes.")
b, e = _item_vector(beta, list(labels), "beta"), _item_vector(eta, list(labels), "eta")
# The complete simulator uses mappings to preserve caller column order.
simulated = simulate_blim(
structure,
len(mask),
beta=dict(zip(labels, b, strict=True)),
eta=dict(zip(labels, e, strict=True)),
pi=pi,
seed=seed,
return_states=True,
)
assert isinstance(simulated, tuple)
complete = simulated[0]
order = [complete.items.index(q) for q in labels]
patterns = complete.patterns[:, order].astype(float)
patterns[~mask] = np.nan
return IncompleteResponseMatrix(labels, patterns)
[docs]
@dataclass(frozen=True)
class IncompleteBootstrap:
"""All parameter/LR replicates, retaining failed and capped refits.
p_value uses add-one Monte Carlo counting only when the observed fits
and EVERY replicate converge and yield finite LR. Otherwise it is NaN:
failures are not silently dropped. A wholly missing generated sample is
failed (NaN parameter rows); a cap is nonconverged. Arrays follow the
observed estimate's items/states. Standard deviations of parameter
samples are descriptive, not automatically valid confidence intervals.
With statistic="X2", the explicitly independent-mask Pearson statistic
replaces G2 for calibration; saturated fitting is unnecessary and is
skipped (observed_gof=None, G2 NaN, saturated iterations zero). Generic
mask callbacks are disallowed for X2. descriptive_tail_fraction retains
finite counting results separately from unresolved calibrated p-values.
"""
estimate: IncompleteBLIMEstimate
observed_gof: IncompleteGOF | None
g2_replicates: np.ndarray
beta_replicates: np.ndarray
eta_replicates: np.ndarray
pi_replicates: np.ndarray
converged_replicates: np.ndarray
saturated_converged: np.ndarray
iterations_replicates: np.ndarray
saturated_iterations: np.ndarray
n_replicates: int
n_failed: int
n_capped: int
p_value: float
seed: int | None
mask_model: str
n_restarts: int | None = None
init_strategy: Literal["uniform", "pks"] = "uniform"
init_range: tuple[float, float] = (0.01, 0.4)
restart_replicates: tuple[tuple[IncompleteBLIMRestart, ...], ...] = ()
replicate_errors: tuple[str | None, ...] = ()
statistic: Literal["G2", "X2"] = "G2"
x2_observed: float | None = None
x2_replicates: np.ndarray | None = None
@property
def statistic_observed(self) -> float:
"""Observed value of the statistic selected for calibration."""
if self.statistic == "G2":
if self.observed_gof is None:
raise ValueError("This result does not contain a likelihood-ratio statistic.")
return self.observed_gof.G2
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 values of the selected statistic, retaining failed rows."""
if self.statistic == "G2":
return self.g2_replicates
if self.x2_replicates is None:
raise ValueError("This result does not contain Pearson replicates.")
return self.x2_replicates
@property
def descriptive_tail_fraction(self) -> float:
"""Finite add-one tail fraction, not calibrated inference for capped fits."""
values, observed = self.statistic_replicates, self.statistic_observed
if not np.isfinite(observed) or not np.isfinite(values).all():
return float("nan")
return float((1 + (values >= observed - 1e-12).sum()) / (1 + self.n_replicates))
[docs]
def bootstrap_blim_incomplete(
structure: KnowledgeStructure,
data: IncompleteResponseMatrix,
*,
mask_model: Literal["fixed", "empirical_independent"] | MaskGenerator,
n_replicates: int = 1000,
seed: int | None = None,
constraints: BLIMConstraints | None = None,
n_restarts: int | None = None,
init_strategy: Literal["uniform", "pks"] = "uniform",
init_range: tuple[float, float] = (0.01, 0.4),
statistic: Literal["G2", "X2"] = "G2",
max_iter: int = 5000,
tol: float = 1e-7,
saturated_max_iter: int = 10000,
saturated_tol: float = 1e-7,
max_patterns: int = 1_048_576,
max_memory_bytes: int = 512_000_000,
) -> IncompleteBootstrap:
"""Fit/simulate/refit BLIM and saturated model with a stated mask model.
'fixed' expands the original masks by integer respondent frequencies;
masks must be exogenous (fixed design/MCAR), not generic MAR. A callback
takes complete individual responses, item labels and the seeded generator
and returns a boolean observed mask. Its statistical assumptions are the
caller's model; the function checks shape/type, not ignorability. The
callback must not alter responses. Sample failures remain visible.
``n_restarts=None`` uses one deterministic standard start for the
observed sample and each replicate. A positive n_restarts instead uses
estimate_blim_incomplete_restarts with the same initialization law,
constraints, cap and tolerance for every fit, drawing from the single
seeded generator. Restart records are retained for the observed fit and
every replicate; all-start failures retain their errors and NaN rows.
Changing this policy also changes the subsequent seeded simulations.
'empirical_independent' resamples whole masks with probabilities from
their weighted empirical frequencies, independently of responses (MCAR).
'fixed' retains each stratum size; the empirical option resamples sizes.
statistic='X2' calibrates independent_mask_pearson for these two explicit
designs only. Generic mask callbacks may use G2 but cannot imply the
independent-mask factorization needed by Pearson. X2 skips saturated
fitting, so observed_gof is None and G2 arrays remain NaN.
"""
if not (
isinstance(mask_model, str) and mask_model in ("fixed", "empirical_independent")
) and not callable(mask_model):
raise ValueError(
"mask_model must be 'fixed', 'empirical_independent' or a mask-generator callable."
)
if statistic not in ("G2", "X2"):
raise ValueError("statistic must be 'G2' or 'X2'.")
if statistic == "X2" and callable(mask_model):
raise ValueError(
"Pearson requires fixed or empirical_independent masks, not an arbitrary MAR callback."
)
n_replicates = _positive_integer(n_replicates, "n_replicates")
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
counts = data.effective_counts
if not np.equal(counts, np.floor(counts)).all():
raise ValueError("Bootstrap requires integer respondent frequencies.")
n = round(data.n_respondents)
if (
8 * (10 * n * data.n_items + n_replicates * (3 * data.n_items + len(structure) + 6))
> max_memory_bytes
):
raise MemoryError("Incomplete bootstrap exceeds max_memory_bytes.")
if n_restarts is not None:
n_restarts = _positive_integer(n_restarts, "n_restarts")
if (
8 * n_restarts * (n_replicates + 1) * (2 * data.n_items + len(structure) + 6)
> max_memory_bytes
):
raise MemoryError("Incomplete bootstrap restart records exceed max_memory_bytes.")
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 refit_sample(sample: IncompleteResponseMatrix) -> IncompleteBLIMEstimate:
if n_restarts is not None:
return estimate_blim_incomplete_restarts(
structure,
sample,
n_restarts=n_restarts,
seed=rng,
init_strategy=init_strategy,
init_range=init_range,
constraints=constraints,
max_iter=max_iter,
tol=tol,
max_memory_bytes=max_memory_bytes,
)
return estimate_blim_incomplete(
structure,
sample,
constraints=constraints,
max_iter=max_iter,
tol=tol,
max_memory_bytes=max_memory_bytes,
)
fit = refit_sample(data)
observed = (
incomplete_gof(
fit,
data,
max_iter=saturated_max_iter,
tol=saturated_tol,
max_patterns=max_patterns,
max_memory_bytes=max_memory_bytes,
)
if statistic == "G2"
else None
)
x2_observed = (
independent_mask_pearson(fit, data, max_memory_bytes=max_memory_bytes)
if statistic == "X2"
else None
)
fixed_mask = np.repeat(data.observed, counts.astype(np.int64), axis=0)
beta, eta = (np.full((n_replicates, data.n_items), np.nan) for _ in range(2))
prior = np.full((n_replicates, len(structure)), np.nan)
g2 = np.full(n_replicates, np.nan)
x2 = np.full(n_replicates, np.nan)
converged, sat_converged = (np.zeros(n_replicates, dtype=bool) for _ in range(2))
iterations, sat_iterations = (np.zeros(n_replicates, dtype=int) for _ in range(2))
failed = 0
restart_records: list[tuple[IncompleteBLIMRestart, ...]] = [()] * n_replicates
errors: list[str | None] = [None] * n_replicates
for i in range(n_replicates):
simulation_mask = np.ones_like(fixed_mask) if callable(mask_model) else fixed_mask
if mask_model == "empirical_independent":
simulation_mask = sample_observation_masks(
data, n, seed=rng, max_memory_bytes=max_memory_bytes
)
simulation = simulate_blim_incomplete(
structure,
simulation_mask,
items=data.items,
beta=fit.beta_dict(),
eta=fit.eta_dict(),
pi=fit.pi_dict(),
seed=rng,
max_memory_bytes=max_memory_bytes,
)
if callable(mask_model):
masks = _mask(
mask_model(simulation.patterns, data.items, rng), simulation.patterns.shape
)
patterns = simulation.patterns.copy()
patterns[~masks] = np.nan
simulation = IncompleteResponseMatrix(data.items, patterns)
if simulation.n_informative == 0:
failed += 1
errors[i] = "No informative responses in generated sample."
continue
with warnings.catch_warnings():
warnings.simplefilter("ignore", ConvergenceWarning)
try:
refit = refit_sample(simulation)
except IncompleteBLIMRestartError as exc:
restart_records[i] = exc.restarts
errors[i] = str(exc)
failed += 1
continue
restart_records[i] = refit.restarts
gof = (
incomplete_gof(
refit,
simulation,
max_iter=saturated_max_iter,
tol=saturated_tol,
max_patterns=max_patterns,
max_memory_bytes=max_memory_bytes,
)
if statistic == "G2"
else None
)
beta[i], eta[i], prior[i] = refit.beta, refit.eta, refit.pi
if gof is not None:
g2[i] = gof.G2
sat_converged[i], sat_iterations[i] = (
gof.saturated.converged,
gof.saturated.n_iterations,
)
if statistic == "X2":
x2[i] = independent_mask_pearson(refit, simulation, max_memory_bytes=max_memory_bytes)
converged[i], iterations[i] = refit.converged, refit.n_iterations
capped = int(
((iterations > 0) & (~converged | (~sat_converged if statistic == "G2" else False))).sum()
)
statistic_observed = observed.G2 if observed is not None else x2_observed
assert statistic_observed is not None
statistic_replicates = g2 if statistic == "G2" else x2
valid = (
fit.converged
and (observed.saturated.converged if observed is not None else True)
and np.isfinite(statistic_observed)
and not failed
and not capped
and np.isfinite(statistic_replicates).all()
)
p_value = (
float((1 + (statistic_replicates >= statistic_observed - 1e-12).sum()) / (1 + n_replicates))
if valid
else float("nan")
)
if not valid:
warnings.warn(
f"Incomplete bootstrap has {failed} failed and {capped} capped replicates, or an unresolved observed fit/statistic; p-value is NaN.",
ConvergenceWarning,
stacklevel=2,
)
for array in (beta, eta, prior, g2, x2, converged, sat_converged, iterations, sat_iterations):
array.flags.writeable = False
return IncompleteBootstrap(
fit,
observed,
g2,
beta,
eta,
prior,
converged,
sat_converged,
iterations,
sat_iterations,
n_replicates,
failed,
capped,
p_value,
seed,
"callable" if callable(mask_model) else mask_model,
n_restarts,
init_strategy,
init_range,
tuple(restart_records),
tuple(errors),
statistic,
x2_observed,
x2 if statistic == "X2" else None,
)