Bootstrap statistics and observation designs

bootstrap_gof fits a BLIM, simulates from it and refits each generated sample with the same estimator, constraints and discrepancy options. Set statistic="G2" (default) for the likelihood-ratio statistic or "X2" for Pearson’s statistic over the full binary response table, including patterns not observed in the sample. Neither option relies on an automatic asymptotic chi-squared reference distribution.

import numpy as np
from knowledgespaces.estimation import (
    BLIMConstraints, ResponseMatrix, bootstrap_gof,
)
from knowledgespaces.structures import KnowledgeStructure

structure = KnowledgeStructure(["a"], [set(), {"a"}])
data = ResponseMatrix(["a"], np.array([[0], [1]]), np.array([12, 28]))
fixed = BLIMConstraints(
    beta_fixed={"a": 0}, eta_fixed={"a": 0},
    pi_fixed={frozenset(): .5, frozenset({"a"}): .5},
)
boot = bootstrap_gof(
    structure, data, constraints=fixed, statistic="X2",
    n_replicates=5, seed=8,
)
assert np.isclose(boot.statistic_observed, 6.4)
assert boot.n_failed == boot.n_capped == 0
assert np.isfinite(boot.p_value)

This deliberately small, fully specified example checks the mechanics; five replicates are insufficient for a precise empirical analysis. Select the number of replications for the required Monte Carlo resolution and inspect that uncertainty when interpreting results.

statistic_observed and statistic_replicates expose the selected statistic. The existing g2_observed/g2_replicates fields always contain G2, even when Pearson determines the p-value. x2_observed/x2_replicates are populated when requested. Pearson enumeration is bounded by max_patterns and max_memory_bytes; large item domains can exceed this limit.

Refit policy and unresolved fits

Without n_restarts, complete and incomplete BLIM bootstrap use deterministic standard starts for the observed data and all refits. A positive integer uses the random-start policy consistently throughout. Supplying a precomputed complete BLIM estimate together with n_restarts is rejected: omit it to ensure the observed search actually uses the same policy. SLM bootstrap preserves the SLM prior restriction and supports its own matching random-start settings. Resampling requires integer respondent frequencies.

Parameter arrays, scores, convergence and iteration counts keep every replicate in order. Numerical failures retain NaN rows and their messages; they are not replaced until a requested number of successful samples is obtained. A capped selected fit remains visible. p_value is NaN when an observed or selected replicate fit is unresolved, or a required statistic is nonfinite. Otherwise it is the add-one tail count (1 + number at least as extreme) / (1 + n_replicates).

descriptive_tail_fraction preserves that count when all statistics are finite but some fits are capped. It is provided for inspection, not as calibrated inference. This corrects the earlier complete BLIM/SLM behavior, which returned a numeric p-value despite an iteration-cap warning. Earlier results with capped fits require reassessment before publication. Increasing the cap and rerunning the whole declared procedure is preferable to dropping the problematic replicates. Convergence still does not prove a global optimum.

Weighted empirical independent masks

For incomplete data, sample_observation_masks(data, n, seed=...) samples whole boolean observation masks, in data.items order. It pools data.effective_counts by mask, so aggregated and individual data represent the same law; zero-weight masks cannot be sampled. Fractional source weights are allowed as probability weights. Missingness dependence between items within each mask is preserved.

simulate_blim_incomplete(structure, masks, ...) then applies these masks to independently generated complete responses. The composition models an empirical MCAR mechanism. It does not infer MAR or MNAR from observed NaN values. bootstrap_blim_incomplete(mask_model="empirical_independent") provides this workflow directly; "fixed" instead preserves original mask stratum sizes, while the empirical option resamples those sizes.

from knowledgespaces.estimation import (
    IncompleteResponseMatrix, sample_observation_masks,
    simulate_blim_incomplete, bootstrap_blim_incomplete,
)

missing = IncompleteResponseMatrix(
    ("a",), np.array([[0.], [1.], [np.nan]]), np.array([12, 28, 10])
)
masks = sample_observation_masks(missing, 50, seed=31)
simulated = simulate_blim_incomplete(structure, masks, beta=0, eta=0, seed=32)
np.testing.assert_array_equal(simulated.observed, masks)
result = bootstrap_blim_incomplete(
    structure, missing, constraints=fixed,
    mask_model="empirical_independent", statistic="X2",
    n_replicates=5, seed=31,
)
assert result.observed_gof is None  # no saturated model needed for Pearson
assert np.isfinite(result.p_value)

Pearson with missing responses: a conditional claim

independent_mask_pearson(fit, data) uses a stronger assumption than the generic observed-response likelihood. For a mask \(O\), its fitted frequency is \(N_O/N\) and a value-and-mask cell has joint probability

\[p(O,r_O) = (N_O/N)P_\theta(r_O).\]

Under independent masks, these are disjoint cells of a joint multinomial distribution. Equivalently, condition on exogenous mask strata of size \(N_O\) and sum their Pearson statistics. In both cases

\[X^2 = \sum_{(O,r_O):n>0}\frac{n_{O,r_O}^2}{N_O P_\theta(r_O)} - N.\]

The normalization identity includes unobserved cells without enumerating them. Entirely missing strata contribute zero. This is the factorization used by the reference MATLAB toolbox’s blim.m (prm = pr.*pm, followed by model.chi); its validity does not extend to arbitrary MAR merely because an observed-data likelihood is ignorable. No chi-squared degrees of freedom are assumed here, and no native MATLAB execution is claimed.

Consequently, bootstrap_blim_incomplete(statistic="X2") accepts only "fixed" and "empirical_independent" masks. It rejects arbitrary missingness callbacks, which can depend on responses and violate this factorization. G2 remains available with an explicitly justified callback. X2 needs no saturated fit: observed_gof=None, G2 replicate entries are NaN and saturated iteration counts are zero. The G2 default retains its observed-likelihood saturated comparison and convergence checks.

The toolbox’s local bootstrap.m mask sampler counts rows of model.miss without the response-pattern frequencies. Python pools those frequencies explicitly, preserving the empirical mask law for aggregated data. It also retains failed/capped refits and uses an add-one tail count rather than discarding runs and interpolating an empirical distribution. These are documented differences, not numerical parity claims across random streams.