Default priors with analytical likelihood

[igrahek]

When fitting regression models without specifying priors (assuming that this leads to default priors) models run when using approx_differentiable, but not analytical. For analytical, the error indicates that the a>= 0 check has failed.

[digicosmos86]

Hi! Could you provide some additional information? Any minimal reproducible code or output will be helpful. Thanks!

[igrahek]

Hi! Sure, sorry about the lack of additional info! The model below will sample fine, but if the likelihood is changed to analytical, it will throw the error pasted below the model. I am getting this error on my own dataset and haven't confirmed that it will be reproduced on an example dataset, but I imagine that it would be.

model1 = hssm.HSSM(
    data=data,
    model = 'ddm',
    loglik_kind="approx_differentiable",
    include=[
        {
            "name": "v",
            "formula": "v ~ 1 + x * y + (1 + x * y | subject) ",
            "link": "identity",
        },
        {            
            "name": "a",
            "formula": "a ~ 1 + x * y + (1 + x * y  | subject) ",
            "link": "identity",
        },
	    {            
            "name": "z",
            "formula": "z ~ 1  + x + (1 + x | subject)",
            "link": "identity",
        },
    ],
)

Error:

Compiling...
Traceback (most recent call last):
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pytensor/compile/function/types.py", line 970, in __call__
    self.vm()
pymc.logprob.utils.ParameterValueError: a >= 0

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/oscar/data/ashenhav/igrahek/Reconfig_Cost_Reward/Experiment7/ddm/bash/Mixed/../../models/model3.py", line 47, in <module>
    modelObject = model.sample(
                  ^^^^^^^^^^^^^
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/hssm/hssm.py", line 554, in sample
    self._inference_obj = self.model.fit(inference_method=sampler, **kwargs)
                          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/bambi/models.py", line 325, in fit
    return self.backend.run(
           ^^^^^^^^^^^^^^^^^
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/bambi/backend/pymc.py", line 96, in run
    result = self._run_mcmc(
             ^^^^^^^^^^^^^^^
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/bambi/backend/pymc.py", line 211, in _run_mcmc
    idata = pymc.sampling_jax.sample_numpyro_nuts(
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pymc/sampling/jax.py", line 623, in sample_numpyro_nuts
    init_params = _get_batched_jittered_initial_points(
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pymc/sampling/jax.py", line 247, in _get_batched_jittered_initial_points
    initial_points = _init_jitter(
                     ^^^^^^^^^^^^^
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pymc/sampling/mcmc.py", line 1221, in _init_jitter
    model.check_start_vals(point)
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pymc/model.py", line 1697, in check_start_vals
    initial_eval = self.point_logps(point=elem)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pymc/model.py", line 1732, in point_logps
    self.compile_fn(factor_logps_fn)(point),
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pymc/pytensorf.py", line 763, in __call__
    return self.f(**state)
           ^^^^^^^^^^^^^^^
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pytensor/compile/function/types.py", line 983, in __call__
    raise_with_op(
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pytensor/link/utils.py", line 535, in raise_with_op
    raise exc_value.with_traceback(exc_trace)
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pytensor/compile/function/types.py", line 970, in __call__
    self.vm()
pymc.logprob.utils.ParameterValueError: a >= 0
Apply node that caused the error: Check{a >= 0}(Switch.0, All{axes=None}.0)
Toposort index: 1639
Inputs types: [TensorType(float64, shape=(32719,)), TensorType(bool, shape=())]
Inputs shapes: [(32719,), ()]
Inputs strides: [(8,), ()]
Inputs values: ['not shown', array(False)]
Outputs clients: [[Check{z >= 0}(Check{a >= 0}.0, All{axes=None}.0)]]

Backtrace when the node is created (use PyTensor flag traceback__limit=N to make it longer):
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pymc/logprob/basic.py", line 541, in conditional_logp
    q_logprob_vars = _logprob(
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/functools.py", line 909, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 132, in logp
    return class_logp(value, *dist_params)
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/hssm/distribution_utils/dist.py", line 381, in logp
    logp = loglik(data, *dist_params)
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/hssm/likelihoods/analytical.py", line 278, in logp_ddm
    checked_logp = check_parameters(logp, a >= 0, msg="a >= 0")
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pymc/distributions/dist_math.py", line 78, in check_parameters
    return CheckParameterValue(msg, can_be_replaced_by_ninf)(expr, all_true_scalar)
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pytensor/graph/op.py", line 295, in __call__
    node = self.make_node(*inputs, **kwargs)
  File "/users/igrahek/anaconda/pyHSSM/lib/python3.11/site-packages/pytensor/raise_op.py", line 90, in make_node
    [value.type()],

HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.

[AlexanderFengler]

Ok this is just a violation of the a >= 0 by (very likely the initial) some sample in the chain. We were using Bambi defaults here which is probably inappropriate with this probably inappropriate in this case.

Thank you for letting us know we will fix that.

To sample in the meanwhile you have multiple options:

1. set the initval parameter on the betas for a to 0, apart from the intercept for which 1 is appropriate.
2. pass a link function that makes a positive
3. pass you own priors of course which will respect the bounds

[igrahek]

Great, thank you!