diff --git a/ssms/config/_modelconfig/__init__.py b/ssms/config/_modelconfig/__init__.py index dde6d3f5..f67a9bd0 100644 --- a/ssms/config/_modelconfig/__init__.py +++ b/ssms/config/_modelconfig/__init__.py @@ -50,6 +50,7 @@ get_ddm_par2_weibull_no_bias_config, ) from .ddm_random import ( + get_ddm_normal_st_config, get_ddm_rayleight_config, get_ddm_sdv_config, get_ddm_st_config, @@ -193,6 +194,7 @@ def get_model_config(): configs = { "ddm": get_ddm_config(), "ddm_st": get_ddm_st_config(), + "ddm_normal_st": get_ddm_normal_st_config(), "ddm_truncnormt": get_ddm_truncnormt_config(), "ddm_rayleight": get_ddm_rayleight_config(), "ddm_sdv": get_ddm_sdv_config(), @@ -319,6 +321,7 @@ def get_model_config(): "get_weibull_config", "get_full_ddm_config", "get_ddm_st_config", + "get_ddm_normal_st_config", "get_ddm_truncnormt_config", "get_ddm_rayleight_config", "get_ddm_sdv_config", diff --git a/ssms/config/_modelconfig/ddm_random.py b/ssms/config/_modelconfig/ddm_random.py index 8a4e41df..aaaf1691 100644 --- a/ssms/config/_modelconfig/ddm_random.py +++ b/ssms/config/_modelconfig/ddm_random.py @@ -56,6 +56,102 @@ def get_ddm_st_config(): } +def get_ddm_normal_st_config(): + """Get configuration for DDM with untruncated-normal non-decision time. + + Description + ----------- + Trial-wise non-decision time is drawn from an untruncated Normal kernel, + ndt ~ Normal(t, st), so ``st`` is the kernel's STANDARD DEVIATION. That is + the contrast with ``ddm_st``, which exposes the same five parameters over + the same bounds but draws ndt ~ Uniform(t - st, t + st), where ``st`` is a + half-width (SD = st / sqrt(3)). The two are not interchangeable at equal + ``st``: equal dispersion needs st_uniform = sqrt(3) * st_normal. + + Parameters + ---------- + v : float + Drift rate. Valid range: [-3.0, 3.0] + a : float + Boundary separation. Valid range: [0.3, 2.5] + z : float + Starting point as a fraction of ``a``. Valid range: [0.3, 0.7] + t : float + Mean non-decision time, in seconds. Valid range: [0.25, 2.25] + st : float + Standard deviation of the non-decision-time kernel, in seconds. + Valid range: [1e-3, 0.25] + + Model Characteristics + --------------------- + - Number of choices: 2 + - Boundary type: constant + - Drift type: constant + - Key assumption: the ndt kernel is Normal and untruncated. + + Notes + ----- + - The kernel's support is unbounded, so ndt has no floor and can fall below + 0 near the lower ``t`` bound, yielding a non-positive RT. Positive RTs + are therefore never guaranteed, anywhere in the box - only made rare, at + a rate that falls off steeply in t / st: measured 8.9% at the corner + t=0.25, st=0.25, a=0.3, 0.013% on the 3-sigma edge 3 * st = t + (st = 0.0833 at t = 0.25), and none in 6e5 draws at st = 0.05. Consumers + that cannot represent a non-positive RT (KDE-based likelihood estimation + in particular) must filter, not assume. ``ddm_st``'s bounded kernel + cannot produce one anywhere in its box. + - ``param_bounds`` are deliberately identical to ``ddm_st``'s so the two + models span the same box; the ``st`` bounds therefore carry different + meanings in the two models. + + See Also + -------- + get_ddm_st_config : same parameters, uniform (bounded) ndt kernel. + get_ddm_truncnormt_config : Normal ndt kernel truncated at 0. + """ + return { + "name": "ddm_normal_st", + "params": ["v", "a", "z", "t", "st"], + "param_bounds": [ + [-3.0, 0.3, 0.3, 0.25, 1e-3], + [3.0, 2.5, 0.7, 2.25, 0.25], + ], + "boundary_name": "constant", + "boundary": bf.constant, + "n_params": 5, + "default_params": [0.0, 1.0, 0.5, 0.25, 1e-3], + "nchoices": 2, + "choices": [-1, 1], + "n_particles": 1, + "simulator": cssm.full_ddm_rv, + "simulator_fixed_params": { + "z_dist": functools.partial(sps.norm.rvs, loc=0, scale=0), + "v_dist": functools.partial(sps.norm.rvs, loc=0, scale=0), + }, + "simulator_param_mappings": { + "t_dist": lambda st: functools.partial(sps.norm.rvs, loc=0, scale=st), + }, + "parameter_transforms": { + "sampling": [], + "simulation": [ + LambdaAdaptation( + lambda theta, cfg, n: ( + theta.update( + { + "z_dist": cfg["simulator_fixed_params"]["z_dist"], + "v_dist": cfg["simulator_fixed_params"]["v_dist"], + } + ) + or theta + ), + name="set_fixed_params", + ), + ApplyMapping("st", "t_dist", "t_dist"), + ], + }, + } + + def get_ddm_truncnormt_config(): """Get configuration for DDM with truncated normal non-decision time.""" return { diff --git a/ssms/config/kde_constants.py b/ssms/config/kde_constants.py index c117fb61..8c55f2f1 100644 --- a/ssms/config/kde_constants.py +++ b/ssms/config/kde_constants.py @@ -2,6 +2,7 @@ "full_ddm_rv", "full_ddm", "ddm_st", + "ddm_normal_st", "ddm_truncnormt", "ddm_rayleight", ] diff --git a/tests/test_ddm_normal_st_config.py b/tests/test_ddm_normal_st_config.py new file mode 100644 index 00000000..b050c025 --- /dev/null +++ b/tests/test_ddm_normal_st_config.py @@ -0,0 +1,158 @@ +"""Tests for the ddm_normal_st model configuration. + +ddm_normal_st draws non-decision time from an untruncated Normal kernel, so `st` +is a standard deviation rather than the uniform half-width `ddm_st` uses. These +tests pin the two consequences: the dispersion semantics of `st`, and the fact +that the kernel has no floor, which makes RT positivity a matter of measured +exposure rather than a guarantee. +""" + +import numpy as np +import pytest +import scipy.stats as sps + +from ssms.basic_simulators.simulator import simulator +from ssms.config import KDE_NO_DISPLACE_T, model_config +from ssms.config.generator_config.data_generator_config import ( + get_default_generator_config, +) +from ssms.dataset_generators.lan_mlp import TrainingDataGenerator +from ssms.support_utils.kde_class import LogKDE + +# Every non-positive-RT rate below is a Monte Carlo estimate: the rate is a +# property of an unbounded kernel, and `random_state` does not pin the ndt +# draws (t_dist samples from scipy's global RNG). Tolerances are set from a +# 3-seed x 200_000-sample measurement, whose seed-to-seed spread at the box +# corner is +-0.06 percentage points. +N_SAMPLES = 100_000 + +# The lower `t` bound with the widest `st`, the worst corner of the box. +BOX_CORNER = {"v": 0.0, "a": 0.3, "z": 0.5, "t": 0.25, "st": 0.25} + +# `st` at which the kernel's 3-sigma floor sits exactly on rt = 0 at t = 0.25. +THREE_SIGMA_ST = BOX_CORNER["t"] / 3.0 + + +def _frac_nonpositive(theta, model="ddm_normal_st", n_samples=N_SAMPLES): + """Fraction of simulated RTs that a KDE cannot take the log of.""" + rts = simulator(model=model, theta=theta, n_samples=n_samples)["rts"] + return float(np.mean(rts <= 0.0)) + + +def test_ddm_normal_st_is_registered_over_ddm_st_s_box(): + """ddm_normal_st exposes ddm_st's parameters over ddm_st's bounds.""" + config = model_config["ddm_normal_st"] + assert config["name"] == "ddm_normal_st" + assert config["params"] == model_config["ddm_st"]["params"] + np.testing.assert_array_equal( + np.asarray(config["param_bounds"], dtype=float), + np.asarray(model_config["ddm_st"]["param_bounds"], dtype=float), + ) + + +def test_ddm_normal_st_simulates(): + """Smoke test at a parameter set the 3-sigma floor keeps clear of zero.""" + result = simulator( + model="ddm_normal_st", + theta={"v": 1.0, "a": 1.5, "z": 0.5, "t": 1.0, "st": 0.1}, + n_samples=100, + ) + + assert result["rts"].shape == (100, 1) + assert result["choices"].shape == (100, 1) + assert set(np.unique(result["choices"])).issubset({-1, 1}) + assert np.all(result["rts"] > 0) + + +def test_st_is_a_standard_deviation_not_a_half_width(): + """`st` is the ndt kernel's SD; ddm_st's `st` is a half-width (SD st/sqrt(3)).""" + st = 0.15 + normal = sps.norm( + **model_config["ddm_normal_st"]["simulator_param_mappings"]["t_dist"]( + st + ).keywords + ) + uniform = sps.uniform( + **model_config["ddm_st"]["simulator_param_mappings"]["t_dist"](st).keywords + ) + + assert normal.std() == pytest.approx(st) + assert uniform.std() == pytest.approx(st / np.sqrt(3.0)) + + # The half-width kernel floors ndt at t - st; the Normal one has no floor, + # which is the whole of the difference in downstream RT handling. + assert uniform.support() == pytest.approx((-st, st)) + assert normal.support()[0] == -np.inf + + +def test_bounded_sibling_never_emits_a_non_positive_rt(): + """ddm_st's floor is t - st >= 0 over its whole box, and rt adds a positive dt.""" + rts = simulator(model="ddm_st", theta=BOX_CORNER, n_samples=N_SAMPLES)["rts"] + assert np.all(rts > 0) + + +def test_non_positive_rts_are_a_declared_exposure_at_the_box_corner(): + """The unbounded kernel puts a large minority of RTs at or below zero here. + + Consumers that cannot represent a non-positive RT must filter, not assume. + Measured 8.92% (8.88-9.00 over 3 seeds x 200_000 samples). + """ + assert 0.05 < _frac_nonpositive(BOX_CORNER) < 0.15 + + +def test_non_positive_rt_exposure_falls_but_never_reaches_zero(): + """At st = t/3 the 3-sigma floor sits on zero, and exposure is ~1e-4, not 0. + + This is why `3 * st <= t` is a rate bound and not a positivity guarantee. + Measured 0.0128% at st = t/3 and 0 in 600_000 draws at st = 0.05. + """ + three_sigma = _frac_nonpositive( + dict(BOX_CORNER, st=THREE_SIGMA_ST), n_samples=200_000 + ) + assert 0.0 < three_sigma < 1e-3 + + assert _frac_nonpositive(dict(BOX_CORNER, st=0.05)) < 1e-4 + + +def test_generator_refuses_displace_t_for_ddm_normal_st(): + """An unbounded ndt kernel has no support edge to displace by, as for ddm_st.""" + assert "ddm_normal_st" in KDE_NO_DISPLACE_T + + generator_config = get_default_generator_config("lan") + generator_config["model"] = "ddm_normal_st" + generator_config["pipeline"]["n_cpus"] = 1 + generator_config["estimator"]["displace_t"] = True + + with pytest.warns(UserWarning, match="displace_t is True"): + generator = TrainingDataGenerator( + config=generator_config, model_config=model_config["ddm_normal_st"] + ) + + assert generator.generator_config["estimator"]["displace_t"] is False + + +def test_logkde_keeps_both_choice_groups_at_the_box_corner(): + """LogKDE drops the non-positive RTs here and keeps both choice groups. + + Unfiltered, one such RT makes its choice's bandwidth NaN, the choice + "no_base_data", and every label for that choice the same constant whatever + the rt. `displace_t=False` is the path the generator takes for this model. + """ + sim = simulator(model="ddm_normal_st", theta=BOX_CORNER, n_samples=2_000) + for choice in (-1, 1): + assert np.any(sim["rts"][sim["choices"] == choice] <= 0) + + kde = LogKDE(simulator_data=sim, displace_t=False) + + assert len(kde.bandwidths) == 2 + for bandwidth in kde.bandwidths: + assert bandwidth != "no_base_data" + assert np.isfinite(bandwidth) and bandwidth > 0 + + # The label LogKDE gives every rt of a choice it has no base data for. + no_base_data_label = -np.log(kde.data["n_trials"] * sim["metadata"]["max_t"]) + rts = np.array([0.1, 0.3, 0.6]) + for choice in (-1, 1): + labels = kde.kde_eval({"rts": rts, "choices": np.full(rts.shape, choice)}) + assert np.all(np.isfinite(labels)) + assert not np.any(np.isclose(labels, no_base_data_label))