"""Saturated observed-data likelihood and BLIM likelihood-ratio comparison.
The saturated model is a common distribution over COMPLETE responses,
observed through each row's mask. No independent mask-distribution factor
or complete-table chi-squared degrees of freedom are assumed.
"""
from __future__ import annotations
import warnings
from dataclasses import dataclass
import numpy as np
from knowledgespaces._patterns import _binary_patterns, _pattern_count, _positive_integer
from knowledgespaces.estimation.blim_em import ConvergenceWarning
from knowledgespaces.estimation.incomplete import IncompleteBLIMEstimate, IncompleteResponseMatrix
[docs]
@dataclass(frozen=True)
class SaturatedIncomplete:
"""Nonparametric ML on a finite response simplex, possibly nonidentified.
complete_patterns is lexicographic in items; probabilities sum to one.
optimality_gap is max_j(dLL/dpi_j)/N - 1 (clamped at zero).
N*gap bounds the likelihood improvement over the current point by
concavity. Converged means gap <= tol, not unique complete probabilities.
Entirely missing rows do not affect the fit. Arrays are read-only.
"""
items: tuple[str, ...]
complete_patterns: np.ndarray
probabilities: np.ndarray
log_likelihood: float
n_iterations: int
converged: bool
optimality_gap: float
objective_history: tuple[float, ...]
[docs]
def saturated_incomplete(
data: IncompleteResponseMatrix,
*,
max_iter: int = 10000,
tol: float = 1e-7,
max_patterns: int = 1_048_576,
max_memory_bytes: int = 512_000_000,
) -> SaturatedIncomplete:
"""Fit a common saturated distribution by multinomial incomplete-data EM.
Uniform positive initialization includes all 2**q completions. Tolerance
controls the concave directional gap, rather than only parameter change.
Enumeration and compatibility arrays are guarded before allocation.
At the iteration cap, return the current fit with a convergence warning.
"""
max_iter = _positive_integer(max_iter, "max_iter")
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
if not np.isfinite(tol) or tol <= 0:
raise ValueError("tol must be finite and positive.")
if data.n_informative == 0:
raise ValueError("No informative responses for the saturated model.")
count = _pattern_count(data.n_items, max_patterns)
if 8 * (4 * data.n_patterns * count + 3 * count * data.n_items + 5 * count) > max_memory_bytes:
raise MemoryError("Saturated incomplete model exceeds max_memory_bytes.")
aggregated = data.aggregate()
keep = aggregated.observed.any(axis=1)
patterns, weights = aggregated.patterns[keep], aggregated.effective_counts[keep]
complete = _binary_patterns(0, count, data.n_items)
compatible = np.ones((len(patterns), count), dtype=bool)
for q in range(data.n_items):
compatible &= np.isnan(patterns[:, q, None]) | (
patterns[:, q, None] == complete[None, :, q]
)
incidence = compatible.astype(float)
prior = np.full(count, 1 / count)
proportions = weights / weights.sum()
history = []
converged = False
for iteration in range(max_iter + 1):
marginal = incidence @ prior
history.append(float(weights @ np.log(marginal)))
score = incidence.T @ (proportions / marginal)
gap = max(float(score.max() - 1), 0.0)
if gap <= tol:
converged = True
break
if iteration < max_iter:
prior *= score
prior /= prior.sum()
if not converged:
warnings.warn(
f"Saturated incomplete fit did not converge in {max_iter} iterations (gap={gap:.3g}).",
ConvergenceWarning,
stacklevel=2,
)
complete.flags.writeable = prior.flags.writeable = False
return SaturatedIncomplete(
data.items, complete, prior, history[-1], iteration, converged, gap, tuple(history)
)
[docs]
@dataclass(frozen=True)
class IncompleteGOF:
"""Observed-data likelihood ratio; no automatic chi-squared reference.
G2 is twice saturated LL minus fitted-model LL on the supplied data.
It is NaN if the approximate saturated solution is worse than the BLIM.
Saturated convergence is retained; its N*optimality_gap bounds the
remaining saturated LL improvement. Use a justified missingness model
when bootstrapping. No mask-frequency Pearson table is constructed.
"""
G2: float
log_likelihood: float
saturated: SaturatedIncomplete
[docs]
def incomplete_gof(
estimate: IncompleteBLIMEstimate,
data: IncompleteResponseMatrix,
*,
max_iter: int = 10000,
tol: float = 1e-7,
max_patterns: int = 1_048_576,
max_memory_bytes: int = 512_000_000,
) -> IncompleteGOF:
"""Compare a supplied fit on these data to their saturated observed model.
Item labels align explicitly; held-out comparisons are descriptive and
are not a training-sample likelihood-ratio test. A nonconverged saturated
fit is retained with a warning, rather than supplying a false exact G2.
"""
sat = saturated_incomplete(
data,
max_iter=max_iter,
tol=tol,
max_patterns=max_patterns,
max_memory_bytes=max_memory_bytes,
)
prediction = estimate.predict(
data.patterns, items=data.items, max_memory_bytes=max_memory_bytes
)
keep = data.effective_counts > 0
ll = float(data.effective_counts[keep] @ prediction.log_probabilities[keep])
g2 = 2 * (sat.log_likelihood - ll)
if g2 < -1e-8:
warnings.warn(
"Saturated likelihood is below the BLIM likelihood; G2 is undefined at this tolerance/cap.",
ConvergenceWarning,
stacklevel=2,
)
g2 = float("nan")
else:
g2 = max(0.0, g2)
return IncompleteGOF(g2, ll, sat)