From ad27b9b2f8e52e095999ef93fec445b671152c7e Mon Sep 17 00:00:00 2001 From: Alexander Fengler Date: Sun, 30 Aug 2026 19:00:54 -0400 Subject: [PATCH 1/4] feat: add bounded group-location prior factory (#1269) --- src/hssm/prior.py | 47 +++++++++++++++++--- tests/unit/test_prior.py | 94 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 136 insertions(+), 5 deletions(-) diff --git a/src/hssm/prior.py b/src/hssm/prior.py index 8b6dcdb0f..2d48cbb3c 100644 --- a/src/hssm/prior.py +++ b/src/hssm/prior.py @@ -255,6 +255,39 @@ def _is_identity_link(link: str | bmb.Link | None) -> bool: return link.name == "identity" +def _has_finite_bounds(bounds: tuple[float, float] | None) -> bool: + """Return whether at least one response-scale bound is finite.""" + return bounds is not None and any(np.isfinite(bound) for bound in bounds) + + +def _make_bounded_group_intercept_prior( + bounds: tuple[float, float], +) -> Prior: + """Build a native hierarchical ``TruncatedNormal`` within ``bounds``. + + The native PyMC distribution keeps every hierarchical argument visible to + Bambi. This is intentionally different from ``Prior(..., bounds=...)``, whose + custom truncation wrapper is suitable for ordinary coefficients but hides the + hyperprior arguments that Bambi requires for group-specific terms. + """ + lower, upper = bounds + mu = mean(bounds) if np.all(np.isfinite(bounds)) else 0.0 + location = Prior( + "TruncatedNormal", + mu=mu, + sigma=0.25, + lower=lower, + upper=upper, + ) + return Prior( + "TruncatedNormal", + mu=location, + sigma=generate_prior("Weibull"), + lower=lower, + upper=upper, + ) + + # AF-TODO: Docstring could benefit from some more details here. def get_default_prior( term_type: str, @@ -268,10 +301,10 @@ def get_default_prior( * common_intercept: Bounded Normal prior (N(mean(bounds), 0.25)). * common: Normal prior (N(0, 0.25)). - * group_intercept: Normal prior N(N(0, 0.25), Weibull(1.5, 0.3)). Under a - transformed link this correctly lives on an unbounded predictor scale. Finite - coefficient bounds for identity-linked group-only intercepts are not yet - supported; see HSSM #1269. + * group_intercept: Under identity with a finite response bound, a hierarchical + TruncatedNormal whose location and group coefficients share that support. + Otherwise, Normal prior N(N(0, 0.25), Weibull(1.5, 0.3)). Transformed-link + coefficients correctly remain unbounded on the predictor scale. * group_specific: Normal prior N(N(0, 0.25), Weibull(1.5, 0.3). This function is taken from bambi.priors.prior.py and modified to handle hssm- @@ -311,7 +344,11 @@ def get_default_prior( "Normal", mu=mean(bounds), sigma=0.25, bounds=bounds ) elif term_type == "group_intercept": - prior = generate_prior("Normal", mu="Normal", sigma="Weibull") + if _is_identity_link(link) and _has_finite_bounds(bounds): + assert bounds is not None + prior = _make_bounded_group_intercept_prior(bounds) + else: + prior = generate_prior("Normal", mu="Normal", sigma="Weibull") elif term_type == "group_specific": prior = generate_prior("Normal", mu="Normal", sigma="Weibull") elif term_type in ["group_intercept_with_common", "group_specific_with_common"]: diff --git a/tests/unit/test_prior.py b/tests/unit/test_prior.py index 15e8584f1..3029db714 100644 --- a/tests/unit/test_prior.py +++ b/tests/unit/test_prior.py @@ -18,6 +18,7 @@ ) from hssm.prior import ( HDDM_SETTINGS_GROUP, + _has_finite_bounds, _is_identity_link, get_default_prior, get_hddm_default_prior, @@ -150,6 +151,20 @@ def test_transformed_link_classification(self, link): """Classify non-identity strings and link objects as transformed.""" assert not _is_identity_link(link) + @pytest.mark.parametrize( + ("bounds", "expected"), + [ + pytest.param(None, False, id="missing"), + pytest.param((-np.inf, np.inf), False, id="unbounded"), + pytest.param((0.0, np.inf), True, id="lower-bounded"), + pytest.param((-np.inf, 1.0), True, id="upper-bounded"), + pytest.param((0.0, 1.0), True, id="finite"), + ], + ) + def test_finite_bound_classification(self, bounds, expected): + """Recognize one- and two-sided finite response bounds.""" + assert _has_finite_bounds(bounds) is expected + @pytest.mark.parametrize("link", IDENTITY_LINKS) @pytest.mark.parametrize( ("bounds", "expected_args"), @@ -178,6 +193,85 @@ def test_generic_common_intercept_transformed_link(self, link): _assert_prior_spec(prior, "Normal", {"mu": 0.0, "sigma": 0.25}, bounds=None) + @pytest.mark.parametrize("link", IDENTITY_LINKS) + @pytest.mark.parametrize( + ("bounds", "location"), + [ + pytest.param((-2.0, 3.0), 0.5, id="finite"), + pytest.param((0.3, np.inf), 0.0, id="lower-bounded"), + pytest.param((-np.inf, 4.0), 0.0, id="upper-bounded"), + ], + ) + def test_generic_group_intercept_uses_native_bounded_hierarchy( + self, link, bounds, location + ): + """Bound both the location hyperprior and identity-scale group values.""" + prior = get_default_prior("group_intercept", "x", bounds, link) + + _assert_prior_tree( + prior, + { + "dist": "TruncatedNormal", + "mu": { + "dist": "TruncatedNormal", + "mu": location, + "sigma": 0.25, + "lower": bounds[0], + "upper": bounds[1], + }, + "sigma": { + "dist": "Weibull", + "alpha": 1.5, + "beta": 0.3, + }, + "lower": bounds[0], + "upper": bounds[1], + }, + ) + assert isinstance(prior, Prior) + assert prior.dist is None + assert not prior.is_truncated + assert prior.bounds is None + + @pytest.mark.parametrize("bounds", [None, (-np.inf, np.inf)]) + @pytest.mark.parametrize("link", IDENTITY_LINKS) + def test_generic_unbounded_group_intercept_retains_normal_hierarchy( + self, link, bounds + ): + """Keep the existing Normal hierarchy when identity has no finite bound.""" + prior = get_default_prior("group_intercept", "x", bounds, link) + + _assert_prior_tree( + prior, + { + "dist": "Normal", + "mu": {"dist": "Normal", "mu": 0.0, "sigma": 0.25}, + "sigma": { + "dist": "Weibull", + "alpha": 1.5, + "beta": 0.3, + }, + }, + ) + + @pytest.mark.parametrize("link", TRANSFORMED_LINKS) + def test_generic_transformed_group_intercept_remains_unbounded(self, link): + """Do not apply response bounds to transformed predictor coefficients.""" + prior = get_default_prior("group_intercept", "x", (0.0, 1.0), link) + + _assert_prior_tree( + prior, + { + "dist": "Normal", + "mu": {"dist": "Normal", "mu": 0.0, "sigma": 0.25}, + "sigma": { + "dist": "Weibull", + "alpha": 1.5, + "beta": 0.3, + }, + }, + ) + @pytest.mark.parametrize("link", IDENTITY_LINKS) @pytest.mark.parametrize( ("param", "bounds", "name", "expected_args"), From 7f4c9ba44c4b253739ab778538847b34a3c79bee Mon Sep 17 00:00:00 2001 From: Alexander Fengler Date: Sun, 30 Aug 2026 19:10:42 -0400 Subject: [PATCH 2/4] fix: bound generated identity group locations (#1269) --- src/hssm/param/regression_param.py | 79 ++++++++--- tests/unit/param/test_params.py | 28 +++- tests/unit/param/test_regression_param.py | 127 +++++++++++++++--- .../param/test_unmatched_group_prior_graph.py | 98 ++++++++++++++ 4 files changed, 287 insertions(+), 45 deletions(-) diff --git a/src/hssm/param/regression_param.py b/src/hssm/param/regression_param.py index 146936d33..6de315f67 100644 --- a/src/hssm/param/regression_param.py +++ b/src/hssm/param/regression_param.py @@ -12,7 +12,12 @@ from formulae.matrices import DesignMatrices from ..link import Link -from ..prior import _is_identity_link, get_default_prior, get_hddm_default_prior +from ..prior import ( + _has_finite_bounds, + _is_identity_link, + get_default_prior, + get_hddm_default_prior, +) from .param import Param from .parameterization import NoncenteredSetting, _resolve_noncentered from .user_param import UserParam @@ -254,6 +259,7 @@ def _make_safe_priors( """ safe_priors = {} generated_unmatched_terms: list[tuple[str, str]] = [] + bounded_group_locations: list[str] = [] get_prior = get_hddm_default_prior if is_ddm else get_default_prior specified_priors = ( @@ -300,34 +306,24 @@ def _make_safe_priors( link=self.link, ) else: - # treat the term as any other group-specific term - if ( - _is_identity_link(self.link) - and self.bounds is not None - and any(np.isfinite(bound) for bound in self.bounds) - ): - _logger.warning( - "The generated group-only intercept for parameter " - "%s is on the response/parameter scale under the " - "identity link, but its coefficient prior does not " - "apply finite HSSM bounds %s due to a current " - "Bambi limitation. Likelihood-level parameter " - "bounds still apply. A support-respecting " - "transformed link instead uses an unconstrained " - "predictor scale; bound-aware identity group " - "priors are " - "tracked in HSSM #1269.", - self.name, - self.bounds, - ) + # Generic identity-scale group locations can use a native + # bounded hierarchy. HDDM-derived families retain their + # calibrated response-scale hierarchies. + group_bounds = None if is_ddm else self.bounds prior = get_prior( "group_intercept", self.name, - bounds=None, + bounds=group_bounds, link=self.link, ) prior.noncentered = False safe_priors[name] = prior + if ( + not is_ddm + and _is_identity_link(self.link) + and _has_finite_bounds(self.bounds) + ): + bounded_group_locations.append(name) generated_unmatched_terms.append( (name, self._group_term_names[name]) ) @@ -349,6 +345,10 @@ def _make_safe_priors( (name, self._group_term_names[name]) ) + self._warn_if_bounded_group_location_is_not_complete_predictor( + dm, bounded_group_locations + ) + if generated_unmatched_terms and _resolve_noncentered(noncentered, self.name): generated_unmatched_terms.sort(key=lambda item: (item[1], item[0])) term_details = [ @@ -384,6 +384,41 @@ def _make_safe_priors( safe_priors.update(self.prior) self.prior = safe_priors + def _warn_if_bounded_group_location_is_not_complete_predictor( + self, + dm: DesignMatrices, + bounded_group_locations: list[str], + ) -> None: + """Warn when a bounded group location is only one additive contribution.""" + if not bounded_group_locations: + return + + common_terms = sorted(dm.common.terms) if dm.common is not None else [] + other_group_terms = ( + sorted( + name for name in dm.group.terms if name not in bounded_group_locations + ) + if dm.group is not None + else [] + ) + additional_terms = common_terms + other_group_terms + if not additional_terms: + return + + _logger.warning( + "The generated identity-scale group-location prior(s) %r for parameter " + "%r constrain those coefficients to HSSM bounds %s, but the formula " + "also contains additive term(s) %r. These coefficient bounds do not " + "constrain the complete linear predictor or final parameter value. " + "Out-of-bounds final values receive HSSM's finite likelihood-floor " + "penalty. Inspect the complete predictor or choose an inverse link whose " + "image matches the required parameter support.", + sorted(bounded_group_locations), + self.name, + self.bounds, + additional_terms, + ) + def _validate_generated_group_locations(self) -> None: """Reject ambiguous safe defaults for repeated group-only expressions. diff --git a/tests/unit/param/test_params.py b/tests/unit/param/test_params.py index 35a72d631..4df26f309 100644 --- a/tests/unit/param/test_params.py +++ b/tests/unit/param/test_params.py @@ -1,6 +1,7 @@ from unittest.mock import Mock import bambi as bmb +import numpy as np import pytest from hssm import HSSM, Link, Prior @@ -10,9 +11,9 @@ from hssm.param.params import ( Params, collect_user_params, - make_params, - make_param_from_user_param, make_param_from_defaults, + make_param_from_user_param, + make_params, ) from hssm.param.regression_param import RegressionParam from hssm.param.simple_param import DefaultParam, SimpleParam @@ -438,6 +439,29 @@ def counted_get_design_matrices(self, data, extra_namespace): assert params["t"]._group_term_names == {} +def test_approximate_ddm_uses_generic_bounded_group_location(data_ddm_reg): + """Use the neural training box rather than HDDM families for safe priors.""" + model = create_mock_model( + "ddm", loglik_kind="approx_differentiable", prior_settings="safe" + ) + model.list_params = ["v"] + model.data = data_ddm_reg.assign(participant_id=np.arange(len(data_ddm_reg)) % 2) + user_params = { + "v": UserParam( + name="v", + formula="v ~ 0 + (1 | participant_id)", + ) + } + + params = make_params(model, user_params, noncentered=True) + + prior = params["v"].prior["1|participant_id"] + assert prior.name == "TruncatedNormal" + assert prior.args["lower"] == -3.0 + assert prior.args["upper"] == 3.0 + assert prior.noncentered is False + + def test_make_params_prepares_formula_without_safe_priors(data_ddm_reg): """Prepare structural metadata even when safe priors are disabled.""" model = create_mock_model("ddm", global_formula="t ~ x + (0 + y | x)") diff --git a/tests/unit/param/test_regression_param.py b/tests/unit/param/test_regression_param.py index bfa23f49a..06ea640cc 100644 --- a/tests/unit/param/test_regression_param.py +++ b/tests/unit/param/test_regression_param.py @@ -908,16 +908,16 @@ def test_hddm_safe_group_only_intercept_uses_preset_identity(cavanagh_test): @pytest.mark.parametrize( - ("link", "expect_bounds_warning"), + ("link", "expected_prior"), [ - pytest.param(None, True, id="omitted-identity"), - pytest.param("identity", True, id="string-identity"), - pytest.param(bmb.Link("identity"), True, id="bambi-identity"), - pytest.param(hssm.Link("identity"), True, id="hssm-identity"), - pytest.param("log", False, id="log"), + pytest.param(None, "TruncatedNormal", id="omitted-identity"), + pytest.param("identity", "TruncatedNormal", id="string-identity"), + pytest.param(bmb.Link("identity"), "TruncatedNormal", id="bambi-identity"), + pytest.param(hssm.Link("identity"), "TruncatedNormal", id="hssm-identity"), + pytest.param("log", "Normal", id="log"), pytest.param( hssm.Link("gen_logit", bounds=(0.0, 1.0)), - False, + "Normal", id="gen-logit", ), pytest.param( @@ -927,15 +927,15 @@ def test_hddm_safe_group_only_intercept_uses_preset_identity(cavanagh_test): linkinv=np.exp, linkinv_backend=pt.exp, ), - False, + "Normal", id="custom", ), ], ) -def test_group_only_bounds_warning_is_identity_specific( - cavanagh_test, caplog, link, expect_bounds_warning +def test_group_only_bounds_follow_effective_link( + cavanagh_test, caplog, link, expected_prior ): - """Do not warn about response bounds on a transformed predictor scale.""" + """Bound only pure identity-scale group locations without a residual warning.""" param = RegressionParam( name="v", formula="v ~ 0 + (1 | participant_id)", @@ -945,13 +945,72 @@ def test_group_only_bounds_warning_is_identity_specific( param.make_safe_priors(cavanagh_test, {}, is_ddm=False, noncentered=False) + prior = param.prior["1|participant_id"] + assert prior.name == expected_prior + assert not any( + "complete linear predictor" in record.message for record in caplog.records + ) + + +@pytest.mark.parametrize( + ("formula", "additional_term"), + [ + pytest.param( + "v ~ 0 + theta + (1 | participant_id)", "theta", id="common-slope" + ), + pytest.param( + "v ~ 0 + (1 + theta | participant_id)", + "theta|participant_id", + id="group-slope", + ), + pytest.param( + "v ~ 0 + hsgp(theta, m=8, c=2) + (1 | participant_id)", + "hsgp(theta, m=8, c=2)", + id="hsgp", + ), + ], +) +def test_bounded_group_location_warns_about_additive_predictor( + cavanagh_test, caplog, formula, additional_term +): + """Distinguish a bounded coefficient from the complete identity predictor.""" + param = RegressionParam( + name="v", + formula=formula, + bounds=(0.0, 1.0), + link="identity", + ) + + param.make_safe_priors(cavanagh_test, {}, is_ddm=False, noncentered=False) + messages = [ - record.message for record in caplog.records if "HSSM #1269" in record.message + record.message + for record in caplog.records + if "complete linear predictor" in record.message ] - assert bool(messages) is expect_bounds_warning - if expect_bounds_warning: - assert len(messages) == 1 - assert "Likelihood-level parameter bounds still apply" in messages[0] + assert len(messages) == 1 + assert "1|participant_id" in messages[0] + assert additional_term in messages[0] + assert "finite likelihood-floor penalty" in messages[0] + + +def test_transformed_mixed_predictor_keeps_unbounded_group_coefficient( + cavanagh_test, caplog +): + """Leave additive transformed-link coefficients on the predictor scale.""" + param = RegressionParam( + name="v", + formula="v ~ 0 + theta + (1 | participant_id)", + bounds=(0.0, np.inf), + link="log", + ) + + param.make_safe_priors(cavanagh_test, {}, is_ddm=False, noncentered=False) + + assert param.prior["1|participant_id"].name == "Normal" + assert not any( + "complete linear predictor" in record.message for record in caplog.records + ) @pytest.mark.parametrize( @@ -968,9 +1027,12 @@ def test_group_only_bounds_warning_requires_finite_bounds( link="identity", ) - param.make_safe_priors(cavanagh_test, {}, is_ddm=True, noncentered=False) + param.make_safe_priors(cavanagh_test, {}, is_ddm=False, noncentered=False) - assert not any("HSSM #1269" in record.message for record in caplog.records) + assert param.prior["1|participant_id"].name == "Normal" + assert not any( + "complete linear predictor" in record.message for record in caplog.records + ) @pytest.mark.parametrize( @@ -1086,12 +1148,12 @@ def test_make_safe_priors(cavanagh_test, caplog, param_name, bounds, is_ddm): param_no_common_intercept.make_safe_priors(cavanagh_test, {}, is_ddm=False) - assert any("limitation" in record.msg for record in caplog.records) + assert any("complete linear predictor" in record.msg for record in caplog.records) assert "Intercept" not in param_no_common_intercept.prior group_intercept_prior = param_no_common_intercept.prior["1|participant_id"] group_slope_prior = param_no_common_intercept.prior["theta|participant_id"] - _check_group_prior(group_intercept_prior) + _check_bounded_group_prior(group_intercept_prior, bounds) _check_group_prior(group_slope_prior) # Change back after testing @@ -1118,6 +1180,27 @@ def _check_group_prior(group_prior): assert sigma.args["beta"] == 0.3 +def _check_bounded_group_prior(group_prior, bounds): + assert isinstance(group_prior, bmb.Prior) + assert group_prior.dist is None + assert group_prior.name == "TruncatedNormal" + assert group_prior.noncentered is False + assert group_prior.args["lower"] == bounds[0] + assert group_prior.args["upper"] == bounds[1] + + mu = group_prior.args["mu"] + assert isinstance(mu, bmb.Prior) + assert mu.name == "TruncatedNormal" + assert mu.args["lower"] == bounds[0] + assert mu.args["upper"] == bounds[1] + + sigma = group_prior.args["sigma"] + assert isinstance(sigma, bmb.Prior) + assert sigma.name == "Weibull" + assert sigma.args["alpha"] == 1.5 + assert sigma.args["beta"] == 0.3 + + def _check_group_prior_with_common(group_prior): assert isinstance(group_prior, bmb.Prior) assert group_prior.dist is None @@ -1267,7 +1350,9 @@ def _check_group_prior_intercept_ddm(group_prior, prior): ) param_no_common_intercept.make_safe_priors(cavanagh_test, {}, is_ddm=True) - assert any("limitation" in record.msg for record in caplog.records) + assert not any( + "complete linear predictor" in record.msg for record in caplog.records + ) assert "Intercept" not in param_no_common_intercept.prior group_intercept_prior = param_no_common_intercept.prior["1|participant_id"] diff --git a/tests/unit/param/test_unmatched_group_prior_graph.py b/tests/unit/param/test_unmatched_group_prior_graph.py index 43f83cf48..61b3af3f6 100644 --- a/tests/unit/param/test_unmatched_group_prior_graph.py +++ b/tests/unit/param/test_unmatched_group_prior_graph.py @@ -3,12 +3,14 @@ import bambi as bmb import numpy as np import pandas as pd +import pymc as pm import pytensor.tensor as pt import pytest from pytensor.graph.traversal import ancestors import hssm from hssm.param.parameterization_check import find_disconnected_free_rvs +from hssm.prior import get_default_prior def _group_only_data() -> pd.DataFrame: @@ -45,6 +47,25 @@ def _build_ddm(parameter: str, formula: str, noncentered=True) -> hssm.HSSM: ) +def _build_bounded_lba(formula: str = "b ~ 0 + (1 | participant_id)") -> hssm.HSSM: + """Build a lower-bounded generic analytical model without sampling.""" + data = _group_only_data().copy() + data["response"] = (data["response"] == 1).astype(int) + return hssm.HSSM( + data=data, + model="lba2", + include=[{"name": "b", "formula": formula}], + A=0.1, + v0=1.0, + v1=1.0, + p_outlier=0.0, + prior_settings="safe", + noncentered=True, + process_initvals=False, + initval_jitter=0.0, + ) + + def _build_group_only_intercept( parameter: str, *, @@ -172,6 +193,83 @@ def test_non_normal_group_only_intercept_builds_centered(): assert find_disconnected_free_rvs(model.pymc_model) == [] +@pytest.mark.parametrize( + "bounds", + [ + pytest.param((0.1, 0.9), id="finite"), + pytest.param((0.1, np.inf), id="lower-bounded"), + pytest.param((-np.inf, 0.9), id="upper-bounded"), + ], +) +def test_native_bounded_group_hierarchy_builds_and_draws_within_support(bounds): + """Exercise the generated hierarchy through Bambi and PyMC directly.""" + data = pd.DataFrame( + { + "y": np.linspace(-1.0, 1.0, 8), + "participant_id": np.repeat(["p0", "p1"], 4), + } + ) + prior = get_default_prior("group_intercept", "x", bounds, "identity") + prior.noncentered = False + model = bmb.Model( + "y ~ 0 + (1 | participant_id)", + data, + priors={"1|participant_id": prior}, + noncentered=True, + ) + model.build() + + pymc_model = model.backend.model + group_name = "1|participant_id" + group_rv = pymc_model.named_vars[group_name] + free_names = {rv.name for rv in pymc_model.free_RVs} + assert type(group_rv.owner.op).__name__ == "TruncatedNormalRV" + assert {group_name, f"{group_name}_mu", f"{group_name}_sigma"} <= free_names + assert f"{group_name}_offset" not in free_names + assert find_disconnected_free_rvs(pymc_model) == [] + + draws = pm.draw(group_rv, draws=64, random_seed=1269, backend="FAST_COMPILE") + lower, upper = bounds + if np.isfinite(lower): + assert np.all(draws > lower) + if np.isfinite(upper): + assert np.all(draws < upper) + + +def test_generated_bounded_group_location_reaches_hssm_graph_and_parameter(): + """Keep a pure identity predictor and its generated location inside support.""" + model = _build_bounded_lba() + + prior = model.params["b"].prior["1|participant_id"] + assert prior.name == "TruncatedNormal" + assert prior.noncentered is False + assert prior.args["lower"] == 0.2 + assert np.isposinf(prior.args["upper"]) + assert prior.args["mu"].name == "TruncatedNormal" + + group_name = "b_1|participant_id" + group_rv = model.pymc_model.named_vars[group_name] + parameter = model.pymc_model.named_vars["b"] + free_names = {rv.name for rv in model.pymc_model.free_RVs} + assert type(group_rv.owner.op).__name__ == "TruncatedNormalRV" + assert { + group_name, + f"{group_name}_mu", + f"{group_name}_sigma", + } <= free_names + assert f"{group_name}_offset" not in free_names + assert find_disconnected_free_rvs(model.pymc_model) == [] + + group_draws, parameter_draws = pm.draw( + [group_rv, parameter], + draws=64, + random_seed=1269, + backend="FAST_COMPILE", + ) + assert np.all(group_draws > 0.2) + assert np.all(parameter_draws > 0.2) + + @pytest.mark.parametrize( ("parameter", "outer_family", "hyperparameters"), [ From f94b9e467d43386bf1e822a9e66d0fb7a9d68055 Mon Sep 17 00:00:00 2001 From: Alexander Fengler Date: Sun, 30 Aug 2026 19:29:25 -0400 Subject: [PATCH 3/4] fix: keep bounded group initvals inside support (#1269) --- src/hssm/base.py | 85 ++++++++++- tests/test_initval_jitter_bounds.py | 211 ++++++++++++++++++++++++++++ 2 files changed, 293 insertions(+), 3 deletions(-) create mode 100644 tests/test_initval_jitter_bounds.py diff --git a/src/hssm/base.py b/src/hssm/base.py index d25488b46..045bfd12e 100644 --- a/src/hssm/base.py +++ b/src/hssm/base.py @@ -2221,14 +2221,93 @@ def __jitter_initvals_vector_only(self, jitter_epsilon: float) -> None: for name_, starting_value in initial_point_dict.items(): name_tmp = name_.replace("_log__", "").replace("_interval__", "") if starting_value.ndim != 0 and starting_value.shape[0] != 1: - starting_value_tmp = starting_value + np.random.uniform( - -jitter_epsilon, jitter_epsilon, starting_value.shape - ).astype(np.float32) + bounds = self._get_group_initval_bounds(name_tmp) + if bounds is None: + starting_value_tmp = starting_value + np.random.uniform( + -jitter_epsilon, jitter_epsilon, starting_value.shape + ).astype(np.float32) + else: + starting_value_tmp = self._jitter_within_bounds( + starting_value, jitter_epsilon, bounds + ) # Note: self._initvals shouldn't be None when this is called dtype = self._initvals[name_tmp].dtype self._initvals[name_tmp] = np.array(starting_value_tmp).astype(dtype) + def _get_group_initval_bounds(self, name: str) -> tuple[float, float] | None: + """Return native ``TruncatedNormal`` bounds for a group term.""" + parameter_name = self._get_prefix(name) + parameter = self.params.get(parameter_name) + if parameter is None or not isinstance(parameter.prior, dict): + return None + + prefix = f"{parameter_name}_" + if not name.startswith(prefix): + return None + term_name = name.removeprefix(prefix) + if term_name not in getattr(parameter, "_group_term_names", {}): + return None + + prior = parameter.prior.get(term_name, parameter.prior.get("group_specific")) + if not isinstance(prior, bmb.Prior) or prior.name != "TruncatedNormal": + return None + lower_arg = prior.args.get("lower", -np.inf) + upper_arg = prior.args.get("upper", np.inf) + lower = np.asarray(-np.inf if lower_arg is None else lower_arg) + upper = np.asarray(np.inf if upper_arg is None else upper_arg) + if ( + lower.ndim != 0 + or upper.ndim != 0 + or not np.issubdtype(lower.dtype, np.number) + or not np.issubdtype(upper.dtype, np.number) + ): + return None + return float(lower), float(upper) + + @staticmethod + def _jitter_within_bounds( + starting_value: np.ndarray, + jitter_epsilon: float, + bounds: tuple[float, float], + ) -> np.ndarray: + """Jitter a constrained vector without crossing an open interval.""" + lower, upper = bounds + dtype = ( + starting_value.dtype + if np.issubdtype(starting_value.dtype, np.floating) + else np.dtype(np.float64) + ) + positive_inf = np.asarray(np.inf, dtype=dtype) + negative_inf = np.asarray(-np.inf, dtype=dtype) + lower_array = np.asarray(lower, dtype=dtype) + upper_array = np.asarray(upper, dtype=dtype) + lower_limit = ( + np.nextafter(lower_array, positive_inf).item() + if np.isfinite(lower) + else -np.inf + ) + upper_limit = ( + np.nextafter(upper_array, negative_inf).item() + if np.isfinite(upper) + else np.inf + ) + if ( + not np.all(np.isfinite(starting_value)) + or np.any(starting_value < lower_limit) + or np.any(starting_value > upper_limit) + ): + # Preserve an already-invalid user value so PyMC can report it instead of + # silently changing the user's specification. + return starting_value.copy() + + jitter_lower = np.maximum(-jitter_epsilon, lower_limit - starting_value) + jitter_upper = np.minimum(jitter_epsilon, upper_limit - starting_value) + jittered = starting_value + np.random.uniform( + jitter_lower, jitter_upper, starting_value.shape + ).astype(dtype) + return np.clip(jittered, lower_limit, upper_limit).astype(dtype, copy=False) + def __jitter_initvals_all(self, jitter_epsilon: float) -> None: # Note: Calling our initial point function here # --> operate on untransformed variables diff --git a/tests/test_initval_jitter_bounds.py b/tests/test_initval_jitter_bounds.py new file mode 100644 index 000000000..4a9620bab --- /dev/null +++ b/tests/test_initval_jitter_bounds.py @@ -0,0 +1,211 @@ +"""Regression tests for support-aware hierarchical initial-value jitter.""" + +import numpy as np +import pytest +from pymc.exceptions import SamplingError +from pymc.initial_point import StartDict, make_initial_point_fns_per_chain + +import hssm +from hssm.defaults import INITVAL_JITTER_SETTINGS +from hssm.likelihoods import logp_ddm + +hssm.set_floatX("float32", update_jax=True) + +GROUP_TERM = "v_1|participant_id" +NARROW_BOUNDS = (0.499, 0.501) + + +def _build_group_model(cavanagh_test, bounds, **kwargs) -> hssm.HSSM: + """Build a tiny custom DDM with one generated group-only location.""" + data = cavanagh_test.groupby("participant_id").head(2).copy() + return hssm.HSSM( + data=data, + model="custom", + model_config={ + "list_params": ["v", "a", "z", "t"], + "choices": [-1, 1], + "bounds": { + "v": bounds, + "a": (0.1, np.inf), + "z": (0.0, 1.0), + "t": (0.0, np.inf), + }, + }, + loglik=logp_ddm, + loglik_kind="analytical", + include=[{"name": "v", "formula": "v ~ 0 + (1 | participant_id)"}], + a=1.5, + z=0.5, + t=0.2, + p_outlier=0.0, + prior_settings="safe", + **kwargs, + ) + + +def _sampler_initial_point(model: hssm.HSSM) -> dict[str, np.ndarray]: + """Compile HSSM's constrained overrides through PyMC's sampler path.""" + overrides: StartDict = { + name: np.asarray(value) for name, value in model._initvals.items() + } + initial_point_fn = make_initial_point_fns_per_chain( + model=model.pymc_model, + overrides=overrides, + jitter_rvs=set(), + chains=1, + )[0] + return initial_point_fn(1269) + + +def _group_value_name(model: hssm.HSSM) -> str: + """Return the transformed PyMC value name for the group random variable.""" + group_rv = model.pymc_model.named_vars[GROUP_TERM] + value_name = model.pymc_model.rvs_to_values[group_rv].name + assert value_name is not None + return value_name + + +def _uniform_endpoint(endpoint: str): + """Return a uniform stub that deterministically selects one endpoint.""" + + def select(low, high, size=None): + selected = low if endpoint == "low" else high + return np.broadcast_to(np.asarray(selected), size).copy() + + return select + + +def test_default_jitter_keeps_narrow_group_initvals_strictly_inside_support( + cavanagh_test, monkeypatch +): + """Default vector jitter cannot cross narrow generated group bounds.""" + monkeypatch.setattr(np.random, "uniform", _uniform_endpoint("high")) + + model = _build_group_model(cavanagh_test, NARROW_BOUNDS) + jittered = model._initvals[GROUP_TERM] + + assert model.initval_jitter == INITVAL_JITTER_SETTINGS["jitter_epsilon"] + assert model.params["v"].prior["1|participant_id"].name == "TruncatedNormal" + assert model._get_group_initval_bounds(GROUP_TERM) == pytest.approx(NARROW_BOUNDS) + assert jittered.dtype == np.dtype("float32") + assert np.all(jittered > NARROW_BOUNDS[0]) + assert np.all(jittered < NARROW_BOUNDS[1]) + assert np.any(jittered != np.float32(0.5)) + + sampler_point = _sampler_initial_point(model) + assert all(np.all(np.isfinite(value)) for value in sampler_point.values()) + + +@pytest.mark.parametrize( + ("bounds", "endpoint"), + [ + pytest.param((0.2, np.inf), "low", id="lower-only"), + pytest.param((-np.inf, 0.2), "high", id="upper-only"), + ], +) +def test_group_jitter_respects_one_sided_support( + cavanagh_test, monkeypatch, bounds, endpoint +): + """One-sided native group priors retain a strict finite boundary.""" + model = _build_group_model(cavanagh_test, bounds, initval_jitter=0.0) + prior = model.params["v"].prior["1|participant_id"] + missing_endpoint = "upper" if np.isfinite(bounds[0]) else "lower" + prior.args[missing_endpoint] = None + detected_bounds = model._get_group_initval_bounds(GROUP_TERM) + assert detected_bounds == pytest.approx(bounds) + + dtype = model._initvals[GROUP_TERM].dtype + if np.isfinite(bounds[0]): + boundary = np.asarray(bounds[0], dtype=dtype) + interior = np.nextafter( + np.nextafter(boundary, np.asarray(np.inf, dtype=dtype)), + np.asarray(np.inf, dtype=dtype), + ).item() + else: + boundary = np.asarray(bounds[1], dtype=dtype) + interior = np.nextafter( + np.nextafter(boundary, np.asarray(-np.inf, dtype=dtype)), + np.asarray(-np.inf, dtype=dtype), + ).item() + model._initvals[GROUP_TERM] = np.full_like(model._initvals[GROUP_TERM], interior) + monkeypatch.setattr(np.random, "uniform", _uniform_endpoint(endpoint)) + + model._jitter_initvals(vector_only=True) + jittered = model._initvals[GROUP_TERM] + + if np.isfinite(bounds[0]): + assert np.all(jittered > bounds[0]) + if np.isfinite(bounds[1]): + assert np.all(jittered < bounds[1]) + sampler_point = _sampler_initial_point(model) + assert np.all(np.isfinite(sampler_point[_group_value_name(model)])) + + +def test_zero_jitter_preserves_group_initial_point(cavanagh_test): + """An explicit zero jitter leaves the generated group vector unchanged.""" + model = _build_group_model(cavanagh_test, NARROW_BOUNDS, initval_jitter=0.0) + + np.testing.assert_array_equal( + model._initvals[GROUP_TERM], model.initial_point()[GROUP_TERM] + ) + sampler_point = _sampler_initial_point(model) + assert np.all(np.isfinite(sampler_point[_group_value_name(model)])) + + +def test_invalid_group_initvals_are_not_silently_repaired(cavanagh_test, monkeypatch): + """Boundary, outside, and NaN starts remain invalid for PyMC to report.""" + model = _build_group_model(cavanagh_test, NARROW_BOUNDS, initval_jitter=0.0) + lower, upper = model._get_group_initval_bounds(GROUP_TERM) + dtype = model._initvals[GROUP_TERM].dtype + lower_value = np.asarray(lower, dtype=dtype) + upper_value = np.asarray(upper, dtype=dtype) + invalid_values = [ + lower_value.item(), + np.nextafter(lower_value, np.asarray(-np.inf, dtype=dtype)).item(), + upper_value.item(), + np.nextafter(upper_value, np.asarray(np.inf, dtype=dtype)).item(), + np.nan, + ] + + def unexpected_uniform(*args, **kwargs): + raise AssertionError("invalid initial values must not be jittered") + + monkeypatch.setattr(np.random, "uniform", unexpected_uniform) + for invalid_value in invalid_values: + supplied = np.full_like(model._initvals[GROUP_TERM], invalid_value) + model._initvals[GROUP_TERM] = supplied.copy() + + model._jitter_initvals(vector_only=True) + + np.testing.assert_equal(model._initvals[GROUP_TERM], supplied) + sampler_point = _sampler_initial_point(model) + assert not np.all(np.isfinite(sampler_point[_group_value_name(model)])) + with pytest.raises(SamplingError): + model.pymc_model.check_start_vals(sampler_point, mode="FAST_COMPILE") + + +def test_unbounded_vector_jitter_matches_legacy_seeded_result(cavanagh_test): + """Unbounded group vectors retain the prior seeded additive-jitter behavior.""" + model = _build_group_model(cavanagh_test, (-np.inf, np.inf), initval_jitter=0.0) + assert model.params["v"].prior["1|participant_id"].name == "Normal" + assert model._get_group_initval_bounds(GROUP_TERM) is None + + starting_value = model._initvals[GROUP_TERM].copy() + jitter_epsilon = INITVAL_JITTER_SETTINGS["jitter_epsilon"] + random_state = np.random.get_state() + try: + np.random.seed(1269) + expected = starting_value + np.random.uniform( + -jitter_epsilon, jitter_epsilon, starting_value.shape + ).astype(np.float32) + np.random.seed(1269) + model._jitter_initvals( + jitter_epsilon=jitter_epsilon, + vector_only=True, + ) + finally: + np.random.set_state(random_state) + + np.testing.assert_array_equal(model._initvals[GROUP_TERM], expected) + sampler_point = _sampler_initial_point(model) + assert np.all(np.isfinite(sampler_point[_group_value_name(model)])) From 9e7002f00e7ff40524a8afc9e83b62723461c98f Mon Sep 17 00:00:00 2001 From: Alexander Fengler Date: Sun, 30 Aug 2026 19:45:32 -0400 Subject: [PATCH 4/4] docs: explain bounded group locations (#1269) --- docs/changelog.md | 4 +- docs/explanations/coming_from_hddm.md | 26 +- docs/how_to/specify_group_priors.md | 27 +- docs/tutorials/link_functions.ipynb | 345 +++++++++++++++++++++----- docs/tutorials/link_functions.py | 132 ++++++++-- 5 files changed, 440 insertions(+), 94 deletions(-) diff --git a/docs/changelog.md b/docs/changelog.md index d02fbaafc..ebf5de1d6 100644 --- a/docs/changelog.md +++ b/docs/changelog.md @@ -16,7 +16,9 @@ 7. **Invalid model-level prior and link presets now fail fast** (#1233). `prior_settings` accepts only `"safe"` or `None`, and `link_settings` accepts only `"log_logit"` or `None`; wrong-case strings, booleans, mappings, and other unsupported values now raise a clear `ValueError` instead of silently changing prior, link, initial-value, or display behavior. The API documentation now clarifies that both presets act on regression parameters: `prior_settings=None` delegates missing regression-term priors to Bambi but leaves HSSM's simple-parameter defaults unchanged. The Poisson-race tutorial has also been migrated to marimo and no longer presents `prior_settings` as a prior dictionary. -8. **Safe priors now preserve a unique unmatched group-only population location and reject ambiguous or unrepresentable specifications** ([#1225](https://github.com/lnccbrown/HSSM/issues/1225)). When HSSM generates a prior for a group term with no exact common Formulae counterpart, the group distribution owns the population location. HSSM therefore preserves its location-bearing hierarchy and centers that generated term even when the model or component requested non-centering, because Bambi's current non-centered construction omits `mu`. For previously generated unmatched Normal terms under default non-centering, this intentionally changes the likelihood by reconnecting the location that Bambi had discarded. Exact common/group matches remain zero-mean deviations and honor the requested parameterization ([#1224](https://github.com/lnccbrown/HSSM/issues/1224)). Repeated unmatched expressions now fail when safe generation has no unique owner. Explicit priors are never rewritten, but HSSM raises before Bambi when a group prior cannot be represented faithfully and warns about buildable centered specifications with multiple free owners of one location. Identity-linked group-only intercepts retain their response-scale hierarchy for every identity spelling ([#1232](https://github.com/lnccbrown/HSSM/issues/1232)); transformed links use a hierarchy on the linear-predictor scale. Follow-ups cover upstream location-aware non-centering ([#1268](https://github.com/lnccbrown/HSSM/issues/1268), [Bambi #1003](https://github.com/bambinos/bambi/issues/1003)), bounded identity group priors ([#1269](https://github.com/lnccbrown/HSSM/issues/1269)), and broader numeric regression-term semantics ([#1271](https://github.com/lnccbrown/HSSM/issues/1271)). +8. **Safe priors now preserve a unique unmatched group-only population location and reject ambiguous or unrepresentable specifications** ([#1225](https://github.com/lnccbrown/HSSM/issues/1225)). When HSSM generates a prior for a group term with no exact common Formulae counterpart, the group distribution owns the population location. HSSM therefore preserves its location-bearing hierarchy and centers that generated term even when the model or component requested non-centering, because Bambi's current non-centered construction omits `mu`. For previously generated unmatched Normal terms under default non-centering, this intentionally changes the likelihood by reconnecting the location that Bambi had discarded. Exact common/group matches remain zero-mean deviations and honor the requested parameterization ([#1224](https://github.com/lnccbrown/HSSM/issues/1224)). Repeated unmatched expressions now fail when safe generation has no unique owner. Explicit priors are never rewritten, but HSSM raises before Bambi when a group prior cannot be represented faithfully and warns about buildable centered specifications with multiple free owners of one location. Identity-linked group-only intercepts retain their response-scale hierarchy for every identity spelling ([#1232](https://github.com/lnccbrown/HSSM/issues/1232)); transformed links use a hierarchy on the linear-predictor scale. Follow-ups cover upstream location-aware non-centering ([#1268](https://github.com/lnccbrown/HSSM/issues/1268), [Bambi #1003](https://github.com/bambinos/bambi/issues/1003)) and broader numeric regression-term semantics ([#1271](https://github.com/lnccbrown/HSSM/issues/1271)). + +9. **Generated identity-linked group locations now honor configured response bounds** ([#1269](https://github.com/lnccbrown/HSSM/issues/1269)). A unique unmatched group-only intercept on a generic parameter receives a centered native hierarchical `TruncatedNormal` whenever at least one response bound is finite. The generated population-location hyperprior and group coefficients retain visible Bambi hierarchy and stay within the configured support; analytical and black-box DDM parameters keep their calibrated HDDM prior families, transformed-link coefficients remain unbounded on the predictor scale, and explicit priors are unchanged. HSSM warns when other additive terms mean that one bounded coefficient cannot constrain the full identity predictor. Vector initial-value jitter is now support-aware for native bounded group terms, including narrow and one-sided intervals, while invalid or boundary starts remain available for PyMC's normal validation. ### 0.4.0 diff --git a/docs/explanations/coming_from_hddm.md b/docs/explanations/coming_from_hddm.md index af4202cf3..c37622496 100644 --- a/docs/explanations/coming_from_hddm.md +++ b/docs/explanations/coming_from_hddm.md @@ -136,17 +136,21 @@ several free centered owners produce a likelihood ridge. [Link functions and safe priors](../tutorials/link_functions.ipynb) explains these scale and ownership changes from first principles. Explicit priors are never rewritten, but incompatible group specifications are rejected before - Bambi can drop or misinterpret them. Finite coefficient bounds are not yet - propagated to generic identity-linked group-only intercept priors; likelihood - bounds still apply, and bounding one coefficient would not constrain a complete - identity-linked predictor after other effects are added. Prefer a - support-respecting transformed link when appropriate; [Specify hierarchical - group priors](../how_to/specify_group_priors.md) explains the remaining - identity-link choices and compatibility rules. The broader rule is the - likelihood: these defaults apply unless you use the neural - (`approx_differentiable`) likelihood, which has its own priors derived from - the network's training bounds. Specifying your own is a different interface - through the same prior controls. + Bambi can drop or misinterpret them. For a generic identity-linked parameter, + a unique generated group-only intercept with a finite configured bound now + receives a centered native `TruncatedNormal` hierarchy. A pure group-intercept + predictor is therefore supported by construction. The bound still applies to + that coefficient, not to a complete additive predictor after slopes or other + effects are added; HSSM warns about that distinction. A transformed link can + instead constrain the assembled predictor when its inverse image matches the + parameter support. Values outside configured likelihood bounds receive HSSM's + finite per-trial log-likelihood floor (`-66.1`), not a hard-support rejection. + [Specify hierarchical group priors](../how_to/specify_group_priors.md) gives + the exact generated and explicit-prior rules. Analytical and black-box DDM + variants keep their calibrated HDDM Gamma, Beta, or Normal hierarchies, while + neural (`approx_differentiable`) variants use generic safe priors derived from + the network's finite training bounds. Specifying your own is a different + interface through the same prior controls. - **Outliers.** HDDM's `p_outlier` exists in HSSM under the same name, and the lapse distribution is configurable rather than fixed. See [Model outliers with lapse probabilities](../tutorials/lapse_prob_and_dist.ipynb). diff --git a/docs/how_to/specify_group_priors.md b/docs/how_to/specify_group_priors.md index 722379769..7d605bcf1 100644 --- a/docs/how_to/specify_group_priors.md +++ b/docs/how_to/specify_group_priors.md @@ -68,6 +68,21 @@ The per-prior `noncentered=False` override takes precedence over a component or model-level `noncentered=True` setting. Generated safe priors apply this centered fallback automatically for a unique group-only location. +When that generated owner is an **intercept on an identity link**, HSSM also +uses the parameter's configured bounds when it can do so without changing an +HDDM-calibrated prior family. A generic parameter with at least one finite bound +receives a native hierarchical `TruncatedNormal`: both its generated population +location and its group coefficients stay inside the configured interval, and +the group term is centered so Bambi retains that hierarchy. A pure formula such +as `b ~ 0 + (1 | participant_id)` therefore keeps the complete predictor inside +the bounds. + +This generated-safe behavior is deliberately narrower than general constraint +propagation. It does not apply to slopes, matching zero-mean deviations, +transformed links, or explicit priors. Analytical and black-box DDM families +keep their calibrated response-scale Gamma, Beta, or Normal hierarchies when +those families already match the built-in parameter support. + If the same unmatched expression occurs under several grouping factors, do not give every group distribution a free location. Add the exact common expression and use zero-mean group deviations, or deliberately choose exactly one group term @@ -86,7 +101,7 @@ The relevant rules are: | Effectively non-centered plain `Normal` | Supported only with hierarchical `sigma`, absent or all-zero `mu`, no truncation or custom distribution, and no extra arguments | | Free or nonzero group `mu` | Use `noncentered=False` so the requested location is retained | | Hierarchical non-Normal or custom outer family | Use `noncentered=False` | -| `hssm.Prior(..., bounds=...)` on a group term | Rejected under either parameterization; HSSM's truncated wrapper cannot satisfy Bambi's group-hyperprior contract | +| Explicit `hssm.Prior(..., bounds=...)` on a group term | Rejected under either parameterization; HSSM's custom truncated wrapper cannot satisfy Bambi's group-hyperprior contract. This is distinct from the native named `TruncatedNormal` hierarchy HSSM generates for the bounded safe case above. | These checks do not rewrite explicit priors. HSSM raises when continuing would either fail in Bambi or silently construct a different prior tree. The same @@ -98,7 +113,8 @@ prior-valued arguments. All common and group coefficients first combine on the linear-predictor scale, and the inverse link is applied afterward. Bounding one identity-linked group intercept therefore does **not** guarantee that the full predictor remains inside -the parameter's support once slopes and other effects are added. +the parameter's support once slopes and other effects are added. HSSM warns when +it generates a bounded group-location coefficient in such a mixed predictor. Use a support-respecting transformed link when it matches the model: @@ -109,8 +125,11 @@ Use a support-respecting transformed link when it matches the model: If an identity link is scientifically required, choose a centered hierarchical family with appropriate natural support when possible, and remember that this -constrains that coefficient rather than the entire predictor. HSSM's likelihood -bounds still apply to the assembled parameter value. +constrains that coefficient rather than the entire predictor. For values outside +configured likelihood bounds, HSSM substitutes a finite per-trial log-likelihood +floor (`-66.1`); this is a penalty with a flat region, not a hard-support prior. +A transformed link is the mechanism that constrains the *complete* additive +predictor when its inverse image matches the parameter support. For the underlying scale and location logic, continue with [Link functions and safe priors](../tutorials/link_functions.ipynb). For the general prior interface, diff --git a/docs/tutorials/link_functions.ipynb b/docs/tutorials/link_functions.ipynb index 95bcb088a..31d946ae3 100644 --- a/docs/tutorials/link_functions.ipynb +++ b/docs/tutorials/link_functions.ipynb @@ -880,6 +880,28 @@ " },\n", "}\n", "group_location_models = {}\n", + "_lba_data = _group_base_kwargs[\"data\"].copy()\n", + "_lba_data[\"response\"] = (_lba_data[\"response\"] == 1).astype(int)\n", + "group_location_models[\"bounded generic identity\"] = build_silent_model(\n", + " data=_lba_data,\n", + " model=\"lba2\",\n", + " loglik_kind=\"analytical\",\n", + " include=[\n", + " {\n", + " \"name\": \"b\",\n", + " \"formula\": \"b ~ 0 + (1 | participant_id)\",\n", + " \"link\": \"identity\",\n", + " }\n", + " ],\n", + " A=0.1,\n", + " v0=1.0,\n", + " v1=1.0,\n", + " p_outlier=0.0,\n", + " prior_settings=\"safe\",\n", + " noncentered=True,\n", + " process_initvals=False,\n", + " initval_jitter=0.0,\n", + ")\n", "for _label, _case in _group_cases.items():\n", " _spec = {\n", " \"name\": _case[\"parameter\"],\n", @@ -936,6 +958,7 @@ " outer family\n", " location\n", " location scale\n", + " group-coefficient support\n", " effective form\n", " direct group RV\n", " offset RV\n", @@ -943,12 +966,25 @@ " \n", " \n", " \n", + " bounded generic identity\n", + " identity\n", + " 1|participant_id\n", + " TruncatedNormal\n", + " TruncatedNormal(mu: 0.0, sigma: 0.25, lower: 0.2, upper: inf)\n", + " bounded response scale\n", + " (0.2, inf)\n", + " centered location owner\n", + " True\n", + " False\n", + " \n", + " \n", " identity HDDM location\n", " identity\n", " 1|participant_id\n", " Gamma\n", " Gamma(mu: 1.5, sigma: 0.75)\n", " response scale\n", + " (0, inf) via Gamma\n", " centered location owner\n", " True\n", " False\n", @@ -960,6 +996,7 @@ " Normal\n", " Normal(mu: 0.0, sigma: 0.25)\n", " log-predictor scale\n", + " real line before exp\n", " centered location owner\n", " True\n", " False\n", @@ -971,6 +1008,7 @@ " Normal\n", " Normal(mu: 0.0, sigma: 0.25)\n", " generalized-log-odds scale\n", + " real line before gen_logit inverse\n", " centered location owner\n", " True\n", " False\n", @@ -982,6 +1020,7 @@ " Normal\n", " array(0.)\n", " identity predictor scale\n", + " real line\n", " non-centered zero-mean deviation\n", " False\n", " True\n", @@ -998,26 +1037,42 @@ "_group_rows = []\n", "_group_specs = [\n", " (\n", + " \"bounded generic identity\",\n", + " \"b\",\n", + " \"1|participant_id\",\n", + " \"bounded response scale\",\n", + " \"(0.2, inf)\",\n", + " ),\n", + " (\n", " \"identity HDDM location\",\n", " \"a\",\n", " \"1|participant_id\",\n", " \"response scale\",\n", + " \"(0, inf) via Gamma\",\n", + " ),\n", + " (\n", + " \"log-scale location\",\n", + " \"a\",\n", + " \"1|participant_id\",\n", + " \"log-predictor scale\",\n", + " \"real line before exp\",\n", " ),\n", - " (\"log-scale location\", \"a\", \"1|participant_id\", \"log-predictor scale\"),\n", " (\n", " \"generalized-logit location\",\n", " \"z\",\n", " \"1|participant_id\",\n", " \"generalized-log-odds scale\",\n", + " \"real line before gen_logit inverse\",\n", " ),\n", " (\n", " \"matched non-centered deviation\",\n", " \"v\",\n", " \"x|participant_id\",\n", " \"identity predictor scale\",\n", + " \"real line\",\n", " ),\n", "]\n", - "for _label, _parameter, _term, _scale in _group_specs:\n", + "for _label, _parameter, _term, _scale, _support in _group_specs:\n", " _model = group_location_models[_label]\n", " _prior = _model.params[_parameter].prior[_term]\n", " _mu = _prior.args.get(\"mu\")\n", @@ -1032,6 +1087,7 @@ " \"outer family\": _prior.name,\n", " \"location\": format_prior(_mu) if _mu_is_free else repr(_mu),\n", " \"location scale\": _scale,\n", + " \"group-coefficient support\": _support,\n", " \"effective form\": (\n", " \"centered location owner\"\n", " if _prior.noncentered is False\n", @@ -1048,6 +1104,11 @@ " .params[\"a\"]\n", " .prior[\"1|participant_id\"]\n", ")\n", + "_bounded_identity_prior = (\n", + " group_location_models[\"bounded generic identity\"]\n", + " .params[\"b\"]\n", + " .prior[\"1|participant_id\"]\n", + ")\n", "_log_prior = (\n", " group_location_models[\"log-scale location\"].params[\"a\"].prior[\"1|participant_id\"]\n", ")\n", @@ -1062,6 +1123,12 @@ " .prior[\"x|participant_id\"]\n", ")\n", "\n", + "assert _bounded_identity_prior.name == \"TruncatedNormal\"\n", + "assert _bounded_identity_prior.noncentered is False\n", + "assert np.isclose(float(_bounded_identity_prior.args[\"lower\"]), 0.2)\n", + "assert np.isposinf(float(_bounded_identity_prior.args[\"upper\"]))\n", + "assert isinstance(_bounded_identity_prior.args[\"mu\"], bmb.Prior)\n", + "assert _bounded_identity_prior.args[\"mu\"].name == \"TruncatedNormal\"\n", "assert _identity_prior.name == \"Gamma\"\n", "assert _log_prior.name == \"Normal\"\n", "assert _gen_logit_prior.name == \"Normal\"\n", @@ -1070,10 +1137,17 @@ "assert _gen_logit_prior.noncentered is False\n", "assert not isinstance(_matched_prior.args[\"mu\"], bmb.Prior)\n", "assert np.all(np.asarray(_matched_prior.args[\"mu\"]) == 0.0)\n", - "assert group_location_prior_table.iloc[:3][\"direct group RV\"].all()\n", - "assert not group_location_prior_table.iloc[:3][\"offset RV\"].any()\n", - "assert not group_location_prior_table.iloc[3][\"direct group RV\"]\n", - "assert group_location_prior_table.iloc[3][\"offset RV\"]\n", + "_rows_by_case = group_location_prior_table.set_index(\"case\")\n", + "_owner_cases = [\n", + " \"bounded generic identity\",\n", + " \"identity HDDM location\",\n", + " \"log-scale location\",\n", + " \"generalized-logit location\",\n", + "]\n", + "assert _rows_by_case.loc[_owner_cases, \"direct group RV\"].all()\n", + "assert not _rows_by_case.loc[_owner_cases, \"offset RV\"].any()\n", + "assert not _rows_by_case.loc[\"matched non-centered deviation\", \"direct group RV\"]\n", + "assert _rows_by_case.loc[\"matched non-centered deviation\", \"offset RV\"]\n", "mo.Html(group_location_prior_table.to_html(index=False, border=0))" ] }, @@ -1086,21 +1160,24 @@ } }, "source": [ - "The first three rows all have one unmatched group intercept and therefore\n", + "The first four rows all have one unmatched group intercept and therefore\n", "one population-location owner. HSSM sets only those generated priors to\n", "`noncentered=False`, retaining a direct group random variable and its\n", "location hyperprior.\n", "\n", "Their links change what that location means:\n", "\n", - "- under **identity**, the `a` hierarchy is the HDDM-derived response-scale\n", - " `Gamma` hierarchy;\n", + "- under **generic identity**, the lower-bounded LBA `b` hierarchy is a native\n", + " response-scale `TruncatedNormal`, so the generated location and every group\n", + " coefficient stay above `0.2`;\n", + "- under **identity for an analytical DDM**, the `a` hierarchy remains the\n", + " calibrated HDDM-derived response-scale `Gamma` hierarchy;\n", "- under **log**, `exp(mu)` is a reference median on the positive parameter\n", " scale, not the mean after integrating over group variation; and\n", "- under **generalized logit**, the inverse-linked `mu` is a bounded reference\n", " value, again not generally an expectation.\n", "\n", - "The fourth row is different: common `x` owns the population slope, so the\n", + "The fifth row is different: common `x` owns the population slope, so the\n", "group term is a zero-mean deviation and can stay non-centered. Under a log\n", "link such a zero deviation would be a neutral multiplicative factor\n", "`exp(0)=1`; under generalized logit it would leave the common predictor\n", @@ -1187,23 +1264,169 @@ "md_prefix": "" } }, + "source": [ + "### Bounded generic identity group location\n", + "\n", + "The pure predictor here is `b ~ 0 + (1 | participant_id)`, so each bounded\n", + "group coefficient is also the complete predictor and the final identity-linked\n", + "parameter. HSSM generates a centered native hierarchical `TruncatedNormal`;\n", + "this keeps Bambi's location and scale hyperpriors visible and constrains both\n", + "the generated population location and each LBA threshold coefficient above\n", + "the configured lower bound `0.2`. LBA2 separately requires `b > A`; this\n", + "construction satisfies that cross-parameter condition because the example\n", + "fixes `A=0.1`, not because a one-coefficient prior can enforce arbitrary\n", + "relationships between parameters.\n", + "\n", + "If a common slope, another group term, an offset, or another additive\n", + "contribution were present, this coefficient bound would protect only the\n", + "baseline—not the complete `eta`. HSSM warns for that mixed safe-generated\n", + "case. A transformed inverse link is what can constrain the complete additive\n", + "predictor when its image matches the parameter support." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "pHFh", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "clusterb_1|participant_id__factor_dim (4)\n", + "\n", + "b_1|participant_id__factor_dim (4)\n", + "\n", + "\n", + "cluster__obs__ (12)\n", + "\n", + "__obs__ (12)\n", + "\n", + "\n", + "cluster__obs__ (12) x rt,response_extra_dim_0 (2)\n", + "\n", + "__obs__ (12) x rt,response_extra_dim_0 (2)\n", + "\n", + "\n", + "\n", + "b_1|participant_id_sigma\n", + "\n", + "b_1|participant_id_sigma\n", + "~\n", + "Weibull\n", + "\n", + "\n", + "\n", + "b_1|participant_id\n", + "\n", + "b_1|participant_id\n", + "~\n", + "Truncated_normal\n", + "\n", + "\n", + "\n", + "b_1|participant_id_sigma->b_1|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "b_1|participant_id_mu\n", + "\n", + "b_1|participant_id_mu\n", + "~\n", + "Truncated_normal\n", + "\n", + "\n", + "\n", + "b_1|participant_id_mu->b_1|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "b\n", + "\n", + "b\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "b_1|participant_id->b\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "rt,response\n", + "\n", + "rt,response\n", + "~\n", + "Lba2_RV\n", + "\n", + "\n", + "\n", + "b->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "bounded_identity_group_location_graph = make_group_location_graph(\n", + " \"bounded generic identity\"\n", + ")\n", + "bounded_identity_group_location_graph" + ] + }, + { + "cell_type": "markdown", + "id": "NCOB", + "metadata": { + "marimo": { + "md_prefix": "" + } + }, "source": [ "### Identity-linked HDDM group location\n", "\n", "The direct `a_1|participant_id` node receives its population location and\n", "scale from the Gamma hierarchy. No offset node replaces that location.\n", - "HSSM's likelihood bounds still apply, but finite coefficient bounds are not\n", - "automatically propagated to generic identity-linked group priors. A bound on\n", - "one coefficient would not constrain the complete predictor after other effects\n", - "are added. Prefer a support-respecting transformed link when appropriate; the\n", - "[group-prior guide](https://lnccbrown.github.io/HSSM/how_to/specify_group_priors/)\n", - "explains the remaining identity-link choices and limitations." + "This analytical DDM path deliberately retains the calibrated HDDM family\n", + "instead of replacing it with the generic bounded Normal hierarchy above.\n", + "Both are response-scale group locations; the distinction is the safe-prior\n", + "family supplied by the likelihood implementation.\n", + "\n", + "Configured likelihood bounds are not themselves a hard-support prior. If a\n", + "complete predictor leaves them, HSSM substitutes a finite per-trial\n", + "log-likelihood floor of `-66.1`. The [group-prior\n", + "guide](https://lnccbrown.github.io/HSSM/how_to/specify_group_priors/)\n", + "separates coefficient support, complete-predictor support, and that finite\n", + "likelihood penalty." ] }, { "cell_type": "code", "execution_count": null, - "id": "pHFh", + "id": "aqbW", "metadata": {}, "outputs": [ { @@ -1313,7 +1536,7 @@ }, { "cell_type": "markdown", - "id": "NCOB", + "id": "TRpd", "metadata": { "marimo": { "md_prefix": "" @@ -1330,7 +1553,7 @@ { "cell_type": "code", "execution_count": null, - "id": "aqbW", + "id": "TXez", "metadata": {}, "outputs": [ { @@ -1440,7 +1663,7 @@ }, { "cell_type": "markdown", - "id": "TRpd", + "id": "dNNg", "metadata": { "marimo": { "md_prefix": "" @@ -1456,7 +1679,7 @@ { "cell_type": "code", "execution_count": null, - "id": "TXez", + "id": "yCnT", "metadata": {}, "outputs": [ { @@ -1566,7 +1789,7 @@ }, { "cell_type": "markdown", - "id": "dNNg", + "id": "wlCL", "metadata": { "marimo": { "md_prefix": "" @@ -1583,7 +1806,7 @@ { "cell_type": "code", "execution_count": null, - "id": "yCnT", + "id": "kqZH", "metadata": {}, "outputs": [ { @@ -1614,11 +1837,11 @@ "\n", "__obs__ (12) x rt,response_extra_dim_0 (2)\n", "\n", - "\n", + "\n", "\n", - "v_x\n", - "\n", - "v_x\n", + "v_Intercept\n", + "\n", + "v_Intercept\n", "~\n", "Normal\n", "\n", @@ -1630,14 +1853,28 @@ "~\n", "Deterministic\n", "\n", + "\n", + "\n", + "v_Intercept->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_x\n", + "\n", + "v_x\n", + "~\n", + "Normal\n", + "\n", "\n", "\n", "v_x->v\n", - "\n", - "\n", + "\n", + "\n", "\n", "\n", - "\n", + "\n", "v_x|participant_id_sigma\n", "\n", "v_x|participant_id_sigma\n", @@ -1645,7 +1882,7 @@ "Weibull\n", "\n", "\n", - "\n", + "\n", "v_x|participant_id\n", "\n", "v_x|participant_id\n", @@ -1653,33 +1890,13 @@ "Deterministic\n", "\n", "\n", - "\n", + "\n", "v_x|participant_id_sigma->v_x|participant_id\n", "\n", "\n", "\n", - "\n", - "\n", - "v_Intercept\n", - "\n", - "v_Intercept\n", - "~\n", - "Normal\n", - "\n", - "\n", - "\n", - "v_Intercept->v\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "v_x|participant_id->v\n", - "\n", - "\n", - "\n", "\n", - "\n", + "\n", "v_x|participant_id_offset\n", "\n", "v_x|participant_id_offset\n", @@ -1687,11 +1904,17 @@ "Normal\n", "\n", "\n", - "\n", + "\n", "v_x|participant_id_offset->v_x|participant_id\n", "\n", "\n", "\n", + "\n", + "\n", + "v_x|participant_id->v\n", + "\n", + "\n", + "\n", "\n", "\n", "rt,response\n", @@ -1723,7 +1946,7 @@ }, { "cell_type": "markdown", - "id": "wlCL", + "id": "wAgl", "metadata": { "marimo": { "md_prefix": "" @@ -1748,7 +1971,7 @@ { "cell_type": "code", "execution_count": null, - "id": "kqZH", + "id": "rEll", "metadata": {}, "outputs": [ { @@ -1815,7 +2038,7 @@ }, { "cell_type": "markdown", - "id": "wAgl", + "id": "dGlV", "metadata": { "marimo": { "md_prefix": "" @@ -1830,7 +2053,7 @@ { "cell_type": "code", "execution_count": null, - "id": "rEll", + "id": "SdmI", "metadata": {}, "outputs": [ { @@ -1871,7 +2094,7 @@ }, { "cell_type": "markdown", - "id": "dGlV", + "id": "lgWD", "metadata": { "marimo": { "md_prefix": "" @@ -1908,7 +2131,7 @@ }, { "cell_type": "markdown", - "id": "SdmI", + "id": "yOPj", "metadata": { "marimo": { "md_prefix": "" @@ -1962,7 +2185,7 @@ "app_config": { "width": "medium" }, - "header": "# /// script\n# requires-python = \">=3.12,<3.15\"\n# dependencies = [\n# \"bambi==0.20.0\",\n# \"graphviz==0.21\",\n# \"hssm @ git+https://github.com/lnccbrown/HSSM.git@a7f6892d387f4b19f35e5db01f648abbb8535910\",\n# \"marimo==0.24.0\",\n# \"matplotlib==3.11.1\",\n# \"numpy==2.4.6\",\n# \"pandas==3.0.5\",\n# \"pymc==6.3.1\",\n# ]\n# ///\n\n\"\"\"Explain link functions, their HSSM role, and link-aware safe priors.\n\nThis construction-only marimo tutorial introduces link functions from first\nprinciples, compares HSSM's identity and ``log_logit`` settings, and verifies\nthe current link-aware safe-prior and group-location behavior. No sampling is\nrequired.\n\nRun the pinned standalone environment locally or in Molab::\n\n uvx marimo edit --sandbox docs/tutorials/link_functions.py\n\nTo exercise an active HSSM checkout instead, ignore the inline environment::\n\n uv run --group notebook --group docs marimo edit --no-sandbox \\\n docs/tutorials/link_functions.py\n uv run --group notebook --group docs marimo check --strict \\\n docs/tutorials/link_functions.py\n uv run --group notebook --group docs marimo export html --no-sandbox \\\n docs/tutorials/link_functions.py \\\n --output /tmp/link-functions.html --force\n uv run --group notebook --group docs marimo export ipynb --no-sandbox \\\n docs/tutorials/link_functions.py \\\n --output docs/tutorials/link_functions.ipynb \\\n --include-outputs --force\n uv run ruff format docs/tutorials/link_functions.ipynb\n\"\"\"\n\n# ruff: noqa: B018, D401, E501, PLR1711 (generated marimo notebook: prose, cell display expressions, and bare returns)\n", + "header": "# /// script\n# requires-python = \">=3.12,<3.15\"\n# dependencies = [\n# \"bambi==0.20.0\",\n# \"graphviz==0.21\",\n# \"hssm @ git+https://github.com/lnccbrown/HSSM.git@f94b9e467d43386bf1e822a9e66d0fb7a9d68055\",\n# \"marimo==0.24.0\",\n# \"matplotlib==3.11.1\",\n# \"numpy==2.4.6\",\n# \"pandas==3.0.5\",\n# \"pymc==6.3.1\",\n# ]\n# ///\n\n\"\"\"Explain link functions, their HSSM role, and link-aware safe priors.\n\nThis construction-only marimo tutorial introduces link functions from first\nprinciples, compares HSSM's identity and ``log_logit`` settings, and verifies\nthe current link-aware safe-prior and group-location behavior. No sampling is\nrequired.\n\nRun the pinned standalone environment locally or in Molab::\n\n uvx marimo edit --sandbox docs/tutorials/link_functions.py\n\nTo exercise an active HSSM checkout instead, ignore the inline environment::\n\n uv run --group notebook --group docs marimo edit --no-sandbox \\\n docs/tutorials/link_functions.py\n uv run --group notebook --group docs marimo check --strict \\\n docs/tutorials/link_functions.py\n uv run --group notebook --group docs marimo export html --no-sandbox \\\n docs/tutorials/link_functions.py \\\n --output /tmp/link-functions.html --force\n uv run --group notebook --group docs marimo export ipynb --no-sandbox \\\n docs/tutorials/link_functions.py \\\n --output docs/tutorials/link_functions.ipynb \\\n --include-outputs --force\n uv run ruff format docs/tutorials/link_functions.ipynb\n\"\"\"\n\n# ruff: noqa: B018, D401, E501, PLR1711 (generated marimo notebook: prose, cell display expressions, and bare returns)\n", "marimo_version": "0.24.0" } }, diff --git a/docs/tutorials/link_functions.py b/docs/tutorials/link_functions.py index 477641129..90377f161 100644 --- a/docs/tutorials/link_functions.py +++ b/docs/tutorials/link_functions.py @@ -3,7 +3,7 @@ # dependencies = [ # "bambi==0.20.0", # "graphviz==0.21", -# "hssm @ git+https://github.com/lnccbrown/HSSM.git@a7f6892d387f4b19f35e5db01f648abbb8535910", +# "hssm @ git+https://github.com/lnccbrown/HSSM.git@f94b9e467d43386bf1e822a9e66d0fb7a9d68055", # "marimo==0.24.0", # "matplotlib==3.11.1", # "numpy==2.4.6", @@ -686,6 +686,28 @@ def _(build_silent_model, hssm, model_kwargs): }, } group_location_models = {} + _lba_data = _group_base_kwargs["data"].copy() + _lba_data["response"] = (_lba_data["response"] == 1).astype(int) + group_location_models["bounded generic identity"] = build_silent_model( + data=_lba_data, + model="lba2", + loglik_kind="analytical", + include=[ + { + "name": "b", + "formula": "b ~ 0 + (1 | participant_id)", + "link": "identity", + } + ], + A=0.1, + v0=1.0, + v1=1.0, + p_outlier=0.0, + prior_settings="safe", + noncentered=True, + process_initvals=False, + initval_jitter=0.0, + ) for _label, _case in _group_cases.items(): _spec = { "name": _case["parameter"], @@ -729,27 +751,43 @@ def _(build_silent_model, hssm, model_kwargs): def _(bmb, format_prior, group_location_models, link_name, mo, np, pd): _group_rows = [] _group_specs = [ + ( + "bounded generic identity", + "b", + "1|participant_id", + "bounded response scale", + "(0.2, inf)", + ), ( "identity HDDM location", "a", "1|participant_id", "response scale", + "(0, inf) via Gamma", + ), + ( + "log-scale location", + "a", + "1|participant_id", + "log-predictor scale", + "real line before exp", ), - ("log-scale location", "a", "1|participant_id", "log-predictor scale"), ( "generalized-logit location", "z", "1|participant_id", "generalized-log-odds scale", + "real line before gen_logit inverse", ), ( "matched non-centered deviation", "v", "x|participant_id", "identity predictor scale", + "real line", ), ] - for _label, _parameter, _term, _scale in _group_specs: + for _label, _parameter, _term, _scale, _support in _group_specs: _model = group_location_models[_label] _prior = _model.params[_parameter].prior[_term] _mu = _prior.args.get("mu") @@ -764,6 +802,7 @@ def _(bmb, format_prior, group_location_models, link_name, mo, np, pd): "outer family": _prior.name, "location": format_prior(_mu) if _mu_is_free else repr(_mu), "location scale": _scale, + "group-coefficient support": _support, "effective form": ( "centered location owner" if _prior.noncentered is False @@ -780,6 +819,11 @@ def _(bmb, format_prior, group_location_models, link_name, mo, np, pd): .params["a"] .prior["1|participant_id"] ) + _bounded_identity_prior = ( + group_location_models["bounded generic identity"] + .params["b"] + .prior["1|participant_id"] + ) _log_prior = ( group_location_models["log-scale location"] .params["a"] @@ -796,6 +840,12 @@ def _(bmb, format_prior, group_location_models, link_name, mo, np, pd): .prior["x|participant_id"] ) + assert _bounded_identity_prior.name == "TruncatedNormal" + assert _bounded_identity_prior.noncentered is False + assert np.isclose(float(_bounded_identity_prior.args["lower"]), 0.2) + assert np.isposinf(float(_bounded_identity_prior.args["upper"])) + assert isinstance(_bounded_identity_prior.args["mu"], bmb.Prior) + assert _bounded_identity_prior.args["mu"].name == "TruncatedNormal" assert _identity_prior.name == "Gamma" assert _log_prior.name == "Normal" assert _gen_logit_prior.name == "Normal" @@ -804,10 +854,17 @@ def _(bmb, format_prior, group_location_models, link_name, mo, np, pd): assert _gen_logit_prior.noncentered is False assert not isinstance(_matched_prior.args["mu"], bmb.Prior) assert np.all(np.asarray(_matched_prior.args["mu"]) == 0.0) - assert group_location_prior_table.iloc[:3]["direct group RV"].all() - assert not group_location_prior_table.iloc[:3]["offset RV"].any() - assert not group_location_prior_table.iloc[3]["direct group RV"] - assert group_location_prior_table.iloc[3]["offset RV"] + _rows_by_case = group_location_prior_table.set_index("case") + _owner_cases = [ + "bounded generic identity", + "identity HDDM location", + "log-scale location", + "generalized-logit location", + ] + assert _rows_by_case.loc[_owner_cases, "direct group RV"].all() + assert not _rows_by_case.loc[_owner_cases, "offset RV"].any() + assert not _rows_by_case.loc["matched non-centered deviation", "direct group RV"] + assert _rows_by_case.loc["matched non-centered deviation", "offset RV"] mo.Html(group_location_prior_table.to_html(index=False, border=0)) return (group_location_prior_table,) @@ -815,21 +872,24 @@ def _(bmb, format_prior, group_location_models, link_name, mo, np, pd): @app.cell def _(mo): mo.md(""" - The first three rows all have one unmatched group intercept and therefore + The first four rows all have one unmatched group intercept and therefore one population-location owner. HSSM sets only those generated priors to `noncentered=False`, retaining a direct group random variable and its location hyperprior. Their links change what that location means: - - under **identity**, the `a` hierarchy is the HDDM-derived response-scale - `Gamma` hierarchy; + - under **generic identity**, the lower-bounded LBA `b` hierarchy is a native + response-scale `TruncatedNormal`, so the generated location and every group + coefficient stay above `0.2`; + - under **identity for an analytical DDM**, the `a` hierarchy remains the + calibrated HDDM-derived response-scale `Gamma` hierarchy; - under **log**, `exp(mu)` is a reference median on the positive parameter scale, not the mean after integrating over group variation; and - under **generalized logit**, the inverse-linked `mu` is a bounded reference value, again not generally an expectation. - The fourth row is different: common `x` owns the population slope, so the + The fifth row is different: common `x` owns the population slope, so the group term is a zero-mean deviation and can stay non-centered. Under a log link such a zero deviation would be a neutral multiplicative factor `exp(0)=1`; under generalized logit it would leave the common predictor @@ -876,6 +936,39 @@ def make_group_location_graph(case): return (make_group_location_graph,) +@app.cell +def _(mo): + mo.md(""" + ### Bounded generic identity group location + + The pure predictor here is `b ~ 0 + (1 | participant_id)`, so each bounded + group coefficient is also the complete predictor and the final identity-linked + parameter. HSSM generates a centered native hierarchical `TruncatedNormal`; + this keeps Bambi's location and scale hyperpriors visible and constrains both + the generated population location and each LBA threshold coefficient above + the configured lower bound `0.2`. LBA2 separately requires `b > A`; this + construction satisfies that cross-parameter condition because the example + fixes `A=0.1`, not because a one-coefficient prior can enforce arbitrary + relationships between parameters. + + If a common slope, another group term, an offset, or another additive + contribution were present, this coefficient bound would protect only the + baseline—not the complete `eta`. HSSM warns for that mixed safe-generated + case. A transformed inverse link is what can constrain the complete additive + predictor when its image matches the parameter support. + """) + return + + +@app.cell +def _(make_group_location_graph): + bounded_identity_group_location_graph = make_group_location_graph( + "bounded generic identity" + ) + bounded_identity_group_location_graph + return + + @app.cell def _(mo): mo.md(""" @@ -883,12 +976,17 @@ def _(mo): The direct `a_1|participant_id` node receives its population location and scale from the Gamma hierarchy. No offset node replaces that location. - HSSM's likelihood bounds still apply, but finite coefficient bounds are not - automatically propagated to generic identity-linked group priors. A bound on - one coefficient would not constrain the complete predictor after other effects - are added. Prefer a support-respecting transformed link when appropriate; the - [group-prior guide](https://lnccbrown.github.io/HSSM/how_to/specify_group_priors/) - explains the remaining identity-link choices and limitations. + This analytical DDM path deliberately retains the calibrated HDDM family + instead of replacing it with the generic bounded Normal hierarchy above. + Both are response-scale group locations; the distinction is the safe-prior + family supplied by the likelihood implementation. + + Configured likelihood bounds are not themselves a hard-support prior. If a + complete predictor leaves them, HSSM substitutes a finite per-trial + log-likelihood floor of `-66.1`. The [group-prior + guide](https://lnccbrown.github.io/HSSM/how_to/specify_group_priors/) + separates coefficient support, complete-predictor support, and that finite + likelihood penalty. """) return