"""Observed-data BLIM likelihood for incomplete binary responses.
Independent marginal-likelihood EM. Ignorability requires MAR and distinct
response/missingness parameters (Rubin, 1976); a NaN input is not evidence
for these assumptions. No missing-not-at-random mechanism is fitted here.
"""
from __future__ import annotations
import warnings
from collections.abc import Sequence
from dataclasses import dataclass
from typing import TYPE_CHECKING
import numpy as np
from scipy.special import logsumexp
from knowledgespaces._patterns import _positive_integer
from knowledgespaces.estimation.blim_em import (
ConvergenceWarning,
ResponseMatrix,
_degenerate_items,
_warn_degenerate,
)
from knowledgespaces.estimation.constraints import BLIMConstraints, _compile_constraints
from knowledgespaces.estimation.prediction import (
BLIMPrediction,
ItemParameter,
StateParameter,
_item_vector,
_state_matrix,
_state_vector,
)
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
if TYPE_CHECKING:
from knowledgespaces.estimation.incomplete_restarts import IncompleteBLIMRestart
[docs]
class UnobservedItemWarning(UserWarning):
"""At least one item has no observed responses; individual errors lack information."""
[docs]
@dataclass(frozen=True, init=False)
class IncompleteResponseMatrix:
"""Immutable 0/1/NaN responses with optional nonnegative frequency weights.
Rows may be individual or aggregated. NaN is the only missing marker;
zero remains an observed incorrect response. Aggregation includes the
missingness mask. Arrays and labels are independent snapshots. Complete
conversion rejects NaN instead of silently deleting or imputing rows.
"""
items: tuple[str, ...]
patterns: np.ndarray
counts: np.ndarray | None = None
def __init__(
self, items: Sequence[str], patterns: np.ndarray, counts: np.ndarray | None = None
) -> None:
items = tuple(items)
patterns = np.array(patterns, dtype=float, copy=True)
if not items or len(items) != len(set(items)):
raise ValueError("items must be nonempty and unique.")
if patterns.ndim != 2 or patterns.shape[1] != len(items) or not len(patterns):
raise ValueError("patterns must be nonempty, 2D, with one column per item.")
if not (np.isin(patterns, [0, 1]) | np.isnan(patterns)).all():
raise ValueError("patterns must contain only 0, 1 and NaN.")
counts = None if counts is None else np.array(counts, dtype=float, copy=True)
if counts is not None:
if (
counts.shape != (len(patterns),)
or not np.isfinite(counts).all()
or (counts < 0).any()
):
raise ValueError("counts must be a finite nonnegative vector matching rows.")
if not 0 < counts.sum() < np.inf:
raise ValueError("Total count must be positive and finite.")
counts.flags.writeable = False
patterns.flags.writeable = False
object.__setattr__(self, "items", items)
object.__setattr__(self, "patterns", patterns)
object.__setattr__(self, "counts", counts)
@property
def observed(self) -> np.ndarray:
"""Boolean mask: True means an observed response."""
return ~np.isnan(self.patterns)
@property
def n_items(self) -> int:
return len(self.items)
@property
def n_patterns(self) -> int:
return len(self.patterns)
@property
def effective_counts(self) -> np.ndarray:
return np.ones(self.n_patterns) if self.counts is None else self.counts
@property
def n_respondents(self) -> float:
return float(self.effective_counts.sum())
@property
def n_informative(self) -> float:
"""Total weight of rows with at least one observed item."""
return float(self.effective_counts[self.observed.any(axis=1)].sum())
[docs]
def aggregate(self) -> IncompleteResponseMatrix:
"""Combine identical values AND masks; remove zero-frequency rows."""
encoded, inverse = np.unique(
np.nan_to_num(self.patterns, nan=2), axis=0, return_inverse=True
)
counts = np.bincount(inverse, weights=self.effective_counts)
encoded[encoded == 2] = np.nan
keep = counts > 0
return IncompleteResponseMatrix(self.items, encoded[keep], counts[keep])
[docs]
def to_complete(self) -> ResponseMatrix:
"""Return an independent complete matrix, only when every cell is observed."""
if np.isnan(self.patterns).any():
raise ValueError("Missing responses remain; select complete rows explicitly first.")
return ResponseMatrix(
list(self.items),
self.patterns.astype(int),
None if self.counts is None else self.counts.copy(),
)
def _log_observed(patterns: np.ndarray, correct: np.ndarray) -> np.ndarray:
yes, no = patterns == 1, patterns == 0
with np.errstate(divide="ignore"):
log_yes, log_no = np.log(correct), np.log1p(-correct)
if np.isfinite(log_yes).all() and np.isfinite(log_no).all():
return yes @ log_yes.T + no @ log_no.T
result = np.zeros((len(patterns), len(correct)))
for q in range(patterns.shape[1]):
result += np.where(yes[:, q, None], log_yes[None, :, q], 0)
result += np.where(no[:, q, None], log_no[None, :, q], 0)
return result
[docs]
def predict_blim_incomplete(
structure: KnowledgeStructure,
patterns: np.ndarray,
*,
items: list[str] | tuple[str, ...] | None = None,
beta: ItemParameter = 0.1,
eta: ItemParameter = 0.1,
pi: StateParameter = None,
max_memory_bytes: int = 512_000_000,
) -> BLIMPrediction:
"""Marginalize missing cells and condition states on the observed cells.
Probabilities are cylinder-event probabilities, not joint probabilities
of response AND missingness. Rows with different masks can overlap and
must not be summed as disjoint cells. All-missing rows have probability
one and the prior posterior. Impossible events have NaN posteriors.
Parameter vectors follow items, state vectors follow canonical order.
"""
data = IncompleteResponseMatrix(
tuple(sorted(structure.domain)) if items is None else tuple(items), patterns
)
if set(data.items) != structure.domain:
raise ValueError("items must match the structure domain.")
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
n, q, m = data.n_patterns, data.n_items, len(structure)
if 8 * (6 * n * m + 5 * n * q + 4 * m * q) > max_memory_bytes:
raise MemoryError("Incomplete prediction exceeds max_memory_bytes.")
labels, states = list(data.items), list(structure)
b, e = _item_vector(beta, labels, "beta"), _item_vector(eta, labels, "eta")
prior, membership = _state_vector(pi, states), _state_matrix(labels, states)
with np.errstate(divide="ignore"):
joint = _log_observed(data.patterns, membership * (1 - b) + (1 - membership) * e) + np.log(
prior
)
logp = logsumexp(joint, axis=1)
posterior = np.full(joint.shape, np.nan)
possible = np.isfinite(logp)
posterior[possible] = np.exp(joint[possible] - logp[possible, None])
# These events carry no information; preserve exact identity, including
# fixed-zero prior masses, instead of rounding log(sum(pi)).
empty = ~data.observed.any(axis=1)
logp[empty], posterior[empty] = 0.0, prior
probability = np.exp(logp)
for array in (probability, logp, posterior):
array.flags.writeable = False
return BLIMPrediction(data.items, tuple(states), probability, logp, posterior)
[docs]
@dataclass(frozen=True)
class IncompleteBLIMEstimate:
"""ML fit to observed cells, without an automatic chi-squared GOF.
npar is the nominal free-coordinate count, not an identified dimension.
n_informative excludes wholly missing rows, which contribute LL=0.
AIC is conventional; BIC uses n_informative and is descriptive under
nonidentifiability/unequal observation designs. For a likelihood-ratio
comparison use incomplete_gof; do not use complete-table residuals.
"""
items: tuple[str, ...]
states: tuple[frozenset[str], ...]
beta: np.ndarray
eta: np.ndarray
pi: np.ndarray
log_likelihood: float
npar: int
n_respondents: float
n_informative: float
n_iterations: int
converged: bool
objective_history: tuple[float, ...]
constraints: BLIMConstraints | None
unobserved_items: tuple[str, ...]
degenerate_items: tuple[str, ...]
restarts: tuple[IncompleteBLIMRestart, ...] = ()
selected_restart: int | None = None
@property
def AIC(self) -> float:
return -2 * self.log_likelihood + 2 * self.npar
@property
def BIC(self) -> float:
return -2 * self.log_likelihood + np.log(self.n_informative) * self.npar
def beta_dict(self) -> dict[str, float]:
return dict(zip(self.items, self.beta.tolist(), strict=True))
def eta_dict(self) -> dict[str, float]:
return dict(zip(self.items, self.eta.tolist(), strict=True))
def pi_dict(self) -> dict[frozenset[str], float]:
return dict(zip(self.states, self.pi.tolist(), strict=True))
def predict(
self,
patterns: np.ndarray,
*,
items: list[str] | tuple[str, ...] | None = None,
max_memory_bytes: int = 512_000_000,
) -> BLIMPrediction:
return predict_blim_incomplete(
KnowledgeStructure(self.items, self.states),
patterns,
items=self.items if items is None else items,
beta=self.beta_dict(),
eta=self.eta_dict(),
pi=self.pi_dict(),
max_memory_bytes=max_memory_bytes,
)
def _validate_incomplete_fit(
structure: KnowledgeStructure,
data: IncompleteResponseMatrix,
max_iter: int,
tol: float,
max_memory_bytes: int,
) -> tuple[int, int]:
"""Validate shared single-fit/search inputs before drawing random starts."""
if not isinstance(data, IncompleteResponseMatrix):
raise TypeError("data must be an IncompleteResponseMatrix.")
if set(data.items) != structure.domain:
raise ValueError("items must match the structure domain.")
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: all positive-frequency rows are missing.")
n, q, m = data.n_patterns, data.n_items, len(structure)
if 8 * (10 * n * m + 8 * n * q + 4 * m * q) > max_memory_bytes:
raise MemoryError("Incomplete BLIM exceeds max_memory_bytes.")
return max_iter, max_memory_bytes
[docs]
def estimate_blim_incomplete(
structure: KnowledgeStructure,
data: IncompleteResponseMatrix,
*,
beta_init: ItemParameter = 0.1,
eta_init: ItemParameter = 0.1,
pi_init: StateParameter = None,
constraints: BLIMConstraints | None = None,
max_iter: int = 5000,
tol: float = 1e-7,
max_memory_bytes: int = 512_000_000,
) -> IncompleteBLIMEstimate:
"""Observed-data ML with fixed/shared errors and optional fixed state prior.
Only observed responses enter error numerators AND denominators. Rows
entirely missing do not alter the fit. Never-observed items are warned
about; unconstrained errors without information retain their starts.
Stopping uses maximum absolute parameter change; a converged EM need
not find a unique or global optimum. Free pi_init must be strictly
positive to avoid locking out states; fixed priors may contain zeros.
No MD/MDML or MNAR model is implied.
"""
max_iter, max_memory_bytes = _validate_incomplete_fit(
structure, data, max_iter, tol, max_memory_bytes
)
m = len(structure)
items, states = list(data.items), list(structure)
compiled = _compile_constraints(constraints or BLIMConstraints(), items, states)
beta, eta = (
_item_vector(beta_init, items, "beta_init"),
_item_vector(eta_init, items, "eta_init"),
)
compiled.beta.project(beta)
compiled.eta.project(eta)
prior = _state_vector(pi_init, states) if compiled.pi is None else compiled.pi.copy()
if compiled.pi is None and np.any(prior == 0):
raise ValueError(
"Free pi_init must be strictly positive; use pi_fixed for structural zeros."
)
# Aggregation by values+mask is sufficient under the common response model.
aggregated = data.aggregate()
informative = aggregated.observed.any(axis=1)
patterns = aggregated.patterns[informative]
counts = aggregated.effective_counts[informative]
proportions = counts / counts.sum()
observed = ~np.isnan(patterns)
responses = np.nan_to_num(patterns, nan=0.0)
membership = _state_matrix(items, states)
unobserved = tuple(q for i, q in enumerate(items) if not observed[:, i].any())
if unobserved:
warnings.warn(
f"No observed responses for items {unobserved}; individual error parameters lack information.",
UnobservedItemWarning,
stacklevel=2,
)
history = []
def expectation() -> np.ndarray:
with np.errstate(divide="ignore"):
joint = _log_observed(
patterns, membership * (1 - beta) + (1 - membership) * eta
) + np.log(prior)
normalizer = logsumexp(joint, axis=1)
if not np.isfinite(normalizer).all():
raise ValueError(
"Observed responses are impossible under the fixed/initial parameters."
)
history.append(float(counts @ normalizer))
return np.exp(joint - normalizer[:, None]) * proportions[:, None]
weights = expectation()
converged = False
for _iteration in range(1, max_iter + 1):
previous = np.concatenate((beta, eta, prior))
compiled.update(
weights,
responses,
membership,
beta,
eta,
eps_on_empty=False,
observed=observed,
minimum_mass=0.0,
)
if compiled.pi is None:
prior = weights.sum(axis=0)
prior /= prior.sum()
weights = expectation()
if np.max(abs(np.concatenate((beta, eta, prior)) - previous)) < tol:
converged = True
break
if not converged:
warnings.warn(
f"Incomplete BLIM did not converge in {max_iter} iterations.",
ConvergenceWarning,
stacklevel=2,
)
degenerate = tuple(_degenerate_items(items, beta, eta))
_warn_degenerate(degenerate)
for array in (beta, eta, prior):
array.flags.writeable = False
return IncompleteBLIMEstimate(
tuple(items),
tuple(states),
beta,
eta,
prior,
history[-1],
compiled.n_parameters(m),
data.n_respondents,
data.n_informative,
_iteration,
converged,
tuple(history),
constraints,
unobserved,
degenerate,
)