Changelog • cogmod

cogmod 0.3.5

Breaking changes

All three change the priors cogmod_priors() returns, so a model refitted under 0.3.5 is not the model it was under 0.3.4.

poutlier’s default prior is centred at 0.25% rather than 0.7%: normal(-6, 1) on the logit link instead of normal(-5, 1), and exponential(280) instead of exponential(100) when poutlier is left out of bf(), the two medians still matching. Every family with an outlier component is affected. Miller (2024) puts the outlier proportion below 0.5%. On the untrimmed Illusion Game data of vignette("outliers") the share of responses at chance speed is -5.3 on the logit scale pooled but -6.2 for the typical participant, which is what the intercept of poutlier ~ 1 + (1 | Participant) estimates; and on 18 such fits of 324,000 trials each the posterior sat at -6.0 to -8.2, below the old centre every time and past its 1st percentile in nine. The advice for trimmed data, normal(-7, 0.5), is unchanged.
cogmod_ddm()’s boundary and bias get priors. brms left boundary flat, intercept and slopes - an improper prior on a softplus link, which brms refuses to sample from with sample_prior = "only" - and the slopes of bias flat. boundary now gets normal(1, 1), a median separation of 1.31 with 95% between 0.32 and 3.0 (lognormal(0.27, 0.5) when left out of bf()), and normal(0, 0.5) on its slopes. bias keeps brms’s logistic(0, 1) intercept, now under 0 + Intercept too, and gets normal(0, 0.5) on its slopes. A group-level SD on either now gets exponential(1), as every other parameter the family names does, instead of brms’s student_t(3, 0, 2.5).
cogmod_invgaussian(), cogmod_exwald() and cogmod_bisa() give boundary a prior, for the same reason: normal(0.5, 1) on the softplus link, a median threshold of 0.97 with 95% between 0.21 and 2.5 (lognormal(0, 0.6) when left out of bf()), normal(0, 0.5) on its slopes and exponential(1) on its group-level SDs.

Bug fixes

cogmod_lba2()’s density could be NaN. The LBA’s closed form needs the difference of two normal densities, phi(z1) - phi(z2), and took it as phi(z1) times -expm1() of their log ratio, whichever of the two was the larger. Where an accumulator’s drift should have carried it past the whole threshold range well before the response - in the decision-making article’s model, the error accumulator with a positive drift and a small SD - phi(z2) is the larger: the expm1() overflowed to infinity while phi(z1) underflowed to zero, and the density was 0 * Inf = NaN. Under the wide driftone prior tried in benchmarks/driftone_priors/, CmdStan met it as log_mix: lambda2 is -nan, and since Stan scores a NaN log density as an infinite potential, each such proposal was a divergent transition. The difference is now taken from the larger of the two densities. Over 36,000 parameter values, from around the article’s fit out to the extremes, the Stan density was NaN at 3,301 and is now NaN at none; through grad_log_prob(), 37 had a finite value with a non-finite gradient and 14 a -Inf value with a non-finite partial - a NaN gradient for the whole model once log_mix() weights that component by zero - and now none do. R’s dcogmod_lba2() had the same formula, and there the NaN reached a mask and errored (“NAs are not allowed in subscripted assignments”); log_lik() and p_outlier() call the same density.
The LBAs’ densities no longer fall to -Inf where they are merely small. The winner’s density was formed linearly and underflowed to zero some 38 standardized units out, so that a 12 s response the race puts at a log density of -777 was handed to the outlier component, at -1864. The loser’s survival is taken given a positive drift, over Phi(drift / sigma), which Stan’s Phi() makes exactly zero below -37.5: the survival came back -Inf where it is near 1, and around there the terms it was built from had lost enough digits to put the race density out by up to 21 log units. Below 1e-280 both are now taken again in log space, from the normal mass on the start-point interval and its mean there, each formed relative to the density at the interval’s near end so that nothing of order z^2 / 2 is subtracted. On the same 36,000 values, 17,964 that were -Inf and 3,293 that were NaN are now finite. Of the 11,310 that were finite under both, the median change is 0 and the 99th percentile 6e-14; of the 56 that moved by more than 1e-12, the 20 that moved most all moved onto a quadrature of the density, to within 1e-5. cogmod_lba1() shares the density and gets the same fix.
The LBAs’ gradients are exact in the tails. The truncation constant log P(v > 0), paid on every observation, was Stan’s std_normal_lcdf(), whose value is right but whose partials are approximate, and cogmod_lba1()’s was log1m_exp(std_normal_lcdf(-x)), the same thing the long way round. Both now go through cogmod_log_Phi(), except where the loser’s survival has just formed Phi(drift / sigma) itself and takes its log directly: exact partials either way, and that reuse is what keeps the fix at the old cost (cogmod_log_Phi() for both accumulators measured 1.09 and 1.12 times the old cost on cogmod_lba2()). Over the gradient check’s grid (benchmarks/gradient_check.R) the worst error against central differences goes from 2.6e-5 to 2.1e-6 for cogmod_lba2(), and 21 of its 22 points from around 1e-6 to around 1e-9; the remaining point, and cogmod_lba1()’s worst at 1.2e-6, are the same to the last digit before and after, so they come from elsewhere in the density. Nothing got dearer: per gradient evaluation cogmod_lba2() costs 0.98 and 1.02 times what it did, and cogmod_lba1(), where one function replaces two, 0.92 and 0.93, over two runs of benchmarks/gradient_cost.R (5,000 trials, median of 21 alternating blocks).
dcogmod_lba1() with a negative drift. Its truncation constant was log1p(-exp(log Phi(-drift / sigma))), which loses digits below a drift ratio of about -6 and is -Inf beyond -8.3, taking the density with it; it is now pnorm(drift / sigma, log.p = TRUE). cogmod_lba1() puts a softplus link on the drift, so fitted models were not affected.

Changes to starting values

cogmod_lnr()’s error accumulator starts slower than the correct one, at nuone = 0 against mu = 0.7, as cogmod_rdm()’s and cogmod_lba2()’s already do. The old equal start was a 50/50 race; the new one implies 16% errors. This is for consistency, not speed: refitted cold on three models, 4 seeds each, the closer start changed the median warmup cost by 8% or less (cheaper on one model, dearer on another), with no divergences or stuck chains under either start. The nuone prior is unchanged, and starting values do not change a posterior, so no fitted model needs revisiting.
cogmod_inits() now covers cogmod_betagate(), which it refused (“nothing to offer”) although cogmod_priors() covered it, so every family one function supports the other now does too. phi starts at 1, the Uniform at mu = 0.5, rather than the U-shaped Beta init = 0 gives it on its softplus link, pex at 0.2 rather than half of all responses on an endpoint, and mu and bex at 0.5. Against simulated data in eight response styles the start’s worst log-likelihood shortfall is 1.23 per observation, against 1.36 for init = 0 and a median of 2.54 for Stan’s random start. Fitted cold, it made no difference to sampling: warmup within 5% of either alternative’s, no divergences, the same posteriors (benchmarks/betagate_start/).

