SLM, ITA and IITA inference¶

Simple learning model¶

The simple learning model (SLM) parametrizes latent state probabilities through item solvability parameters \(g_q\). For a state \(K\) and its outer fringe \(K^\mathcal O\),

\[ \pi_g(K)=\prod_{q\in K}g_q\prod_{q\in K^\mathcal O}(1-g_q). \]

These products sum to one on a finite learning space, including boundary values \(g_q=0,1\): Falmagne and Doignon (2011), Eq. 11.18 and Theorem 11.5.4. The API checks this sufficient structural hypothesis. It rejects arbitrary knowledge structures instead of renormalizing their products. Well-gradedness alone or union closure alone is insufficient (Examples 11.5.1–11.5.2). Some other families may normalize for particular parameters; they are outside this API’s guaranteed domain. Solvability is generally not marginal mastery.

from knowledgespaces.datasets import load_doignon_falmagne_7
from knowledgespaces.estimation import (
    BLIMConstraints, estimate_slm, slm_state_probabilities,
)

example = load_doignon_falmagne_7()
fit = estimate_slm(example.structure, example.data, max_iter=10000)
print(fit.converged, fit.g_dict(), fit.gof.G2)
assert abs(slm_state_probabilities(example.structure, fit.g_dict()).sum() - 1) < 1e-12
prediction = fit.predict(example.data.patterns, items=example.data.items)
residuals = fit.residuals(example.data)

The response model uses the same slip beta and guess eta convention as BLIM. Given expected state counts \(M_K\), the M-step updates

\[ g_q=\frac{\sum_{K:q\in K}M_K} {\sum_{K:q\in K}M_K+\sum_{K:q\in K^\mathcal O}M_K}. \]

ML uses all states; MD makes one discrepancy assignment; MDML restricts its E-step to the selected states. BLIMConstraints supports fixed/equal error parameters. pi_fixed is rejected because the prior is determined by \(g\); fixed or shared \(g\) is not currently exposed. Free errors and \(g\) are clipped to \([10^{-6},1-10^{-6}]\) during estimation. Fixed errors may be exactly zero. ML/MDML stop on maximum absolute change in errors and \(g\), as in pks::slm. objective_history retains the initial and updated EM objectives (empty for MD); MDML’s objective is the restricted joint sum. Reported likelihood and G2 always refer to the full response model. Convergence is not a guarantee of a global maximum.

gof.npar counts free error groups plus \(|Q|\), not free state masses. fit.jacobian() or slm_jacobian differentiates this SLM parameterization; columns are free beta groups, free eta groups, then \(g\) in sorted item order. Its rank is a pointwise numerical diagnostic. A free-prior BLIM rank report cannot establish SLM identifiability.

bootstrap_slm fits, simulates and refits the SLM with identical method, constraints and discrepancy settings. It retains g_replicates as well as beta/eta/pi, G2, convergence and iterations. State probabilities remain SLM-derived in every replicate. bootstrap_gof rejects an SLM estimate to prevent accidentally refitting a different model. Inspect convergence and identifiability before interpreting parameter dispersions; these samples are not automatically valid confidence intervals. The SLM estimator guards working-array memory; retained bootstrap arrays consume additional storage.

SLM here specifies a cross-sectional distribution. It does not estimate individual learning transitions or a general longitudinal Markov process. The probability dataset’s K1/K2 families are not learning spaces: use BLIM for these candidates, not this SLM.

Random SLM starts and search diagnostics¶

estimate_slm_restarts supports ML and MDML with reproducible random starts for beta, eta and g. n_restarts=1, init_strategy="pks" supplies the random initialization law of pks::slm(randinit=TRUE): errors from U(0,1), sequential reflection of beta then eta, followed by independent U(0,1) solvabilities. R and NumPy seeds produce different numbers. Equality groups are projected before Python’s first E-step; pks initially leaves them unequal. Fixed values override the draws, and free starts use the estimator’s numerical box.

