Parameter Estimation¶

The knowledgespaces.estimation module estimates BLIM parameters from observed response data — by maximum likelihood via Expectation-Maximization (default), or by minimum discrepancy (MD) and minimum-discrepancy maximum likelihood (MDML; Heller & Wickelmaier, 2013). It also provides BLIM data simulation and identifiability diagnostics.

Predicting states after fitting¶

For a complete-data BLIM or SLM estimate fit, provide a binary matrix responses with columns in training item order, or give items= explicitly:

bayes = fit.predict(responses)              # default, regardless of fit.method
md = fit.predict(responses, method="MD")
mdml = fit.predict(responses, method="MDML")
state_weights = mdml.state_probabilities    # response rows × canonical states
states = mdml.classify(ties="min")          # tuple of frozensets (or None)
random_states = md.classify(ties="random", seed=42)

The assignment methods follow the formulas in the pks prediction manual:

ML uses Bayes’ formula with fitted parameters: $w_{RK}=P(R\mid K)\pi_K / P(R)$. This is a posterior over states; the name does not mean maximizing the conditional likelihood over states.
MD gives equal weight to states at minimum Hamming distance from the response, or within the requested excess radius.
MDML restricts $P(R\mid K)\pi_K$ to the included states and normalizes over that selection.

For every method, .probabilities, .log_probabilities, and .posteriors retain the full BLIM $P(R)$, log $P(R)$, and Bayesian $P(K\mid R)$. Use .state_probabilities for the chosen assignment method. An MD assignment is a distance-based rule, not the posterior of the original BLIM. SLM prediction uses its derived state prior and the same response model.

classify selects a state with the largest computed assignment probability. Exact ties use smallest cardinality (min), largest cardinality (max), or a uniform draw over all tied states (random). Remaining cardinality ties use canonical order: size, then lexicographic item labels. Random ties use a local NumPy generator; seed can also be a numpy.random.Generator. R and NumPy seeds do not imply identical draws. Each input row receives one classification; compressed frequency rows do not simulate individual draws.

from knowledgespaces.estimation import MDOptions

wider = fit.predict(responses, method="MDML", discrepancy=MDOptions(radius=1))

For MD/MDML, omitted options inherit the fit’s discrepancy rule and radius; otherwise MDOptions() supplies minimum distance with radius zero. The radius is an excess distance from the response, not a distance between states. Hyperbolic MDOptions(rule="hypblc1"/"hypblc2", exponent=...) are also available for MD prediction using the estimation weights; they are an extension beyond pks’ default minimum-rule prediction. A hyperbolic fit requires an explicit minimum-rule override to predict by MDML.

Alternatively, inclusion= supplies a binary row-by-canonical-state mask (pks i.RK). It overrides the distance rule and fixed-zero exclusions for MD; it is mutually exclusive with explicitly supplied discrepancy and is invalid for ML. Unlike R’s permissive matrix arithmetic, negative, fractional, missing, or misdimensioned masks are rejected. Numeric hyperbolic weights are specified through MDOptions, not disguised as an indicator matrix.

Computed discrepancy assignments exclude errors forbidden by declared fixed-zero beta/eta, including equality groups. In low-level predict_blim, pass constraints= explicitly; restrictions must agree with supplied parameters. Merely supplying a zero parameter without declaring it fixed does not modify the MD distance rule. A zero state prior does not exclude that state from MD, but does give it zero joint weight in MDML.

An empty inclusion set or zero total selected joint mass gives an all-NaN assignment row and None from classification. The empty knowledge state remains a legitimate frozenset(), distinct from None. Bayesian posterior rows are undefined only when the model assigns zero probability to the response; calculations in log space preserve well-defined conditionals even when the floating-point value of $P(R)$ underflows to zero. Missing responses remain handled by the separate incomplete-data Bayesian interface.

See executable cookbook/16_state_prediction.py for a complete example.

What it estimates¶

Given a knowledge structure and a matrix of student responses, the EM algorithm estimates:

$\beta_q$ (slip per item): $P(\text{incorrect} \mid q \text{ mastered})$
$\eta_q$ (guess per item): $P(\text{correct} \mid q \text{ not mastered})$
$\pi_K$ (state prior): $P(\text{student is in state } K)$

High-level API¶

import knowledgespaces as ks

structure = ks.space_from_prerequisites(
    ["add", "sub", "mul"],
    [("add", "sub"), ("sub", "mul")],
)