Performance

The classic cogmod_ddm() costs 0.77x what it did per gradient in CmdStan, and 0.83x in stanli, with the same log density and gradient to the last bit. “Classic” means sigmadrift, sigmabias and sigmandt all fixed at 0 in bf(). On every trial the Stan code chooses between two of Stan’s Wiener densities from an estimate of the log density. The estimate only decides a branch, yet autodiff recorded it like any other expression, and it was a fifth of the gradient. A few products now prove, for nearly every trial, that the estimate would have chosen the faster density, and the estimate runs only where they cannot. So every trial takes the branch it took before. Checked at 200,000 parameter values, and at 120 points of a 10-participant model in CmdStan and stanli. The ratios are from that model, on two runs; benchmarks/gradient_cost.R gives 0.77 too (benchmarks/ddm_router/). The cost of estimating sigmadrift, relative to the classic model, goes from about 2.8x to about 3.5x in ?cogmod_ddm, because the classic model got cheaper.

Documentation

New article, Comparison with Other Packages, setting cogmod beside brms’s native families and bmm.
?cogmod_priors said the RDM’s drift rates were left alone, which stopped being true when driftone got its prior; it now documents that prior, and its table of omitted parameters has the new rows. ?rcogmod_rdm now explains the driftone prior that the decision-making article sends readers to. ?rcogmod_exwald quoted the ndt prior as normal(-1.2, 0.2); it is normal(-1.2, 0.5).
Why the race families’ driftone priors are unchanged. Under the dec(Error) coding every example uses, driftone is the error accumulator, and wherever errors are rare its posterior sits in the lower tail of its prior - on 324,000 Illusion Game trials at prior quantiles down to 3e-6. ?cogmod_priors, ?rcogmod_rdm and ?rcogmod_lba2 now say that this is the fence working: below a drift of a few hundredths the likelihood is flat, and nothing observable depends on where on the plateau the drift sits. Priors centred on the error accumulator were measured before being turned down (benchmarks/driftone_priors/). On the decision-making article’s data they predicted the same error rates and RTs, to within 0.003 and 5 ms. On the RDM they cost 9-26% more gradient evaluations per draw. On the LBA2 they raised the divergent transitions from 1 in 12 chains to between 5 in 12 and 10 in 11.
Why cogmod_exwald()’s chains can disagree. The ex-Wald can make a right tail in two ways: from the Wald itself, at a low drift with a short tau, or from a long exponential stage behind a nearly Gaussian Wald, an ex-Gaussian with a shift. On the RT article’s data the two predict response times within 8 ms of each other from the 1st to the 95th percentile. Each is a separate mode of the posterior, more data separates them further rather than merging them, and a chain settles in whichever its warmup reaches. ?rcogmod_exwald and the RT article now say so.

The measurements are in benchmarks/exwald_modes/. The split the driftone_priors refits showed is one chain in the second solution, held there by a valley of 12 log-density units. The density is exact across the seam where it changes branch. On data simulated from the ex-Wald itself the second mode was present in every data set from 2,000 to 32,000 trials, and its valley grew with the data.
?cogmod_priors and ?cogmod_inits are rewritten to be read. ?cogmod_priors now has one table of every prior it sets - the parameter, its families, the prior in bf(), the prior when left out, and the prior on its slopes, which were documented for only some parameters - in place of four overlapping tables. It sorts the reasons into five kinds rather than a paragraph per family, each family’s own page keeping the full account. It also stops saying five things that were not true: it left cogmod_lba1() out of the families covered, sent cogmod_geg() readers to a section that did not exist, said shape existed only in cogmod_loggamma(), and said an omitted shape was left flat, where brms gives it gamma(0.01, 0.01), which cogmod_priors() replaces. With ?rcogmod_exgaussian, it also said brms’s default intercept for the ex-Gaussian’s mu was centred on zero seconds. It is centred on the median response; the trouble is its scale of 2.5 s, and under 0 + Intercept there is no default at all. Neither function returns anything different for a family it already covered.

cogmod 0.3.4

Breaking changes

cogmod_priors() now states the SD of a smooth term (sds) in link units, whatever the basis. One unit of sds moves the linear predictor by an amount set by the basis mgcv builds: about 0.45 link units for s(x), 6 for s(x, bs = "cr") and 0.1 for t2(x, z, k = c(5, 5), bs = c("cr", "cr")), almost regardless of the number of rows or the units of the covariates. The fixed exponential(1) of 0.3.3 therefore allowed a median of 4 link units of wiggle on the second of those and 0.07 on the third, and it applied only to the parameters a family names, the rest keeping brms’s equally basis-blind student_t(3, 0, 2.5). Each smooth term now gets exponential(rate) with the rate set to its own scale, read off brms::make_standata(), which is exponential(1) on the wiggle measured in link units. It applies to every smooth in the model, the response’s own included; the response’s slopes are still left to brms. On the tensor smooths of a 324,000-trial study the posterior sds ran from 0.4 to 37 (median 7) across four families, where exponential(1) puts 0.1% of its mass above 7. Any model with a smooth gets different priors from 0.3.3. cogmod_priors() now builds the basis of each distinct smooth term once (1.8 s for a t2() on 324,000 rows). See the Smooth terms section of ?cogmod_priors.
The unpenalised part of a smooth on a dpar the family names gets the dpar’s slope prior stated in link units. brms reports those coefficients (bs_*: the linear trend of an s(x), the marginal trends of a t2()) under class b, so until now they took the dpar’s blanket slope prior, normal(0, 0.2) on ndt. mgcv sets the scale of those columns, not the user: their root mean square is about 0.16 for s(x) whatever the units of x, 0.8 to 1 for a t2(), 2 for s(x, bs = "cr"). On the default s(x) that prior was normal(0, 0.03) on the trend at a typical row, which pinned ndt’s linear trend to about 0.1 link units across its range while the penalty left that part of the smooth free, and on s(x, bs = "cr") it allowed 0.4. Each such column now gets normal(0, s / rms), with s the dpar’s slope SD, as its own b row. The response’s keep brms’s default, as its slopes do. Any model with a smooth on ndt, poutlier or another family-named dpar gets different priors from 0.3.3.
In cogmod_lnr(), cogmod_rdm() and cogmod_lba2(), mu’s group-level SD gets the same prior as the other accumulator’s. In a race mu is the rate of accumulator 0, and the dec() coding decides which accumulator that is, so brms’s student_t(3, 0, 2.5) on one accumulator’s SD against exponential(1) on the other’s changed whenever the coding did. Both now get exponential(1). mu’s intercept and slopes are still left to brms.

Changes to starting values

