New aDDM estimation function

[addm_recovery_issue_mre.py](https://github.com/user-attachments/files/30268131/addm_recovery_issue_mre.py)

## Summary

While validating the `hssm.addm.aDDM` model with a parameter-recovery test, I found
that `eta` and `x0` recover well, but `a` (boundary), `b` (collapse) and `t`
(non-decision time) come out **systematically biased** — `a↑`, `b↑`, `t↓` — with
**tight posteriors that exclude the true values**, even though the chains converge
cleanly (`r_hat ≤ 1.01`).

I'm not assuming this is a bug. I'd like to ask whether it reflects (a) expected
partial identifiability of the collapsing-boundary parameterization, (b) a
convention mismatch on my side in how I construct the per-trial covariates, or
(c) an actual discrepancy between the `ssms` simulator and the JAX likelihood
kernel.

## Environment

- **hssm 0.4.0**, installed from source at commit `1823172`. I installed from my
own fork, but it is in sync with upstream `main` and carries **no local
modifications**, so this should be reproducible on upstream as-is.
- ssms 0.13.2, jax 0.11.0, numpyro 0.21.0, pymc 6.1.0, bambi 0.19.0,
arviz 1.2.0, pytensor 3.1.3, numpy 2.4.6
- Python 3.12.13, Linux, CUDA 12 (RTX 3090)

## Reproduction

The attached script (`addm_recovery_issue_mre.py`) is self-contained. It:

1. Draws a known `theta = (eta, kappa, a, b, x0, t)`.
2. Generates per-trial covariates `r1`, `r2`, `flag`, `sigma = 1`, plus a long
alternating fixation schedule (`sacc_array` onsets, first entry exactly `0.0`,
~0.3 s fixations; the first fixation is forced to be longer than `t` so the
NDT gate is satisfied).
3. Simulates `(rt, response)` with the `ssms` `'addm'` simulator, conditioning on
those fixations through `extra_fields`.
4. Truncates the schedule to onsets `< rt`, giving the *observed* `sacc_array`
and `d`.
5. Fits with `aDDM(data=df)` (six scalar parameters, default priors) using
`sampler="numpyro"`, `draws=750`, `tune=750`, `chains=2`.

## Results — 3 random parameter sets, 700 trials each

| param | set 1 (true → est) | set 2 (true → est) | set 3 (true → est) | in 94% HDI |
|-------|--------------------|--------------------|--------------------|:----------:|
| eta | 0.40 → 0.44 | 0.87 → 0.90 | 0.46 → 0.46 | 2/3 |
| x0 | 0.11 → 0.09 | −0.12 → −0.09 | 0.17 → 0.14 | 3/3 |
| kappa | 1.21 → 0.85 | 1.03 → 0.98 | 0.50 → 0.49 | 1/3 |
| a | 1.74 → **2.18** | 2.09 → **2.19** | 2.12 → **2.24** | 2/3 |
| b | 0.96 → **1.47** | 0.47 → **0.59** | 0.16 → **0.20** | 2/3 |
| t | 0.12 → **0.00** | 0.06 → **0.04** | 0.07 → **0.04** | 2/3 |

Overall 94% HDI coverage: **12/18 (67%)**, `max r_hat = 1.010`.

The direction of the bias is consistent rather than random: `a` is overestimated
in all three sets, `b` overestimated in all three, and `t` underestimated in all
three.

## Posterior correlations of (a, b, t)

| | set 1 | set 2 | set 3 |
|-------|------:|------:|------:|
| a ↔ b | +0.87 | +0.89 | +0.78 |
| a ↔ t | −0.47 | −0.83 | −0.68 |
| b ↔ t | −0.41 | −0.73 | −0.52 |

So `(a, b, t)` clearly trade off along a ridge, which by itself is unsurprising for
a collapsing boundary. What puzzles me is that the **marginal HDIs are narrow yet
still exclude the truth** — e.g. in set 1, `t` has a true value of 0.117 but a 94%
HDI of `[0.000, 0.012]`. Pure non-identifiability would give wide intervals that
still cover the true value; a narrow interval sitting confidently away from the
truth looks more like the likelihood optimum being displaced relative to the
generating process.

## Questions

1. Is the `a` / `b` / `t` trade-off expected for this parameterization, and are
these three intended to be jointly identifiable in practice?
2. Are my covariate conventions correct? Specifically: should `sacc_array` be
fixation onsets **relative to stimulus onset, with the first entry exactly
`0.0`**, and is `t` meant to fall *inside* the first fixation on that same
clock? If onsets are instead expected on the post-NDT decision clock, that
alone could explain the collapse of `t` toward zero.
3. Should the `ssms` `'addm'` simulator and the JAX likelihood kernel agree
exactly (same boundary-collapse formula and discretisation)? In other words,
is simulate-then-fit a valid recovery check for this model, or is some
mismatch between the two implementations expected?

Happy to run further diagnostics or share the posterior traces if that helps.

Start with the attached addm_recovery_issue_mre.py script and reproduce the parameter-recovery results using the stated environment and settings. Compare the ssms 'addm' simulator inputs, fixation-onset conventions, and the JAX likelihood kernel for the a, b, and t parameters. Done means determining whether the bias is expected, caused by covariate conventions, or reflects a simulator-likelihood discrepancy.

Written by the indexing model from the issue text.