Clamp default initial values into the model's declared bounds

docs/changelog.md

@@ -18,6 +18,8 @@
 
 8. **`hssm.load_data` now returns a `pd.DataFrame` unconditionally** (#1146). Its `dataset` argument is required, and the return type is no longer `pd.DataFrame | str`, which forced type-checkers (and users) to narrow away a `str` branch that existed only to print the dataset listing. Use the new **`hssm.list_data()`** to get the names of the built-in datasets as a `tuple[str, ...]` (mirroring `hssm.list_models()`). Breaking: `hssm.load_data()` with no argument now raises `TypeError` instead of returning a listing string.
 
+9. **Default initial values now respect a parameter's declared bounds** (#1293). The starting values in `INITVAL_SETTINGS` are shared across models, so a model whose bounds exclude one of them — for example `t` declared as `(0.25, 2.0)` against the shared default of `0.025` — started sampling at a point with `-inf` log-probability, which no sampler can move away from: every chain froze with 100% divergences and no error was raised. `process_initvals=True` now moves such a default to a point 5% of the bound width inside the violated endpoint and logs a warning naming the parameter and the substituted value. Bounds are read from the parameter itself, so both routes work: `include=[{"name": "t", "bounds": (0.25, 2.0)}]` and `model_config=ModelConfig(bounds=...)`. Initial values you supply yourself are never touched, and defaults already inside their bounds are unchanged. Note that `initval_jitter` (default `0.01`) is applied after this clamp and does not itself consult the bounds, so a bound narrower than about `0.2` can still be jittered outside.
+
 ### 0.4.0
 
 This version contains major breaking updates for HSSM. Please read the release notes below to migrate to HSSM 0.4.0.

src/hssm/base.py

@@ -68,6 +68,7 @@
     emit_parameterization_warnings,
     find_disconnected_free_rvs,
 )
+from .param.utils import _clamp_default_initval_to_bounds
 
 _logger = logging.getLogger("hssm")
 
@@ -204,6 +205,9 @@ class HSSMBase(ABC, DataValidatorMixin, MissingDataMixin):
         either `missing_data` or `deadline` is not `False`. Defaults to `None`.
     process_initvals : optional
         If `True`, the model will process the initial values. Defaults to `True`.
+        Processing also clamps a default initial value that falls outside its
+        parameter's declared `bounds` to a point just inside them; initial values
+        you supply yourself are used as given.
     initval_jitter : optional
         The jitter value for the initial values. Defaults to `0.01`.
     noncentered : optional