cogmod_inits() starts a smooth’s SD at 0.1 link units of wiggle instead of sds = 0.05, for the reason above: 0.05 was 0.3 link units on s(x, bs = "cr") and 0.004 on the t2(), where the posterior sds sat five or six log units above the start. The smooth coefficients still start at zero, so the starting smooth stays flat.
cogmod_inits() scales the jitter on regression coefficients by how far their design columns reach. The noise that disperses the chains is 0.25 on every population-level parameter, argued as “a quarter of a link unit”. For a slope it reaches the linear predictor multiplied by its column, so on a column that spans +-3 it was three quarters, and the start tilted the whole surface. The unpenalised part of a t2(..., bs = c("cr", "cr")) smooth (bs_*) is such a block: on the Illusion Game’s Muller-Lyer models it tilted every smoothed parameter by a median of 2.3 link units across the rows, and 30% of chains started with more than 1% of their responses below their own ndt. Each block’s jitter SD is now divided by the largest row norm of its design wherever that exceeds 1, which keeps the SD of what reaches the predictor at or under 0.25 link units at every row (the noise is Normal, so a single draw can still exceed that): a median tilt of 0.5 there, and 0.7% of chains past 1%. A dummy column or a standardised covariate keeps its 0.25; a slope on an unstandardised one (age in years) no longer has an SD of six link units at the extremes.
cogmod_lba2()’s error accumulator starts slower than the correct one, at driftone = 1 against mu = 3, as cogmod_rdm()’s has since 0.3.2. The old equal start was a 50/50 race. On simulated data with error rates from 3% to 43% it was the costliest of the starts 3, 2, 1 and 0 every time, and moving it to 1 took 2,400-3,800 log-likelihood units per 2,000 trials off the shortfall of the whole start. Starting values do not change a posterior, so no fitted model needs revisiting for this.

Documentation

?rcogmod_lba2 now opens its Negative drift rates section with the advice the rest of it explains: wherever one response is rare in some condition, keep the drift SD of its accumulator free of predictors (sigmaone ~ 1, or a participant intercept only). A condition-varying SD gives every such condition a direction of its own along which drift and SD slide together without changing the fit.
The Performance article has a new, explicitly experimental section on sampling a cogmod model with stanli, a Stan interpreter that needs no C++ compilation, and turning the result into a brmsfit through brms’s "mock" backend. On the LNR it gave CmdStan’s posterior in 0.74x CmdStan’s wall time on Linux and under half of it on Windows (benchmarks/stanli/).

cogmod 0.3.3

CRAN release: 2026-09-25

New features

A distributional parameter pinned at a boundary is now checked against the response. The three bounded-scale families are mixtures of a continuous part and one or more point masses, and every weight involved can be fixed in bf(). Fixing one at a boundary switches its component off, so any observation belonging to that component has zero probability and the log-likelihood is -Inf everywhere - which CmdStan reports only as Initialization failed after 100 attempts, naming nothing. cogmod_priors() now errors first, saying how many rows are involved and how to proceed. Eleven combinations are covered: for cogmod_choco(), pmid = 0 with a response exactly at the midpoint, pmid = 1 with one anywhere else, pex = 0 with an exact 0 or 1, bex = 0 with a 1 and bex = 1 with a 0; for cogmod_betagate() the same pex = 0 and bex cases plus pex = 1 with any interior response; and for cogmod_betadiscrete(), pzero = 0 with a 0 and pzero = 1 with any rating. The first and the last are the ones the documentation invites - ?rcogmod_betadiscrete offers pzero = 0 as the way to say a scale has no zero category, and ?cogmod_priors suggests fixing pmid or pzero at 0 to switch the parameter off.
cogmod_priors() now covers cogmod_choco(), cogmod_betagate() and cogmod_betadiscrete(). Every distributional parameter of the three bounded-scale rating families previously arrived either flat - improper - or with a brms default aimed at a different parameterization. Two of those were worse than unhelpful. The point-mass probabilities pmid and pzero were flat on a logit link, which is improper in the posterior as well as the prior whenever the event is absent from the data: with no exact midpoints anywhere, the likelihood in pmid increases monotonically all the way to zero and nothing stops the logit running to minus infinity. That is poutlier’s failure, and they now get poutlier’s treatment, including an omitted form whose mode is at zero. Measured on 400 slider responses with no midpoints and no extremes: the flat defaults put pex at -1.1e14 and pmid at -3.4e13, Rhat 2.9, ESS 5, with 12% divergent transitions; these priors give -4.70 and -5.15, Rhat 1.00, ESS 3000-4000, no divergences. phi arrived with student_t(3, 0, 2.5), which brms supplies because it recognises the name from its own beta family - the same trap cogmod_exgaussian()’s sigma falls into; on cogmod_betadiscrete()’s log link that prior’s 95% interval runs to phi = 2853, where the Beta has collapsed onto a single rating category, so it is overridden rather than filled. The Beta precisions are fenced away from both ends: below about 1 the underlying Beta is U-shaped and unbounded, and a large one narrows it sharply. mu is left to brms in all three, since on a logit link its student_t(3, 0, 2.5) is the standard weakly informative choice. See ?cogmod_priors for the full table and the reasoning behind each number.
cogmod_inits() now covers cogmod_choco() and cogmod_betadiscrete(). Both are bounded-scale families for subjective ratings, and neither has a flat region of the kind that makes an init mandatory for the RT families - every start is a proper density and a chain begun at the default does move. What it has to move away from is a description of a rating scale that no data set matches. The logit origin puts pmid and pzero at 0.5, i.e. half of every response exactly on the midpoint or outside the scale altogether; they now start at 0.05. cogmod_choco()’s two Beta precisions are behind a softplus link, so softplus(0) = 0.69 puts both Beta shapes below 1 - a U-shaped rating distribution, unbounded at both ends of each half of the scale; they now start at 2, in the middle of a plateau where the fit is insensitive to the exact value. pex starts at 0.1 rather than 0.5. On a 2000-trial slider data set the old start sat about 1900 log-likelihood units from the new one, 983 of them pmid alone. For cogmod_betadiscrete(), phi = 1 with mu = 0.5 is the discrete Uniform, which is where the origin already is - it is named so that init = "random", which draws phi anywhere from 0.14 to 7.4, no longer decides it. So this buys warmup rather than a fit that would otherwise fail.
cogmod_inits() lists cogmod_geg() among the families it supports. It had targets for it but was not naming it, so the error message for an unsupported family said cogmod_geg() was one. The supported list is now derived from the targets rather than kept beside them.
cogmod_inits() jitters the hierarchical blocks a fifth as much as the population-level ones, and starts smooths flat. One jitter (default 0.25) now applies to intercepts and slopes; the standardized group-level effects z_*, the smooth coefficients zs_* and their scales sd_* and sds_* get a fifth of it (0.05), and sds_* starts at 0.05 rather than the generic 0.25. Two numbers set the two tiers directly. The reason is that a unit of noise on a group effect or a spline coefficient is not a unit on the linear predictor: it is multiplied by a scale and a design column - tensor basis values reach tens - once per participant or basis function. On a production model with a tensor smooth on five distributional parameters and a participant intercept on six, the old jitter moved linear predictors by one to two and a half link units at some rows, which started one participant’s non-decision time above 98 of their 128 trials and a sigma at 0.07 s where 0.5 was intended, all of it then explained by the outlier component; at 0.05 on those blocks the same model started where the targets say. Chains still start apart where it matters for Rhat, on the intercepts and slopes. Warm starts (cogmod_warmstart(), warmstart =) use the same rule at their smaller default.
The 7-parameter cogmod_ddm() is 1.7x cheaper per gradient when both sigmabias and sigmandt are estimated. Stan’s 7-parameter wiener_lpdf() integrates the start-point and non-decision-time ranges out by adaptive cubature - eight integrals per observation, one for the density and one per partial derivative - and takes a tolerance for them that the Stan code never passed, so it ran at Stan’s default of 1e-4. It now passes 1e-3 (.DDM_WIENER_PRECISION). Measured per observation and gradient on CmdStan 2.38.0: with both ranges open, 274 us falls to 163; with one open the cost does not move (102 to 100 us), because a one-dimensional integral is already settled by cubature’s first pass, and nothing much happens past 1e-3 either (153 us at 1e-2), which is why it stops there. On the intercept-only program benchmarks/gradient_cost.R times, 5000 trials with all seven parameters estimated, a gradient went from 3.25 s to 1.50 s (ratio 0.46, median of 21 alternating blocks). The price: the Stan density still agrees with the R one to 1e-5 over the test grid, but its gradient sits about 2e-5 relative from central differences at a typical point where it sat 2e-6 before, and 1.2e-4 with ndt pushed 2.5 log units above its start - a point no chain visits after warmup, and small next to the leapfrog error a sampler already carries, but a change of an order of magnitude, so it is recorded with the constant. Any value between 3e-4 and 1e-3 measures the same gradient error and 1e-3 is the cheaper. A paired fit settles it (benchmarks/ddm_precision/): the 7-parameter model on 100 trials, both tolerances from one start, metric and seed, end to end and again with adaptation held fixed - the same 24 leapfrog steps per iteration, the same step size, acceptance 0.90 against 0.92, fewer divergences (7 against 10 and 11), min bulk ESS 248 against 206 and 261 against 224, at half the cost per gradient (445 against 883 ms): 2.2 to 2.4 times the effective samples per CPU second. That fit also puts the standalone figures above in perspective: at its posterior a gradient cost 8.3 ms per observation at Stan’s default tolerance, thirty times the 274 us of the benign benchmark point. The same cost measurements are now in ?cogmod_ddm, under the between-trial variability section: estimating sigmadrift costs about 2.8x the classic model per gradient, estimating one of sigmabias / sigmandt about 18x and both about 30x, and the fast path is a test for exactly zero, so a tight prior does not buy it back - only sigmandt = 0 in bf() does. Nothing user-facing said so before, and it is the difference between a day and weeks on a large data set.

