From cd893e1baed134155a48d1ffab7b68675ab1054a Mon Sep 17 00:00:00 2001 From: "eitan_perlin@brown.edu" Date: Mon, 31 Aug 2026 13:01:58 -0400 Subject: [PATCH] Clamp default initial values into the parameter's declared bounds - HSSM's default initvals are fixed values in INITVAL_SETTINGS. When a parameter's declared bounds exclude its default, sampling starts outside the support: measured 100% divergences and rhat nan for a model declaring t >= 0.25 against the default t = 0.025. - The bounds are read from the Param, not from model_config. This is the route users actually take: bounds passed as include=[{"name": "t", "bounds": (0.25, 2.0)}] land on the Param, and only the model_config=ModelConfig(bounds=...) route merges into model_config. Reading model_config.bounds would leave the common case unfixed. The Param is already built in this loop, so the lookup is a hoist rather than an addition, and _get_prefix does the name resolution instead of a local replace("_Intercept", ""). - A default already inside its bounds is returned unchanged, so every shipped model keeps byte-identical initvals; verified by exact comparison, not tolerance. Out-of-bounds defaults move to 5% of the bound width inside, and the parameter is logged. - Verified across ddm, ddm_sdv and angle, over finite two-sided bounds, one-sided (0.3, inf), a regression _Intercept, and an already-inside control: the new end-to-end test fails 5 of 6 cases before this change and passes all 6 after, with the control byte-identical throughout. - The helper lives in param/utils.py beside validate_bounds and _make_default_prior, which already handles the same four infinite-endpoint cases. - Known follow-up, not addressed here: initval_jitter is applied after this clamp and does not consult bounds, so for a bound narrower than 0.2 a clamped value can be jittered back out. Co-Authored-By: Claude Fable 5 - The default's scale is chosen from the parameter's own link, not from the model-wide link_settings, which a regression may override. Identity means a natural-scale default, where the declared bounds apply; HSSM's own log and gen_logit links mean the link-space default; any other custom link keeps the model-wide setting, since its scale is not known here. Previously a log_logit model with an explicit identity override kept the link-space default - measured t_Intercept = -4.0 against bounds (0.25, 2.0) - and a link_settings=None model with an explicit log link got the natural-scale default on a log scale, a_Intercept = 1.5, i.e. a = 4.48. Both now resolve from the effective link; every other case is byte-identical. --- docs/changelog.md | 2 + src/hssm/base.py | 34 +++++++-- src/hssm/hssm.py | 3 + src/hssm/param/utils.py | 66 +++++++++++++++++ tests/test_initvals.py | 155 ++++++++++++++++++++++++++++++++++++++++ 5 files changed, 253 insertions(+), 7 deletions(-) diff --git a/docs/changelog.md b/docs/changelog.md index 08eda7404..d9f0d162e 100644 --- a/docs/changelog.md +++ b/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. diff --git a/src/hssm/base.py b/src/hssm/base.py index 59572b9d9..02cb5e16b 100644 --- a/src/hssm/base.py +++ b/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. diff --git a/src/hssm/hssm.py b/src/hssm/hssm.py index aad7dcf35..106851c3d 100644 --- a/src/hssm/hssm.py +++ b/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 diff --git a/src/hssm/param/utils.py b/src/hssm/param/utils.py index 6d630f673..718836f54 100644 --- a/src/hssm/param/utils.py +++ b/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 diff --git a/tests/test_initvals.py b/tests/test_initvals.py index b485c269b..c277ed8fc 100644 --- a/tests/test_initvals.py +++ b/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