Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add 5-parameter diffusion drift; fixes #2584 #2585

Conversation

kendalfoster
Copy link

Summary

This pull request implements the 5-parameter variant of the DDM by adding two new functions: ddm_lpdf() and ddm_lcdf().

This pull request addresses issue #2584: Adding 5-Parameter Drift Diffusion Model (DDM) PDF and CDF to Stan Math. This pull request has been discussed in this Discourse thread.

Tests

As recommended in the guidelines, I added the following tests to mirror that of the related wiener_lpdf() function that already exists in Stan Math:

  • test/prob/ddm/ddm_cdf_log_test.hpp
  • test/prob/ddm/ddm_test.hpp
  • test/unit/math/prim/prob/ddm_cdf_log_test.cpp
  • test/unit/math/prim/prob/ddm_test.cpp

Side Effects

There should be no side effects.

Release notes

The 5-parameter variant of the DDM will become available in Stan Math.

Checklist

@bob-carpenter bob-carpenter changed the title fixes #2584 add 5-parameter diffusion drift; fixes #2584 Sep 21, 2021
@bob-carpenter
Copy link
Contributor

Thanks for the PR! I don't understand these densities well enough to review myself. I'm afraid that everyone who might review this is scrambling to fix outstanding issues for the 2.28 code freeze this week. If someone doesn't review in the next few days, please ping again.

@kendalfoster
Copy link
Author

Hi, thanks for the comment! As you suggested, I'm adding this comment to serve as a reminder for this PR.

In the failed check, it looks as though there are 3 main types of error in the dist.log file:

  1. test/prob/ddm/ddm_00059_generated_vv_test.cpp:9:163: error: type name does not allow storage class to be specified (example in line 1686 of dist.log file)
  2. ./test/prob/test_fixture_distr.hpp:83:26: error: no matching member function for call to 'log_prob' (example in line 1771 of dist.log file)
  3. ./test/prob/ddm/ddm_test.hpp:179:68: note: candidate function template not viable: requires 8 arguments, but 6 were provided (example in line 1764 of dist.log file)

The third type of error seems the most straightforward to me, but I can't figure out how only 6 arguments were provided to the function in question. If you could advise me on how to fix any of these errors, I'd really appreciate it! Thanks!

Copy link
Member

@syclik syclik left a comment

Choose a reason for hiding this comment

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

Lots of things to change. I'm gonna stop here for now.

I haven't looked at the math yet, but cleaning this up will help get to the math.

#include <stan/math/prim/err/throw_domain_error.hpp>
#include <stan/math/prim/err/invalid_argument.hpp>

// Open the Namespace
Copy link
Member

Choose a reason for hiding this comment

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

Please remove this comment; this information is assumed for C++.


/**
* The log of the first passage time distribution function for a
* (Ratcliff, 1978) drift diffusion model with intrinsic trial-trial variability
Copy link
Member

Choose a reason for hiding this comment

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

Could you include a full reference for the "(Ratcliff, 1978)" below? Maybe line 24.

#ifndef STAN_MATH_PRIM_PROB_DDM_LPDF_HPP
#define STAN_MATH_PRIM_PROB_DDM_LPDF_HPP

#define _USE_MATH_DEFINES
Copy link
Member

Choose a reason for hiding this comment

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

@kendalfoster What is this used for? This looks like it could cause unintended side effects.

static const int max_terms_large = 1; // heuristic for switching mechanism
static const double ERR_TOL = 0.000001; // error tolerance for PDF approx
static const double SV_THRESH = 0; // threshold for using variable drift rate
static const double LOG_PI = log(M_PI);
Copy link
Member

Choose a reason for hiding this comment

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

There is already a declaration of LOG_PI in the math library. Please:

  1. include this header at the top: stan/math/prim/fun/constants.hpp
  2. remove this line. LOG_PI should be defined in this namespace and it should resolve ok.

static const double ERR_TOL = 0.000001; // error tolerance for PDF approx
static const double SV_THRESH = 0; // threshold for using variable drift rate
static const double LOG_PI = log(M_PI);
static const double LOG_2PI_2 = 0.5 * log(2 * M_PI);
Copy link
Member

Choose a reason for hiding this comment

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

Please replace LOG_2PI_2 with the constant HALF_LOG_TWO_PI found in the same constants.hpp above.

if (Nt0 < 1) { // t0
return 0;
} else {
for (size_t i = 0; i < Nt0; i++) {
Copy link
Member

Choose a reason for hiding this comment

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

Please use check_positive_finite().

You'll have to include stan/math/prim/err/check_positive_finite.hpp at the top.

if (Nw < 1) { // w
return 0;
} else {
for (size_t i = 0; i < Nw; i++) {
Copy link
Member

Choose a reason for hiding this comment

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

Please use check_bounded() with 0, 1.

if (Nsv < 1) { // sv
return 0;
} else {
for (size_t i = 0; i < Nsv; i++) {
Copy link
Member

Choose a reason for hiding this comment

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

Please use check_positive().

}
}

if (!include_summand<propto, T_rt, T_response, T_a, T_v, T_t0, T_w,
Copy link
Member

Choose a reason for hiding this comment

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

After all the checks, this is where you'd put something like:

  if (size_zero(rt, response, a, v, t0, w, ...)) {
    return 0.0;
  }

double t, a_i, v_i, w_i, sv_i;
for (size_t i = 0; i < Nmax; i++) {
// Check Parameter Values
t = rt_vec[i % Nrt]
Copy link
Member

Choose a reason for hiding this comment

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

The pattern of indexing into rt_vec isn't right.

The purpose of the view (scalar_seq_view<T_rt_ref> rt_vec(rt_ref)) is to treat rt_vec as if it were a vector even if it's a scalar.

So... this can be simplified to something like:

t = rt_vec[i] - t0_vec[i];

But... when you see that, you'll see the types are wrong now. t isn't a double, right? It could either be double or var (or fvar<T>). I don't think this should work if you've tested this with rt as var and t0 as double.

Copy link
Member

@syclik syclik left a comment

Choose a reason for hiding this comment

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

Lots of things to change. I'm gonna stop here for now.

I haven't looked at the math yet, but cleaning this up will help get to the math.

@syclik syclik marked this pull request as draft May 23, 2023 13:26
@andrjohns
Copy link
Collaborator

Closed by #2822

@andrjohns andrjohns closed this Mar 21, 2024
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.

None yet

5 participants