result = ks.fit_blim(
    structure,
    items=["add", "sub", "mul"],
    responses=[[1,1,1], [1,1,0], [1,0,0], [0,0,0]],
    counts=[45, 30, 20, 5],  # optional: pattern frequencies
)

print(result["converged"])           # True
print(result["n_iterations"])        # number of EM iterations
print(result["beta"])                # slip per item (dict)
print(result["eta"])                 # guess per item (dict)
print(result["log_likelihood"])      # final log-likelihood
print(result["degenerate_items"])    # items with beta + eta >= 1 (see below)
print(result["gof"]["BIC"])          # conventional Schwarz BIC; see limitations below
print(result["gof"]["BIC_npatterns"])# pks-compatible variant (see below)
print(result["gof"]["AIC"])          # Akaike information criterion
print(result["gof"]["G2"])           # likelihood ratio statistic

Low-level API¶

For full control:

from knowledgespaces.estimation import estimate_blim, ResponseMatrix
import numpy as np

data = ResponseMatrix(
    items=["add", "sub", "mul"],
    patterns=np.array([[1,1,1], [1,1,0], [1,0,0], [0,0,0]]),
    counts=np.array([45, 30, 20, 5]),
)

result = estimate_blim(
    structure, data,
    max_iter=500,
    tol=1e-6,
    beta_init=np.array([0.05, 0.1, 0.15]),  # per-item initialization
    eta_init=0.1,                          # or global scalar
)

print(result.beta_for("add"))
print(result.eta_for("mul"))
print(result.pi)  # state prior distribution

How it works¶

The EM algorithm alternates:

E-step: For each response pattern, compute the posterior probability of each knowledge state.
M-step: Re-estimate $\beta_q$, $\eta_q$, and $\pi_K$ from the weighted sufficient statistics.

An exact EM step cannot decrease the log-likelihood. Numerical clipping and finite precision can cause small departures from this property. Convergence is declared when the change in log-likelihood falls below tol.

With no fixed/shared constraints, the M-step clips each $\beta_q$ and $\eta_q$ independently into $(0, 1)$, as in pks::blim(). The joint condition $\beta_q + \eta_q < 1$ is the informative-item condition — an informativeness (positive-discrimination) requirement, not part of the parameter space (Falmagne & Doignon, 2011, §11) — and is deliberately left unenforced, matching pks::blim(). Enforcing it in the unrestricted item-error model would keep the closed form and EM monotonicity: the order-restricted M-step pools the two success probabilities of an offending item. But the pooled update parks the item exactly on the uninformative boundary $\beta_q + \eta_q = 1$, so the estimator instead leaves the violation visible and reports such items as degenerate (see below).

The same parameter domain applies to initial values: beta_init and eta_init must each be finite and in [0, 1), but their sum may equal or exceed one. Free starting errors are clipped individually to [1e-6, 1 - 1e-6]; fixed constraints keep their specified values. This allows a valid fitted model to initialize a new fit, including the embedded null in a bootstrap comparison, without rejecting or reflecting reversed discrimination. The final diagnostic still reports uninformative or reversed items; accepting a start does not guarantee identifiability or a global optimum.

Estimation methods: ML, MD, MDML¶

Three estimation methods are available, matching pks::blim() in R (Heller & Wickelmaier, 2013):

ml   = estimate_blim(structure, data)                 # default: EM / ML
md   = estimate_blim(structure, data, method="MD")    # minimum discrepancy
mdml = estimate_blim(structure, data, method="MDML")  # ML among MD solutions
print(md.method, md.n_iterations)   # "MD", 1 — non-iterative, deterministic

ML maximizes the multinomial likelihood via EM. Convergence is declared on the absolute log-likelihood change (tol).
MD assigns each observed pattern uniformly to the knowledge states requiring the fewest response errors to explain it, then reads the parameters off a single M-step. It is deterministic given the data (initialization plays no role), which makes it ideal for exact cross-package comparisons.
MDML runs EM with the posterior restricted to each pattern’s minimum-discrepancy states — maximum likelihood within the minimum-discrepancy solutions. Convergence is declared on the maximum absolute parameter change, mirroring pks.

The returned log-likelihood always uses the unrestricted BLIM marginal, including for MDML. This is not the restricted objective used by MDML, so do not use its monotonicity as a convergence test or assume a universal MD-to-MDML ordering for the reported likelihood. A global ML optimum dominates feasible parameter points under the same constraints, but single-start EM need not find it. Restarts mitigate local optima without certifying global optimality. All three methods are numerically compared with pks::blim() in tests/test_cross_validation.py.

Fixed and shared parameters¶

