Question about Wiener Diffusion Model Distribution in rstan_2.7.0-2

[Andrew Ellis]

Has anyone successfully used this?

I tried to implement the example given here: http://www.cidlab.com/supp.php?o=wienerstan

The 'wiener.hpp' file states that "y should contain reaction times in seconds, with upper-boundary responses strictly positive and lower-boundary response times coded as strictly negative numbers." 

My example code:

library(rstan)
library(RWiener)


# generate data using RWiener package
set.seed(0)
dat <- rwiener(n = 1000, alpha = 1, tau = 0.2, beta = 0.5, delta = 0.5)

standata <- list(y = ifelse(dat$resp == "upper", dat$q, -1 * dat$q),
                 N = length(dat$resp))


wiener_model <- "
data {
    int<lower=0> N;
    real y[N];
}
parameters {
    real<lower=0> a;
    real<lower=0> t;
    real<lower=0,upper=1> b;
    real v;
}
model {
    v ~ normal(0, 2);


    for (n in 1:N) {
        y[n] ~ wiener(1, 0.2, 0.5, v);
    }
}"


fit <- stan(model_code = wiener_model, data = standata)

I get the error message: 

Rejecting initial value:
  Error evaluating the log probability at the initial value.


Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception thrown at line 16:
stan::math::wiener_log(%1%): Random variable is -0.395369, but must be > 0!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

According to the manual (p. 396), Stan returns the first passage time of the accumulation process over the upper boundary only. How do I deal with both decision boundaries, if I cannot pass in lower-boundary response times as negative numbers?

cheers,

Andrew

> sessionInfo()
R version 3.2.1 (2015-06-18)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.4 (Yosemite)


locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8


attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base    


other attached packages:
[1] RWiener_1.2-0  rstan_2.7.0-2  Rcpp_0.12.0    magrittr_1.5   devtools_1.8.0


loaded via a namespace (and not attached):
 [1] codetools_0.2-14   digest_0.6.8       R6_2.1.0           stats4_3.2.1       git2r_0.10.1       httr_1.0.0        
 [7] stringi_0.5-5      curl_0.9.1         xml2_0.1.1         tools_3.2.1        stringr_1.0.0.9000 inline_0.3.14    
[13] rversions_1.0.2    memoise_0.2.1    

[Andrew Ellis]

this seems to work:

lower boundary responses should no longer be coded as negative response times; instead the response times and responses should be passed to stan as data:

library(rstan)

library(RWiener)

set.seed(0)
dat <- rwiener(n = 5000, alpha = 1, tau = 0.2, beta = 0.75, delta = 0.5)
standata <- list(y = dat$q,
                 resp = ifelse(dat$resp == "upper", 1, 0),
                 N = length(dat$resp))

and in the stan model, y ~ wiener(a, t, b, v) if the "upper" response was selected, and y ~ wiener(a, t, 1-b, -v) if the "lower" response was selected.

data {
  int<lower=0> N;
  real y[N];

  int<lower=0, upper=1> resp[N];

}
parameters {
  real<lower=0> a;
  real<lower=0> t;

  real<lower=0, upper=1> b;

  real v;
}
model {
  v ~ normal(0, 2);

  a ~ gamma(3, 3);
  t ~ gamma(4, 10);
  b ~ beta(5, 5);

  for (n in 1:N) {

      if (resp[n] == 1) {
          y[n] ~ wiener(a, t, b, v);
      }
      else {
          y[n] ~ wiener(a, t, 1-b, -v);
      }


  }
}

Full code example:

library(rstan)
library(RWiener)


set.seed(0)
dat <- rwiener(n = 5000, alpha = 1, tau = 0.2, beta = 0.75, delta = 0.5)

## stan model ----
wiener_model <- "

data {
    int<lower=0> N;
    real y[N];

    int<lower=0, upper=1> resp[N];

}
parameters {
    real<lower=0> a;
    real<lower=0> t;

    real<lower=0, upper=1> b;

    real v;
}
model {
    v ~ normal(0, 2);

    a ~ gamma(3, 3);
    t ~ gamma(4, 10);
    b ~ beta(5, 5);

    for (n in 1:N) {

        if (resp[n] == 1) {
            y[n] ~ wiener(a, t, b, v);
        }
        else {
            y[n] ~ wiener(a, t, 1-b, -v);
        }


    }
}"