The default uniform strategy uses errors from U(.01,.4); with a wider init_range, each pair is halved until its sum is below .95. Solvabilities still range over U(0,1). Neither strategy constrains final beta+eta, proves identifiability, or guarantees finding the global optimum.

from knowledgespaces.estimation import estimate_slm_restarts

searched = estimate_slm_restarts(
    example.structure, example.data, n_restarts=5, seed=42,
    method="ML", max_iter=10000,
)
print(searched.selected_restart, searched.converged)
for run in searched.restarts:
    print(run.log_likelihood, run.objective, run.converged)

The result is an SLMEstimate with every run’s effective starts, criterion, full likelihood, iterations and convergence, plus the zero-based selected index. Replaying a recorded start uses estimate_slm with numpy arrays for beta_init, eta_init and g_init. Selection defaults to the optimized objective: full likelihood for ML, restricted joint sum for MDML. With selection="likelihood", MDML runs are instead ranked by full likelihood, matching the selection convention of estimate_blim_restarts. These are different criteria. First run wins exact ties. Unconverged runs remain eligible; a selected unconverged fit warns and must be inspected.

bootstrap_slm(..., n_restarts=5, seed=42) applies the same search policy to observed and simulated data, using fresh draws from a shared local generator. Replicate arrays retain each search’s selected fit. Default n_restarts=None preserves the deterministic 0.1 starts. MD is independent of starts and rejects random-restart requests.

Discrepancy inclusion options¶

from knowledgespaces.estimation import MDOptions, estimate_blim

radius_fit = estimate_blim(
    example.structure, example.data, method="MDML",
    discrepancy=MDOptions(radius=1),
)
weighted_fit = estimate_blim(
    example.structure, example.data, method="MD",
    discrepancy=MDOptions(rule="hypblc1", exponent=2),
)

Let \(d(R,K)\) be Hamming distance and \(d_{\min}(R)\) the minimum over feasible states for response \(R\). These are the pks::blim/blimMD conventions:

Rule	Inclusion weight before row normalization	Methods
`minimum`	\(1\{d(R,K)\le d_{\min}(R)+r\}\)	MD, MDML
`hypblc1`	\((1+d(R,K)-d_{\min}(R))^{-m}\)	MD
`hypblc2`	\((1+d(R,K))^{-m}\)	MD

Radius \(r\) is a nonnegative integer excess distance from the response. It is not the distance between a candidate state and a nearest state. Exponent \(m\) must be positive and finite. Hyperbolic rules include all feasible states and cannot be combined with a radius. MDML uses minimum-rule inclusion as a restriction on its posterior, not a uniform assignment after each iteration. Explicit options are rejected for unrestricted ML.

Fixed-zero errors and fixed state priors can make states infeasible. Their weights remain zero, even under a large radius. An observed response with no feasible state raises an error. Zero-frequency impossible rows do not affect the fit. Hyperbolic weights use log scaling to avoid total underflow for large exponents. Free BLIM errors retain the library’s numerical box; unlike pks::blimMD, this can replace a free zero estimate by \(10^{-6}\).

The options are stored on fitted objects and forwarded by restarts and bootstrap. A supplied bootstrap estimate passes its options to refits; explicitly conflicting options are rejected. Default None retains the ordinary MD/MDML behavior.

Threshold item tree analysis¶

derivation.ita is the threshold method, distinct from inductive derivation.iita. Pair (a,b) enters at threshold \(L\) if its weighted number of observations \(a=0,b=1\) is at most \(L\). Thus mastery of b implies mastery of a. A threshold relation must already be transitive; taking its closure would change the method.

from knowledgespaces.derivation import ita

result = ita(example.data, search="global")
print(result.threshold, result.discrepancy.fit, result.discrepancy.complexity)

The criterion is the sum of mean response-to-structure Hamming distance (weighted by respondent counts) and mean state-to-observed-response distance (uniform over states). Zero-frequency rows are excluded from both terms. This descriptive tradeoff is not a hypothesis test. Fractional analysis weights are supported, in which case thresholds are weighted counts.