Use BLIMConstraints to express model hypotheses independently of starting values. A scalar beta_init=0.1 initializes every item at the same value; it does not make their fitted errors equal.

from knowledgespaces.datasets import load_doignon_falmagne_7
from knowledgespaces.estimation import BLIMConstraints, estimate_blim, check_identifiability

ds = load_doignon_falmagne_7()
constraints = BLIMConstraints(
    beta_equal=[ds.data.items],  # one shared slip parameter
    eta_equal=[ds.data.items],   # one shared guess parameter
)
fit = estimate_blim(ds.structure, ds.data, constraints=constraints, max_iter=2000)
print(fit.gof.npar)  # |K| - 1 + 2
diagnostic = check_identifiability(ds.structure, constraints=constraints)
print(diagnostic.rank, diagnostic.npar)

beta_fixed={"a": 0.05} and eta_fixed={"e": 0} fix labelled parameters. Fixed values are in [0, 1); zero is exact. Overlapping equality groups merge transitively. A fixed member fixes its entire equality group; conflicting declarations raise an error. The M-step pools expected error counts and denominators within each free group: the constrained Bernoulli maximizer. Free values retain the estimator’s interior clipping.

pi_fixed supplies a complete labelled state distribution, including explicit zero masses. Partial fixed priors are not supported. ML and MDML support fixed priors; assignment-based MD rejects them. MD and MDML exclude states made impossible by fixed zeros before computing minimum discrepancy. Positive-frequency observations impossible under the declared model raise ValueError; zero-frequency rows do not alter the fit.

estimate_blim_restarts uses the constraints at every start. bootstrap_gof(..., estimate=fit) inherits them for every refit, and checks that the supplied parameters satisfy the declared constraints. If passing constraints explicitly too, supply the same specification. blim_jacobian and check_identifiability accept the same object, sum derivative columns for shared parameters and omit fixed coordinates. Pass constraints=fit.constraints when diagnosing a fitted model. Fixed coordinates are imposed on the diagnostic point; starting values within a free equality group are replaced by their mean.

GOF parameter counts reflect these constraints; they remain conventional counts, not a proof of observable model dimension. Shared/fixed errors need substantive justification and are distinct from imposing beta + eta < 1.

Simulating data from a BLIM¶

simulate_blim generates response patterns from a known BLIM — the counterpart of pks::simulate.blim() and the building block for parameter-recovery studies:

from knowledgespaces.estimation import simulate_blim

data = simulate_blim(structure, 1000, beta=0.1, eta=0.15, seed=42)
est = estimate_blim(structure, data)  # recover the generating parameters

# For assessment simulations, keep each respondent's true state:
data, true_states = simulate_blim(structure, 500, seed=1, return_states=True)

beta/eta accept a scalar or a per-item dict; pi an optional state prior (uniform by default).

Identifiability diagnostics¶

The BLIM is not identifiable in general: different parameter vectors can produce the same distribution over response patterns, so two correct programs may return different $(\beta, \eta, \pi)$ with equal likelihood. Check before interpreting (or comparing) parameter estimates:

from knowledgespaces.estimation import check_identifiability

report = check_identifiability(structure)   # beta = eta = 0.1, uniform pi, + random points
print(report.rank, report.npar)             # Jacobian rank vs. free parameters
print(report.locally_identifiable)          # rank == npar?
print(report.forward_graded)                # items with an eta ↔ pi trade-off
print(report.backward_graded)               # items with a beta ↔ pi trade-off

The rank criterion follows Stefanutti, Heller, Anselmi & Robusto (2012); the gradedness conditions follow Spoto, Stefanutti & Vidotto (2013). When the structure is not identifiable, compare fits on likelihood-based quantities only (log-likelihood, $G^2$, AIC, BIC, predicted pattern probabilities). The classic DoignonFalmagne7 structure, for instance, has rank 14 < npar 18.

The rank at a single parameter point only lower-bounds the generic rank of the model: on some symmetric structures the default point $\beta = \eta = 0.1$ with uniform $\pi$ falls in the exceptional null set where the rank drops. check_identifiability therefore also evaluates the Jacobian at up to n_points (default 5) random interior points drawn from a fixed-seed generator (seed=0), reports the maximum rank, and stops early once full rank is found — so repeated calls are deterministic. The verdict is one-sided: exact full rank at an interior point establishes generic local identifiability. The implementation estimates rank with a floating-point singular-value threshold; inspect the reported tolerance and singular-value margin. Deficient numerical ranks at all sampled points provide evidence, not a proof, of generic rank deficiency. Neither result establishes global identifiability. Pass n_points=0 to check exactly one point, e.g. fitted estimates: check_identifiability(s, beta=est.beta_dict(), eta=est.eta_dict(), pi=est.pi_dict(), n_points=0). The dictionaries align parameters by item and state labels, even when the response columns were not supplied in sorted order. A deficient rank at a singular point alone does not prove local nonidentifiability.

