"""Simple learning model: Falmagne & Doignon (2011), Eq. 11.18, Thm. 11.5.4.
The state mass is the product of g on the state and (1-g) on its outer
fringe. On a learning space these masses sum to one without normalization.
Independent EM implementation using Bernoulli sufficient statistics.
"""
from __future__ import annotations
import warnings
from dataclasses import dataclass, field, replace
from typing import Literal
import numpy as np
from scipy.special import logsumexp
from knowledgespaces._patterns import _positive_integer
from knowledgespaces.estimation._data_signature import _data_signature
from knowledgespaces.estimation._starts import _random_error_starts, _validate_start_options
from knowledgespaces.estimation.blim_em import (
BLIMEstimate,
ConvergenceWarning,
ResponseMatrix,
SparseGOFWarning,
_compute_gof,
_degenerate_items,
_warn_degenerate,
)
from knowledgespaces.estimation.bootstrap import BootstrapGOF
from knowledgespaces.estimation.constraints import BLIMConstraints, _compile_constraints
from knowledgespaces.estimation.discrepancy import (
MDOptions,
_discrepancy_weights,
_validate_options,
)
from knowledgespaces.estimation.identifiability import blim_jacobian
from knowledgespaces.estimation.prediction import (
ItemParameter,
_item_vector,
_log_conditional,
_state_matrix,
)
from knowledgespaces.estimation.simulate import simulate_blim
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
def _slm_matrices(structure: KnowledgeStructure, items: list[str]) -> tuple[np.ndarray, np.ndarray]:
if not structure.is_learning_space:
raise ValueError("The SLM requires a learning space (well-graded and union closed).")
if len(items) != len(set(items)) or set(items) != structure.domain:
raise ValueError("items must match the structure's domain exactly once.")
states = list(structure)
return _state_matrix(items, states), _state_matrix(
items, [structure.outer_fringe(s) for s in states]
)
def _log_prior(g: np.ndarray, membership: np.ndarray, fringe: np.ndarray) -> np.ndarray:
with np.errstate(divide="ignore"):
return np.where(membership == 1, np.log(g), 0).sum(axis=1) + np.where(
fringe == 1, np.log1p(-g), 0
).sum(axis=1)
[docs]
def slm_state_probabilities(
structure: KnowledgeStructure,
g: ItemParameter = 0.1,
*,
items: list[str] | None = None,
) -> np.ndarray:
"""SLM state distribution in canonical state order, including g=0 or 1.
Arrays of g follow ``items`` (sorted domain by default); mappings align
by label. The learning-space hypothesis is checked, not bypassed by
renormalizing arbitrary families. The returned array is read-only.
Solvability g is not generally the marginal probability of mastering q.
"""
items = sorted(structure.domain) if items is None else list(items)
membership, fringe = _slm_matrices(structure, items)
result = np.exp(_log_prior(_item_vector(g, items, "g"), membership, fringe))
result.flags.writeable = False
return result
[docs]
@dataclass(frozen=True)
class SLMRestart:
"""One search run, with effective starts in the estimate's item order.
Starts include fixed/equality projection and the numerical box. The
objective is full log likelihood for ML and the restricted joint sum
for MDML; ``log_likelihood`` always evaluates the full response model.
"""
beta_init: tuple[float, ...]
eta_init: tuple[float, ...]
g_init: tuple[float, ...]
log_likelihood: float
objective: float
n_iterations: int
converged: bool
[docs]
@dataclass(frozen=True)
class SLMEstimate(BLIMEstimate):
"""SLM fit, with BLIM response prediction/residuals and an SLM prior.
``g`` follows ``items``. ``gof.npar`` counts free error groups plus ``|Q|``
solvability parameters; state masses are derived, not freely estimated.
``objective_history`` records the EM criterion before and after updates:
ordinary log likelihood for ML, restricted joint sum for MDML, empty
for single-step MD. Reported log_likelihood/G2 always use the full model.
Use slm_jacobian (or .jacobian()) for this parameterization. A BLIM rank
report with free state masses does not test SLM identifiability.
"""
g: np.ndarray = field(default_factory=lambda: np.empty(0))
objective_history: tuple[float, ...] = ()
restarts: tuple[SLMRestart, ...] = ()
selected_restart: int | None = None
restart_selection: Literal["objective", "likelihood"] | None = None
[docs]
def g_dict(self) -> dict[str, float]:
"""Solvability estimates by item label."""
return dict(zip(self.items, self.g.tolist(), strict=True))
[docs]
def jacobian(self, *, max_memory_bytes: int = 512_000_000) -> np.ndarray:
"""Response-map Jacobian for the fitted SLM, in sorted item order."""
return slm_jacobian(
KnowledgeStructure(self.items, self.states),
g=self.g_dict(),
beta=self.beta_dict(),
eta=self.eta_dict(),
constraints=self.constraints,
max_memory_bytes=max_memory_bytes,
)
[docs]
def estimate_slm(
structure: KnowledgeStructure,
data: ResponseMatrix,
*,
method: Literal["ML", "MD", "MDML"] = "ML",
g_init: ItemParameter = 0.1,
beta_init: ItemParameter = 0.1,
eta_init: ItemParameter = 0.1,
constraints: BLIMConstraints | None = None,
discrepancy: MDOptions | None = None,
max_iter: int = 5000,
tol: float = 1e-7,
max_memory_bytes: int = 512_000_000,
) -> SLMEstimate:
"""Fit the simple learning model by ML, MD or MDML.
The E-step uses BLIM errors and the SLM state distribution. The g update
is expected membership / (expected membership + expected outer-fringe
membership). Error updates pool the same sufficient statistics as BLIM.
ML/MDML stop on max absolute change in beta, eta and g, as pks::slm.
Convergence does not guarantee a global optimum or identifiability.
MD uses one discrepancy assignment and M-step, independent of starts.
Optional MDOptions adds a radius (MD/MDML) or hyperbolic weights (MD).
Fixed/equal beta/eta are supported; pi_fixed is rejected because pi is
determined by g. All g and free errors use the numerical box
[1e-6,1-1e-6]; only fixed errors can be exactly zero. No g constraints
are imposed and no latent-state distribution is silently normalized.
"""
if method not in ("ML", "MD", "MDML"):
raise ValueError("method must be 'ML', 'MD', or 'MDML'.")
options = _validate_options(discrepancy, method)
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.")
items, states = list(data.items), list(structure)
n, q, m = data.n_patterns, data.n_items, len(states)
if 8 * (10 * n * m + 5 * n * q + 6 * m * q) > max_memory_bytes:
raise MemoryError("SLM exceeds max_memory_bytes; reduce the data or structure size.")
membership, fringe = _slm_matrices(structure, items)
compiled = _compile_constraints(constraints or BLIMConstraints(), items, states)
if compiled.pi is not None:
raise ValueError("SLM state probabilities are determined by g; pi_fixed is not supported.")
g = np.clip(_item_vector(g_init, items, "g_init"), 1e-6, 1 - 1e-6)
beta, eta = (
_item_vector(beta_init, items, "beta_init"),
_item_vector(eta_init, items, "eta_init"),
)
compiled.beta.project(beta)
compiled.eta.project(eta)
responses, counts = data.patterns.astype(float), data.effective_counts
inclusion = None
if method != "ML":
correct = membership * (1 - beta) + (1 - membership) * eta
inclusion = _discrepancy_weights(
responses,
membership,
counts,
options,
np.isfinite(_log_conditional(responses, correct)),
)
history: list[float] = []
def expectation() -> np.ndarray:
correct = membership * (1 - beta) + (1 - membership) * eta
joint = _log_conditional(responses, correct) + _log_prior(g, membership, fringe)
if inclusion is not None:
with np.errstate(divide="ignore"):
joint += np.log(inclusion)
normalizer = logsumexp(joint, axis=1)
possible = np.isfinite(normalizer)
if np.any((counts > 0) & ~possible):
raise ValueError("Observed responses are impossible under the fixed parameters.")
history.append(float(counts[counts > 0] @ normalizer[counts > 0]))
weights = np.zeros_like(joint)
weights[possible] = (
np.exp(joint[possible] - normalizer[possible, None]) * counts[possible, None]
)
return weights
if method == "MD":
assert inclusion is not None
weights = counts[:, None] * inclusion / inclusion.sum(axis=1, keepdims=True)
else:
weights = expectation()
converged = False
for _iteration in range(1, (1 if method == "MD" else max_iter) + 1):
previous = np.concatenate((beta, eta, g))
compiled.update(weights, responses, membership, beta, eta, eps_on_empty=True)
mass = weights.sum(axis=0)
mastered, available = mass @ membership, mass @ fringe
g = np.clip(
np.divide(
mastered, mastered + available, out=np.zeros(q), where=mastered + available > 0
),
1e-6,
1 - 1e-6,
)
if method != "MD":
weights = expectation()
if method == "MD" or np.max(np.abs(np.concatenate((beta, eta, g)) - previous)) < tol:
converged = True
break
if not converged:
warnings.warn(
f"SLM {method} did not converge in {max_iter} iterations.",
ConvergenceWarning,
stacklevel=2,
)
prior = np.exp(_log_prior(g, membership, fringe))
log_prob = logsumexp(
_log_conditional(responses, membership * (1 - beta) + (1 - membership) * eta)
+ _log_prior(g, membership, fringe),
axis=1,
)
ll = float(counts[counts > 0] @ log_prob[counts > 0])
npar = compiled.beta.n_free + compiled.eta.n_free + q
gof = _compute_gof(
ll,
m,
q,
data.n_respondents,
beta,
eta,
prior,
membership,
responses,
counts,
npar_override=npar,
)
degenerate = _degenerate_items(items, beta, eta)
_warn_degenerate(degenerate)
for array in (beta, eta, prior, g):
array.flags.writeable = False
return SLMEstimate(
beta=beta,
eta=eta,
pi=prior,
log_likelihood=ll,
n_iterations=_iteration,
converged=converged,
items=items,
states=states,
gof=gof,
degenerate_items=degenerate,
method=method,
constraints=constraints,
discrepancy=discrepancy,
g=g,
objective_history=tuple(history),
data_signature=_data_signature(data),
)
[docs]
def estimate_slm_restarts(
structure: KnowledgeStructure,
data: ResponseMatrix,
*,
method: Literal["ML", "MDML"] = "ML",
n_restarts: int = 10,
seed: int | np.random.Generator | None = None,
init_range: tuple[float, float] = (0.01, 0.4),
init_strategy: Literal["uniform", "pks"] = "uniform",
selection: Literal["objective", "likelihood"] = "objective",
constraints: BLIMConstraints | None = None,
discrepancy: MDOptions | None = None,
max_iter: int = 5000,
tol: float = 1e-7,
max_memory_bytes: int = 512_000_000,
) -> SLMEstimate:
"""Search SLM fits from random error and solvability starts.
``n_restarts=1, init_strategy="pks"`` supplies the initialization law
of pks::slm(randinit=TRUE): two U(0,1) error draws with sequential
reflection, then independent U(0,1) solvabilities. NumPy and R seeds
do not produce identical draws. ``uniform`` uses BLIM's narrowed error
starts (U(low, high), halved together until their sum is below .95);
g still uses U(0,1). init_range is validated but unused for pks.
Fixed/equal errors are projected before the first E-step, as in
estimate_slm; pks does not initially project its equality groups.
Constraints can override the informative starting condition, which
is not imposed on subsequent estimates. Free starts are clipped to
[1e-6,1-1e-6]. All fit settings are forwarded unchanged to each run.
Selection maximizes the optimized criterion by default: likelihood
for ML, restricted joint sum for MDML. ``selection="likelihood"``
instead ranks full response likelihoods (the BLIM restart convention).
First run wins exact ties; unconverged runs are retained and eligible.
The result records every start's diagnostics and the zero-based chosen
index. A selected unconverged fit warns. Restarts offer no global-optimum
or identifiability guarantee. MD is deterministic and is rejected.
"""
if method not in ("ML", "MDML"):
raise ValueError("SLM restarts require ML or MDML; use estimate_slm for deterministic MD.")
n_restarts = _positive_integer(n_restarts, "n_restarts")
_validate_start_options(init_range, init_strategy)
if selection not in ("objective", "likelihood"):
raise ValueError("selection must be 'objective' or 'likelihood'.")
compiled = _compile_constraints(constraints or BLIMConstraints(), data.items, list(structure))
rng = np.random.default_rng(seed)
records: list[SLMRestart] = []
best: SLMEstimate | None = None
best_score = -np.inf
selected = 0
for index in range(n_restarts):
beta, eta = _random_error_starts(rng, data.n_items, init_range, init_strategy)
g = np.clip(rng.uniform(size=data.n_items), 1e-6, 1 - 1e-6)
compiled.beta.project(beta)
compiled.eta.project(eta)
with warnings.catch_warnings():
warnings.simplefilter("ignore", ConvergenceWarning)
if index:
warnings.simplefilter("ignore", SparseGOFWarning)
fit = estimate_slm(
structure,
data,
method=method,
beta_init=beta,
eta_init=eta,
g_init=g,
constraints=constraints,
discrepancy=discrepancy,
max_iter=max_iter,
tol=tol,
max_memory_bytes=max_memory_bytes,
)
objective = fit.objective_history[-1]
records.append(
SLMRestart(
tuple(beta.tolist()),
tuple(eta.tolist()),
tuple(g.tolist()),
fit.log_likelihood,
objective,
fit.n_iterations,
fit.converged,
)
)
score = objective if selection == "objective" else fit.log_likelihood
if best is None or score > best_score:
best, best_score, selected = fit, score, index
assert best is not None
if not best.converged:
failed = sum(not record.converged for record in records)
warnings.warn(
f"Selected SLM restart {selected} has not converged; "
f"{failed}/{n_restarts} runs did not converge. Increase max_iter.",
ConvergenceWarning,
stacklevel=2,
)
_warn_degenerate(best.degenerate_items)
return replace(
best, restarts=tuple(records), selected_restart=selected, restart_selection=selection
)
[docs]
def slm_jacobian(
structure: KnowledgeStructure,
*,
g: ItemParameter = 0.1,
beta: ItemParameter = 0.1,
eta: ItemParameter = 0.1,
constraints: BLIMConstraints | None = None,
max_items: int | None = None,
max_memory_bytes: int = 512_000_000,
) -> np.ndarray:
"""Analytic response Jacobian: free beta groups, eta groups, then all g.
Items are sorted; groups follow blim_jacobian's constraint convention.
The chain rule uses exact polynomial derivatives of Eq. 11.18, including
boundary g values. Rank is a pointwise numerical diagnostic, not a
generic/global proof. Enumeration and memory guards match blim_jacobian.
"""
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
items, states = sorted(structure.domain), list(structure)
extra_bytes = 8 * (3 * len(states) * len(items) + (2 ** len(items)) * 3 * len(items))
if extra_bytes >= max_memory_bytes:
raise MemoryError("SLM Jacobian exceeds max_memory_bytes.")
membership, fringe = _slm_matrices(structure, items)
g_vec = _item_vector(g, items, "g")
prior = np.exp(_log_prior(g_vec, membership, fringe))
compiled = _compile_constraints(constraints or BLIMConstraints(), items, states)
if compiled.pi is not None:
raise ValueError("SLM does not support pi_fixed.")
jac = blim_jacobian(
structure,
beta=beta,
eta=eta,
pi=prior,
constraints=constraints or BLIMConstraints(),
max_items=max_items,
max_memory_bytes=max_memory_bytes - extra_bytes,
)
derivative = np.empty((len(states), len(items)))
for i in range(len(items)):
other = np.arange(len(items)) != i
derivative[:, i] = (membership[:, i] - fringe[:, i]) * np.exp(
_log_prior(g_vec[other], membership[:, other], fringe[:, other])
)
n_errors = compiled.beta.n_free + compiled.eta.n_free
return np.column_stack((jac[:, :n_errors], jac[:, n_errors:] @ derivative[1:]))
[docs]
@dataclass(frozen=True)
class SLMBootstrap(BootstrapGOF):
"""SLM bootstrap with g replicates in addition to beta/eta/pi and G2.
g_replicates follows estimate.items. The inherited parameter_summary
describes beta/eta/derived pi; use g_replicates.std(axis=0, ddof=1) for
solvability dispersion. All refits are retained, including unconverged.
"""
g_replicates: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
[docs]
def bootstrap_slm(
structure: KnowledgeStructure,
data: ResponseMatrix,
*,
method: Literal["ML", "MD", "MDML"] = "ML",
n_replicates: int = 1000,
seed: int | None = None,
n_restarts: int | None = None,
init_range: tuple[float, float] = (0.01, 0.4),
init_strategy: Literal["uniform", "pks"] = "uniform",
selection: Literal["objective", "likelihood"] = "objective",
constraints: BLIMConstraints | None = None,
discrepancy: MDOptions | None = None,
max_iter: int = 5000,
tol: float = 1e-7,
max_memory_bytes: int = 512_000_000,
) -> SLMBootstrap:
"""Fit, simulate and refit the same SLM, preserving its prior restriction.
Integer respondent counts are required. Same settings are used for
observed and bootstrap fits; beta/eta/g starts are the default 0.1.
With n_restarts supplied, each fit instead uses estimate_slm_restarts
with identical search settings and fresh draws from the seeded generator.
Random restarts are unavailable for deterministic MD. Replicate arrays
retain the selected fit of each search, including unconverged selections.
The add-one Monte Carlo p-value and parameter dispersions have the same
model/convergence qualifications as bootstrap_gof. No BLIM with free
state masses is substituted for the SLM during refitting.
max_memory_bytes guards each fit's working arrays; retained replicate
arrays additionally require ``O(n_replicates * (|Q| + |K|))`` storage.
"""
n_replicates = _positive_integer(n_replicates, "n_replicates")
if not np.equal(data.effective_counts, np.floor(data.effective_counts)).all():
raise ValueError("Bootstrap requires integer response frequencies.")
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
if 8 * n_replicates * (3 * data.n_items + len(structure) + 4) > max_memory_bytes:
raise MemoryError("SLM bootstrap replicate arrays exceed max_memory_bytes.")
rng = np.random.default_rng(seed)
def fit_data(sample: ResponseMatrix) -> SLMEstimate:
if n_restarts is not None:
if method not in ("ML", "MDML"):
raise ValueError("SLM bootstrap restarts require ML or MDML.")
return estimate_slm_restarts(
structure,
sample,
method=method,
n_restarts=n_restarts,
seed=rng,
init_range=init_range,
init_strategy=init_strategy,
selection=selection,
constraints=constraints,
discrepancy=discrepancy,
max_iter=max_iter,
tol=tol,
max_memory_bytes=max_memory_bytes,
)
return estimate_slm(
structure,
sample,
method=method,
constraints=constraints,
discrepancy=discrepancy,
max_iter=max_iter,
tol=tol,
max_memory_bytes=max_memory_bytes,
)
fit = fit_data(data)
beta, eta, g = (np.full((n_replicates, data.n_items), np.nan) for _ in range(3))
prior = np.full((n_replicates, len(structure)), np.nan)
g2, iterations = np.full(n_replicates, np.nan), np.zeros(n_replicates, dtype=int)
converged = np.zeros(n_replicates, dtype=bool)
errors: list[str | None] = [None] * n_replicates
with warnings.catch_warnings():
warnings.simplefilter("ignore")
for i in range(n_replicates):
sim = simulate_blim(
structure,
round(data.n_respondents),
beta=fit.beta_dict(),
eta=fit.eta_dict(),
pi=fit.pi_dict(),
seed=rng,
)
assert isinstance(sim, ResponseMatrix)
order = [sim.items.index(q) for q in data.items]
sim = ResponseMatrix(list(data.items), sim.patterns[:, order], sim.counts)
try:
refit = fit_data(sim)
except ValueError as exc:
if str(exc) != "Observed responses are impossible under the fixed parameters.":
raise
errors[i] = str(exc)
continue
except FloatingPointError as exc:
errors[i] = str(exc)
continue
if not np.isfinite(refit.log_likelihood):
errors[i] = "Non-finite fitted log likelihood."
continue
beta[i], eta[i], g[i], prior[i] = refit.beta, refit.eta, refit.g, refit.pi
g2[i], iterations[i], converged[i] = refit.gof.G2, refit.n_iterations, refit.converged
capped = int((~converged & (iterations > 0)).sum())
failed = sum(error is not None for error in errors)
valid = (
fit.converged
and not capped
and not failed
and np.isfinite(fit.gof.G2)
and np.isfinite(g2).all()
)
if not valid:
warnings.warn(
f"SLM bootstrap: {capped}/{n_replicates} refits did not converge, {failed} failed, "
"or the observed fit/statistic is unresolved; p-value is NaN.",
ConvergenceWarning,
stacklevel=2,
)
extreme = int((g2 >= fit.gof.G2 - 1e-12).sum())
for array in (beta, eta, g, prior, g2, iterations, converged):
array.flags.writeable = False
return SLMBootstrap(
(1 + extreme) / (1 + n_replicates) if valid else float("nan"),
fit.gof.G2,
g2,
n_replicates,
extreme,
capped,
method,
seed,
fit,
beta,
eta,
prior,
converged,
iterations,
g,
n_failed=failed,
replicate_errors=tuple(errors),
n_restarts=n_restarts,
init_strategy=init_strategy,
init_range=init_range,
)