Bug fixes

cogmod_priors() now sets the wiggliness prior (sds) of a smooth on a distributional parameter. The table in ?cogmod_priors has always listed sds alongside sd at exponential(1), but brms fills a smooth’s blanket sds row itself and leaves the per-term rows empty, so a function that fills only what arrives empty never reached it: ndt ~ s(x) kept student_t(3, 0, 2.5) on its smooth’s scale. On a log or logit link that is a half-t with median 1.9 in link units - loose enough for the smooth alone to move a sigmandt by a factor of seven or walk a sigmabias across its range, behind an intercept the family had deliberately fenced. The blanket row is now replaced, for every family and every dpar cogmod_priors() handles; the response’s own smooth is left to brms, as its slopes are. Group-level sd rows were never affected.
cogmod_rdm()’s gradient is now exact. Every normal tail in the RDM’s Stan code went through Stan’s std_normal_lcdf(), whose value is right but whose partial derivatives are an approximation, and in a race those partials are the gradient of the drifts and the boundary. Measured against central differences of the log probability, the gradient sat 2e-4 relative from the truth at a typical start and as far as 7e-2 where the drift is small, while the log probability itself was smooth to 1e-7. HMC stays exact under an inexact gradient - the accept step corrects for it - but pays in step size and acceptance. All eleven calls now go through cogmod_log_Phi(), the function introduced below for the LNR, which brings the worst error over the same grid to 2e-7 at about 14% more per gradient evaluation. The function moved into a prelude of its own (.LOG_PHI_STAN_PRELUDE) so that any family can take it; the LogNormal’s and the RDM’s preludes both start with it.
cogmod_invgaussian()’s gradient is now exact. The same defect as the RDM’s, in the family’s own Stan code: every normal tail of the Wald - the two terms of its CDF and survival, the two integrated ones a non-decision range needs, and the truncation factors that make sigmadrift a truncated normal drift - went through std_normal_lcdf(). sigmadrift is differentiated almost entirely through those truncation factors, and its partial sat 1.2e-3 relative from central differences of the log probability, mu’s 2.7e-4, against 1e-8 for the families already on cogmod_log_Phi(). Six of the eighteen points the gradient check walks failed on it. All of them now go through cogmod_log_Phi(), which brings the worst error over the same grid to 5e-8; the truncation factors lose a round trip through log1m_exp() in the bargain, since log1m_exp(std_normal_lcdf(-z)) is cogmod_log_Phi(z) written the long way. cogmod_exwald() builds on the same prelude and gets the same fix. Values are unchanged.
cogmod_geg()’s CDF term is no longer Stan’s exp_mod_normal_lcdf(). The alpha-power construction puts (shape - 1) * log F_EG inside the density, so the ex-Gaussian CDF is differentiated on every evaluation, not only when a response is censored. The built-in agrees with the R side to 1e-8 over the range the tests walk, but outside it - at sigma / tau around 9, where the two terms of F_EG cancel hardest - both its value and its partials go wrong: Stan’s own central differences put d/d sigma at 28.5 where its autodiff said 77.5, and d/d tau at 113.6 against 53.0. The Stan side now subtracts the two terms in log space itself, through cogmod_log_Phi(), exactly as .lcdf_exgaussian() always has in R; over the same grid the worst gradient error is 2e-9, and the log CDF matches quadrature of the density to 5e-15. It costs about 10% per gradient evaluation (3.1 ms to 3.5 ms over 5000 trials, median of 21 alternating blocks). Two things kept it to that. The ex-Gaussian density is the CDF’s second term over tau, so the GEG writes both out from one pair of normal tails rather than calling for each; and the pair is written inline rather than returned from a helper - a vector[2] return measured half again the cost of the inline form, through the autodiff-stack allocation it makes on every call.
cogmod_exgaussian()’s Stan density, CDF and survival go through cogmod_log_Phi(). The density is the same expression as before - exp_mod_normal_lpdf() and dcogmod_exgaussian() in R evaluate it too, and log-probabilities are unchanged to every digit that matters - but it is now written out over cogmod_log_Phi(), so it stays finite about 38 standardized units into the left tail where the built-in’s bare erfc underflows, and it shares its normal tail with the CDF that cogmod_geg() is built on. The cens() path gets the accurate CDF with it: cogmod_exgaussian_lcdf() was exp_mod_normal_lcdf() and had the same wrong partials as the GEG’s, and cogmod_exgaussian_lccdf() was on std_normal_lcdf(). About 10% per gradient evaluation, the cost of cogmod_log_Phi() over the built-in.
cogmod_lnr() and cogmod_lognormal() no longer hand Stan a non-finite gradient in the tails. cogmod_lognormal_ldiff_Phi() formed log(Phi(y + c) - Phi(y)) from erfc as log(u1) + log1m(u2 / u1), and the A == 0 branch of cogmod_lognormal_acc_ltails() called lognormal_lcdf() and lognormal_lccdf(), which are erfc alone. erfc underflows near x = -38 - about 38 standardized log units from an accumulator’s median finishing time - and past that the value is log(0) and the partials are inf or 0 / 0. The outlier mixture then hides it: log_mix() stays finite with one component at -inf, but reverse mode multiplies the zero adjoint into the stored partial and 0 * inf is NaN, so one response in a data set turned the gradient of the whole model to NaN. That reads as Gradient evaluated at the initial value is not finite at the start of a fit and as divergent transitions afterwards. Both now go through one function, cogmod_log_Phi(): erfc in the body of the distribution and, below x = -25, the asymptotic expansion of the tail, whose leading term is the exponent itself, so nothing underflows and the result stays finite and differentiable as far as x = -1e150; its six terms agree with R’s pnorm(log.p = TRUE) to 4e-16 relative, so the two branches meet with no step in the density. cogmod_lognormal_ldiff_Phi() now takes the difference in logs with log1m_exp() rather than as a quotient of two minute numbers. Measured over a grid of decision times from 1 ms to 300 s and sigmas from 0.02 to 1.2, with and without a start-point range: 6 of 72 gradients were non-finite before and none are now, and the densities agree with the R kernels everywhere, with no value newly truncated to -inf.

