Skip to content

initial poisson race implementation#840

Open
hayden-johnson wants to merge 13 commits intomainfrom
add-poisson-race
Open

initial poisson race implementation#840
hayden-johnson wants to merge 13 commits intomainfrom
add-poisson-race

Conversation

@hayden-johnson
Copy link
Collaborator

@hayden-johnson hayden-johnson commented Oct 28, 2025

Issue #801

This PR introduces a preliminary implementation of the Poisson race model as described in Appendix A of:

A comparison of two response time models applied to perceptual matching

Key Differences from the original formulation

Stimulus-dependent rate parameters are not implemented. This would require trial-by-trial stimulus identifiers for log-probability evaluation.

Shape parameter generalization: we do not restrict the shape parameter to be an integer. This simplifies the implementation but requires evaluation of the incomplete gamma function using scipy, rather than using factorials.

@codecov
Copy link

codecov bot commented Oct 28, 2025

Codecov Report

✅ All modified and coverable lines are covered by tests.

Files with missing lines Coverage Δ
src/hssm/_types.py 100.00% <ø> (ø)
src/hssm/likelihoods/__init__.py 100.00% <ø> (ø)
src/hssm/likelihoods/analytical.py 100.00% <100.00%> (ø)
src/hssm/modelconfig/poisson_race_config.py 100.00% <100.00%> (ø)
tests/test_likelihoods_poisson_race.py 100.00% <100.00%> (ø)

... and 4 files with indirect coverage changes

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Copy link
Member

@AlexanderFengler AlexanderFengler left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @hayden-johnson this looks reasonable as far as the likelihood is concerned. We might want to check in on the trial-wise covariate issues if it remains unresolved.

Apart from that a few more high level comments to get the PR over the finish line.

  1. We need some tests for the basic likelihood
    2 We need the basic simulator implemented in ssm-simulators (separate small PR)
  2. We need a small notebook to that illustrates the model and shows a model fit (goal is to start having these for all incoming models and backfill for existing ones, so people have a much easier time getting starting with using specific models in the toolbox).

@AlexanderFengler AlexanderFengler marked this pull request as ready for review November 17, 2025 22:32
@cpaniaguam cpaniaguam requested a review from Copilot November 18, 2025 21:04
Copy link
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull Request Overview

This PR implements the Poisson race model for two-choice response time modeling, based on a comparison of response time models. The implementation uses continuous shape parameters and evaluates log-likelihoods analytically using gamma distributions.

Key changes:

  • Adds analytical log-likelihood computation for the Poisson race model with two accumulators
  • Includes comprehensive unit tests validating vectorization, exponential edge cases, and parameter validation
  • Registers the new model in HSSM's type system and likelihood exports

Reviewed Changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated 1 comment.

Show a summary per file
File Description
src/hssm/likelihoods/analytical.py Core implementation of logp_poisson_race function with parameter bounds and distribution registration
src/hssm/modelconfig/poisson_race_config.py Configuration defining model parameters, priors, and choices for the Poisson race model
tests/test_likelihoods_poisson_race.py Unit tests covering vectorization, exponential special case, parameter validation, and edge cases
src/hssm/likelihoods/__init__.py Exports the new logp_poisson_race function
src/hssm/_types.py Registers "poisson_race" as a valid model type
Comments suppressed due to low confidence (2)

src/hssm/likelihoods/analytical.py:1

  • Corrected grammar: 'We do not condition on the underlying stimulus condition.'
"""pytensor implementation of the Wiener First Passage Time Distribution.

src/hssm/likelihoods/analytical.py:1

  • The documentation states 'Response > 0 indicates accumulator 1', but the code at line 618 uses r_c = pt.switch(flip, r2, r1) where flip = response > 0. This means positive responses select r2/k2 (accumulator 2), not accumulator 1. The documentation should be corrected to match the implementation, or clarify the accumulator indexing convention.
"""pytensor implementation of the Wiener First Passage Time Distribution.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copy link
Collaborator

@cpaniaguam cpaniaguam left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The tests look great! I added a suggestion for adopting the more declarative and pythonic list comprehension pattern with a function to do the number crunching.

@review-notebook-app
Copy link

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

@hayden-johnson
Copy link
Collaborator Author

@AlexanderFengler

Added an implementation of the Poisson race model in ssm-simulators (see PR 246), along with a minimal example notebook in HSSM/docs/poisson_race.ipynb.

Lmk what you think! Feel free to suggest additions to the tutorial notebook.

hayden-johnson and others added 3 commits December 29, 2025 19:55
fix shape bounds

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
add pythonic syntax

Co-authored-by: Carlos Paniagua <cpaniaguam@gmail.com>
Copy link
Collaborator

@cpaniaguam cpaniaguam left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looking real good! Left a few more observations. I haven't looked at the notebook yet.

Comment on lines 575 to 579
Each accumulator time follows a Gamma distribution (with continuous
shape parameter) with shape parameter k and rate r. The per-trial
likelihood decomposes into the density of the winning accumulator
evaluated at the observed decision time and the survival function of
the losing accumulator at the same time.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be okay to update the first sentence from Each accumulator time follows a Gamma distribution (with continuous shape parameter) with shape parameter k and rate r. to Each accumulator time follows a Gamma distribution with continuous shape parameter k and rate r.?

Comment on lines 602 to 603
pytensor.tensor.TensorVariable
Per-trial log-likelihoods compatible with PyTensor graphs.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The return type in the function's signature says it's np.ndarray. Consider revising this.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Resolved, and added some more detail in the docstring.

k_l = pt.switch(flip, k1, k2)

