"""
Simulation of response patterns from a BLIM.
Generates synthetic response data given a knowledge structure and BLIM
parameters (state prior π, careless-error β, lucky-guess η). This is the
counterpart of ``pks::simulate.blim()`` in R and the building block for
parameter-recovery studies.
References:
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
from collections.abc import Mapping
import numpy as np
from numpy.typing import NDArray
from knowledgespaces.estimation.blim_em import ResponseMatrix
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
def _param_vector(
value: float | Mapping[str, float],
items: list[str],
name: str,
) -> np.ndarray:
"""Expand a scalar or per-item mapping into an array aligned with items."""
if isinstance(value, Mapping):
missing = [q for q in items if q not in value]
if missing:
raise ValueError(f"{name} is missing items {missing}.")
extra = set(value) - set(items)
if extra:
raise ValueError(f"{name} has unknown items {sorted(extra)}.")
vec = np.array([float(value[q]) for q in items], dtype=np.float64)
else:
vec = np.full(len(items), float(value), dtype=np.float64)
# Finiteness first: NaN passes every comparison below (NaN < x is
# always False) and would silently produce all-incorrect responses
# for the affected item instead of an error.
if not np.all(np.isfinite(vec)) or np.any(vec < 0) or np.any(vec >= 1):
raise ValueError(f"All {name} values must be finite and in [0, 1).")
return vec
[docs]
def simulate_blim(
structure: KnowledgeStructure,
n_respondents: int,
*,
beta: float | Mapping[str, float] = 0.1,
eta: float | Mapping[str, float] = 0.1,
pi: Mapping[frozenset[str], float] | np.ndarray | None = None,
seed: int | np.random.Generator | None = None,
return_states: bool = False,
) -> ResponseMatrix | tuple[ResponseMatrix, list[frozenset[str]]]:
"""Simulate response patterns from a BLIM.
Each simulated respondent is assigned a knowledge state drawn from
``pi``, then answers every item independently: an item in the state is
answered correctly with probability ``1 - beta[q]`` (careless error
``beta[q]``), an item outside the state is answered correctly with
probability ``eta[q]`` (lucky guess). This is the data-generating
process of the basic local independence model (Falmagne & Doignon
2011, ch. 11) and mirrors ``pks::simulate.blim()``.
Parameters
----------
structure : KnowledgeStructure
The knowledge structure whose states generate the data.
n_respondents : int
Number of respondents to simulate.
beta : float or Mapping[str, float]
Careless-error probability, scalar (homogeneous) or per-item.
Values must be in [0, 1). Default 0.1.
eta : float or Mapping[str, float]
Lucky-guess probability, scalar or per-item, in [0, 1).
Default 0.1. The informative-item condition ``beta + eta < 1`` is
*not* enforced: simulating degenerate items is a legitimate use
(e.g., to exercise degenerate-item diagnostics).
pi : Mapping[frozenset[str], float] or np.ndarray or None
State prior. Either a mapping from state to probability, or an
array in the canonical state order (states sorted by size, then
lexicographically — the same order used by
:func:`~knowledgespaces.estimation.estimate_blim`), or None for
the uniform prior. Must be non-negative and sum to 1 (tolerance
1e-8).
seed : int, np.random.Generator, or None
Seed or generator for reproducibility.
return_states : bool
If False (default), return an aggregated :class:`ResponseMatrix`
(unique patterns with counts). If True, return the tuple
``(data, states)`` where ``data`` has one *non-aggregated* row per
respondent and ``states[r]`` is respondent ``r``'s true knowledge
state — the form needed for recovery and assessment simulations.
Returns
-------
ResponseMatrix or (ResponseMatrix, list[frozenset[str]])
Simulated response data; with ``return_states=True``, also the
true state of each respondent (row-aligned).
Raises
------
ValueError
If ``n_respondents < 1``, parameters are out of range, or ``pi``
does not match the structure's states.
Examples
--------
>>> from knowledgespaces import space_from_prerequisites
>>> s = space_from_prerequisites(["a", "b"], [("a", "b")])
>>> data = simulate_blim(s, 500, beta=0.1, eta=0.05, seed=42)
>>> int(data.effective_counts.sum())
500
"""
if n_respondents < 1:
raise ValueError(f"n_respondents must be >= 1, got {n_respondents}.")
items = sorted(structure.domain)
states = sorted(structure.states, key=lambda s: (len(s), sorted(s)))
n_items = len(items)
n_states = len(states)
beta_vec = _param_vector(beta, items, "beta")
eta_vec = _param_vector(eta, items, "eta")
if pi is None:
pi_vec: NDArray[np.float64] = np.full(n_states, 1.0 / n_states)
elif isinstance(pi, Mapping):
missing = [s for s in states if s not in pi]
if missing:
raise ValueError(f"pi is missing states {sorted(missing, key=sorted)}.")
extra = [s for s in pi if frozenset(s) not in set(states)]
if extra:
raise ValueError(f"pi has unknown states {extra}.")
pi_vec = np.array([float(pi[s]) for s in states], dtype=np.float64)
else:
pi_vec = np.asarray(pi, dtype=np.float64)
if pi_vec.shape != (n_states,):
raise ValueError(
f"pi array has shape {pi_vec.shape}, expected ({n_states},) "
f"in canonical state order (by size, then lexicographically)."
)
if not np.all(np.isfinite(pi_vec)):
raise ValueError("pi must be finite.")
if np.any(pi_vec < 0):
raise ValueError("pi must be non-negative.")
if abs(pi_vec.sum() - 1.0) > 1e-8:
raise ValueError(f"pi must sum to 1, got {pi_vec.sum():.10f}.")
rng = seed if isinstance(seed, np.random.Generator) else np.random.default_rng(seed)
# State membership matrix S[k, q]
item_idx = {q: i for i, q in enumerate(items)}
S = np.zeros((n_states, n_items), dtype=np.float64)
for k, state in enumerate(states):
for q in state:
S[k, item_idx[q]] = 1.0
state_idx = rng.choice(n_states, size=n_respondents, p=pi_vec)
mastered = S[state_idx] # (n_respondents, n_items)
u = rng.random((n_respondents, n_items))
# mastered: correct unless careless error; not mastered: lucky guess
responses = np.where(mastered == 1.0, u >= beta_vec, u < eta_vec).astype(int)
if return_states:
data = ResponseMatrix(items=items, patterns=responses)
true_states = [states[k] for k in state_idx]
return data, true_states
unique, counts = np.unique(responses, axis=0, return_counts=True)
return ResponseMatrix(items=items, patterns=unique, counts=counts.astype(np.float64))