std_normal_lcdf() is not used for this, although it has the range - its value is exact against pnorm(log.p = TRUE) as far as x = -1e7. Its analytic partials are not: on 20000 responses they sat 1.7e-3 from central differences of the log probability where the erfc route sat 4e-6, and in a race those partials are the gradient of nu and sigma.

Where the old code’s gradient was finite it was not always right. On 20000 responses with a start-point range it sat 8.3e-4 from central differences and the new one sits 9.6e-7, which is the finite differences’ own noise; without a start-point range both sit at 4.1e-6, so nothing there was given up for it.

Sampling is not slower in the case most models are in. cogmod_lognormal_acc_logcdf() and cogmod_lognormal_acc_logsurv() now take the single tail they were asked for when sigmabias = 0, instead of building the pair and discarding one: the two share no work there, and cogmod_lnr() reads the survival alone. On 20000 responses a gradient of the plain LNR came out about 20% cheaper than before and one with a start-point range about 15% dearer, the latter buying the corrected gradient above.

Documentation

Help pages that sent the reader to vignette("rt_models") or vignette("performance") now link to those articles on the package website. Both are website articles, not installed vignettes, so the call failed.
The cogmod_warmstart() example is now self-contained and runs (with CmdStan): it simulates its data, fits the pilot and the warm-started model, and writes the warm-start table to a temporary file rather than to the working directory.

cogmod 0.3.2

New features

cogmod_lognormal() gains sigmabias, the same between-trial start-point range as the LNR’s below: the decision time is the LogNormal multiplied by a Uniform(1, 1 + sigmabias) distance, which is the single-accumulator LBA with a LogNormal drift rate and its threshold offset pinned at 1. At sigmabias = 0 - the default of rcogmod_lognormal(), dcogmod_lognormal() and pcogmod_lognormal(), placed after ndt and before poutlier as the Wald’s sigmadrift is, and the value to fix in the formula unless the design speaks to start-point variability - the family is the shifted LogNormal exactly as before, bit for bit and at the same cost. The CDF and survival cens() needs, the mean posterior_epred() reports (which gains a factor 1 + sigmabias / 2) and the Stan functions all carry the range, and agree with quadrature over the start point to 1e-7 and with each other to 1e-10. The kernels are shared with cogmod_lnr(), which is now a race of two of these accumulators in code as well as in theory. As for the LNR, a formula that omits sigmabias estimates it, cogmod_priors() fences the flat direction at zero, and fits made with earlier versions have to be refit; the two vignette LogNormal models were. cogmod_logstudent() and cogmod_loggamma() do not get the parameter: the density needs a partial first moment of the rate distribution, which does not exist for a Student-t on the log scale and needs incomplete gamma functions for the log-Gamma.
cogmod_lnr() gains sigmabias, a between-trial start-point range: each accumulator now starts at Uniform(0, sigmabias) and runs to a threshold 1 + sigmabias at its LogNormal rate, so its finishing time is the distance divided by the rate rather than the reciprocal of the rate alone. At sigmabias = 0 - the default of rcogmod_lnr() and dcogmod_lnr(), and the value to fix in the formula unless the design can identify a start-point range - the family is the LNR exactly as before, bit for bit and at the same cost, since both the R and the Stan kernels take the plain lognormal branch there. Above zero it is the LBA with LogNormal drift rates (Heathcote & Love, 2012), the model the LNR was introduced as a limit of: with a LogNormal distance as well as a LogNormal rate the two fold into one sigma, which is why the LNR never had a start-point parameter, whereas a Uniform distance leaves a shape the rate alone cannot produce. The threshold offset is pinned at 1 rather than sigma at 1 because rescaling the evidence axis shifts nu and scales the range and the threshold but leaves a LogNormal rate’s sigma untouched, so sigma cannot pin the scale; sigmabias is therefore read in units of the threshold offset. The density is a difference of two normal CDFs and the survival one more, evaluated from a series below a start-point range of about 1e-4 * sigma and from ratios of log-CDFs in the tails, and both agree with one-dimensional quadrature over the start point to 1e-7. cogmod_priors() fences the flat direction at zero the way it does for cogmod_lba1()’s sigmabias, and cogmod_inits() starts it at 0.5. A formula that omits sigmabias now estimates it, as brms does with any dpar; fits made with earlier versions cannot be post-processed, because their family carries no sigmabias, and have to be refit. In rcogmod_lnr() and dcogmod_lnr() the argument sits between ndt (or response) and poutlier, where the DDM keeps its between-trial variabilities, so a call that passed poutlier by position needs it named. The vignette LNR model pins it at zero and was refit.
cogmod_priors() gains warmstart, which re-centres the priors on a previous fit: every population-level intercept and coefficient, every group-level SD, and every dpar left out of the formula gets normal(median, prior_scale * sd) from that fit’s posterior, matched parameter by parameter on the class, dpar, coef and group a get_prior() row carries. No transformation is involved - the parameter a prior row is about is the parameter the source sampled, the Intercept prior being stated on the centred intercept in both - so the scales line up by construction. cogmod_warmstart() extracts the posterior median and SD alongside the means for this, as.data.frame() carries them in two new columns, and a table written before they existed still reads (the priors are then left alone, with a message). The correlations keep their LKJ and the standardized effects have no stated prior to change. Unlike the rest of the warm start, this changes the posterior, and if the source was fitted to data the new model also contains it double-counts it - a pilot on half the participants used to centre the priors for the fit on all of them uses that half twice. prior_scale, 3 by default, is what stands between a prior that only says roughly where the parameter lives and one that is the source’s posterior outright; ?cogmod_priors says when the argument is and is not legitimate.

