Repository navigation
Add the Linear Ballistic Accumulator (lba) on the shared race core (#349) - #500
Open
GidonFrischkorn wants to merge 35 commits into
Open
GidonFrischkorn wants to merge 35 commits into
GidonFrischkorn wants to merge 35 commits into
Conversation
…frastructure (#349) Rebuilt from 349-lba-model (2c0a8c6) on race-model-core: the model, its distribution functions (now R/distributions-lba.R), Stan chunks and tests, with the coding, validation, revert, observables, log-space helpers and epred integrator taken from the shared race code. The model registers its shared-test cases in helper-model-lba.R and its column roles and response annotations in its own file; the benchmark moves from internal-model-tests.Rmd to tests/internal/benchmark_lba.Rmd. Changes against the branch: posterior_epred() integrates with race_expected_time() over the range race_time_range() finds instead of the 20-points-per-decade trapezoid; category names follow the shared checks (case-sensitive Stan keywords, the model's own parameters, Intercept, and the identifiers of the generated likelihood such as b, n, t), so log and exp are now allowed and b, A, drift and the like refused; accumulators must name exactly the formula categories and constant counts must be positive integers; a response time of exactly 0 is refused.
brms passes a fixed dpar as one number, and s is fixed by default, so redrawing the non-positive normal drifts indexed s[neg] into NA: the draw errored or looped forever whenever a drift past the first draw was not positive. Present on 349-lba-model.
… gamma and Frechet kernels dlba() errored for the gamma distribution at rt == ndt and for an NA response, and for the Frechet distribution at rt <= ndt and for an NA response; plba() errored at q <= ndt for both, and qlba(0) errored for the Frechet distribution. Times at or below the non-decision time and NA responses are now masked before the kernels run, so the density is 0 or NA and the probability is 0 as for the normal and lognormal distributions.
The help page said the distribution functions accept sp = 0, but dlba() returned NaN for every distribution and plba() returned 0 for three of them. sp must now be positive, as in the lba() model, and the help page says so.
…r draw when ndt is a single number When a prep stored ndt as one number, the number of draws was taken from ndt and the start points collapsed to a single draw, so the predicted RTs were far too concentrated (SD 0.0009 instead of 0.04 in a fast race). The number of draws now comes from the prep.
…rift gradient is exact The normal-drift LBA density took log(Phi(hi) - Phi(lo)) from std_normal_lcdf, whose Stan gradient is a piecewise polynomial approximation. The log-density itself was correct, but its partial derivative with respect to the drift rate was off by up to 17 % for slow responses and negative drifts, which biases the sampler's trajectories rather than the posterior's shape. The Phi differences now come from erfc, whose derivative is exact, with the lcdf kept only beyond z = -37 where erfc underflows. Values are unchanged (2e-8 against the R twin). A cmdstan gradient diagnostic now checks the drift partial against finite differences at two positive-drift points the old kernel got wrong by 2.2 % and 0.17 %.
…race in posterior_epred The category-name test checked one identifier of the generated Stan code, so removing another from the list went unnoticed. A second test compares posterior_epred() of a slow normal race with direct quadrature, which a fixed integration range of 1e-3 to 10 seconds misses by about 2 percent.
…or as well After the Phi differences moved to erfc, the drift partial of the normal-drift LBA density was still off by 0.2 % at slow negative-drift responses (v / s = -6.7): the truncation denominator log Phi(v / s) was still taken from std_normal_lcdf, whose Stan gradient is approximate. The denominator of the density and of the survivor now use the same erfc-based log Phi. Values are unchanged. The gradient test now covers the two negative-drift points as well.
…mma shape partial, and what plba()/qlba() describe
…or as well After the Phi differences moved to erfc, the drift partial of the normal-drift LBA density was still off by 0.2 % at slow negative-drift responses (v / s = -6.7): the truncation denominator log Phi(v / s) was still taken from std_normal_lcdf, whose Stan gradient is approximate. The denominator of the density and of the survivor now use the same erfc-based log Phi, which takes its upper tail through log1m so that the survivor keeps 1 - Phi(v / s) at full relative precision (a plain erfc there drove the survivor to its floor for v / s >= 8). Values are unchanged. The gradient test now covers the two negative-drift points as well.
The posterior_epred() paragraph no longer quotes measured values, the gamma note names where the shape partial is worst without the test-grid shares, and ?lba_dist states which arguments qlba() takes as single values. The early return in qlba() is spelled as an if/else.
# Conflicts: # NAMESPACE # NEWS.md
The macOS and Windows pull-request jobs ship no CmdStan and set
BMM_TEST_STAN=false. With only skip_if_not_installed("cmdstanr"), the
seven tests that compile the kernel program errored there with
"CmdStan path has not been set yet"; they now skip.
…s in test-conventions.R test-model_lba-stan.R compiles two Stan programs and samples one; test-model_lba.R passes backend = "cmdstanr" to stancode() six times.
The test draws 20 trials without a seed. One draw in fifty has no error response, and the custom version then warns that the category 'wrong' is not in the data, so the test depended on where the RNG stream stood.
…ollow-up PR The gamma, Frechet and lognormal drift distributions carry most of the hard numerics of the LBA and cannot be checked against one source, so they move to their own pull request. Removed: their entries in the spec table, their R kernels and Stan chunks, their tests and benchmark fits, and their text in ?lba, ?lba_dist and NEWS. The distribution argument of lba() and the d/p/q/rlba() functions stays, with "normal" as the only value, and the normal model's generated Stan code is unchanged.
# Conflicts: # tests/testthat/test-conventions.R
…ires Rule R4 now checks that every provider has them. The sd rates and the settable links are the values default_prior() and settable_links() give the simple and custom versions, written out.
…te mean posterior_epred() integrates the race survivor up to t_max = 1e7 s. On a trial where only one normal-drift accumulator races the mean is infinite and the returned value is the mean of the truncated race. Such trials now trigger one warning per call with their count; values are unchanged. race_time_range() marks a bracket truncated whenever S(t) t is above its 1e-12 tolerance at t_max, which every two-accumulator normal race is, although the mass it leaves out is below 1e-8 s. The warning therefore also requires the survivor to decay no faster than t^-1.5 between t_max/10 and t_max, which singles out a lone accumulator (a normal-drift race decays like t^-k for k racing accumulators).
A custom formula with one category describes a race without a choice, and with normal drift a single accumulator has no finite mean response time. check_model() now stops with a message naming the category. The help page says the custom version needs two categories, and its example category names no longer end in a number, which the version refuses.
The generated custom likelihood declared b, A, t, n, i, j, drift, win, lp and reps, so response categories with those names were refused. The locals are now lba_thr, lba_sp, lba_dt, lba_cnt, lba_i, lba_j, lba_drift, lba_win, lba_lp and lba_reps. A category cannot contain an underscore, so none of them can collide, and the refused names shrink to the arguments rt and response. The generated Stan code is identical to the old one once the new names are mapped back, for both versions, and the compiled likelihood returns identical values on 60 observations with zero, one and two accumulators per category. stanc accepts categories named b, A, t, n, j, win, lp and delta (a local of a helper function).
The Default priors section called the rate-1/rate-2 split a package-wide rule, which the circular and SDT models do not follow. It now names ddm() and ezdm(), whose rates lba() shares. The rates themselves are unchanged.
The help page now says that posterior_epred() integrates the race per draw and observation, at about half a millisecond each (measured on 200-1000 draws by 10-50 observations), and that ndraws limits the cost. The paragraph is reflowed.
The warning said the mean was infinite, which holds for the normal-drift race it fires on (a single accumulator, survivor ~ 1/t) but not for a race that decays between t^-1 and t^-1.5, whose finite mean rests on times beyond 1e7 s. It now states the decay order it tests and covers both cases.
The core now reserves n, which brms uses as the index of its random-effects loop, so the LBA test of the freed former locals keeps b, A and t.
A race in which a single normal-drift accumulator runs has no finite mean response time. posterior_epred() returned the mean of the race truncated at 1e7 s for it; it now returns Inf for every draw whose survivor decays no faster than t^-1.5, skips the integral for those draws, and warns once per call with the number of observations affected, pointing to the median or quantiles of posterior_predict(). The other draws keep the range searched over all draws, so their values are unchanged.
The help page pointed only to dlba() and rlba(); it now lists all four exported distribution functions, as ?lnr does.
…tor counts A lone normal-drift accumulator has an infinite mean in every draw, but the local decay test only saw it while the 1/t tail stayed above the range search's tolerance: from a mean drift of about 8 (s = 1) such trials returned a finite value without a warning, and one low-drift draw turned the other draws of the observation to Inf. An observation with a single normal-drift racer now returns Inf in every draw from its counts alone; every other observation is integrated over its own range as before.
This was referenced Oct 18, 2026
…o race-lba # Conflicts: # tests/testthat/test-conventions.R
check_data() builds .lba_cat and .lba_n1..n<k> through the race core, but lba had no built_data_columns() method, so the declaration contract test failed and standata(fit, newdata = user) treated user-layout newdata as the model frame and asked brms for the coding columns. The simple version declares the two pooled count columns, the custom version one per formula category.
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Add the Linear Ballistic Accumulator (lba) on the shared race core (#349)
The Linear Ballistic Accumulator (Brown & Heathcote, 2008) for choice response times with two or more options. Each option has an accumulator that starts at a uniform point between 0 and
spand rises linearly at a drift rate drawn on every trial from a normal distribution truncated at zero. The first to reach the thresholdb = gap + spgives the response, and a non-decision time is added.lba()has two versions:simple(one mean drift for the correct response and one shared by the errors) andcustom(one mean drift per named response category, with optional per-categoryaccumulators). Also addsdlba(),plba(),qlba()andrlba().This PR adds the normal-drift LBA only. The gamma, Fréchet and lognormal drift distributions are in #565, which builds on this branch. The PR now targets
dev-1.6.0, since the race core (#499) is merged.What changed since the last description
posterior_epred()with an infinite mean. A race in which only one normal-drift accumulator runs has no finite mean response time. That happens on a trial where anaccumulatorscolumn leaves one category.posterior_epred()used to return a finite number for such a trial: the mean of the race cut off at 1e7 s, which looks like a real mean. It now returnsInffor such a trial, decided from the accumulator counts so that it holds whatever the drift, and warns once per call with the number of observations affected, pointing to the median or quantiles ofposterior_predict(). All other values are unchanged.lba_prefix, so categories such asb,Aortwork. The model itself reserves onlyrtandresponse. The names brms declares (nin the random-effects loop,seq/start/end/ptarget/nnunder threading) are refused for every race model by the race core. That needs the small core PR that addsn,ptargetandnn, which is merged into this branch: merge that PR first. The generated Stan code is identical to the old one after mapping the names back, and old and new likelihoods give identical values.?lbastates the cost ofposterior_epred()(437 to 537 µs per draw and observation here) and thatndraws =limits it. It names the models whose random-effect SD rateslba()follows (ddm(),ezdm()), and See Also lists all four distribution functions. Makingposterior_epred()faster for all race models is a separate issue.ndtandsalready had theexponential(2)SD rate thatlnr()andrdm()now use too.Checks (local, macOS, R 4.6.1)
NOT_CRAN=true, Stan tests on, 4 workers): 2001 tests, 0 failed, 0 skipped.R CMD check --no-manual: 0 errors, 0 warnings, 1 NOTE. The NOTE ("Cannot extract version info" for NEWS.md) comes from the two development headings that ADR 0002 puts ondev-*branches. Everydev-*branch has it.Merging next to #501 and #502
The three model PRs merge in any order. The second and third conflict in
NEWS.md(keep both bullets),NAMESPACEandman/bmm.Rd(rundevtools::document()), and, with #501,tests/testthat/test-conventions.R(keep both files' entries). After this PR merges, #565 needs a retarget todev-1.6.0.For the reviewer
Most of the code was written with AI assistance and has had AI-assisted reviews. What it needs from a person is judgement on the model: the likelihood against Brown and Heathcote (2008), the parameterisation and defaults, and whether the documentation makes sense to someone who fits LBAs. About 4 hours.
Pass 1: likelihood and Stan kernel (about 1.5 h).
inst/stan_chunks/lba_normal_functions.stanandlba_shared_functions.stan(the single-accumulator density and survivor, with the drift truncated at zero), then the generated likelihood (.lba_stan_code()) andconfigure_model()inR/model_lba.R. Check the density and CDF against the paper, and the truncation of the drift distribution.Pass 2: R distribution functions and postprocessing (about 1.5 h).
R/distributions-lba.R(dlba(),plba(),qlba(),rlba()), thenlog_lik,posterior_predictandposterior_epredinR/model_lba.R. Tests:tests/testthat/test-model_lba.R,test-model_lba-stan.R(Stan against R),test-model_lba-stan-accuracy.R.Pass 3: docs (about 1 h).
?lbaand?lba_dist. Would someone who knows the LBA from the paper or fromrtdistsknow how to fit it from these pages?Skim or skip:
man/,NAMESPACE,tests/testthat/helper-model-lba.R(cases for the cross-model tests),tests/internal/benchmark_lba.Rmd.Questions where an outside view helps most:
gap(threshold minus start-point range) andsp(start-point range) on the log link, so thatb > spholds by construction. Papers andrtdistsusebandA. Is the mapping clear in?lba, and is it the right choice?sfixed to 1. Only ratios of drift,gap,spandsare identified, sosis fixed by default andcheck_model()warns if a formula frees it. Do you fix the scale this way in your own work, or per condition?normal(3, 1)(correct) andnormal(1, 1)(errors);gap,spandndtare on the log link. Do they fit the tasks you would use the LBA for?loo()by a constant per error trial. Is that the right default?posterior_epred()returningInffor trials where one accumulator races (see above). Is that what you would expect from the model's expected value?Already checked, so no need to repeat:
log_likwithin 3.2e-14 relative; seeded random draws are identical.sp = 0is refused rather than returningNaN.S, refused becausesis a parameter).