Add the ddm_normal_st model config (config part of #361, stacked on #349)

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",

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 {

ssms/config/kde_constants.py

@@ -2,6 +2,7 @@
     "full_ddm_rv",
     "full_ddm",
     "ddm_st",
+    "ddm_normal_st",
     "ddm_truncnormt",
     "ddm_rayleight",
 ]

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))