Bug fixes

cogmod_rdm() no longer freezes one chain in four on a cold start. cogmod_inits() started both drifts at 3; the error accumulator now starts at 1. The failure looked like a stuck chain - every transition at the maximum treedepth, step size a thousand times smaller than the other chains’, Rhat 1.5 to 2.9 - and it was traced to the very first warmup transition. A Wald density is thin on the fast side and flat on the slow side, so on data whose error drift is about 0.2 (the speed-accuracy data of vignette("performance")) the old start sat 400 log-density units above the posterior. The first trajectory converted that into momentum along the flat driftone direction (the plateau cogmod_priors() fences, where the likelihood no longer changes) and carried the chain from a link value of +3 to -20 in eight leapfrog steps; the step size then collapsed to 1e-5 within a dozen iterations and the metric windows that followed were estimated from a chain that no longer moved. Stan’s model methods found nothing numerical at the frozen position - the log-density and its gradient are smooth across the driftless branch of the survival function, and finite differences agree with autodiff - so the fix is where the chain starts, not what it computes. On an 800-trial mixed model with a 200-iteration warmup the old start froze a chain in 1 run in 4 with Stan’s default metric and in 5 of 6 when handed a metric adapted to the bulk (4 of 4 when the benchmark cell was rerun); the new start has done so in 0 of 6 under the latter, the harder case, and each such fit ran in 3 minutes instead of 30. When the error accumulator really is as fast as the correct one the start is off by a factor of three on the cheap side, which costs a few dozen units and changes nothing.
cogmod_inits() now starts ndt at half the first percentile of the observed response times instead of a fixed 0.1 s. Same mechanism as the previous item, other cold start. The fixed value was a third of the prior median and safely below any ordinary data, which was the whole argument for it; on data whose non-decision time is 0.6 s it sat half a second low, every decision time looked far too long, a driftless race then fit better than a fast one, and the first trajectory of a cold cogmod_rdm() chain threw both drifts onto their flat regions - three chains in four on the benchmark’s shifted target when handed a metric adapted to the bulk. Half the first percentile is still below essentially every response, so the gradient the small start was protecting is intact, and it follows the scale of the data. Every ndt + poutlier family gets it. Starting values do not change a posterior, so no fitted model needs revisiting for this.

Breaking changes

The default ndt prior is wider: normal(-1.2, 0.5) on the log scale, lognormal(-1.2, 0.5) for an omitted ndt, in place of the 0.2 SD. The centre is unchanged at 0.30 s; 95% of the mass now sits between about 0.11 and 0.80 s instead of 0.20 to 0.44 s. The old SD put a 0.6 s non-decision time - not unusual for older participants or more demanding responses - 3.5 SDs from the centre, and on the warm-start benchmark’s shifted target (vignette("performance")), whose non-decision time is about 0.6 s, cogmod_rdm() produced divergent transitions in every run with that prior, the 600-warmup reference included; the wider prior removed them (Rhat 1.06 and a minimum ESS of 67 became 1.01 and 392 on the same seed). The prior’s job is to fence the ndt -> 0 direction, where the likelihood goes flat and a flat prior would make the posterior improper; at 0.5 it still does (0.01 s is 6.8 SDs out) without telling the data where in 0.1 to 0.8 s the non-decision time is. Fits with the default priors will move slightly, most where the data put ndt far from 0.30 s, which is where they should have been free to move. Every ndt + poutlier family is affected.

cogmod 0.3.1

New features

New cogmod_warmstart() turns what a previous fit’s warmup produced - the adapted inverse metric, the step size and the posterior means - into the init, inv_metric and step_size arguments of a new brm() call, so that a refit, or the same model on more participants, can run a much shorter warmup. Stan adapts one variance per unconstrained parameter, and the function labels each with its Stan name (read off the generated program with the parser cogmod_inits() already has) and joins the two models on those names: population-level entries carry over, a pilot participant’s standardized effects follow it by level name to its position in the bigger model, new participants take their effect’s average variance and start at zero, and anything without a counterpart gets Stan’s default variance and a generic start, with a count in print(). The standardized effects and Cholesky factors that brms drops from a saved fit are rebuilt from the r_, sd_ and cor_ it keeps. as.data.frame() gives a table of a few kilobytes that survives write.csv() and can be passed back as a file path, so a pilot fitted on a laptop can warm-start an array job on a cluster. On the other side, cogmod_inits() gains a warmstart argument and the new cogmod_inv_metric() and cogmod_step_size() share its signature - the model’s formula and data, then the source, be it a fit, a table or a file - so that each argument of brm() has one helper and the table is mapped onto the model on the way. Whatever of formula and data is not given is taken from the source fit. Works for any brms model fitted with the cmdstanr backend and the diagonal metric. On a mixed LNR and a mixed DDM, a pilot on 4 of 8 participants warm-started the full fit to about twice the effective draws per second of a cold start with the full warmup, and four to six times those of a cold start with the same short warmup; the starting values alone bought nothing, so the metric and step size are the product (vignette("performance")).
cogmod_invgaussian() gains sigmandt, the between-trial range of the non-decision time (st0): each trial’s non-decision time is drawn from Uniform(ndt, ndt + sigmandt), so ndt becomes its lower bound, exactly as in cogmod_ddm(). Spreading the shift turns the Wald density into a difference of two CDFs and the CDF into a difference of two integrated CDFs, both closed form at a fixed drift, so the parameter costs a few normal CDFs per observation, works with cens() unchanged, and rides the existing drift quadrature when sigmadrift > 0 too. It is on a log link with cogmod_ddm()’s prior for the same quantity. It is hard to estimate and should be fixed at zero for most applications (sigmandt = 0 in bf()); it shares the leading edge of the distribution with ndt and poutlier, and should only be freed with a lot of data, a strong prior, or both. As with sigmadrift, leaving it out of bf() estimates it, so existing Wald formulas that do not mention it now fit one more parameter unless they add sigmandt = 0, and fits made before this version cannot be post-processed with it; the vignette models were refit. rcogmod_invgaussian(), dcogmod_invgaussian() and pcogmod_invgaussian() take sigmandt right after sigmadrift, so a poutlier passed by position moves along one.
Censored reaction times: brms’s cens() works on the RT-only families. bf(rt | cens(error) ~ ...) scores an error trial as a right-censored correct response: its RT is a lower bound on when the correct process would have finished, so it contributes that process’s survival rather than its density. On cogmod_invgaussian() this is the simple censored shifted Wald of Miller et al. (2018), their Eq. 4, the version = "simple" of the cswald model in bmm. It is not their competing-risks variant (Eq. 5, a race of two Wald accumulators with drifts v and -v, as implemented in rtdists and bmm’s version = "crisk"), which is a choice model rather than a censoring construction; cogmod_ddm() with bias fixed at 0.5 covers that ground. Here censoring is not a family but a construction, so the same formula works on every RT-only family with a closed-form CDF: cogmod_lognormal(), cogmod_logstudent(), cogmod_gamma(), cogmod_invgamma(), cogmod_weibull(), cogmod_invweibull(), cogmod_logweibull(), cogmod_bisa(), cogmod_exgaussian() and cogmod_geg(). Left- and interval-censoring come with it, and log_lik() - hence loo() - honours all three, which brms leaves to a custom family’s own method. posterior_predict() predicts the latent, uncensored RT, as brms does for its own families. cogmod_priors() and cogmod_stanvars() refuse cens() on the families that cannot take it, and cogmod_priors() warns above 20% censored trials - a threshold that is exact for timeouts and omissions at any rate but lenient for commission errors, where the construction is biased well before it (see ?rcogmod_invgaussian). What the model is for, what it assumes and the one check to run before using it are in ?rcogmod_invgaussian and the Censored Shifted Wald section of vignette("rt_models").