Classic datasets¶

from knowledgespaces.datasets import load_doignon_falmagne_7

ds = load_doignon_falmagne_7()   # Doignon & Falmagne (1999, ch. 7)
est = estimate_blim(ds.structure, ds.data)

Random restarts¶

EM can converge to local optima. To mitigate this, use estimate_blim_restarts, which runs EM several times from random initial values and returns the best fit:

from knowledgespaces.estimation import estimate_blim_restarts

result = estimate_blim_restarts(
    structure, data,
    n_restarts=10,
    seed=42,              # reproducibility
    init_strategy="uniform",  # or "pks" to mirror pks::blim(..., randinit=TRUE)
)

Degenerate items¶

Items whose final $\beta_q + \eta_q \geq 1 - 10^{-3}$ are flagged by a numerical tolerance around the informative-item boundary. At exactly $\beta_q + \eta_q = 1$, success does not depend on mastery; above one, discrimination is reversed. Values just below one still discriminate positively, but weakly. The estimator surfaces these items via result["degenerate_items"] and emits a ConvergenceWarning. Inspect coding, fit stability and the assumed structure before revising or removing an item; this flag alone does not identify the cause.

Goodness of fit and model selection¶

result["gof"] collects the conventional summaries:

G2, df, p_value — likelihood-ratio $G^2$ against the saturated multinomial, with reported degrees of freedom $\max(\min(2^Q - 1, N) - \mathrm{npar}, 0)$ (pks convention; caps at $N$ when the sample is smaller than the full pattern space). The p_value is asymptotic and descriptive: whenever the cap binds ($N < 2^Q - 1$) the estimator emits a SparseGOFWarning. This screens for sparsity; it does not prove failure of every chi-squared approximation. Expected counts and regularity matter (Koehler & Larntz, 1980). For a calibrated test use the parametric bootstrap below. The cap is a software convention, not a change in the population model’s dimension; p_value is NaN at zero reported degrees of freedom.
`AIC = -2,\mathrm{LL} + 2,\mathrm{npar}$.
BIC = -2\,\mathrm{LL} + \log(N)\,\mathrm{npar} — the Schwarz (1978) form based on the number of respondents $N$. Its consistency argument holds in regular models; BLIM selection problems are often nonregular (estimates on the boundary, rank-deficient structures), so read AIC and BIC as conventional summaries, and prefer predictive comparisons (e.g., cross-validated log-scores) when a selection matters.
BIC_npatterns = -2\,\mathrm{LL} + \log(n_\text{patterns})\, \mathrm{npar} — the variant returned by R pks::blim() via its nobs.blim override. Exposed for cross-package replication only; $n_\text{patterns}$ is bounded above by $2^Q$, so the Schwarz argument does not apply to it at all.

Calibration by parametric bootstrap¶

bootstrap_gof simulates n_replicates data sets of the observed size from the fitted model, refits each with the same estimator settings, and locates the observed $G^2$ in the bootstrap distribution (Monte Carlo p-value with the add-one convention):

from knowledgespaces.estimation import bootstrap_gof

boot = bootstrap_gof(structure, data, n_replicates=1000, seed=42)
boot.p_value       # calibrated p-value
boot.g2_observed   # G2 of the fit to the observed data
boot.n_capped      # refits stopped by the iteration cap (should be 0)

The run is sequential by design (one process, one core). Refits that hit the iteration cap are counted in n_capped and reported with a ConvergenceWarning: unfinished optimization can inflate replicated $G^2$ values. Bootstrap validity still depends on the fitted model, estimator and Monte Carlo precision; it is not an exact finite-sample guarantee. Respondent frequencies must be integers. A supplied estimate is checked against the structure, method and observed-data likelihood.

The result also retains parameter replicates and refit convergence status. See diagnostics for their dispersion summaries, full-table residuals and Jacobian trade-offs by parameter family.

Additional classical estimators¶

See the classical methods guide for SLM estimation and its Jacobian/bootstrap, plus MDOptions radius and hyperbolic inclusion rules.

For incomplete 0/1/NaN responses, use the dedicated observed-data ML API described in the incomplete-data guide. Complete-data estimators keep their original binary-only contract.