## stan data ----
standata <- list(y = dat$q,
                 resp = ifelse(dat$resp == "upper", 1, 0),
                 N = length(dat$resp))

fit <- stan(model_code = wiener_model, data = standata)

[Bob Carpenter]

Our doc needs some work --- I certainly couldn't answer your
question from looking at it. I'm going to add an issue.

https://github.com/stan-dev/stan/issues/1576

Thanks for bringing it up. We'll cite you in the manual for
the great comments when we update our reference manual entry and
hopefully provide an example. Do you mind if we use yours if
we cite it? It all gets released under BSD for code and CC-BY
for text.

- Bob

[Marco Inacio]

I think the problem here was that Andrew looked at the documentation in
wiener_log.hpp file which wrongly describes the possibility of lower
boundary by negating y (random variable) which was originally what the
patch of Joachim Vandekerckhove did.

I changed this (thus forcing y to be strictly positive) because there
was no computational advantage in doing so and thought it would be
potentially dangerous for users. I realize I should have discussed such
change with stan-dev at the time.

The Stan manual, however, correctly describes the current behavior
(although it does deserves a better explanation), and I didn't put it
there that y must positive because it is already in the category
"Non-negative Continuous Distributions".

[Bob Carpenter]

> On Aug 1, 2015, at 12:17 AM, Marco Inacio <marcoig...@gmail.com> wrote:
>
> I think the problem here was that Andrew looked at the documentation in wiener_log.hpp file which wrongly describes the possibility of lower boundary by negating y (random variable) which was originally what the patch of Joachim Vandekerckhove did.
>
> I changed this (thus forcing y to be strictly positive) because there was no computational advantage in doing so and thought it would be potentially dangerous for users. I realize I should have discussed such change with stan-dev at the time.
>
> The Stan manual, however, correctly describes the current behavior (although it does deserves a better explanation), and I didn't put it there that y must positive because it is already in the category "Non-negative Continuous Distributions".

I'd prefer to follow the pattern in the other distribution doc
which clearly indicates both the support for the outcome and any
constraints on parameters. Would you mind adding that? I'm still not
sure exactly what the behavior is, or I'd do it myself. I have no
idea what this density does!

- Bob

[Marco Inacio]

> I'd prefer to follow the pattern in the other distribution doc
> which clearly indicates both the support for the outcome and any
> constraints on parameters. Would you mind adding that?

OK. Don't worry, I'll take care of it.

[Andrew Ellis]

hi Bob, of course I don't mind, but my example is similar to the one provided by Joachim ( http://www.cidlab.com/supp.php?o=wienerstan). But feel free to use my Stan code.

[Andrew Ellis]

hi all, 

sorry for the confusion.

On Saturday, August 1, 2015 at 6:17:52 AM UTC+2, Marco Inacio wrote:

I think the problem here was that Andrew looked at the documentation in
wiener_log.hpp file which wrongly describes the possibility of lower
boundary by negating y (random variable) which was originally what the
patch of Joachim Vandekerckhove did.

I had previously used this, so I assumed my example (using negative response times) would  work.

I changed this (thus forcing y to be strictly positive) because there
was no computational advantage in doing so and thought it would be
potentially dangerous for users. I realize I should have discussed such
change with stan-dev at the time.

The Stan manual, however, correctly describes the current behavior

it does 

(although it does deserves a better explanation), and I didn't put it
there that y must positive because it is already in the category
"Non-negative Continuous Distributions".

This is true, but I found the Wiener distribution in the manual by searching, rather than in the TOC, so I missed the bit about it being a non-negative continuous distribution. Sorry

[Matt]

Hi Marco,

On Saturday, 1 August 2015 00:17:52 UTC-4, Marco Inacio wrote:

I changed this (thus forcing y to be strictly positive) because there
was no computational advantage in doing so and thought it would be
potentially dangerous for users.

 
I'm curious as to the potential dangers of not having y as strictly positive for the Wiener distribution?

Matt

[Bob Carpenter]

I don't know anything about this distribution and there's
no definition in the manual. Is it symmetric around 0?

If so, I don't see any reason to constrain it to be positive.
You might want to do so in a model to avoid multimodality if
+theta and -theta have the same density.

- Bob

[Mike Lawrence]

Yes, it'll be strictly positive. It represents the distribution of times it takes a quasi-random walker to get from a given starting point to one of a pair of given ending points on a 1-D plane. The walker is quasi-random because it may have a slight preference for toward one boundary or the other.

[Marco Inacio]

Suppose for example the user does some manipulation in his/her data that should be strictly positive, but something goes wrong, and it gets negative, then the distribution should give him an error since negative time does not indeed exist.

If you check the implementation of this distribution in R (same author as this one), the function has an argument resp = "upper"/"lower"/"both", even though, internally the function uses negative/positive responses.

So I think that's the correct way of doing it, and in the future an extra argument for lower/upper/both responses can be added (with default to "upper"), but then this wouldn't be backward compatible if the function was accepting negative responses right now.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Mike Lawrence]

A little off-topic, but possibly of interest to those wondering what this distribution is used for:

Elaborations of this model have been used extensively in the study of speeded human decisions (ex. simple perceptual decisions like "is this square black or white?" or more abstract decisions like "is this word English or French?") since the 1970s to connect the observed data of response speed and response accuracy to latent variables that can be more informative on the structure of the mind. For example, a given experimental manipulation might yield fast-but-error-prone responses in one condition and slow-but-accurate responses in another. The manipulation obviously causes a trade-off of speed and accuracy, but is there any additional changes to information processing as well? Where the diffusion model explicitly operationalizes information processing efficiency (bias of the walker; more biased=faster information processing) as independent of the speed-accuracy trade-off criterion (separation of the bounds; closer bounds = more emphasis on speed, less on accuracy), one can answer this question. Other parameters of the more elaborated models can address other interesting questions and there are also several competing models (ex. linear ballistic accumulator, log-normal race, etc) that generally seek to address similar questions. I'll actually be coding up an example of both the diffusion and log-normal race (which is easy with the distributions already in stan) in the next few days, and I'll post it to the examples for anyone interested.

[Bob Carpenter]

> On Aug 4, 2015, at 11:34 AM, Marco Inacio <marcoig...@gmail.com> wrote:
>
> Suppose for example the user does some manipulation in his/her data that should be strictly positive, but something goes wrong, and it gets negative, then the distribution should give him an error since negative time does not indeed exist.

Usually it just means the past and most of these random walks
are reversible. That's how we use negative time in NUTS, for
example.

> If you check the implementation of this distribution in R (same author as this one), the function has an argument resp = "upper"/"lower"/"both", even though, internally the function uses negative/positive responses.

I'm still confused about the function, so I'll have to wait
for the doc.

Just to state the obvious, if there are constraints on variables,
either outcomes or parameters, then they should be checked.

> So I think that's the correct way of doing it, and in the future an extra argument for lower/upper/both responses can be added (with default to "upper"), but then this wouldn't be backward compatible if the function was accepting negative responses right now.

The way to maintain backward compatibility is to just add
a new function.

- Bob

[Marco Inacio]

>> If you check the implementation of this distribution in R (same author as this one), the function has an argument resp = "upper"/"lower"/"both", even though, internally the function uses negative/positive responses.
> I'm still confused about the function, so I'll have to wait
> for the doc.
>
> Just to state the obvious, if there are constraints on variables,
> either outcomes or parameters, then they should be checked.

In this distribution, we model the probability (actually probability
density) of the an specific process reach an upper boundary (alpha) or a
lower boundary (zero) (or both) in an specific time (time is the
response variable).


Examples:

> require(RWiener)

> # Probability of reaching the upper boundary alpha=2 in at most 100
units of time
> pwiener(100, alpha=2,tau=0.3, beta=0.45, delta=.4, "upper")
[1] 0.6430842

> # Probability of reaching the lower boundary 0 in at most 100 units
of time
> pwiener(100, alpha=2,tau=0.3,beta=0.45,delta=.4, "lower")
[1] 0.3569158

> # Probability of reaching any of the two boundaries in at most 100
units of time
> pwiener(100, alpha=2,tau=0.3,beta=0.45,delta=.4, "both")
[1] 1

(100 is so large for this values of the parameters that the pdf is
integrating to 1 already)

We are discussing which is the best way to select "upper"/"lower"/"both"
in Stan. I think the best way is to have an extra argument at the end
wiener_log just like the R examples above, but the author originally had
negative response variable meaning "lower", positive response variables
meaning "upper" (although the response variable of the process is in
fact time elapsed, so it's strictly positive)

[Guido Biele]

For what it's worth I believe the current implementation of the wiener is fine.

Of course, a short example in the documentation that explains how the parameters are associated with different aspects of decision making would be helpful.

(I can help with that if needed.)

I think more interesting than having a wiener function that allows to specify if the likelihood for the lower or upper boundary should be returned would be to see if the current implementation can be made faster as described here: http://cvc.psy.ku.dk/gondan-2014.pdf.

From my short experience i got the impression that the wiener function is relatively slow (but it converges extremely fast to the right parameter values!)

anyhow, thanks for making the function available!

cheers  - guido

[Bob Carpenter]

Thanks --- that's very helpful. I don't understand why diffusions and
crossing time are a good model of trading speed for accuracy.

If you create the models and simulate data (or even better have
a program you can share in Python or R to generate data),
we can put them in stan-dev/example-models if you're willing
to release them under BSD.

Are there plots of what the Wiener density looks like for various
parameters somewhere? I suppose I could make them with Stan now.

- Bob

[Bob Carpenter]

> On Aug 4, 2015, at 1:59 PM, Guido Biele <guido...@gmail.com> wrote:
>
> For what it's worth I believe the current implementation of the wiener is fine.
> Of course, a short example in the documentation that explains how the parameters are associated with different aspects of decision making would be helpful.
> (I can help with that if needed.)

Thanks, but I'm pretty sure Marco's most recent pull
requests have already done that.

> I think more interesting than having a wiener function that allows to specify if the likelihood for the lower or upper boundary should be returned would be to see if the current implementation can be made faster as described here: http://cvc.psy.ku.dk/gondan-2014.pdf.

Going faster is always good (assuming we don't have to
sacrifice accuracy, if I can be self-referential to Mike
Lawrence's application). It's always easier to optimize
once there's a working version with tests in place.

I don't know if any of the Stan developers will have time to
do this, but we do take pull requests!

- Bob

[Marco Inacio]

On 2015-08-04 11:25 PM, Bob Carpenter wrote:
>> On Aug 4, 2015, at 1:59 PM, Guido Biele <guido...@gmail.com> wrote:
>>
>> For what it's worth I believe the current implementation of the wiener is fine.
>> Of course, a short example in the documentation that explains how the parameters are associated with different aspects of decision making would be helpful.
>> (I can help with that if needed.)
> Thanks, but I'm pretty sure Marco's most recent pull
> requests have already done that.

Actually it didn't (how disappointing it must be :-( ), since I never
used this model, so help would be appreciated indeed.

[Joachim Vandekerckhove]

Hi all,

Sorry I'm a little late to this party, I had a long summer of travel.  The reason we originally implemented the lower-bound hits as negative reaction times was to have a cheap way of encoding what is essentially a bivariate distribution (strictly positive RT and binary response).  Because both the RT and the outcome are random variables, the current implementation (with only one bound) isn't entirely satisfactory if at any point we want to generate random numbers from the distribution.

I always liked the negative-RT solution even though it's a little counter-intuitive.  I'd probably prefer it in Stan since it's how it works in JAGS... but the more conventional solution would be to make it actually bivariate, and let it operate on a [RT response] pair, where the first is a strictly positive real and the second is a boolean.

Joachim

[Bob Carpenter]

Thanks --- we'll try to make that clear in the manual. Here's the
manual issue comment for it:

https://github.com/stan-dev/stan/issues/1617#issuecomment-144792409

- Bob

> To post to this group, send email to stan-...@googlegroups.com.

[Royce Anders]

in this implementation of the two-boundary Wiener process in STAN/JAGS, is it missing a piece of information: particularly when the process terminates due to hitting the opposite boundary than the sign of the drift rate? 

e.g., :    if (resp[n] == 1) {

          y[n] ~ wiener(a, t, b, v);

      }

      else {

          y[n] ~ wiener(a, t, 1-b, -v);

      }

Is it not a completely true two-boundary Wiener process then? 

I understand that the probability is quite low that the opposite boundary (e.g., the lower boundary will be hit when the drift v is positive), however it becomes increasingly probable as the threshold values (a) are small, or the drifts (v) are slow. In the case of (a) it may be quite important to capture speed-accuracy manipulations in experiments.  It's my understanding that this is an important/defining aspect of the diffusion decision model (DDM), hence am I missing something here or is this mirrored application (e.g., as in two conjugate inverse gaussians)  capturing this important aspect of the data when we fit by this implementation? 

[Mike Lawrence]

My recollection is that the distribution is defined so that negative values represent the times of opposite-boundary crossings. So you code your data with negative signs on the rts for error trials, and the distribution knows to treat these as opposite boundary crossings.

Mike

[Royce Anders]

ah right.. I think the trick is in what I missed before, the 1-b, with the -v, and the Wiener dist pdf here is only calculating the probability that the upper bound is hit, as acting like the lower bound due to inversing these two parameters together. Great :)

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/-6wJfA-t2cQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Maciej Swat]

Hi all

I have a related comment.
There are three potential typos in the STAN manual (2.9.0) in the
PDF for 'Wiener Diffusion Model Distributions’, but I’m not sure if these
bugs 2&3 proliferated into the code.

1. y is defined as 'y \in (0, \tau)’ but it should be 'y \in (\tau, \infty)’
otherwise the term \sqrt{y - \tau} is undefined.


Two others based on reference paper (Blurton et at. 2012)

2. the exponent of the denominator (y - \tau) should be 3/2 instead of 3
See attached PDF plots for
alphaB=3.25;tauB=0.25;betaB=0.4;deltaB=0.65; N = 150;

3. the numerator in standard normal distribution function \phi reads
‘(2k + \beta) \alpha’ instead of '2k + \alpha \beta’.
The difference between these two results in very minor changes.

Best Regards, Maciej

[Bob Carpenter]

[Cross-posted to stan-dev]

Thanks for the detailed comments.

Could someone who understands this Wiener thing please respond
on stan-users?

I've never understood the function or its doc and we get a steady
stream of these messages saying it's wrong, which in the past
have been due to poor doc, not buggy code.

The next step would be to create reproducible errors.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <Rplot03.pdf><Rplot05.pdf>

>
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <Screen Shot 2016-02-19 at 10.44.08.png>

>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <Screen Shot 2016-02-19 at 10.44.19.png>

[Guido Biele]

Hi,

I agree with Maciej that there is a typo in the manual (which I hadn't seen myself, but i double-checked after Maciej's message).

regarding the code: my impression when doing a simple parameter recovery test was that the wiener distribution (which should really be called wiener first passage time distribution) in stan recovers parameters from data generated with the rwiener package. 

I don't habe time to redo this now, but maybe Maciej's can double-check my parameter recovery exercise here: https://github.com/gbiele/stan_wiener_test.

best guido

[Guido Biele]

Hello,

I am attaching a figure for a parameter recovery exercise I did.

It looks as if the wiener function in stan does a good job recovering most parameters except the non-decision time tau 

(where the parameter tau is not really needed to calculate the wiener first passage time distribution, but models non-decision related 

aspects of a task including time for perceiving a stimulus and executing a motor-response).

The figure shows for 4 parameters how 5 different parameter value were recovered. I did each recovery experiment 5 times to get an idea of the variability of the effect.

I think the results are kind of a mixed bag. on the one hand it is good that one can generally recover parameter values well.

on the other hand I am not sure how much I should be concerned about the fact that one can recover the order, but not all values for the parameter tau.

Still, before concluding that something could be wrong with the code itself, one should run a larger parameter recovery experiment.

Cheers - Guido

[Eric-Jan Wagenmakers]

I support Joachim's statement. I have extensive experience with this model, and the way it was defined before (with negative RTs signifying lower-bound responses) is more straightforward, and, importantly, greatly reduces the probability of coding the model incorrectly. It also yields cleaner code. This may not be immediate obvious to those who have not worked with the model, but it is nevertheless true.

Cheers,