Under the hood every censorable family gets a <family>_lcdf and a <family>_lccdf beside its _lpdf, generated from two new registry slots so a family cannot drift out of step with itself. The survivals are written as survivals - never as log(1 - exp(lcdf)) - and the half Normal outlier’s through std_normal_lcdf(-z) rather than std_normal_lccdf(z), which is -inf from 1.66 s on: the two places bmm’s implementation broke. With sigmadrift > 0 the Wald CDF has no closed form and is taken by 64-point Gauss-Legendre quadrature over the drift, in R and Stan alike off one node table.
pcogmod_*() for every censorable family. pcogmod_lognormal(), pcogmod_logstudent(), pcogmod_gamma(), pcogmod_invgamma(), pcogmod_weibull(), pcogmod_invweibull(), pcogmod_logweibull(), pcogmod_bisa() and pcogmod_exgaussian() join pcogmod_invgaussian(), with lower.tail and log.p. The upper tail is computed as the upper tail rather than as 1 - CDF; these are the R side of the Stan _lcdf/_lccdf pair, and the tests hold the two to each other.
cogmod_priors() now checks the response before returning. Everything in this package is stated in seconds and none of it is unit-equivariant - the ndt prior means 170-300 ms, .POUTLIER_SCALE is 0.2 s - but brms fills its own defaults from the data, so a column of milliseconds produces a model whose two halves silently describe different quantities. It compiles, it samples, it converges, and the estimates are meaningless. The check catches that and the handful of other mistakes with the same character.

It stops where the offending rows would make the fit impossible or wrong in a way Stan cannot report: a non-positive reaction time under a family that places no density below ndt; a response outside [0, 1] for cogmod_choco() or cogmod_betagate(); a non-integer rating for cogmod_betadiscrete(); a non-numeric response; and a third level in dec(), which the choice families would otherwise fold silently into option 1, since their Stan code tests dec == 0 and takes the else branch for everything else.

It warns about the rest: a median implying milliseconds, NAs, and either tail running past what poutlier can absorb. The tails are judged as proportions rather than counts, because the outlier component is supposed to produce the occasional fast response - rcogmod_lognormal(200, ndt = 0.2, poutlier = 0.02) puts one at 81 ms - and a count-based test fires on the package’s own generator. Over 20000 draws the component sends 0.8% of responses below 0.1 s at poutlier = 0.02 and 1.9% at 0.05, the top of the default prior, so the warning sits at 5%.

Families with neither ndt nor poutlier - cogmod_exgaussian() and cogmod_geg() - are exempted from the tail checks, and a non-positive response is a warning rather than an error there, their support being the whole real line. A formula or family the check cannot read is passed through untouched, so brms’s own error is what the user sees.

Bug fixes

dcogmod_ddm() and everything built on it (log_lik(), loo(), p_outlier(), the other R-side post-processing of cogmod_ddm() fits) are now accurate when the sigmandt range reaches down to fast decision times. The R density integrates the non-decision time out with a fixed 25-node Gauss-Legendre rule, and when that range covers decision times from about zero up to the response - a fast response, or a wide sigmandt - the integrand holds the whole early peak of the first-passage density inside a sliver of it, which 25 nodes on the plain time scale cannot resolve. It is the defect reported against rtdists::ddiffusion() in rtdists issue 28, and it was here too: on the issue’s own example (boundary = 0.5, drift = 0.5, bias = 0.3, sigmandt = 0.16, decision time 0.16) the density was out by 5e-4, by 3% with bias = 0.1 and sigmandt = 0.2, and by 50% with the start point almost on the responding boundary. The rule now runs over log decision time, from the point where the density is dead rather than from zero, so the peak is about one log unit wide wherever it sits and the same 25 nodes resolve it at any time scale, and the rule takes more nodes only when the log range is wide enough to need them: the worst error over the issue’s sweep is now below 1e-9, and below 4e-12 over a much broader grid, against a converged 1600-node rule. Densities with sigmandt = 0 are unchanged to the last bit, and the common case - sigmandt well inside the response time, where the old rule was already accurate - gives the same values at the same cost. The Stan likelihood uses Stan’s own adaptive wiener_lpdf() and was never affected, so fitted models are unchanged; only their R-side post-processing moves, and only in that regime. The R-versus-Stan density test is tightened from a relative 1e-4 to 1e-5 accordingly, which is Stan’s own tolerance.
cogmod_rdm() no longer produces divergent transitions by the hundred on healthy posteriors. The Stan log-survival of the losing accumulator formed its reflection term as log_diff_exp(log R(b), log R(k)), and for a response less than about half a millisecond above the non-decision time both normal CDFs in R round to exactly 1, so it evaluated log_diff_exp(0, 0). The value is fine (-Inf for a term that really is negligible there, which is why the R-versus-Stan density tests never caught it), but its reverse-mode adjoint is 0 / 0, and that NaN propagated into the gradient of every parameter. Stan reports a NaN gradient as a divergent transition, and because ndt is estimated a few milliseconds below the fastest responses, most trajectories crossed one of those windows: on the lexical decision data of the decision-making article, 900 trials from 6 participants, between a third and two thirds of the transitions were divergent - with population-level effects only or with participant intercepts, under diag_e or dense_e - while Rhat and the effective sample sizes said the posterior was fine, because it was. The difference is now assembled from three pieces that each stay away from the saturated end of the normal CDF, at the cost of two extra normal CDFs on early responses only. Values are unchanged to 1e-12 on the log scale; the same fits now run without a divergence (population-level) or with the handful the other families also show under dense_e with random effects. The other race families are unaffected. A gradient regression test guards it, gated behind COGMOD_TEST_SLOW like the other tests that compile a model of their own.
rcogmod_ddm() no longer returns the odd 10-40 s response in place of a fast one. The sampler inverts the CDF with Newton’s method, and took a step under 1e-10 in log time as convergence. A pass that lands far in the tail finds the survival and the density both denormal, and their ratio makes the step look tiny while the residual is still hundreds of log units off, so the draw was accepted where it stood. It hit about one draw in 4,000 to 20,000 at short boundaries or strong drifts (boundary = 0.3, bias = 0.3, drift = -5 is one such cell; boundary = 2, bias = 0.7, drift = 6 another), which is rare in rcogmod_ddm() and a handful of absurd draws per observation in posterior_predict(), where every observation gets thousands. Convergence now also requires a small residual, and a test pushes every draw back through the CDF to check it lands on its own quantile.