@@ -2052,10 +2056,20 @@ def _postprocess_initvals_deterministic(
             # If the user actively supplies a link function, the user
             # should also have supplied an initial value insofar it matters.
 
-            if self.params[self._get_prefix(name_tmp)].is_regression:
-                param_link_setting = self.link_settings
-            else:
+            param = self.params[self._get_prefix(name_tmp)]
+            # The parameter's own link decides which default scale applies. A
+            # regression may override the model-wide link_settings: under an
+            # identity link the default is natural-scale, where the declared
+            # bounds apply; under HSSM's own log or gen_logit links it is the
+            # link-space default. Any other link keeps the model-wide setting,
+            # since its scale is not known here.
+            link_name = getattr(param.link, "name", param.link)
+            if not param.is_regression or link_name == "identity":
                 param_link_setting = None
+            elif link_name in ("log", "gen_logit"):
+                param_link_setting = "log_logit"
+            else:
+                param_link_setting = self.link_settings
             if name_tmp in initval_settings[param_link_setting].keys():
                 if self._check_if_initval_user_supplied(name_tmp):
                     _logger.info(
@@ -2065,11 +2079,17 @@ def _postprocess_initvals_deterministic(
                     )
                     continue
 
-                # Apply specific settings from initval_settings dictionary
+                # Apply specific settings from initval_settings dictionary,
+                # clamped into the parameter's declared bounds (natural-scale
+                # defaults only; log_logit defaults are link-space and any
+                # user-supplied value was already skipped above).
+                value = initval_settings[param_link_setting][name_tmp]
+                if param_link_setting is None:
+                    value = _clamp_default_initval_to_bounds(
+                        value, name_tmp, param.bounds
+                    )
                 dtype = self._initvals[name_tmp].dtype
-                self._initvals[name_tmp] = np.array(
-                    initval_settings[param_link_setting][name_tmp]
-                ).astype(dtype)
+                self._initvals[name_tmp] = np.array(value).astype(dtype)
 
     def _get_prefix(self, name_str: str) -> str:
         """Resolve parameter prefix, handling underscore-containing RL param names.

src/hssm/hssm.py

@@ -204,6 +204,9 @@ class HSSM(HSSMBase):
         either `missing_data` or `deadline` is not `False`. Defaults to `None`.
     process_initvals : optional
         If `True`, the model will process the initial values. Defaults to `True`.
+        Processing also clamps a default initial value that falls outside its
+        parameter's declared `bounds` to a point just inside them; initial values
+        you supply yourself are used as given.
     initval_jitter : optional
         The jitter value for the initial values. Defaults to `0.01`.
     **kwargs

src/hssm/param/utils.py

@@ -1,8 +1,12 @@
 """Utility functions for the parameter classes."""
 
+import logging
+
 import bambi as bmb
 import numpy as np
 
+_logger = logging.getLogger("hssm")
+
 
 def validate_bounds(bounds: tuple[float, float]) -> None:
     """Validate the bounds."""
@@ -41,3 +45,65 @@ def _make_default_prior(bounds: tuple[float, float] | None) -> bmb.Prior:
         prior = bmb.Prior(name="Uniform", lower=lower, upper=upper)
 
     return prior
+
+
+def _clamp_default_initval_to_bounds(
+    value: float, name: str, bounds: tuple[float, float] | None
+) -> float:
+    """Clamp a default initial value into a parameter's declared bounds.
+
+    The value is returned unchanged when it lies strictly inside ``bounds``;
+    otherwise it is moved to the point 5% of the bound width inside the violated
+    endpoint. The defaults in ``INITVAL_SETTINGS`` are shared across models, so a
+    model's declared bounds may exclude them, and a start outside the bounds has
+    -inf log-probability from which sampling cannot move. This applies only to
+    defaults on the natural scale (the ``None``-link branch); user-supplied
+    initial values are never touched.
+
+    Parameters
+    ----------
+    value
+        The default initial value, on the natural scale of the parameter.
+    name
+        The parameter name as it appears in the model's initial point. Used only
+        in the warning emitted when the value is moved.
+    bounds
+        The parameter's ``(lower, upper)`` bounds, or ``None`` if the parameter
+        declares none.
+
+    Returns
+    -------
+        ``value`` itself when it lies strictly inside ``bounds``, or when
+        ``bounds`` is ``None``; otherwise a finite value strictly inside
+        ``bounds``, a distance of 5% of the bound width from the violated
+        endpoint.
+    """
+    if bounds is None:
+        return value
+    lower, upper = bounds
+    if lower < value < upper:
+        return value
+    # A one-sided bound - a: (0, inf), st: (0, inf) - makes the width infinite,
+    # and a margin proportional to it would place the result at +/-inf. Scale
+    # the margin off whichever endpoint is finite instead, so the clamped value
+    # is always finite and strictly inside the bounds.
+    width = upper - lower
+    if np.isfinite(width):
+        margin = 0.05 * width
+    else:
+        # The 1.0 floor keeps the margin strictly positive when the only finite
+        # endpoint is 0.0, as in (0, inf); a margin of 0 there would clamp onto
+        # the excluded boundary itself.
+        finite_endpoints = [abs(b) for b in (lower, upper) if np.isfinite(b)]
+        margin = 0.05 * max(finite_endpoints + [1.0])
+    clamped = float(np.clip(value, lower + margin, upper - margin))
+    _logger.warning(
+        "Default initial value %s for %s lies outside the declared bounds "
+        "(%s, %s); using %s instead. Pass an explicit initval to override.",
+        value,
+        name,
+        lower,
+        upper,
+        clamped,
+    )
+    return clamped

tests/test_initvals.py

@@ -3,6 +3,8 @@
 import hssm
 import logging
 
+from hssm.param.utils import _clamp_default_initval_to_bounds
+
 hssm.set_floatX("float32", update_jax=True)
 logger = logging.getLogger("hssm")
 
@@ -309,3 +311,156 @@ def test_process_no_process(caplog):
         model_on.initvals != model_off.initvals
     ), """Initial values should not be the same when
     initval processing is turned off vs. turned on."""
+
+
+def test_default_initval_clamped_into_bounds():
+    """A default initval outside the declared bounds is moved inside them."""
+    # Below the lower bound -> moved inside by 5% of the bound width.
+    clamped_low = _clamp_default_initval_to_bounds(0.025, "t", (0.25, 2.25))
+    assert clamped_low == 0.25 + 0.05 * 2.0
+
+    # Above the upper bound -> moved inside.
+    clamped_high = _clamp_default_initval_to_bounds(5.0, "a", (0.3, 2.5))
+    assert 0.3 < clamped_high < 2.5
+
+    # Inside the bounds -> returned unchanged.
+    assert _clamp_default_initval_to_bounds(0.4, "t", (0.25, 2.25)) == 0.4
+
+    # No bounds declared -> returned unchanged.
+    assert _clamp_default_initval_to_bounds(0.025, "t", None) == 0.025
+
+
+@pytest.mark.parametrize(
+    ("bounds", "value"),
+    [
+        ((0.0, np.inf), -1.0),
+        ((0.3, np.inf), 0.025),
+        ((-np.inf, 1.0), 5.0),
+    ],
+)
+def test_default_initval_clamped_into_one_sided_bounds(bounds, value):
+    """Bounds with an infinite endpoint still yield a finite interior value.
+
+    Shipped configs declare one-sided bounds - a: (0, inf) and t: (0, inf) in
+    the analytical likelihoods, sz/st: (0, inf) in full_ddm - and a user may
+    merge their own. A margin proportional to an infinite width would return
+    +/-inf, which is a worse starting value than the unclamped default.
+    """
+    lower, upper = bounds
+    result = _clamp_default_initval_to_bounds(value, "t", bounds)
+    assert np.isfinite(result)
+    assert lower < result < upper
+
+
+def test_doubly_infinite_bounds_leave_finite_default_untouched():
+    """A doubly-infinite bound already contains every finite default."""
+    assert _clamp_default_initval_to_bounds(0.0, "v", (-np.inf, np.inf)) == 0.0
+
+
+def test_clamp_warns_naming_the_parameter_and_replacement(caplog):
+    """Moving a default is announced, so a surprising start is traceable."""
+    caplog.set_level(logging.WARNING, logger="hssm")
+
+    _clamp_default_initval_to_bounds(0.025, "t", (0.25, 2.25))
+
+    assert len(caplog.records) == 1
+    message = caplog.records[0].getMessage()
+    assert "Default initial value 0.025 for t" in message
+    assert "outside the declared bounds (0.25, 2.25)" in message
+    assert "using 0.35 instead" in message
+    assert "Pass an explicit initval to override." in message
+
+
+@pytest.mark.parametrize(
+    ("model", "name", "bounds"),
+    [
+        # A finite two-sided bound excluding the shared default t = 0.025.
+        ("ddm", "t", (0.25, 2.0)),
+        ("ddm_sdv", "t", (0.25, 2.0)),
+        ("angle", "t", (0.25, 2.0)),
+        # A one-sided bound, whose infinite width still yields a finite start.
+        ("ddm", "t", (0.3, np.inf)),
+        # A bound that already contains the default, which stays untouched.
+        ("ddm", "t", (0.0, 2.0)),
+        # An _Intercept name resolves to its base parameter's bounds.
+        ("ddm", "t_Intercept", (0.25, 2.0)),
+    ],
+)
+def test_declared_bounds_reach_the_initial_value(cavanagh_test, model, name, bounds):
+    """Bounds passed through ``include=`` land on the resulting initial value.
+
+    ``include=[{"name": ..., "bounds": ...}]`` stores bounds on the ``Param``,
+    not on ``model_config``, and it is the only route by which a user of a
+    shipped model can declare a bound that excludes a default initval.
+    """
+    lower, upper = bounds
+    spec: dict = {"name": "t", "bounds": bounds}
+    if name.endswith("_Intercept"):
+        spec["formula"] = "t ~ 1"
+
+    fitted = hssm.HSSM(
+        data=cavanagh_test.iloc[:12],
+        model=model,
+        include=[spec],
+        p_outlier=0.0,
+        prior_settings=None,
+        link_settings=None,
+        process_initvals=True,
+        initval_jitter=0.0,
+    )
+
+    assert fitted.params["t"].bounds == bounds
+    initval = fitted._initvals[name]
+    assert np.isfinite(initval)
+    assert lower < initval < upper
+
+    # A default that is already inside its bounds is passed through as-is.
+    default = hssm.defaults.INITVAL_SETTINGS[None][name]
+    if lower < default < upper:
+        assert initval == np.array(default).astype(initval.dtype)
+
+
+def test_identity_link_override_clamps_into_bounds(caplog, cavanagh_test):
+    """An explicit identity link puts the default on the natural scale.
+
+    ``link_settings="log_logit"`` is model-wide, but a regression may override
+    the link for one parameter. Under identity the default is natural-scale, so
+    the declared bounds apply to it and it is clamped into them.
+    """
+    model = hssm.HSSM(
+        data=cavanagh_test,
+        model="ddm",
+        link_settings="log_logit",
+        initval_jitter=0,
+        include=[
+            {
+                "name": "t",
+                "formula": "t ~ 1 + stim",
+                "link": "identity",
+                "bounds": (0.25, 2.0),
+            }
+        ],
+    )
+    assert getattr(model.params["t"].link, "name", model.params["t"].link) == "identity"
+    initval = float(np.asarray(model._initvals["t_Intercept"]))
+    # the link-space default of -4.0 would be far outside these bounds
+    assert 0.25 < initval < 2.0
+    assert initval == pytest.approx(0.25 + 0.05 * (2.0 - 0.25))
+
+
+def test_user_log_link_gets_link_space_default(cavanagh_test):
+    """A regression's own log link selects the link-space default.
+
+    The model-wide ``link_settings`` is ``None`` here, so the natural-scale
+    table would have applied before; the parameter's effective link decides.
+    """
+    model = hssm.HSSM(
+        data=cavanagh_test,
+        model="ddm",
+        link_settings=None,
+        initval_jitter=0,
+        include=[{"name": "a", "formula": "a ~ 1 + stim", "link": "log"}],
+    )
+    assert getattr(model.params["a"].link, "name", model.params["a"].link) == "log"
+    # the log-space default, i.e. a = exp(0) = 1, not the natural-scale 1.5
+    assert float(np.asarray(model._initvals["a_Intercept"])) == 0.0