Search evaluates the distinct observed counterexample levels whose relations are transitive, largest first. local stops at the first increase in the criterion; global evaluates all levels. Ties select the largest threshold. The result exposes all transitive thresholds and the evaluated fit/complexity components. Supplying threshold=L bypasses search and checks transitivity. With an explicit threshold, make_structure=False avoids state enumeration; automatic search still needs structures to calculate its criterion. max_states bounds enumeration and max_memory_bytes preflights arrays.

Quasi-ordinal simulation and supplied IITA candidates¶

from knowledgespaces.estimation import random_surmise_relation, simulate_quasiordinal
from knowledgespaces.derivation import counterexamples, inductive_generation, iita

simulation = simulate_quasiordinal(
    ["a", "b", "c", "d"], 1000, delta=.15,
    beta=.08, eta=.12, seed=42,
)
candidates = inductive_generation(
    counterexamples(simulation.data), items=simulation.data.items,
)
comparison = iita(simulation.data, selection_set=candidates, version="minimized")
print(comparison.diff, comparison.error_rates)

random_surmise_relation(items, delta) independently includes each ordered nonreflexive pair with probability delta, then takes transitive closure. Reflexivity is implicit; (a,b) means a is a prerequisite for b. Cycles produce equivalent items. This is the sampling scheme documented by DAKS 2.1-3, not uniform sampling among quasi orders. Delta controls inclusion before closure, not final edge density. Pair draws follow sorted labels for reproducibility.

simulate_quasiordinal accepts exactly one of delta or a supplied, already transitive relation. It returns the relation, its full downset structure, one response row per respondent and row-aligned true_states. Equivalent items are supported: a quasi-ordinal space need not be a learning space. viz.plot_hasse(simulation.structure) displays inclusion covers, including covers spanning multiple items when items are equivalent. The default state prior is uniform; explicit pi and heterogeneous beta/eta use simulate_blim conventions. Enumeration is bounded by max_states, relation pairs by max_pairs, and numerical arrays by max_memory_bytes. Python set/object overhead is additional. Large sparse relations can still have exponentially many states.

Responses follow the canonical BLIM: P(correct | mastered)=1-beta and P(correct | unmastered)=eta. DAKS 2.1-3’s implementation differs from its stated model: its two sequential mutations allow a careless error to be reversed by a lucky guess, giving effective slip ce*(1-lg). That discrepancy is not reproduced here. With uniform latent 00/11, ce=.4 and lg=.3, DAKS produces marginal success .51, while the stated BLIM gives .45. The native R fixture reproduces this example; conditional-rate tests verify Python’s model. This distinction is about response generation, not IITA’s gamma, which is a separate criterion parameter and is not a BLIM slip estimate.

inductive_generation exposes candidate generation without fitting or state enumeration. Its counterexample matrix follows the supplied item order. iita(selection_set=...) evaluates expert-specified or other prespecified relations using any of the three criteria, corresponding to DAKS’s orig_iita, corr_iita and mini_iita A argument. Candidates must already be transitive, have matching labels and contain a nonreflexive implication; otherwise gamma is undefined. Their order and duplicates are retained, with the first minimum selected. error_rates reports all candidate gammas; error_rate is the selected one. These descriptive outputs carry no post-selection inferential guarantee.

Population IITA and fixed-relation inference¶

population_iita computes response probabilities under a BLIM and evaluates all three IITA criteria without generating a sample. It accepts a relation (expanded to compatible states) or a general structure, scalar/item-specific errors and uniform/nonuniform state masses. Candidate relations come from population counterexamples, an explicit selection_set, or a separate data argument. Computed counterexample near-ties may change inductive candidate generation under rounding; explicit candidates make comparisons reproducible. All response patterns are enumerated, with size/memory guards.

