Adding 7-Parameter Drift Diffusion Model (DDM) PDF with partial derivatives to Stan Math

[Franzi2114]

Description

The Drift Diffusion Model (DDM) introduced by Ratcliff provides a framework for analyzing data of fast binary decisions jointly modelling accuracy and response times. There are several parameterizations of the DDM. Currently Stan Math contains the PDF of the 4-parameter variant (i.e., the parameters for threshold separation, non-decision time, relative starting point, and drift rate) in the function wiener_lpdf(). Furthermore, @kendalfoster is working on a 5-parameter variant of the DDM (i.e., including the parameter for inter-trial variability in the drift rate) (see issue #2584)

I am working on a pull request to add a 7-parameter variant with partial derivatives of the DDM (i.e., additionaly including the parameters for inter-trial variabilities in relative starting point and non-decision time). For this goal, I first want to implement the lpdf and lcdf of the 4-parameter model wiener4_lpdf and wiener4_lcdf containing partial derivatives. In the next step, I want to use these functions to expand the model to a 7-parameter model containing partial derivatives wiener7_lpdf and wiener7_lcdf. These wiener7 functions will be built in the way that they call the wiener4 functions for the case that all three variability parameters are set to 0. Almost all computations are done on the log scale which makes computations more stable. All functions will be based on the implementation of the R-package WienR by @RaphaelHartmann. (Therefore, checks whether the correct values are computed can be done via the WienR package in R).

The 4-parameter functions will accept the following parameters: response time, threshold separation, drift rate, relative starting point, non-decision time, precision, response.

The 7-parameter functions will accept the following parameters: response time, threshold separation, drift rate, relative starting point, non-decision time, inter-trial variability drift rate, inter-trial variability relative starting point, inter-trial variability non-decision time, precision, response.

(There may occur minor changes in the parameter lists).

The idea is to implement the 7-parameter model in the way that it also can be used as a 4-, 5-, or 6-parameter model. For this, the user only has to set the corresponding variabilities to 0.

Example

4-parameter model:
An example call to the lpdf is wiener4_lpdf(1, 1, -1, .5, .1, .01, 1) and would return the log of the PDF.
An example call to the lcdf is wiener4_lcdf(1, 1, -1, .5, .1, .01, 1) and would return the log of the CDF.

7-parameter model
An example call to the lpdf is wiener7_lpdf(1, 1, -1, .5, .1, .5, .6, .7, .01, 1) and would return the log of the PDF.
An example call to the lcdf is wiener7_lcdf(1, 1, -1, .5, .1, .5, .6, .7, .01, 1) and would return the log of the CDF.

Expected Output

We should expect the example call to wiener4_lpdf(1, 1, -1, .5, .1, .01, 1) to return -4.24659 .
We should expect the example call to wiener4_lcdf(1, 1, -1, .5, .1, .01, 1) to return -1.323102.

We should expect the example call to wiener7_lpdf(1, 1, -1, .5, .1, .5, .6, .7, .01, 1) to return -1.90939.
We should expect the example call to wiener7_lcdf(1, 1, -1, .5, .1, .5, .6, .7, .01, 1) to return -1.289187.

Edit

In a first step, I want to implement wiener_full_lpdf. This function will be able to compute a 4-, 5-, 6-, and 7-parameter model. So: wiener4_lpdf and wiener7_lpdf will be parts of wiener_full_lpdf. The cdf function will be postponed.

Second Edit

The 7-parameter model wiener_full_lpdf is now implemented. We do not implement the 4-parameter model seperately and we also do not implement the CDF functions. The 7-parameter PDF functions are now in the testing phase and will soon be availabe.

Open Questions

I already started the implementation. There are still some open questions:

• How to handle the precision? (see this Discourse thread)

• Set the precision to a fixed value?
• Or let the user set a precision and set a default for the case the user does not give a precision?
• Then: How to implement a default value? Or if-statement whether the user set 0 and then internally set the precision?

• How to compute a multidimensional integral? (see this Discourse thread)

• At the moment I am working on a solution to use hcubature or a similar routine.
• Edit To incorporate the other two inter-trial variabilities I have to integrate over two variables. Therefore, I need to compute two bounded integrals to obtain the 7-parameter PDF. Is there a way to compute these integrals in Stan?
• Is there a way to translate the GPLv2+ hcubature into Stan-style and use it in BSD3 Stan-math?
• Can I somehow combine Julia with Stan?

• How to organise the file? (see this Discourse thread)

• As the computation of the derivatives is quite long, I compute the derivatives outside the lpdf function in an extra function. Shall I let the functions for the derivatives in the same file as the lpdf or save them in additional files? Then: Would stan/math/prim/prob be the correct spot to save the derivatives?

• When do I need the bool<propto>?
• How to deal best with responses?

• In the DDM there are two response options: Upper bound = 1, lower bound = 0. The current wiener_lpdf only computes the lpdf for upper bound responses. For the lower bound responses the user has to change the input data. But it is also possible to give the response as a parameter to the function and then choose the upper/lower bound model internally. Do you have suggestions/pros/cons?

• Edit Now, the functions for wiener_full_lpdf are almost ready for the pull request. Which tests and steps do I have to do before I send the pull request? Can anybody give me a list or a link where I find everything?

Current Version:

v4.3.0

[bob-carpenter]

Thanks, @Franzi2114. I can help with some of the questions:

• if there's a variable precision algorithm that's really costly as precision goes up, then the most friendly API makes its user-controllable with reasonable defaults; usually we try to hit 5 or 6 digits at least, as we've found less than that won't work with our ODE solvers (of course, this depends if the error control is stepwise or global); the best thing to do is run some experiments with typical values

• I have no idea bout hquadrature, but the usual bottleneck for ay quadrature algorithm is dimensionality; is this something that could be evaluated with nested MCMC (we don't have a way to do that with HMC at the moment)

• helper functions can go in the same file and should be placed in namespace stan::math::internal

• when propto is true (1), you can drop constants in the log density that don't depend on any of the autodiff variables (all parameters will be autodiff variables), as they don't matter to Stan's HMC algorithm for sampling

• as for the response options, I'd suggest two different functions rather than one with a flag; then you can have a reasonable default with no extra args and then one marked with a suffix like _upper. Or, if the user can workaround with data transforms, it makes sense to only implement one

[Franzi2114]

Hi, I have two more questions:

with reasonable defaults

@bob-carpenter, how is the way to define defaults for our new function? Is it related to the Stan_math_signatures.ml or is it in the definition of the function itself?

The function is now nearly finished and ready for a pre-publication. I first want to add a new branch with this feature in Stan/math and not make a Pull Request to the develop branch at the moment. This shall come a bit later. How can I add a new branch to Stan/math without writing access?

Best,
Franzi

[bob-carpenter]

Sorry for not answering sooner, but I've been traveling.

how is the way to define defaults for our new function?

Stan doesn't have a notion of default arguments or of named arguments. Every argument to a Stan function is required. What we've done with things like the ODE solver is provide two functions (see https://mc-stan.org/docs/functions-reference/functions-ode-solver.html for examples), one with defaults and one where you specify all the tolerances.

I first want to add a new branch with this feature in Stan/math

You can't do that until you are a Stan developer with permissions. Which we don't grant to people until they've submitted at least one substantial PR. Instead, you will need to fork the project, branch from there, then submit a PR from your branch. When you do that, we'll be able to import your new branch for code review.

[Franzi2114]

Hi @bob-carpenter , thanks for your reply!

then submit a PR from your branch. When you do that, we'll be able to import your new branch for code review.

Is this PR from my fork then the official PR such that my new branch will be merged into the develop branch when it passes code review and all tests OR will it end up in a new branch in Stan/math and another PR is needed for the final merge into the develop branch? I would prefer the latter one, since my plan is that with the pre-publication people who are interested can use the function, test it and give feedback if something should be changed.

[bob-carpenter]

@Franzi2114: The former, but PRs can be modified as much as you want before they're merged. You just push new commits to your branch. I'd suggest reading some introductions to how PRs work on GitHub, but most of them are so complicated and contain so much extraneous detail that I find them all very confusing.

[Franzi2114]

The PR is in progress! (#2822)

The code can be found in my fork.

An instruction how to integrate this stan version into a local environment can be found here.

[Franzi2114]

Hey @bob-carpenter,

now the PR for this issue is open (#2822) and all tests run successfully. I have some questions to the next steps: Who starts the workflows and who reviews the code? Shall I ask specific people for this who are firm with this type of function or does the developer team announce reviewers?

Best, Franzi

[bob-carpenter]

Thanks for submitting. Anyone with push permission on the project can kick off the continuous integration tests. So I started them off for you. I think after that, they will continue to run after updates, but if not, let me know.

Any of our developers are allowed to review anything---we expect people to use their judgement in what they review. And we only have a single reviewer deciding, so as not to get conflicting reviews, but anyone can comment on a PR. If you have someone in mind to review, that'd be great. I'm not really at all familiar with these models, so I shouldn't be the one to to do it.

[Franzi2114]

Hey @bob-carpenter,
I already tried to find a person but unfortunately this person was familiar with the function but not with the Stan source code. I checked the Stan developers list but found no other person that had sufficient contact with diffusion models, on the first glance.

What is the idea of this review process? I ask this because what we implemented strongly follows the WienR package of R, which is implemented in C++. In principle, there is nothing new in the code we provide. It is just adapted to the Stan language. Also the hcubature implementation exists in C++ and in Julia. Therefore, we are pretty sure that the single parts of the function do the things we claim, since we can check 1:1 against the R package. This means that the reviewing person need not necessarily be very familiar with diffusion models. What do you think?

Is it maybe possible to combine knowledge of different people: We have the ones familiar with the function, but not with the Stan language itself, and the ones familiar with the Stan language but not with the function, and the ones maybe with knowledge of both but with no write access.

It seems that in the meantime, no other person did anything with the PR. Should I continue searching for a suitable person? Shall I open a discourse thread on this? Or is there something like an editorial board that asks people for a review? Or do the people still exist that reviewed the wiener_lpdf() function when it had its PR? Is it possible to find out who this was?

Best, Franzi

[bob-carpenter]

Sorry about that @Franzi2114. I know how frustrating it is to have written a function and not get it reviewed. We can definitely have multiple people review.

Someone familiar with the function should check that the definition makes sense and that there is a sufficient testing plan. @charlesm93 may be able to do that.

Then someone who knows the current state of uur C++ library should be able to review the code. The most active people on the C++ side these days who are reviewing distribution code are @spinkney and @rok-cesnovar and @andrjohns. @SteveBronder may also be able to do it.

You can do a search for "wiener" on the closed PRs: https://github.com/stan-dev/math/pulls?q=is%3Apr+is%3Aclosed

But that's basically broken, not on the GitHub side, but because we do continuous integration regression tests, one of which is the Wiener distribution code, so there are way too many false positives

Scrolling down to ones with "wiener" in the title, it looks like @t4c1 and @bbbales2 were the most recently active, but they graduated and finished their postdocs respectively, and aren't really active on the project any more.