"""
EM algorithm for BLIM parameter estimation.
Estimates slip (β) and guess (η) parameters from observed response
patterns using Expectation-Maximization. Per-item parameters are the
default; BLIMConstraints supplies fixed values and homogeneous or grouped
parameterizations.
The algorithm iterates between:
- E-step: compute posterior P(state | response pattern) for each pattern.
- M-step: re-estimate β, η, and state prior π from sufficient statistics.
References:
Doignon, J.-P., & Falmagne, J.-C. (1999).
Knowledge Spaces, Chapter 7. Springer-Verlag.
Falmagne, J.-C., & Doignon, J.-P. (2011).
Learning Spaces, Chapter 11. Springer-Verlag.
Heller, J., & Wickelmaier, F. (2013).
Minimum discrepancy estimation in probabilistic knowledge structures.
Electronic Notes in Discrete Mathematics, 42, 49-56.
"""
from __future__ import annotations
import warnings
from dataclasses import dataclass
from typing import TYPE_CHECKING, Literal
import numpy as np
from numpy.typing import NDArray
from scipy import stats
from scipy.special import logsumexp
from knowledgespaces.estimation._data_signature import _data_signature
from knowledgespaces.estimation._starts import (
_random_error_starts,
_random_prior_start,
_validate_start_options,
)
from knowledgespaces.estimation.constraints import BLIMConstraints, _compile_constraints
from knowledgespaces.estimation.discrepancy import (
MDOptions,
_discrepancy_weights,
_validate_options,
)
from knowledgespaces.estimation.prediction import (
BLIMPrediction,
StateParameter,
_log_conditional,
_positive_integer,
_state_vector,
predict_blim,
)
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
if TYPE_CHECKING:
from knowledgespaces.estimation.diagnostics import BLIMResiduals
from knowledgespaces.metrics.reliability import BLIMReliability
[docs]
class ConvergenceWarning(UserWarning):
"""Emitted when an iterative estimator fails to meet its convergence criterion."""
[docs]
class SparseGOFWarning(UserWarning):
"""The response table is sparse; assess the chi-squared approximation.
Emitted when the sample size is smaller than the number of possible
response patterns minus one (``N < 2^|Q| - 1``), the regime in which
the ``pks`` degrees-of-freedom convention caps the saturated
dimension at ``N``. This is a screening rule, not a theorem that every
sparse-table approximation fails. Expected cell counts and model
regularity matter (Koehler & Larntz, 1980). Treat ``p_value`` as
descriptive and consider :func:`bootstrap_gof` for calibration.
"""
[docs]
@dataclass
class ResponseMatrix:
"""Observed response patterns from a group of respondents.
Parameters
----------
items : list[str]
Item labels (columns), must match the structure's domain.
patterns : np.ndarray
Binary matrix of shape (n_respondents, n_items).
patterns[r, q] = 1 if respondent r answered item q correctly.
counts : np.ndarray or None
Optional frequency for each unique pattern. If None, each
row in patterns is one respondent (count=1 each).
"""
items: list[str]
patterns: np.ndarray
counts: np.ndarray | None = None
def __post_init__(self) -> None:
# Item uniqueness
seen: set[str] = set()
for item in self.items:
if item in seen:
raise ValueError(f"Duplicate item label: '{item}'.")
seen.add(item)
if self.patterns.ndim != 2:
raise ValueError(f"patterns must be 2D, got {self.patterns.ndim}D.")
if self.patterns.shape[0] == 0:
raise ValueError("patterns must contain at least one response pattern.")
if self.patterns.shape[1] != len(self.items):
raise ValueError(
f"patterns has {self.patterns.shape[1]} columns but {len(self.items)} items."
)
if not np.isin(self.patterns, [0, 1]).all():
raise ValueError("patterns must contain only 0 and 1.")
if self.counts is not None:
if self.counts.ndim != 1:
raise ValueError("counts must be one-dimensional.")
if not np.isfinite(self.counts).all():
raise ValueError("counts must contain only finite values.")
if len(self.counts) != self.patterns.shape[0]:
raise ValueError(
f"counts length {len(self.counts)} != patterns rows {self.patterns.shape[0]}."
)
if np.any(self.counts < 0):
raise ValueError("counts must be non-negative.")
if self.counts.sum() == 0:
raise ValueError("Total count must be positive (at least one respondent).")
if not np.isfinite(self.counts.sum()):
raise ValueError("Total count must be finite.")
@property
def n_patterns(self) -> int:
return self.patterns.shape[0]
@property
def n_items(self) -> int:
return len(self.items)
@property
def effective_counts(self) -> np.ndarray:
"""Counts for each pattern (ones if not provided)."""
if self.counts is not None:
return self.counts
return np.ones(self.n_patterns)
@property
def n_respondents(self) -> float:
return float(self.effective_counts.sum())
[docs]
@dataclass(frozen=True)
class GoodnessOfFit:
"""Goodness-of-fit statistics for a BLIM estimate.
Follows the approach of Heller & Wickelmaier (2013, Electronic
Notes in Discrete Mathematics) and the conventions of the R
``pks`` package (Wickelmaier, Heller & Mollenhauer).
The primary statistic is the likelihood ratio G2 (deviance),
tested against a chi-squared distribution.
Attributes
----------
G2 : float
Likelihood ratio statistic: 2 * sum_r N_r ln(N_r / E_r).
df : int
Degrees of freedom: ``max(min(2^Q - 1, N) - npar, 0)``.
The ``min(2^Q - 1, N)`` cap follows the ``pks`` convention:
it is a software convention, not the dimension of the saturated
multinomial model, which remains 2^Q - 1. Sparse observations
do not change population identifiability.
p_value : float
P-value from the chi-squared reference for G2. Asymptotic and
descriptive: when the pattern table is sparse relative to the
sample size (``N < 2^Q - 1``), the reference distribution is
unreliable and a :class:`SparseGOFWarning` is emitted;
:func:`bootstrap_gof` provides a simulated reference under the
fitted model and refit policy, without a general calibration
guarantee for singular models.
npar : int
Number of free parameters: ``|K| - 1 + 2 * Q``.
AIC : float
Akaike Information Criterion: ``-2*LL + 2*npar``.
BIC : float
Bayesian Information Criterion using the total sample size:
``-2*LL + ln(N)*npar``. This is the standard BIC definition of
Schwarz (1978, p. 461), where ``N`` is the number of independent
observations contributing to the likelihood. For BLIM each of
the ``N`` respondents supplies one i.i.d. draw from the pattern
distribution, so the Laplace-approximation derivation of the
``log(N)·npar`` penalty applies *in regular models*. BLIM
selection problems are often nonregular — estimates on the
boundary of the parameter space, rank-deficient structures —
and the standard consistency argument then needs qualification,
so read AIC and BIC as conventional summaries rather than
oracles; predictive comparisons (e.g., cross-validated
log-scores) are the more defensible primary criterion.
BIC_npatterns : float
Variant Bayesian Information Criterion using the number of
distinct observed response patterns: ``-2*LL + ln(n_patterns)*npar``.
This matches what R ``pks::blim()`` returns: pks does not define
an explicit ``BIC`` method, instead overriding ``nobs.blim`` to
return the count of distinct patterns and delegating to
``stats::BIC`` (see ``cran/pks/R/blim.R``, ``logLik.blim`` /
``nobs.blim``). Provided for cross-package replication; not
recommended as a primary selection criterion because the count
of distinct patterns is bounded above by ``2^Q`` and therefore
does not satisfy the asymptotic-consistency conditions of
Schwarz (1978).
"""
G2: float
df: int
p_value: float
npar: int
AIC: float
BIC: float
BIC_npatterns: float
[docs]
@dataclass(frozen=True)
class BLIMEstimate:
"""Result of BLIM parameter estimation via EM.
Attributes
----------
beta : np.ndarray
Slip parameters, shape (n_items,). beta[q] = P(incorrect | q mastered).
eta : np.ndarray
Guess parameters, shape (n_items,). eta[q] = P(correct | q not mastered).
pi : np.ndarray
State prior probabilities, shape (n_states,).
log_likelihood : float
Final log-likelihood of the data.
n_iterations : int
Number of EM iterations until convergence.
converged : bool
True if converged within max_iter.
items : list[str]
Item labels corresponding to beta/eta indices.
states : list[frozenset[str]]
Knowledge states corresponding to pi indices (same order).
gof : GoodnessOfFit
Goodness-of-fit statistics (G2, df, p-value, AIC, BIC).
degenerate_items : tuple[str, ...]
Items whose final ``beta[q] + eta[q] >= 1 - 1e-3``. Such items are
flagged near, at, or beyond the non-discriminating boundary.
Equality beta + eta = 1 means independence of mastery; larger
sums reverse discrimination. The tolerance also flags small
positive discrimination. Inspect coding and fit stability before
changing the structure or removing items.
method : str
Estimation method that produced this estimate: ``"ML"`` (maximum
likelihood via EM), ``"MD"`` (minimum discrepancy), or ``"MDML"``
(maximum likelihood restricted to minimum-discrepancy state
assignments); see Heller & Wickelmaier (2013).
"""
beta: np.ndarray
eta: np.ndarray
pi: np.ndarray
log_likelihood: float
n_iterations: int
converged: bool
items: list[str]
states: list[frozenset[str]]
gof: GoodnessOfFit
degenerate_items: tuple[str, ...]
method: str = "ML"
constraints: BLIMConstraints | None = None
discrepancy: MDOptions | None = None
data_signature: str | None = None
"""Fingerprint of the training response measure; None for older/manual fits."""
[docs]
def residuals(
self,
data: ResponseMatrix,
*,
chunk_size: int = 1024,
max_patterns: int = 1_048_576,
max_memory_bytes: int = 512_000_000,
) -> BLIMResiduals:
"""Full-table Pearson/deviance residuals on training or held-out data."""
from knowledgespaces.estimation.diagnostics import blim_residuals
return blim_residuals(
self,
data,
chunk_size=chunk_size,
max_patterns=max_patterns,
max_memory_bytes=max_memory_bytes,
)
[docs]
def predict(
self,
patterns: np.ndarray,
*,
items: list[str] | None = None,
method: Literal["ML", "MD", "MDML"] = "ML",
discrepancy: MDOptions | None = None,
inclusion: np.ndarray | None = None,
max_memory_bytes: int = 512_000_000,
) -> BLIMPrediction:
"""Predict responses and states; columns default to training item order.
Prediction defaults to ML (Bayes), independently of the estimation
method. MD/MDML inherit fitted discrepancy options unless overridden
or an explicit inclusion mask is supplied. Hyperbolic MD fits require
a minimum-rule override for MDML prediction. Fixed-zero constraints
are preserved. See :func:`predict_blim` for assignment semantics.
"""
return predict_blim(
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(),
method=method,
discrepancy=(
self.discrepancy
if method != "ML" and discrepancy is None and inclusion is None
else discrepancy
),
inclusion=inclusion,
constraints=self.constraints,
max_memory_bytes=max_memory_bytes,
)
[docs]
def reliability(
self,
*,
chunk_size: int = 1024,
max_patterns: int = 1_048_576,
max_memory_bytes: int = 512_000_000,
tie_tolerance: float = 1e-12,
) -> BLIMReliability:
"""Full-test reliability at this fit; parameter uncertainty is omitted."""
from knowledgespaces.metrics.reliability import blim_reliability
return blim_reliability(
KnowledgeStructure(self.items, self.states),
beta=self.beta_dict(),
eta=self.eta_dict(),
pi=self.pi_dict(),
chunk_size=chunk_size,
max_patterns=max_patterns,
max_memory_bytes=max_memory_bytes,
tie_tolerance=tie_tolerance,
)
[docs]
def beta_for(self, item: str) -> float:
"""Get beta (slip) for a specific item."""
return float(self.beta[self.items.index(item)])
[docs]
def eta_for(self, item: str) -> float:
"""Get eta (guess) for a specific item."""
return float(self.eta[self.items.index(item)])
[docs]
def beta_dict(self) -> dict[str, float]:
"""Return beta as {item: value} dict."""
return dict(zip(self.items, self.beta.tolist(), strict=True))
[docs]
def eta_dict(self) -> dict[str, float]:
"""Return eta as {item: value} dict."""
return dict(zip(self.items, self.eta.tolist(), strict=True))
[docs]
def pi_dict(self) -> dict[frozenset[str], float]:
"""Return pi as {state: probability} dict."""
return dict(zip(self.states, self.pi.tolist(), strict=True))
def _compute_gof(
log_likelihood: float,
n_states: int,
n_items: int,
n_respondents: float,
beta: np.ndarray,
eta: np.ndarray,
pi: np.ndarray,
S: np.ndarray,
R: np.ndarray,
counts: np.ndarray,
*,
npar_override: int | None = None,
) -> GoodnessOfFit:
"""Compute goodness-of-fit statistics for a BLIM estimate.
Follows the R ``pks`` package conventions (Heller & Wickelmaier, 2013).
G2 is computed over *unique* response patterns with their observed
frequencies, compared against expected frequencies from the model.
"""
N = n_respondents
# Aggregate to unique patterns with observed frequencies
unique_patterns, inverse = np.unique(R, axis=0, return_inverse=True)
observed_freq = np.zeros(len(unique_patterns))
for i, idx in enumerate(inverse):
observed_freq[idx] += counts[i]
n_unique = int(np.count_nonzero(observed_freq > 0))
# Predicted probability of each unique pattern: P(R) = sum_K P(R|K)*P(K)
p_correct = S * (1 - beta) + (1 - S) * eta # (n_states, n_items)
log_pc = np.log(np.clip(p_correct, 1e-300, None))
log_pinc = np.log(np.clip(1 - p_correct, 1e-300, None))
log_lik_rk = unique_patterns @ log_pc.T + (1 - unique_patterns) @ log_pinc.T
log_prior = np.log(np.clip(pi, 1e-300, None))
log_joint = log_lik_rk + log_prior
max_log = log_joint.max(axis=1, keepdims=True)
log_P_R = (max_log + np.log(np.exp(log_joint - max_log).sum(axis=1, keepdims=True))).ravel()
if np.any(p_correct == 0) or np.any(p_correct == 1) or np.any(pi == 0):
with np.errstate(divide="ignore"):
log_joint = _log_conditional(unique_patterns, p_correct) + np.log(pi)
log_P_R = logsumexp(log_joint, axis=1)
# G2 = 2 * sum_r N_r * ln(N_r / E_r), skipping zero-count patterns
# Keep expected counts in log space: exponentiation can underflow for
# very unlikely patterns even when their log probability is finite.
mask = observed_freq > 0
fractions = observed_freq[mask] / N
log_fractions = np.log(observed_freq[mask]) - np.log(N)
# Taking log(f/N) directly avoids subtracting two nearly equal logs.
# Keep the difference-of-logs fallback for subnormal analysis weights.
np.log(fractions, out=log_fractions, where=fractions > 0)
G2 = float(2 * np.sum(observed_freq[mask] * (log_fractions - log_P_R[mask])))
# Free parameters: |K|-1 (pi) + Q (beta) + Q (eta)
npar = (n_states - 1) + 2 * n_items if npar_override is None else npar_override
# Degrees of freedom: the saturated model has one free parameter per
# observable cell minus the sum-to-1 constraint, i.e. 2^Q - 1. When
# the sample size N is smaller, pks caps this at N. This is a software
# convention, not a change in the saturated model's identifiability.
n_possible = 2**n_items
n_saturated = min(n_possible - 1, round(N))
if round(N) < n_possible - 1:
warnings.warn(
SparseGOFWarning(
f"G2 p-value: the sample size N={round(N)} is below the "
f"2^{n_items} - 1 = {n_possible - 1} possible response "
f"patterns, so the saturated dimension is capped at N (pks "
f"convention). Assess the chi-squared approximation "
f"(Koehler & Larntz, 1980). Treat p_value as descriptive; "
f"consider bootstrap_gof for model-based calibration."
),
stacklevel=3,
)
if n_items > 20:
warnings.warn(
f"Domain has {n_items} items — the full 2^Q pattern space "
f"({n_possible}) is very large. GOF df may be unreliable.",
stacklevel=3,
)
raw_df = int(n_saturated - npar)
df = max(raw_df, 0)
# P-value. When raw_df <= 0 the model is saturated or over-parametrized
# (more free parameters than identifiable cells), so the chi-squared
# reference distribution is undefined and no finite p-value exists; the
# function emits a warning and returns NaN.
if raw_df > 0:
p_value = float(stats.chi2.sf(G2, df))
else:
warnings.warn(
f"BLIM goodness-of-fit: npar={npar} free parameters for only "
f"{n_saturated} dimensions under the pks sample-size convention "
f"(df={raw_df} <= 0); the chi-squared p-value is undefined "
f"(returning NaN). This count alone does not establish population "
f"saturation or nonidentifiability.",
stacklevel=3,
)
p_value = float("nan")
# Information criteria.
# Primary BIC follows Schwarz (1978): N = number of independent
# observations (here, total respondents). pks::blim() instead returns
# log(n_patterns)*npar via nobs.blim override; that variant is exposed
# as BIC_npatterns for cross-package replication but is not the
# asymptotically consistent BIC because n_patterns is bounded by 2^Q.
AIC = -2 * log_likelihood + 2 * npar
BIC = -2 * log_likelihood + np.log(N) * npar # Schwarz (1978)
BIC_npatterns = -2 * log_likelihood + np.log(n_unique) * npar # pks variant
return GoodnessOfFit(
G2=G2,
df=df,
p_value=p_value,
npar=npar,
AIC=float(AIC),
BIC=float(BIC),
BIC_npatterns=float(BIC_npatterns),
)
def _log_likelihood_at_params(
structure: KnowledgeStructure,
data: ResponseMatrix,
beta: np.ndarray,
eta: np.ndarray,
pi: np.ndarray,
) -> float:
"""Compute the BLIM log-likelihood at the given parameters.
Mirrors the E-step log-marginal used inside :func:`estimate_blim`.
Used to keep ``BLIMEstimate.log_likelihood`` consistent with the
returned ``beta/eta/pi`` in the non-convergence branch, and exposed
for regression testing.
"""
items = data.items
states = sorted(structure.states, key=lambda s: (len(s), sorted(s)))
n_items = len(items)
n_states = len(states)
item_idx = {item: i for i, item in enumerate(items)}
S = np.zeros((n_states, n_items), dtype=np.float64)
for k, state in enumerate(states):
for item in state:
S[k, item_idx[item]] = 1.0
R = data.patterns.astype(np.float64)
counts = data.effective_counts
p_correct = S * (1 - beta) + (1 - S) * eta
if np.any(p_correct == 0) or np.any(p_correct == 1) or np.any(pi == 0):
with np.errstate(divide="ignore"):
log_joint = _log_conditional(R, p_correct) + np.log(pi)
log_marginal = logsumexp(log_joint, axis=1)
positive = counts > 0
return float((counts[positive] * log_marginal[positive]).sum())
log_p_correct = np.log(np.clip(p_correct, 1e-300, None))
log_p_incorrect = np.log(np.clip(1 - p_correct, 1e-300, None))
log_lik_rk = R @ log_p_correct.T + (1 - R) @ log_p_incorrect.T
log_joint = log_lik_rk + np.log(np.clip(pi, 1e-300, None))
max_log = log_joint.max(axis=1, keepdims=True)
log_marginal = max_log + np.log(np.exp(log_joint - max_log).sum(axis=1, keepdims=True))
return float((counts * log_marginal.ravel()).sum())
_DEGENERATE_TOL = 1e-3
"""Slack on the ``beta + eta < 1`` informativeness condition. Items whose
final ``beta[q] + eta[q] >= 1 - _DEGENERATE_TOL`` are flagged as degenerate."""
def _degenerate_items(items: list[str], beta: np.ndarray, eta: np.ndarray) -> tuple[str, ...]:
"""Return items whose ``beta + eta`` is at the degenerate boundary."""
threshold = 1.0 - _DEGENERATE_TOL
return tuple(items[q] for q in range(len(items)) if beta[q] + eta[q] >= threshold)
def _update_item_params(
W: np.ndarray,
R: np.ndarray,
S: np.ndarray,
beta: np.ndarray,
eta: np.ndarray,
*,
eps_on_empty: bool,
) -> None:
"""M-step re-estimation of ``beta``/``eta`` (in place) from weights ``W``.
``W[r, k]`` is the (expected) number of respondents with pattern ``r``
assigned to state ``k``. Per-coordinate clip into ``[eps, 1 - eps]``
mirrors ``pks::blim()``; the canonical BLIM parameter space is the open
box (0, 1) per item independently — the joint constraint
``beta + eta < 1`` is the "informative item" condition discussed in the
identifiability literature (Stefanutti, Heller, Anselmi & Robusto 2012;
Spoto, Stefanutti & Vidotto 2013), not part of the parameter space, so
it is *not* enforced here. Enforcing it would keep both the closed
form and the EM ascent property: the order-restricted M-step pools
the two success probabilities of an offending item into their common
success rate, an exact constrained maximization. The pooled update,
however, parks the item exactly on the uninformative boundary
``beta + eta = 1``, so the violation is left visible instead. Items
violating it at convergence are surfaced
via :attr:`BLIMEstimate.degenerate_items`.
When a denominator is (numerically) zero the coordinate is left
unchanged if ``eps_on_empty`` is False (ML behaviour, keeps the current
iterate), or set to ``eps`` if True (MD/MDML behaviour, mirroring the
``NA -> eps`` rule of ``pks``).
"""
eps = 1e-6
n_items = S.shape[1]
N_total = W.sum()
for q in range(n_items):
states_with_q = S[:, q] # (n_states,) binary
mastered_mass = W @ states_with_q # (n_patterns,)
N_mastered_q = mastered_mass.sum()
N_not_mastered_q = N_total - N_mastered_q
C_mastered_q = (R[:, q] * mastered_mass).sum()
C_not_mastered_q = (R[:, q] * (W @ (1 - states_with_q))).sum()
if N_mastered_q > 1e-10:
beta[q] = np.clip(1 - C_mastered_q / N_mastered_q, eps, 1 - eps)
elif eps_on_empty:
beta[q] = eps
if N_not_mastered_q > 1e-10:
eta[q] = np.clip(C_not_mastered_q / N_not_mastered_q, eps, 1 - eps)
elif eps_on_empty:
eta[q] = eps
def _warn_degenerate(degenerate: tuple[str, ...]) -> None:
"""Emit a ConvergenceWarning if any item is at the degenerate boundary."""
if not degenerate:
return
warnings.warn(
f"BLIM fit has degenerate items {list(degenerate)} "
f"(beta + eta >= 1 - {_DEGENERATE_TOL:g}). Such items are "
f"near, at, or beyond the non-discriminating boundary. Inspect "
f"response coding, fit stability, and the assumed structure; the "
f"flag alone does not establish why discrimination is weak or reversed.",
category=ConvergenceWarning,
stacklevel=3,
)
[docs]
def estimate_blim(
structure: KnowledgeStructure,
data: ResponseMatrix,
*,
method: Literal["ML", "MD", "MDML"] = "ML",
max_iter: int = 500,
tol: float = 1e-6,
beta_init: float | np.ndarray = 0.1,
eta_init: float | np.ndarray = 0.1,
pi_init: StateParameter = None,
constraints: BLIMConstraints | None = None,
discrepancy: MDOptions | None = None,
max_memory_bytes: int = 8_000_000_000,
) -> BLIMEstimate:
"""Estimate BLIM parameters by ML (EM), MD, or MDML.
Parameters
----------
structure : KnowledgeStructure
The knowledge structure defining valid states.
data : ResponseMatrix
Observed response patterns.
method : {"ML", "MD", "MDML"}
Estimation method (Heller & Wickelmaier 2013; same options as
``pks::blim()``). ``"ML"`` (default) maximizes the likelihood via
EM. ``"MD"`` is the non-iterative minimum-discrepancy estimator:
each observed pattern is assigned (uniformly) to the knowledge
states requiring the fewest response errors to explain it, and the
parameters are read off that assignment in a single M-step —
deterministic given the data, independent of ``beta_init``/
``eta_init``/``tol``. ``"MDML"`` runs EM with the posterior
restricted to the minimum-discrepancy states of each pattern
(maximum likelihood among minimum-discrepancy solutions).
max_iter : int
Maximum number of EM iterations (ML/MDML). Default 500.
tol : float
Convergence tolerance. For ``"ML"`` this is the absolute change in
log-likelihood between iterations; for ``"MDML"`` it is the maximum
absolute change in the parameters ``(pi, beta, eta)``, mirroring
the stopping rule of ``pks::blim()``. Ignored for ``"MD"``.
Default 1e-6.
beta_init : float or np.ndarray
Initial slip values. A scalar initializes every item equally;
it does not impose equality during estimation. Free initial errors
are clipped to [1e-6, 1 - 1e-6], as in the subsequent M-steps.
Values must be finite and in [0, 1); beta_init + eta_init may
equal or exceed 1, consistently with the fitted parameter space.
eta_init : float or np.ndarray
Initial guess values, finite and in [0, 1). Equality during fitting
requires constraints. No joint informative-item inequality is imposed.
pi_init : Mapping or np.ndarray or None
Initial state distribution (canonical order); None uses uniform.
Ignored by MD. Zeros are absorbing in EM; use a positive start to
allow every state positive fitted mass.
constraints : BLIMConstraints or None
Fixed and shared error parameters, optionally a fully fixed pi.
Grouped updates pool sufficient statistics. The free parameter
count excludes fixed classes and accounts for equality classes.
A fixed pi is supported by ML/MDML, not the MD assignment rule.
Zero fixed error rates exclude incompatible assignments exactly.
discrepancy : MDOptions or None
Optional radius for MD/MDML, or a hyperbolic inclusion rule for MD.
None uses ordinary minimum discrepancy. Not applicable to ML.
max_memory_bytes : int
Preflight limit on estimated peak array allocation, including
response-by-state and response-by-item intermediates. This is an
allocation estimate, not an OS-level memory cap. Default 8 GB.
Returns
-------
BLIMEstimate
Estimated parameters, log-likelihood, and convergence info.
``log_likelihood`` is always the *unrestricted* BLIM log-likelihood
evaluated at the returned parameters (also for MD/MDML, matching
``pks``).
Raises
------
ValueError
If data items don't match structure domain, or if init
parameters are out of range.
MemoryError
If the estimated posterior allocation would exceed
``max_memory_bytes``.
Notes
-----
The M-step independently clips ``beta[q]`` and ``eta[q]`` into
``[1e-6, 1 - 1e-6]``, mirroring R ``pks::blim()``. The canonical BLIM
parameter space is the open box per item; the joint condition
``beta[q] + eta[q] < 1`` is the *informative item* condition of the
identifiability literature (Stefanutti, Heller, Anselmi & Robusto
2012; Spoto, Stefanutti & Vidotto 2013), not part of the parameter
space, and is therefore *not* enforced inside the loop. This is a
design choice, not a necessity: the order-restricted M-step retains
a closed form (the two success probabilities of an offending item
pool into their common success rate) and, being an exact constrained
maximization, preserves the EM ascent property. The pooled update,
however, parks the item exactly on the uninformative boundary
``beta + eta = 1``, and enforcing the strict inequality would
require an arbitrary margin, so the violation is left visible.
Items whose final ``beta[q] + eta[q] >= 1 - 1e-3`` are surfaced via
:attr:`BLIMEstimate.degenerate_items` and a :class:`ConvergenceWarning`
is emitted. This is a diagnostic of weak, absent or reversed
discrimination, not a sufficient reason to remove an item.
References
----------
Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum
likelihood from incomplete data via the EM algorithm. *J. R. Stat.
Soc. B*, 39(1), 1-38.
Heller, J., & Wickelmaier, F. (2013). Minimum discrepancy estimation
in probabilistic knowledge structures. *ENDM*, 42, 49-56.
Stefanutti, L., Heller, J., Anselmi, P., & Robusto, E. (2012).
Assessing the local identifiability of probabilistic knowledge
structures. *Behavior Research Methods*, 44(4), 1197-1211.
Spoto, A., Stefanutti, L., & Vidotto, G. (2013). Considerations about
the identification of forward- and backward-graded knowledge
structures. *Journal of Mathematical Psychology*, 57(5), 249-254.
"""
md_options = _validate_options(discrepancy, method)
# Validate hyperparameters
if method not in ("ML", "MD", "MDML"):
raise ValueError(f"method must be 'ML', 'MD', or 'MDML', got {method!r}.")
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(f"tol must be > 0, got {tol}.")
# Validate domain match
if set(data.items) != structure.domain:
raise ValueError(
f"ResponseMatrix items {set(data.items)} don't match "
f"structure domain {set(structure.domain)}."
)
items = list(data.items)
n_items = len(items)
counts = data.effective_counts
# Build state matrix: S[k, q] = 1 if item q is in state k
states = sorted(structure.states, key=lambda s: (len(s), sorted(s)))
n_states = len(states)
compiled = None if constraints is None else _compile_constraints(constraints, items, states)
if method == "MD" and compiled is not None and compiled.pi is not None:
raise ValueError("pi_fixed is supported by ML/MDML, not the assignment-based MD estimator.")
npar = None if compiled is None else compiled.n_parameters(n_states)
# Preflight: at its peak the E/M-step holds up to nine
# (n_patterns × n_states) float64 matrices at once (log_lik_rk,
# log_joint, posterior, W, the goodness-of-fit likelihood matrix, and
# broadcasting temporaries — several names stay bound for the whole
# iteration), on top of about three (n_patterns × n_items) arrays (the
# float64 pattern matrix, its int64 source copy, and the 1 - R
# temporary). For powerset structures over large domains this explodes
# to tens of GB and the OS kills the process without a useful message.
# The multipliers are calibrated as an upper bound on tracemalloc peaks
# across sparse (|K| ~ |Q|) and powerset (|K| = 2^|Q|) shapes (see
# tests/test_estimation.py::TestEMMemoryPreflight); fail fast above a
# hard cap (default 8 GB), warn above 1 GB.
estimated_bytes = data.n_patterns * (9 * n_states + 3 * n_items) * 8
if estimated_bytes > max_memory_bytes:
raise MemoryError(
f"BLIM EM would allocate ~{estimated_bytes / 1e9:.2f} GB for the "
f"E-step matrices ({data.n_patterns} patterns x [9 x {n_states} "
f"states + 3 x {n_items} items] x 8 bytes). Exceeds "
f"max_memory_bytes={max_memory_bytes / 1e9:.2f} GB. Reduce domain "
f"size, use a sparser knowledge structure, or pass "
f"max_memory_bytes=... to override."
)
if estimated_bytes > 1_000_000_000:
warnings.warn(
f"BLIM EM will allocate ~{estimated_bytes / 1e9:.2f} GB for the "
f"E-step matrices ({data.n_patterns} patterns x {n_states} "
f"states).",
category=ResourceWarning,
stacklevel=2,
)
item_idx = {item: i for i, item in enumerate(items)}
S = np.zeros((n_states, n_items), dtype=np.float64)
for k, state in enumerate(states):
for item in state:
S[k, item_idx[item]] = 1.0
# Initialize and validate parameters
if isinstance(beta_init, (int, float)):
beta: NDArray[np.float64] = np.full(n_items, float(beta_init))
else:
beta = np.array(beta_init, dtype=np.float64).copy()
if beta.shape != (n_items,):
raise ValueError(f"beta_init array has shape {beta.shape}, expected ({n_items},).")
if isinstance(eta_init, (int, float)):
eta: NDArray[np.float64] = np.full(n_items, float(eta_init))
else:
eta = np.array(eta_init, dtype=np.float64).copy()
if eta.shape != (n_items,):
raise ValueError(f"eta_init array has shape {eta.shape}, expected ({n_items},).")
if not np.all(np.isfinite(beta)) or np.any(beta < 0) or np.any(beta >= 1):
raise ValueError("All beta_init values must be finite and in [0, 1).")
if not np.all(np.isfinite(eta)) or np.any(eta < 0) or np.any(eta >= 1):
raise ValueError("All eta_init values must be finite and in [0, 1).")
pi = _state_vector(pi_init, states)
if compiled is None:
# A zero starting guess/slip is not a fixed-zero model hypothesis.
# Keep free starts inside the same numerical box as free M-steps.
np.clip(beta, 1e-6, 1 - 1e-6, out=beta)
np.clip(eta, 1e-6, 1 - 1e-6, out=eta)
if compiled is not None:
compiled.beta.project(beta)
compiled.eta.project(eta)
if compiled.pi is not None:
pi = compiled.pi.copy()
R = data.patterns.astype(np.float64) # (n_patterns, n_items)
# Fixed-zero parameters restrict feasibility before distance minimization.
i_rk = None
if method in ("MD", "MDML"):
feasible = None
if compiled is not None:
correct = S * (1 - beta) + (1 - S) * eta
feasible = np.isfinite(_log_conditional(R, correct))
if method == "MDML" and compiled.pi is not None:
feasible &= compiled.pi > 0
i_rk = _discrepancy_weights(R, S, counts, md_options, feasible)
if method == "MD":
# Non-iterative (pks::blimMD): assign each pattern uniformly to its
# minimum-discrepancy states and read the parameters off a single
# M-step. Deterministic given the data; beta_init/eta_init/tol play
# no role in the result.
assert i_rk is not None
m_rk = i_rk / i_rk.sum(axis=1, keepdims=True)
W = counts[:, np.newaxis] * m_rk
pi = W.sum(axis=0) / counts.sum()
if compiled is None:
_update_item_params(W, R, S, beta, eta, eps_on_empty=True)
else:
compiled.update(W, R, S, beta, eta, eps_on_empty=True)
ll = _log_likelihood_at_params(structure, data, beta, eta, pi)
gof = _compute_gof(
ll,
n_states,
n_items,
data.n_respondents,
beta,
eta,
pi,
S,
R,
counts,
npar_override=npar,
)
degenerate = _degenerate_items(items, beta, eta)
_warn_degenerate(degenerate)
return BLIMEstimate(
beta=beta,
eta=eta,
pi=pi,
log_likelihood=ll,
n_iterations=1,
converged=True,
items=items,
states=states,
gof=gof,
degenerate_items=degenerate,
method="MD",
constraints=constraints,
discrepancy=discrepancy,
data_signature=_data_signature(data),
)
prev_ll = -np.inf
for iteration in range(1, max_iter + 1):
# ============================================================
# E-step: compute P(pattern | state) and posterior P(state | pattern)
# ============================================================
# P(R_r | K_k) = prod_q [ P(R_rq | K_k, q) ]
# where P(correct | q in K) = 1-beta_q, P(correct | q not in K) = eta_q
# Likelihood per (pattern, state, item)
# P(R_rq=1 | K_k) = S[k,q]*(1-beta[q]) + (1-S[k,q])*eta[q]
p_correct = S * (1 - beta) + (1 - S) * eta # (n_states, n_items)
p_incorrect = 1 - p_correct # (n_states, n_items)
# P(R_r | K_k) for each (pattern, state)
# log for numerical stability
log_p_correct = np.log(np.clip(p_correct, 1e-300, None))
log_p_incorrect = np.log(np.clip(p_incorrect, 1e-300, None))
# log P(R_r | K_k) = sum_q [ R_rq * log_p_correct[k,q] + (1-R_rq) * log_p_incorrect[k,q] ]
log_lik_rk = R @ log_p_correct.T + (1 - R) @ log_p_incorrect.T # (n_patterns, n_states)
# Add log prior
log_joint = log_lik_rk + np.log(np.clip(pi, 1e-300, None)) # (n_patterns, n_states)
# Log-sum-exp for numerical stability
max_log = log_joint.max(axis=1, keepdims=True)
log_marginal = max_log + np.log(np.exp(log_joint - max_log).sum(axis=1, keepdims=True))
# Posterior: P(K_k | R_r) = exp(log_joint - log_marginal)
posterior = np.exp(log_joint - log_marginal) # (n_patterns, n_states)
if np.any(p_correct == 0) or np.any(p_correct == 1) or np.any(pi == 0):
with np.errstate(divide="ignore"):
log_joint = _log_conditional(R, p_correct) + np.log(pi)
log_marginal = logsumexp(log_joint, axis=1, keepdims=True)
possible = np.isfinite(log_marginal[:, 0])
if np.any((counts > 0) & ~possible):
raise ValueError("Observed responses are impossible under the fixed parameters.")
posterior = np.zeros(log_joint.shape)
posterior[possible] = np.exp(log_joint[possible] - log_marginal[possible])
# MDML: restrict the posterior to the minimum-discrepancy states of
# each pattern and renormalize (pks::blimEM with md = em = 1).
if i_rk is not None:
posterior *= i_rk
restricted_mass = posterior.sum(axis=1, keepdims=True)
posterior /= np.maximum(restricted_mass, 1e-300)
# A nearest state can be arbitrarily unlikely under the current
# prior/errors. Normalize its restricted likelihood in log space
# before exponentiation when the ordinary posterior underflows.
tiny = restricted_mass[:, 0] < 1e-200
if np.any(tiny):
restricted_log = np.where(i_rk[tiny] > 0, log_joint[tiny], -np.inf)
normalizer = logsumexp(restricted_log, axis=1, keepdims=True)
possible = np.isfinite(normalizer[:, 0])
if np.any((counts[tiny] > 0) & ~possible):
raise ValueError(
"No possible minimum-discrepancy state for an observed response."
)
normalized = np.zeros(restricted_log.shape)
normalized[possible] = np.exp(restricted_log[possible] - normalizer[possible])
posterior[tiny] = normalized
# Log-likelihood (unrestricted BLIM marginal, also under MDML)
ll = float((counts[counts > 0] * log_marginal.ravel()[counts > 0]).sum())
# Check convergence (ML: absolute log-likelihood change; MDML
# converges on parameter change after the M-step below)
if method == "ML" and abs(ll - prev_ll) < tol:
gof = _compute_gof(
ll,
n_states,
n_items,
data.n_respondents,
beta,
eta,
pi,
S,
R,
counts,
npar_override=npar,
)
degenerate = _degenerate_items(items, beta, eta)
_warn_degenerate(degenerate)
return BLIMEstimate(
beta=beta,
eta=eta,
pi=pi,
log_likelihood=ll,
n_iterations=iteration,
converged=True,
items=items,
states=states,
gof=gof,
degenerate_items=degenerate,
method="ML",
constraints=constraints,
discrepancy=discrepancy,
data_signature=_data_signature(data),
)
prev_ll = ll
# ============================================================
# M-step: re-estimate beta, eta, pi from sufficient statistics
# (clipping policy documented in _update_item_params)
# ============================================================
# Weighted posterior: w[r, k] = counts[r] * posterior[r, k]
W = counts[:, np.newaxis] * posterior # (n_patterns, n_states)
N_total = counts.sum()
pi_old, beta_old, eta_old = pi, beta.copy(), eta.copy()
# Pi: state prior
pi = W.sum(axis=0) / N_total # (n_states,)
if compiled is None:
_update_item_params(W, R, S, beta, eta, eps_on_empty=(method == "MDML"))
else:
if compiled.pi is not None:
pi = compiled.pi.copy()
compiled.update(W, R, S, beta, eta, eps_on_empty=(method == "MDML"))
# MDML converges on the maximum absolute parameter change, mirroring
# pks::blimEM (the unrestricted log-likelihood is not the objective
# being maximized under the minimum-discrepancy restriction).
if method == "MDML":
max_diff = max(
float(np.abs(pi - pi_old).max()),
float(np.abs(beta - beta_old).max()),
float(np.abs(eta - eta_old).max()),
)
if max_diff < tol:
final_ll = _log_likelihood_at_params(structure, data, beta, eta, pi)
gof = _compute_gof(
final_ll,
n_states,
n_items,
data.n_respondents,
beta,
eta,
pi,
S,
R,
counts,
npar_override=npar,
)
degenerate = _degenerate_items(items, beta, eta)
_warn_degenerate(degenerate)
return BLIMEstimate(
beta=beta,
eta=eta,
pi=pi,
log_likelihood=final_ll,
n_iterations=iteration,
converged=True,
items=items,
states=states,
gof=gof,
degenerate_items=degenerate,
method="MDML",
constraints=constraints,
discrepancy=discrepancy,
data_signature=_data_signature(data),
)
criterion = "|Δ log-likelihood|" if method == "ML" else "max |Δ (pi, beta, eta)|"
warnings.warn(
f"BLIM {method} did not converge in {max_iter} iterations "
f"({criterion} > tol = {tol:.3e}). "
f"Consider increasing max_iter, using estimate_blim_restarts, "
f"or inspecting the knowledge structure for identifiability issues.",
category=ConvergenceWarning,
stacklevel=2,
)
# Recompute log-likelihood on the post-M-step parameters so that the
# returned `log_likelihood` is consistent with the returned beta/eta/pi.
# Inside the loop, `ll` is computed at the E-step (before that
# iteration's M-step), so at loop exit `prev_ll` reflects the params
# from one M-step ago — off-by-one vs. the returned parameters.
final_ll = _log_likelihood_at_params(structure, data, beta, eta, pi)
gof = _compute_gof(
final_ll,
n_states,
n_items,
data.n_respondents,
beta,
eta,
pi,
S,
R,
counts,
npar_override=npar,
)
degenerate = _degenerate_items(items, beta, eta)
_warn_degenerate(degenerate)
return BLIMEstimate(
beta=beta,
eta=eta,
pi=pi,
log_likelihood=final_ll,
n_iterations=max_iter,
converged=False,
items=items,
states=states,
gof=gof,
degenerate_items=degenerate,
method=method,
constraints=constraints,
discrepancy=discrepancy,
data_signature=_data_signature(data),
)
[docs]
def estimate_blim_restarts(
structure: KnowledgeStructure,
data: ResponseMatrix,
*,
method: Literal["ML", "MDML"] = "ML",
n_restarts: int = 10,
max_iter: int = 500,
tol: float = 1e-6,
seed: int | np.random.Generator | None = None,
max_memory_bytes: int = 8_000_000_000,
init_range: tuple[float, float] = (0.01, 0.4),
init_strategy: Literal["uniform", "pks"] = "uniform",
constraints: BLIMConstraints | None = None,
discrepancy: MDOptions | None = None,
) -> BLIMEstimate:
"""Estimate BLIM parameters with multiple random restarts.
Runs :func:`estimate_blim` ``n_restarts`` times with random
initial values for beta, eta (and pi for ``init_strategy="pks"``),
and selects the result with the
highest log-likelihood. This helps avoid local optima.
The R ``pks`` package does not provide this natively — users
must loop manually with ``randinit=TRUE``.
Parameters
----------
structure : KnowledgeStructure
The knowledge structure defining valid states.
data : ResponseMatrix
Observed response patterns.
method : {"ML", "MDML"}
Estimation method forwarded to :func:`estimate_blim`. ``"MD"`` is
rejected here: the minimum-discrepancy estimator is deterministic
given the data, so restarting it is pointless — call
:func:`estimate_blim` with ``method="MD"`` directly instead.
n_restarts : int
Number of random restarts. Default 10.
max_iter : int
Maximum EM iterations per restart.
tol : float
Convergence tolerance per restart.
seed : int, numpy.random.Generator or None
Random seed for reproducibility.
max_memory_bytes : int
Forwarded to :func:`estimate_blim`. Default 8 GB.
discrepancy : MDOptions or None
Forwarded unchanged to every fit. MDML supports a minimum-distance
radius; options are not applicable to ML.
init_range : tuple[float, float]
Lower/upper bounds for the random ``U(low, high)`` draw of
``beta_init`` and ``eta_init`` when ``init_strategy="uniform"``.
Default ``(0.01, 0.4)``. Ignored when ``init_strategy="pks"``.
init_strategy : {"uniform", "pks"}
Random initialization strategy. ``"uniform"`` (default) draws
both parameters from ``U(*init_range)`` and rescales in-place
until ``beta[q] + eta[q] < 0.95`` on each item. ``"pks"`` mirrors
``pks::blim(..., randinit=TRUE)`` (R source:
``cran/pks/R/blim.R``): each parameter is drawn from ``U(0, 1)``,
then reflected as ``1 - x`` sequentially — ``beta`` first on
items where ``beta[q] + eta[q] >= 1``, then ``eta`` on items
where the constraint still fails with the updated ``beta`` —
to restore the informative-item condition. The state prior uses
Dirichlet(1, ..., 1), uniform on the simplex, independently of the
errors. This has the same law as pks's uniform spacings. NumPy and
R random streams differ. The default ``"uniform"`` strategy keeps
the prior at equal state masses. A fixed prior overrides the draw.
Returns
-------
BLIMEstimate
The best result (highest log-likelihood) across all restarts.
Notes
-----
The default ``init_range=(0.01, 0.4)`` is a narrowed basin that
avoids near-boundary draws at the uninformative frontier
``beta + eta = 1``, where EM can stall in a degenerate-item
attractor (cf. Stefanutti, Heller, Anselmi & Robusto 2012). The ``"pks"``
strategy is provided for reproducibility with the R ``pks`` package:
note that pks uses ``runif(nitems)`` = ``U(0, 1)`` with reflection,
**not** ``U(0, 0.5)`` as sometimes reported.
References
----------
Heller, J., & Wickelmaier, F. (2013). Minimum discrepancy estimation
in probabilistic knowledge structures. *ENDM*, 42, 49-56.
"""
if method not in ("ML", "MDML"):
raise ValueError(
f"method must be 'ML' or 'MDML' for restarts, got {method!r}. "
f"The MD estimator is deterministic given the data — call "
f"estimate_blim(method='MD') directly."
)
n_restarts = _positive_integer(n_restarts, "n_restarts")
_validate_start_options(init_range, init_strategy)
rng = np.random.default_rng(seed)
n_items = len(data.items)
best: BLIMEstimate | None = None
n_failed = 0
for restart_i in range(n_restarts):
beta_init, eta_init = _random_error_starts(rng, n_items, init_range, init_strategy)
pi_init = _random_prior_start(rng, len(structure), init_strategy)
# Suppress per-restart ConvergenceWarning: a single aggregate warning
# is emitted below when all restarts fail. Per-restart non-convergence
# is expected during the search and only relevant in aggregate.
# ResourceWarning (memory preflight) is deduplicated to fire only on
# the first restart — the estimate is identical across all restarts.
with warnings.catch_warnings():
warnings.simplefilter("ignore", ConvergenceWarning)
if restart_i > 0:
warnings.simplefilter("ignore", ResourceWarning)
result = estimate_blim(
structure,
data,
method=method,
max_iter=max_iter,
tol=tol,
beta_init=beta_init,
eta_init=eta_init,
pi_init=pi_init,
constraints=constraints,
discrepancy=discrepancy,
max_memory_bytes=max_memory_bytes,
)
if not result.converged:
n_failed += 1
if best is None or result.log_likelihood > best.log_likelihood:
best = result
assert best is not None # n_restarts is positive
if n_failed == n_restarts:
warnings.warn(
f"BLIM EM did not converge in any of the {n_restarts} restarts "
f"(max_iter={max_iter}, tol={tol:.3e}). "
f"The returned estimate is the best (highest log-likelihood) "
f"non-converged fit. Consider increasing max_iter or n_restarts.",
category=ConvergenceWarning,
stacklevel=2,
)
elif not best.converged:
warnings.warn(
"The highest-likelihood restart has not converged, although some other "
"restarts did. Increase max_iter before interpreting the selected fit.",
category=ConvergenceWarning,
stacklevel=2,
)
# Re-emit degenerate-item warning on the *selected* best fit. Per-restart
# warnings are suppressed above to avoid n_restarts duplicates; surfacing
# the diagnostic on the chosen estimate keeps the user informed.
_warn_degenerate(best.degenerate_items)
return best