For inference, iita_variance uses integer respondent frequencies and population_iita_variance uses a specified probability distribution. Both treat the supplied, transitively closed relation as fixed and recompute the IITA error rate from the distribution. Every item must have positive success probability and at least one nonreflexive implication is required.

Write \(D=\mathrm{diff}/N^2\) for the normalized sample criterion, or its population analogue. For multinomial response probabilities \(\rho\), the first-order delta coefficient is

\[ V=\nabla D^\top[\operatorname{diag}(\rho)-\rho\rho^\top]\nabla D. \]

asymptotic_variance is \(V\), the coefficient for \(\sqrt N(\widehat D-D)\); standard_error is \(\sqrt{\widehat V/N}\) in a sample and None for a population. Analytic derivatives include the estimated error rate. The read-only influence array is the centered derivative in input-row order. Zero-probability cells contribute zero, so sample inference does not need powerset enumeration. Original, corrected and minimized variants are supported.

iita_z_test compares one fixed relation’s \(D\) against a prespecified null or compares \(D_1-D_2\). With other_relation only, the criteria use the same respondents and their covariance is included. Supplying other_data assumes independent samples. greater uses the upper normal tail; less uses the lower tail and an interval \((-\infty,U)\). Intervals are untruncated normal intervals and may extend beyond the parameter range.

A single-relation null \(D=0\) is nonregular. The criterion is a sum of squares, so its gradient vanishes at that null. A positive sample plug-in variance does not repair the null calibration. The function warns and returns NaN p-value/interval; any finite Z is descriptive. Zero plug-in variance also returns undefined Z. regular_null=False flags these detected failures; True is not a proof of remaining regularity assumptions. A prespecified positive single-relation null or a comparison of fixed relations may use the ordinary approximation under regularity. No second-order test for exact fit or post-selection adjustment is currently implemented.

Relations discovered and selected on these same respondents cannot be interpreted as prespecified hypotheses. These functions do not account for that selection uncertainty. Population evaluation and descriptive criteria remain available independently of normal calibration.

Deliberate differences from DAKS¶

The criteria and population probabilities are compared with DAKS 2.1-3. The delta derivatives are instead checked against central finite differences of the criterion and the full multinomial covariance. On the documented three-item fork fixture, minimized sample \(V\) is \(8.68195355867706\times10^{-5}\) in Python versus \(7.18199712441497\times10^{-5}\) in DAKS; the minimized chain case agrees. DAKS derivative expressions omit/change contributions on some configurations. The Python implementation does not reproduce those discrepancies. Its same-respondent covariance, one-sided tails and zero-null handling also differ from DAKS z_test. Numerical agreement is compatibility evidence on specific cases, not a proof of universal correctness.

See tests/reference_data/generate_classical.R, classical.json and tests/test_iita_inference.py for the reproducible reference and derivative checks. The executable cookbook/06_classical_methods.py illustrates SLM and fixed-relation comparisons on complete binary data.

References¶

Falmagne, J.-C., & Doignon, J.-P. (2011). Learning Spaces, Chapter 11, Eq. 11.18, Examples 11.5.1–11.5.2, Theorem 11.5.4.
pks SLM documentation and ITA documentation, with the archived pks 0.7-0 source and manuals for inclusion/EM conventions.
Heller, J., & Wickelmaier, F. (2013). Minimum discrepancy estimation in probabilistic knowledge structures. Electronic Notes in Discrete Mathematics, 42, 49–56.
Schrepp, M. (1999). On the empirical construction of implications between bi-valued test items. Mathematical Social Sciences, 38, 361–375.
Ünlü, A., & Sargin, A. (2010). DAKS: An R Package for Data Analysis Methods in Knowledge Space Theory. Journal of Statistical Software, 37(2), Section 3.2, footnote 4 for the delta coefficient. Normal inference additionally requires nondegeneracy.

Run python cookbook/12_randomized_methods.py for an offline example of SLM restart diagnostics, quasi-ordinal simulation and supplied IITA candidates.