Performance

rcogmod_ddm() is 3x faster, and posterior_predict() on a cogmod_ddm() model up to 16x. The sampler inverts a series whose length is set by the fastest response it could be asked for, and used that length for every draw: 41 terms at the default start point where the median draw needs 6, 205 at a start point of 0.1. It now runs in stages - 16 terms settle the bulk of the draws, and only the responses too fast for that many go round again with four times as many, until the full series is reached. Nothing is approximated: a draw is only accepted from a stage whose series is exact at its root, and the draws agree with the CDF to 1e-12. The gain is largest where the parameters vary across draws, as they do in posterior_predict(), because one extreme posterior draw used to set the series length for all of them. Converged draws now also drop out of the Newton iteration, and a fast response starts from the single-barrier small-time approximation rather than from the floor of the bracket.
The choice families’ posterior_predict_*() methods take a vector of observations, returning the draws stacked with those for i[1] first. brms::posterior_predict() calls the method once per observation, and with a few dozen draws per call about half of each call is fixed cost and the loop adds as much again; predicting in chunks of ~50 observations from a prepared brmsprep instead runs a posterior predictive check on 2,500 DDM trials in about a third of the time. The recipe is in ?posterior_predict_cogmod_ddm. The DDM sampler’s own fixed cost per call is also down by about 15%, from indexed assignment in place of ifelse() and no column copies while every draw is still active; the draws are bit-identical.
The R-side DDM density no longer goes through brms::dwiener(). dcogmod_ddm(), and with it log_lik(), loo(), p_outlier() and every other post-processing method that evaluates the likelihood in R, now use a vectorised Navarro and Fuss (2009) series written in log space. The 4-parameter density is about eight times cheaper per element and agrees with brms::dwiener() to 1e-12 on the log scale. The 7-parameter density, which evaluates that series 625 times per observation under Gauss-Legendre quadrature, goes from about 5 ms to about 0.6 ms per draw-observation - a LOO over 4000 draws of 500 trials drops from close to three hours to about twenty minutes. Both now return a finite log-density in the far tails where dwiener() returns log(0). The Stan likelihood is unchanged. RWiener is still needed by the test suite, which uses dwiener() as the reference.

Documentation

The performance article is reorganised from the suggestions with no downside to the ones that need judgement, and gains four sections. Compiler optimizations: stanc’s O1 and CmdStan’s STAN_CPP_OPTIMS and STAN_NO_RANGE_CHECKS, passed through stan_model_args, and what each one does. Mass matrix adaptation (metric = "dense_e"): why the boundary/ndt and drift/boundary trade-offs of evidence accumulation models make the default diagonal metric a poor fit, what the dense metric costs as the number of parameters grows, and how to pass it through either backend. Warm starts: reusing the adapted metric and step size that brms keeps in a fit’s metadata to shorten the warmup of a refit; carrying a pilot fit’s metric, step size and posterior means over to the same model on more participants, which has more parameters, by mapping the metric across by parameter name (about twice the effective draws per second of a cold start on a mixed LNR and a mixed DDM, where the pilot’s initial values alone bought nothing); and a cmdstanr-level pipeline that initializes MCMC from Pathfinder draws and their unconstrained covariance, then wraps the result back into a brmsfit, with the reasons never to fix the metric to a variational approximation. The approximation section now also covers the Laplace approximation (algorithm = "laplace") and how it compares with Pathfinder. Each section reports what the option bought on the DDM, LBA, LNR and RDM in a local benchmark; the scripts behind those numbers live in benchmarks/ (not part of the installed package) and can be rerun on any model.

Breaking changes

cogmod_lba2() now truncates each drift rate at zero, the convention of rtdists (posdrift = TRUE), DMC, EMC2 and ggdmc. Previously the pair of drifts was conditioned on at least one being positive and a losing accumulator was allowed a negative rate, which it kept forever. The two are different models of the same race wherever a drift is small relative to its SD: densities up to about 40-50% apart in the tails, choice probabilities a few percentage points apart. The change makes cogmod LBA estimates directly comparable with those packages and with the literature built on them, and dcogmod_lba2() now reproduces rtdists::dLBA() at the same parameter values. It also simplifies the sampler, which draws each drift from its truncated Normal rather than splitting the conditional law into cases. Fits made with earlier versions cannot be post-processed with this one, and their estimates are not comparable with new ones at low drift rates; the vignette model was refit. The loser’s survival is computed as P(v > 0, unfinished) / P(v > 0) from whichever tail keeps its digits, so the density stays accurate for a loser with a strongly negative drift, where both quantities are tiny, and in the far tail, where 1 - CDF would cancel. cogmod_lba1() is unaffected: with one accumulator the two conventions coincide.

The truncation has one cost, and cogmod_priors() now covers it: once an accumulator rarely wins, its drift and sigma are identified only through |drift| / sigma^2 (the truncated Normal converges to an Exponential along that ray), so a flat prior lets the drift run off - the vignette’s error accumulator sat at -12 with an interval of -23 to -6.5. cogmod_priors() therefore puts normal(1, 2) on driftone and normal(0, 1.5) on its slopes, the treatment cogmod_lnr()’s nuone already had. Existing formulas that left driftone to brms get this prior on their next cogmod_priors() call.

cogmod 0.3.0

CRAN release: 2026-09-12

CRAN Publication.