log_pdf = (
k_c * pt.log(pt.maximum(r_c, epsilon))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not define r_c_safe as well?

Copy link
Collaborator Author

@hayden-johnson hayden-johnson Jan 22, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was originally only using the safe option for problematic computations (i.e., logarithms). For other computations, I was using unsafe variables; in either case, the evaluation is invalid and should be caught by the parameter checks.

I was inconsistent with this, however, in that I used rt_safe outside of the log function later in the pdf - I changed this now.

log_pdf = (
k_c * pt.log(pt.maximum(r_c, epsilon))
+ (k_c - 1.0) * pt.log(rt_safe)
- r_c * rt_safe
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just curious: why use safe versions for a standard product for one variable but not the other? It seems like it should be fine to use the untreated versions.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

True. I meant to only use safe values inside the log function. I changed this for consistency.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Notebooks looks good to me!

Copy link
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

Copilot reviewed 5 out of 6 changed files in this pull request and generated 8 comments.


💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Comment on lines +136 to +139
" \"r1\": ssms_params[\"r1\"],\n",
" \"r2\": ssms_params[\"r0\"],\n",
" \"k1\": ssms_params[\"k1\"],\n",
" \"k2\": ssms_params[\"k0\"],\n",
Copy link

Copilot AI Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The notebook creates a mapping between ssms_params (with r0/r1/k0/k1 naming) and true_params (with r1/r2/k1/k2 naming), where ssms r1 maps to HSSM r1 and ssms r0 maps to HSSM r2. However, there's an inconsistency: ssms k1 maps to HSSM k1, but ssms k0 maps to HSSM k2. This creates confusion because r1 from ssms goes to r1 in HSSM, but k1 from ssms goes to k1 in HSSM while k0 goes to k2. For clarity and consistency, the mapping should be: ssms r0->HSSM r1, ssms r1->HSSM r2, ssms k0->HSSM k1, ssms k1->HSSM k2 (consistently mapping lower accumulator to 1 and upper to 2, or vice versa).

Suggested change
" \"r1\": ssms_params[\"r1\"],\n",
" \"r2\": ssms_params[\"r0\"],\n",
" \"k1\": ssms_params[\"k1\"],\n",
" \"k2\": ssms_params[\"k0\"],\n",
" \"r1\": ssms_params[\"r0\"],\n",
" \"r2\": ssms_params[\"r1\"],\n",
" \"k1\": ssms_params[\"k0\"],\n",
" \"k2\": ssms_params[\"k1\"],\n",

Copilot uses AI. Check for mistakes.
"logp_ddm_bbox",
"logp_ddm_sdv_bbox",
"logp_full_ddm",
"logp_poisson_race",
Copy link

Copilot AI Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The POISSON_RACE distribution object is created but not exported in the init.py file. The init.py only exports logp_poisson_race function but not the POISSON_RACE distribution object itself. Looking at similar models like LBA2 and LBA3, both the logp function and the distribution object should be exported for consistency.

Copilot uses AI. Check for mistakes.
log_pdf = (
k_c * pt.log(r_c_safe)
+ (k_c - 1.0) * pt.log(rt_safe)
- r_c * rt
Copy link

Copilot AI Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the log_pdf computation, the term "- r_c * rt" should use "rt_safe" instead of "rt" for consistency with the rest of the computation. When rt is negative (rt <= 0), rt_safe is set to epsilon, but the original rt (which is negative) is still used in this calculation. This could lead to incorrect log-likelihood values for edge cases where rt is barely positive after subtracting t.

Suggested change
- r_c * rt
- r_c * rt_safe

Copilot uses AI. Check for mistakes.

# set bounds
poisson_race_params = ["r1", "r2", "k1", "k2", "t"]
poisson_race_bounds = {param: (0.0, np.inf) for param in poisson_race_params}
Copy link

Copilot AI Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The bounds for all parameters are set to (0.0, np.inf), which means the lower bound is 0.0 (inclusive). However, the parameter validation checks r1 > 0, r2 > 0, k1 > 0, and k2 > 0 (strictly greater than), meaning 0 is not actually a valid value. The bounds should be defined with an exclusive lower bound, or the implementation should use epsilon as the lower bound to be consistent with the parameter checks.

Suggested change
poisson_race_bounds = {param: (0.0, np.inf) for param in poisson_race_params}
poisson_race_bounds = {
"r1": (np.finfo(float).eps, np.inf),
"r2": (np.finfo(float).eps, np.inf),
"k1": (np.finfo(float).eps, np.inf),
"k2": (np.finfo(float).eps, np.inf),
"t": (0.0, np.inf),
}

Copilot uses AI. Check for mistakes.
Comment on lines +90 to +104
@pytest.mark.parametrize(
"param,bad_value",
[
("r1", 0.0),
("r2", -0.1),
("k1", 0.0),
("k2", -1.0),
("t", -0.2),
],
)
def test_poisson_race_parameter_validation(poisson_race_data, param, bad_value):
"""Invalid parameter values should produce a ParameterValueError."""
assert_parameter_value_error(
poisson_race_data, theta_poisson_race, **{param: bad_value}
)
Copy link

Copilot AI Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The test only validates that 0.0 and negative values are rejected, but doesn't test boundary values like very small positive numbers (e.g., 1e-10). Additionally, there's no test case for when t equals 0, which is a valid edge case according to the bounds (t >= 0) but might behave differently than other values.

Copilot uses AI. Check for mistakes.
Comment on lines 13 to +14
"logp_full_ddm",
"logp_poisson_race",
Copy link

Copilot AI Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The name 'logp_poisson_race' is exported by all but is not defined.

Suggested change
"logp_full_ddm",
"logp_poisson_race",
"logp_full_ddm",

Copilot uses AI. Check for mistakes.
@AlexanderFengler
Copy link
Member

Overall looks like this is basically in shape, apart from some of the CoPilot comments (some of which seem worth addressing). thanks @hayden-johnson

fix doctring typo

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants