"""Population IITA and multinomial delta-method inference.
Independent moment/gradient implementation of Ünlü & Sargin (2010), JSS
37(2), Section 3.2, footnote 4. The normalized criterion is diff/N**2.
No GPL derivative code is ported. Same-sample comparisons use covariance;
normal-tail conventions follow the stated alternative hypothesis.
"""
from __future__ import annotations
import warnings
from collections.abc import Sequence
from dataclasses import dataclass
from typing import Literal
import numpy as np
from scipy.stats import norm
from knowledgespaces._patterns import _binary_patterns, _pattern_count, _positive_integer
from knowledgespaces.derivation.iita import (
_diff,
_error_rate,
_error_rate_minimized,
_inductive_generation,
_pairs,
counterexamples,
)
from knowledgespaces.estimation.blim_em import ResponseMatrix
from knowledgespaces.estimation.prediction import ItemParameter, StateParameter, predict_blim
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
from knowledgespaces.structures.relations import SurmiseRelation
Version = Literal["original", "corrected", "minimized"]
[docs]
class IITAInferenceWarning(UserWarning):
"""First-order normal calibration is undefined or the null is nonregular."""
def _moments(responses: np.ndarray, probabilities: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
p = probabilities @ responses
if np.any(p <= 0):
raise ValueError("Every item must have positive response success probability.")
return p, ((1 - responses) * probabilities[:, None]).T @ responses
def _criterion_gradient(
responses: np.ndarray,
probabilities: np.ndarray,
pairs: frozenset[tuple[int, int]],
version: Version,
) -> tuple[float, float, np.ndarray]:
if version not in ("original", "corrected", "minimized"):
raise ValueError("version must be 'original', 'corrected', or 'minimized'.")
p, b = _moments(responses, probabilities)
n, q = responses.shape
if version == "minimized":
numerator = denominator = 0.0
dnum, dden = np.zeros(n), np.zeros(n)
for i in range(q):
for j in range(q):
db = (1 - responses[:, i]) * responses[:, j]
if (i, j) in pairs:
value, slope, delta, ddelta = p[j], responses[:, j], b[i, j], db
elif (j, i) in pairs:
value, slope = p[i], responses[:, i]
delta, ddelta = b[i, j] - p[j] + p[i], db - responses[:, j] + responses[:, i]
else:
continue
numerator += delta * value
denominator += value * value
dnum += ddelta * value + delta * slope
dden += 2 * value * slope
gamma = numerator / denominator
dgamma = (dnum - gamma * dden) / denominator
else:
gamma = _error_rate(b, p, pairs)
dgamma = np.zeros(n)
for i, j in pairs:
db = (1 - responses[:, i]) * responses[:, j]
dgamma += (db * p[j] - b[i, j] * responses[:, j]) / p[j] ** 2
dgamma /= len(pairs)
gradient = np.zeros(n)
for i in range(q):
for j in range(q):
if i == j:
continue
di, dj = responses[:, i], responses[:, j]
db = (1 - di) * dj
if (i, j) in pairs:
expected = gamma * p[j]
derivative = dgamma * p[j] + gamma * dj
elif version == "original":
expected = (1 - p[i]) * p[j] * (1 - gamma)
derivative = (-di * p[j] + (1 - p[i]) * dj) * (1 - gamma) - (1 - p[i]) * p[
j
] * dgamma
elif (j, i) in pairs:
expected = p[j] - p[i] + gamma * p[i]
derivative = dj - di + dgamma * p[i] + gamma * di
else:
expected = (1 - p[i]) * p[j]
derivative = -di * p[j] + (1 - p[i]) * dj
gradient += 2 * (b[i, j] - expected) * (db - derivative)
gradient /= q * (q - 1)
return _diff(b, p, 1.0, pairs, gamma, version), float(gamma), gradient
[docs]
@dataclass(frozen=True)
class IITAVariance:
"""Normalized diff and first-order multinomial uncertainty.
asymptotic_variance is V for sqrt(N)*(diff_hat-diff), not V/N.
standard_error is sqrt(V/N) for a sample and None for a population.
influence is the mean-centered derivative in input response-row order;
its probability-weighted squared mean is V. It is read-only. Zero V
invalidates the first-order normal approximation (e.g. exact fit).
Relations are treated as fixed; selection on the same data is not
accounted for. No claim of valid inference after exploratory selection.
"""
normalized_diff: float
error_rate: float
asymptotic_variance: float
standard_error: float | None
n_respondents: float | None
influence: np.ndarray
def _variance(
data: ResponseMatrix,
relation: SurmiseRelation,
version: Version,
*,
population: bool,
) -> IITAVariance:
pairs = _pairs(relation, data.items)
probabilities = data.effective_counts / data.n_respondents
diff, gamma, gradient = _criterion_gradient(
data.patterns.astype(float), probabilities, pairs, version
)
influence = gradient - probabilities @ gradient
coefficient = float(probabilities @ (influence * influence))
influence.flags.writeable = False
return IITAVariance(
diff,
gamma,
coefficient,
None if population else float(np.sqrt(coefficient / data.n_respondents)),
None if population else data.n_respondents,
influence,
)
[docs]
def iita_variance(
data: ResponseMatrix,
relation: SurmiseRelation,
*,
version: Version = "minimized",
) -> IITAVariance:
"""Sample plug-in delta variance; counts represent integer respondents.
Only observed rows are needed: zero-probability multinomial cells
contribute nothing. Duplicate rows/aggregated frequencies give the same
result. Analytic differentiation includes the estimated gamma rate;
it is not treated as fixed. The original variant is supported too,
beyond the corrected/minimized scope of DAKS 2.1-3 variances.
"""
if not np.equal(data.effective_counts, np.floor(data.effective_counts)).all():
raise ValueError("Sample inference requires integer respondent frequencies.")
return _variance(data, relation, version, population=False)
[docs]
def population_iita_variance(
patterns: np.ndarray,
probabilities: np.ndarray,
relation: SurmiseRelation,
*,
items: list[str] | None = None,
version: Version = "minimized",
) -> IITAVariance:
"""Delta variance at a specified population distribution on responses.
Rows may describe only the positive support; probabilities must be
finite, nonnegative, and sum to one. Columns follow items (sorted
relation domain by default). Gamma is recomputed from the distribution,
so an inconsistent independently supplied error rate cannot enter.
"""
data = ResponseMatrix(
sorted(relation.items) if items is None else list(items),
np.asarray(patterns),
np.asarray(probabilities, dtype=float),
)
if not np.isclose(data.n_respondents, 1, rtol=0, atol=1e-10):
raise ValueError("Population probabilities must sum to one.")
return _variance(data, relation, version, population=True)
[docs]
@dataclass(frozen=True)
class PopulationIITAResult:
"""Population response probabilities and all candidate diff/error rates.
diff uses probabilities (sample diff divided by N squared), and arrays
use input items and lexicographic binary response order. Candidate
generation uses exact computed counterexample levels; near ties can
change the selection set under rounding. Supply a fixed selection_set
for controlled comparisons. No Monte Carlo sample is generated.
"""
items: tuple[str, ...]
patterns: np.ndarray
probabilities: np.ndarray
selection_set: tuple[SurmiseRelation, ...]
diff: tuple[float, ...]
error_rates: tuple[float, ...]
version: str
@property
def selected_index(self) -> int:
"""First minimum, with candidate order retained."""
return int(np.argmin(self.diff))
@property
def relation(self) -> SurmiseRelation:
"""Selected relation."""
return self.selection_set[self.selected_index]
[docs]
def population_iita(
structure: KnowledgeStructure | SurmiseRelation,
*,
beta: ItemParameter = 0.1,
eta: ItemParameter = 0.1,
pi: StateParameter = None,
version: Version = "minimized",
selection_set: Sequence[SurmiseRelation] | None = None,
data: ResponseMatrix | None = None,
max_patterns: int = 1_048_576,
max_memory_bytes: int = 512_000_000,
) -> PopulationIITAResult:
"""Evaluate IITA under a BLIM population, with optional sample candidates.
A relation is expanded to its compatible states. Defaults are uniform
state masses and homogeneous errors, as in DAKS::pop_iita. General
structures, item-specific errors and nonuniform pi are also supported.
Candidate relations come from population moments, supplied selection_set,
or data (mutually exclusive). All evaluated relations must be transitive
and have at least one nonreflexive pair. Enumeration is exponential.
"""
if version not in ("original", "corrected", "minimized"):
raise ValueError("Unknown IITA version.")
if data is not None and selection_set is not None:
raise ValueError("Specify either data or selection_set, not both.")
if isinstance(structure, SurmiseRelation):
_pattern_count(len(structure.items), max_patterns)
structure = KnowledgeStructure.from_surmise_relation(structure)
items = sorted(structure.domain)
count = _pattern_count(len(items), max_patterns)
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
if (
8 * (count * (6 * len(structure) + 4 * len(items)) + 4 * len(structure) * len(items))
> max_memory_bytes
):
raise MemoryError("Population IITA exceeds max_memory_bytes.")
patterns = _binary_patterns(0, count, len(items))
probabilities = predict_blim(
structure, patterns, beta=beta, eta=eta, pi=pi, max_memory_bytes=max_memory_bytes
).probabilities
p, b = _moments(patterns.astype(float), probabilities)
if data is not None:
if set(data.items) != structure.domain:
raise ValueError("data must match the population item domain.")
order = [data.items.index(q) for q in items]
b_generate = counterexamples(ResponseMatrix(items, data.patterns[:, order], data.counts))
else:
b_generate = b
if selection_set is None:
selection_set = [
SurmiseRelation(items, [(items[i], items[j]) for i, j in pairs])
for pairs in _inductive_generation(b_generate)
]
if not selection_set:
raise ValueError("selection_set must contain at least one nonempty relation.")
diffs, gammas = [], []
for relation in selection_set:
pairs = _pairs(relation, items)
gamma = (
_error_rate_minimized(b, p, pairs)
if version == "minimized"
else _error_rate(b, p, pairs)
)
diffs.append(_diff(b, p, 1.0, pairs, gamma, version))
gammas.append(gamma)
patterns.flags.writeable = False
return PopulationIITAResult(
tuple(items),
patterns,
probabilities,
tuple(selection_set),
tuple(diffs),
tuple(gammas),
version,
)
[docs]
@dataclass(frozen=True)
class IITAZTest:
"""Normal test of a fixed relation's diff or a difference of two diffs.
estimate is diff_1 (one relation) or diff_1-diff_2. dependence is 'single',
'paired' (same response data, including covariance), or 'independent'.
A zero standard error returns NaN z/p/interval with IITAInferenceWarning.
A single-relation null diff=0 has zero gradient in the population:
p/interval are NaN even if the plug-in variance is positive (z is then
descriptive). regular_null only flags these detected failures; True
does not establish the remaining asymptotic/model assumptions.
Confidence intervals are untruncated normal intervals, potentially
outside the parameter range. No post-selection correction is applied.
"""
estimate: float
standard_error: float
z: float
p_value: float
confidence_interval: tuple[float, float]
diff_values: tuple[float, ...]
alternative: str
null_difference: float
confidence: float
dependence: str
regular_null: bool
[docs]
def iita_z_test(
data: ResponseMatrix,
relation: SurmiseRelation,
*,
other_relation: SurmiseRelation | None = None,
other_data: ResponseMatrix | None = None,
version: Version = "minimized",
alternative: Literal["two-sided", "less", "greater"] = "two-sided",
null_difference: float = 0,
confidence: float = 0.95,
) -> IITAZTest:
"""Delta-method Z test, with correct pairing and normal tails.
With other_relation only, both criteria use the same respondents and
covariance is included. With other_data too, samples are assumed
independent. 'greater' tests estimate > null_difference using the upper
normal tail; 'less' uses the lower tail. Applicable to fixed relations,
not automatically to winners selected on these data. At exact fit the
first-order variance can vanish; normal calibration is then undefined.
In particular, a single relation tested against zero has a zero-gradient
null because diff is a sum of squares. This returns a warning and NaN
p-value/interval even with positive sample variance. A prespecified
positive single-relation null, or a difference of two fixed relations,
may admit the ordinary approximation under regularity assumptions.
"""
if alternative not in ("two-sided", "less", "greater"):
raise ValueError("Unknown alternative.")
if not np.isfinite(null_difference) or not 0 < confidence < 1:
raise ValueError("Require finite null_difference and confidence in (0,1).")
if other_data is not None and other_relation is None:
raise ValueError("other_data requires other_relation.")
if other_relation is None and null_difference < 0:
raise ValueError("A single relation's diff is nonnegative; the null cannot be negative.")
first = iita_variance(data, relation, version=version)
assert first.standard_error is not None
estimate, variance = first.normalized_diff, first.standard_error**2
dependence = "single"
values: tuple[float, ...] = (first.normalized_diff,)
if other_relation is not None:
second = iita_variance(
data if other_data is None else other_data, other_relation, version=version
)
estimate -= second.normalized_diff
values = (first.normalized_diff, second.normalized_diff)
if other_data is None:
influence = first.influence - second.influence
variance = float(
(data.effective_counts / data.n_respondents) @ influence**2 / data.n_respondents
)
dependence = "paired"
else:
assert second.standard_error is not None
variance += second.standard_error**2
dependence = "independent"
se = float(np.sqrt(variance))
regular_null = se > 0 and (other_relation is not None or null_difference > 0)
if se == 0:
warnings.warn(
"Zero delta variance: the first-order normal test is undefined.",
IITAInferenceWarning,
stacklevel=2,
)
z = p_value = float("nan")
interval = (float("nan"), float("nan"))
else:
z = (estimate - null_difference) / se
if alternative == "two-sided":
p_value = float(2 * norm.sf(abs(z)))
width = float(norm.ppf((1 + confidence) / 2)) * se
interval = (estimate - width, estimate + width)
elif alternative == "greater":
p_value = float(norm.sf(z))
interval = (estimate - float(norm.ppf(confidence)) * se, float("inf"))
else:
p_value = float(norm.cdf(z))
interval = (float("-inf"), estimate + float(norm.ppf(confidence)) * se)
if not regular_null:
warnings.warn(
"The single-relation null diff=0 has zero gradient: first-order normal "
"calibration is invalid even when the sample plug-in variance is positive. "
"Returning NaN p-value and interval; z is descriptive only.",
IITAInferenceWarning,
stacklevel=2,
)
p_value = float("nan")
interval = (float("nan"), float("nan"))
return IITAZTest(
estimate,
se,
z,
p_value,
interval,
values,
alternative,
null_difference,
confidence,
dependence,
regular_null,
)