E.J.

[Bob Carpenter]

I've never understood any of the discussion around the Wiener
distribution. Did its behavior change at some point? If so,
that got snuck in without my realizing, as I've been trying really
hard to enforce backward compatibility in the modeling language.

- Bob

> To post to this group, send email to stan-...@googlegroups.com.

[Eric-Jan Wagenmakers]

The behavior is OK. It was just that the JAGS version used a more
intuitive coding that was less error-prone. The Wiener process
produces two variables. One indicates response choice and is binary
(correct/error, or left/right, or whatever), the other indicates
response duration (or response time: rt) and is continuous. In the
JAGS code, these bivariate outcomes were captured in a single number
by using the sign. So .313 indicates an rt of 313 ms that ended up at
the upper bound, and -.400 indicates an rt of 400 ms that ended up at
the lower bound.

In the Stan code, the negative sign was eliminated and only positive
numbers can be used. Confronted with lower bound responses, the user
now needs to invert both the starting point and the drift rate. This
doubles the Stan code length and increases the probability of making
an error.

I can understand that the original idea was "what the hell are these
negative times good for?", but they were there for a very good reason.
The code still works, but reverting to the old input procedure (the
one with the negative rts) would be a good idea, imo.

Cheers,
E.J.
********************************************
Eric-Jan Wagenmakers
Department of Psychological Methods, room G 0.29
University of Amsterdam, Nieuwe Prinsengracht 129B
Letter: PO Box 15906, 1001 NK Amsterdam
Parcel: Valckenierstraat 59, 1018 XE Amsterdam

Web: ejwagenmakers.com
Book: bayesmodels.com
Stats: jasp-stats.org
Email: EJ.Wage...@gmail.com
Phone: (+31) 20 525 6420

“Man follows only phantoms.”
Pierre-Simon Laplace, last words
********************************************

> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/-6wJfA-t2cQ/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

Thanks for the clarification. We'd probably add a new function
and deprecate the old one --- I don't like breaking backward
compatibility!

If someone wants to do this, it should be a relatively straightforward
project (if you can figure out what the whole thing does, that is).

- Bob

[matus.s...@privatdemail.net]

On Wednesday, April 27, 2016 at 8:35:39 PM UTC+2, Eric-Jan Wagenmakers wrote:

In the JAGS code, these bivariate outcomes were captured in a single number
by using the sign. So .313 indicates an rt of 313 ms that ended up at
the upper bound, and -.400 indicates an rt of 400 ms that ended up at
the lower bound.

How do you represent missing values? How do you represent situation in which only one of the two variables is missing? How do you represent censored values? How do I know whether a value is an observed score or an upper or lower threshold of a censored score?

I use the negative response time values to flag the special cases and hence I'm in favor of the two-vector representation. In addition, it's consistent with the representation used by generalized linear models of response time-accuracy trade-off in IRT (e.g. van der Linden, W. J. (2007). A hierarchical framework for modeling speed and accuracy on test items. Psychometrika, 72(3), 287-308.) Sure, I and the IRT folks can translate our code to fit the specification used by response-time modelers. Sure I can add a vector or two that encodes thresholds. But then I have to complain that  this "doubles the Stan code length and increases the probability of making an error".

[Eric-Jan Wagenmakers]

Well, the model is very popular in RT research, from which field it originated. I would use the coding that most users are comfortable with.
Cheers,
EJ

--

[Joachim Vandekerckhove]

There really is a single canonical solution; at least as far as the diffusion model is concerned. It's a bivariate model so the likelihood should be bivariate. I certainly understand that people are uncomfortable with the negative-RTs hack, and it's not good coding practice. But neither is implementing a bivariate distribution as two conditional univariate distributions.

Since there is a concern with breaking backward compatibility, I second Bob's suggestion of introducing a new distribution, suggest we explicitly label it bivariate and call it the "double Wiener distribution" and I can provide code to implement ~dwiener() if someone can help me translate it to the correct C prototypes and whatnot.

I think that will satisfy everyone involved.

Joachim

[Bob Carpenter]

That'd be great. Thanks.

You should just be able to follow the existing C++
code. Also see:

https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan

And the rest of the wiki has info on the full developer process.

- Bob

> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.