Translating a JAGS module for Stan?

[Mike Lawrence]

Hi folks,

For decades it has been known that human response data from experiments involving simple perceptual discriminations are very well fit by a "diffusion model" consisting of a Wiener process unfolding through time between two bounds (reaching a bound causes the process to terminate, determining a response and the response time). I'd like to use Stan to perform inference on a hierarchical version of such a model (aggregating data from multiple human subjects, each measured hundreds of times). Apparently there already exists code in the form of a JAGS module (http://sourceforge.net/projects/jags-wiener/) that implements the sort of model I'd like to implement. Any advice on how I should go about translating this for use with Stan? Or are there tools already in Stan that could be used to implement this sort of model?

Mike

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[Andrew Gelman]

Mike,

Is there a pdf describing the model that is being fit in Jags?  I clicked through and found code but no description of the probability model that is being fit.  If there is an algorithm for computing the probability model, it should be possible to fit it in Stan, I'd think.

A

--
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/groups/opt_out.
 
 

[Mike Lawrence]

Hm, after some further reading (http://ppw.kuleuven.be/okp/_pdf/Tuerlinckx2004TECOT.pdf‎), it seems that the PDF is rather treacherous (see Eq 3), with parts that need numerical estimation on the fly. So possibly not a good project for the mathematically naive like myself :Op 

[Ben Goodrich]

Some of Francis' other papers on the diffusion model may be more directly applicable to Stan. Can you write down the model in .stan code and just use the (temporarily nonexistent) weiner() function to indicate where / how the log-likelihood comes into play? We might be able to do something reasonable to approximate the infinite sums with a finite number of terms, but the oscillation is going to be a challenge.

Ben

On Sunday, June 23, 2013 9:18:23 AM UTC-4, Mike Lawrence wrote:

Hm, after some further reading (http://ppw.kuleuven.be/okp/_pdf/Tuerlinckx2004TECOT.pdf‎), it seems that the PDF is rather treacherous (see Eq 3), with parts that need numerical estimation on the fly. So possibly not a good project for the mathematically naive like myself :Op 

On Saturday, June 22, 2013 10:25:32 PM UTC-3, Andrew Gelman wrote:

[Bob Carpenter]

1. If it's easier to define the joint density p(y,t), then
you could always just model both y and t and then forget
about y. Basically, MCMC would do the integration for you.

2. To implement Tuerlinckx's approach, we'd have to port the MATLAB
or Fortran iterative code to templated C++. The derivatives would
be slow if many iterations are required. If there were
a way to compute the derivatives directly (even if iterative),
that would be much faster.

3. Next year, we'll be working on building in diff eq solvers
for implicit equations, but I'm not sure that'll help here.

- Bob

On 6/23/13 9:18 AM, Mike Lawrence wrote:
> Hm, after some further reading (http://ppw.kuleuven.be/okp/_pdf/Tuerlinckx2004TECOT.pdf), it seems that the PDF is

> rather treacherous (see Eq 3), with parts that need numerical estimation on the fly. So possibly not a good project for
> the mathematically naive like myself :Op


> On Saturday, June 22, 2013 10:25:32 PM UTC-3, Andrew Gelman wrote:
>

> Mike,
> Is there a pdf describing the model that is being fit in Jags? I clicked through and found code but no description
> of the probability model that is being fit. If there is an algorithm for computing the probability model, it should
> be possible to fit it in Stan, I'd think.
> A
>
> On Jun 22, 2013, at 8:07 PM, Mike Lawrence wrote:
>
>> Hi folks,
>>
>> For decades it has been known that human response data from experiments involving simple perceptual
>> discriminations are very well fit by a "diffusion model" consisting of a Wiener process unfolding through time
>> between two bounds (reaching a bound causes the process to terminate, determining a response and the response
>> time). I'd like to use Stan to perform inference on a hierarchical version of such a model (aggregating data from
>> multiple human subjects, each measured hundreds of times). Apparently there already exists code in the form of a

>> JAGS module (http://sourceforge.net/projects/jags-wiener/ <http://sourceforge.net/projects/jags-wiener/>) that

[Andrew Gelman]

F

See below.  It looks like the Stan people are talking about your model!

A

On Jun 23, 2013, at 9:18 AM, Mike Lawrence wrote:

Hm, after some further reading (http://ppw.kuleuven.be/okp/_pdf/Tuerlinckx2004TECOT.pdf), it seems that the PDF is rather treacherous (see Eq 3), with parts that need numerical estimation on the fly. So possibly not a good project for the mathematically naive like myself :Op

On Saturday, June 22, 2013 10:25:32 PM UTC-3, Andrew Gelman wrote:

[Mike Lawrence]

Thanks Andrew and Bob for your responses. Apologies also: I should have looked at this a little more deeply before I bothered the list, as I now realize that the complexity of the diffusion model is beyond my rather meagre mathematical ability. I think I will instead focus my efforts on a recent alternative and increasingly substantiated model of timed discrimination responses that has a closed-form PDF (which, in turn relies only on the gaussian pdf and cdf). I'll probably need help in getting that one implemented, but I'll post such queries to a different thread...

[Peter Meilstrup]

I've also been puzzled by trying to implement this model in Stan. If it helps the discussion along, I have some R code to simulate and illustrate this kind of model:

https://gist.github.com/crowding/5846850

That Tuerlinckx paper is pretty heady. Here's a simpler introduction to the topic: http://www.shadlenlab.columbia.edu/~mike/papers/mine/okinawa_chapter2006.pdf

For most of the questions the mathematical psychologists have about such models (e.g. are bounds time-varying, is a single decision variable enough, how to deal with N>2 choices), finding closed form solutions for the likelihood function is a fool's errand. In practice, what people do is discretize the response time and the decision variable(s) over a grid, for each distinct type of trial, and then numerically sum out the Fokker-Planck equation from there to get likelihoods.

That's also the best idea I have for implementing this model in Stan, but I'm not sure if the gradient calculations can cope with it.

(Other possibilities are making the noise an explicit vector parameter for each trial, which would lead to an explosion of parameters.)

Peter Meilstrup

Physiology and Biophysics

University of Washington

On Sunday, June 23, 2013 6:18:23 AM UTC-7, Mike Lawrence wrote:

Hm, after some further reading (http://ppw.kuleuven.be/okp/_pdf/Tuerlinckx2004TECOT.pdf‎), it seems that the PDF is rather treacherous (see Eq 3), with parts that need numerical estimation on the fly. So possibly not a good project for the mathematically naive like myself :Op 

On Saturday, June 22, 2013 10:25:32 PM UTC-3, Andrew Gelman wrote:

[Andrew Gelman]

Hi, Joachim.  If/when you have this working in Stan, please let us all know.  This would be worth writing up.

Thanks

A

On Jun 24, 2013, at 3:58 AM, Joachim Vandekerckhove wrote:

I'm already talking to Mike Lawrence. I haven't tried to implement in Stan yet and I've been swamped with other stuff lately but will happily work with him (or others) to get it running because I'm pretty enthusiastic about Stan.
Given that pretty efficient/stable C++ code for the diffusion model exists I think it should be relatively easy.

Cheers,
Joachim

--
Joachim Vandekerckhove
Assistant Professor
Department of Cognitive Sciences
University of California, Irvine

Francis Tuerlinckx <Francis.T...@ppw.kuleuven.be> wrote:

hi andrew,

 

tell them that they can contact Joachim Vandekerckhove (in cc). he has c-coded PDFs for the simple wiener diffusion and maybe he already tried to link them to stan.

 

gr,

francis

 

From: Andrew Gelman [mailto:gel...@stat.columbia.edu]
Sent: maandag 24 juni 2013 1:58
To: Francis Tuerlinckx
Subject: Re: [stan-users] Translating a JAGS module for Stan?

 

F

More…

A

 

On Jun 23, 2013, at 7:57 PM, Peter Meilstrup wrote:

I've also been puzzled by trying to implement this model in Stan. If it helps the discussion along, I have some R code to simulate and illustrate this kind of model:

 

https://gist.github.com/crowding/5846850

 

That Tuerlinckx paper is pretty heady. Here's a simpler introduction to the topic: http://www.shadlenlab.columbia.edu/~mike/papers/mine/okinawa_chapter2006.pdf

 

For most of the questions the mathematical psychologists have about such models (e.g. are bounds time-varying, is a single decision variable enough, how to deal with N>2 choices), finding closed form solutions for the likelihood function is a fool's errand. In practice, what people do is discretize the response time and the decision variable(s) over a grid, for each distinct type of trial, and then numerically sum out the Fokker-Planck equation from there to get likelihoods.

 

That's also the best idea I have for implementing this model in Stan, but I'm not sure if the gradient calculations can cope with it.

 

(Other possibilities are making the noise an explicit vector parameter for each trial, which would lead to an explosion of parameters.)

 

Peter Meilstrup

Physiology and Biophysics

University of Washington

On Sunday, June 23, 2013 6:18:23 AM UTC-7, Mike Lawrence wrote:

Hm, after some further reading (http://ppw.kuleuven.be/okp/_pdf/Tuerlinckx2004TECOT.pdf), it seems that the PDF is rather treacherous (see Eq 3), with parts that need numerical estimation on the fly. So possibly not a good project for the mathematically naive like myself :Op 

On Saturday, June 22, 2013 10:25:32 PM UTC-3, Andrew Gelman wrote:

Mike,

Is there a pdf describing the model that is being fit in Jags?  I clicked through and found code but no description of the probability model that is being fit.  If there is an algorithm for computing the probability model, it should be possible to fit it in Stan, I'd think.

A

 

On Jun 22, 2013, at 8:07 PM, Mike Lawrence wrote:

Hi folks,

 

For decades it has been known that human response data from experiments involving simple perceptual discriminations are very well fit by a "diffusion model" consisting of a Wiener process unfolding through time between two bounds (reaching a bound causes the process to terminate, determining a response and the response time). I'd like to use Stan to perform inference on a hierarchical version of such a model (aggregating data from multiple human subjects, each measured hundreds of times). Apparently there already exists code in the form of a JAGS module (http://sourceforge.net/projects/jags-wiener/) that implements the sort of model I'd like to implement. Any advice on how I should go about translating this for use with Stan? Or are there tools already in Stan that could be used to implement this sort of model?

 

Mike

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

 

--
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/groups/opt_out.
 
 

 

[Joachim Vandekerckhove]

Hi all,

I had some time to look into it, and I now have a diffusion model working in Stan. I was able to use much of the code in the JAGS module, which was nice. I've barely tested this at all, though, so... caveat emptor. The implementation requires three files, all available here: http://www.cidlab.com/supp.php?o=wienerstan

Put the files in the right location, overwrite where needed, remake Stan, now you have  ~wiener(a,t,b,v).

If anyone needs help getting this to run, let me know.

I certainly wouldn't mind if one of the Stan developers would care to have a look at the wiener.hpp file I wrote. Suggestions to make it more style compliant are welcome.

Joachim

[Mike Lawrence]

Fantastic!

Am I correct in my understanding that this models the distribution of accurate response times "only" (quotes to indicate that this alone is still a great help!)? 

If so, and if we have the error RTs as well, would it be reasonable to let them help constrain inference by doing:

    accRTs ~ wiener(a,t,b,v)

    errRTs ~ wiener(a,t,-b,-v)

?

Mike

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

--

[Joachim Vandekerckhove]

No, it works just like the JAGS and WinBUGS ones, with negative values for error RTs (or lower bound hits). Very concise documentation is inside the wiener.hpp file.

What you wrote would probably also work, if you replace the -b with 1-b.

Joachim

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/29IcLbQ__uQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Mike Lawrence]

Oh, sorry, I should have considered looking at the docs for the jags version again to refresh my memory. Beautiful :)

[Bob Carpenter]

On 9/28/13 3:31 AM, Joachim Vandekerckhove wrote:

> I had some time to look into it, and I now have a diffusion model working in Stan. I was able to use much of the code in
> the JAGS module, which was nice.

Thanks!

> I've barely tested this at all, though, so... caveat emptor. The implementation
> requires three files, all available here: http://www.cidlab.com/supp.php?o=wienerstan
> Put the files in the right location, overwrite where needed, remake Stan, now you have ~wiener(a,t,b,v).
>
> If anyone needs help getting this to run, let me know.
>
> I certainly wouldn't mind if one of the Stan developers would care to have a look at the wiener.hpp file I wrote.
> Suggestions to make it more style compliant are welcome.

I'll absolutely take a look at it. Right now, I'm
having to prioritize finishing Stan 2.0, but we should
be able to get this incorporated for Stan 2.1, assuming
we can figure out the test cases for it.

I added an issue to remind us:

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

- Bob

[Joachim Vandekerckhove]

Awesome. I'd be happy to contribute to the testing/writing of doc.

In the meantime, do you think you could sneak in a "void add_quinary(const::std::string& name);" into ast.hpp and ast_def.hpp? The Wiener CDF function needs four input parameters.

- Bob

--
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/29IcLbQ__uQ/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Mike Lawrence]

Has anyone had any luck using this with either v2.0.0 or v2.0.1? I've tried both here on my mac but first get the message: 

/Library/Frameworks/R.framework/Versions/3.0/Resources/library/rstan/include//stansrc/stan/prob/distributions/univariate/continuous/wiener.hpp:126:68: error: no member named 'default_policy' in namespace 'stan::math'; did you mean 'boost::math::policies::default_policy'?

Then, after I modify wiener.hpp to replace all instances of "stan::math::default_policy" to "boost::math::policies::default_policy" as the error message suggests, I get the next error message:

/Library/Frameworks/R.framework/Versions/3.0/Resources/library/rstan/include//stansrc/stan/prob/distributions/univariate/continuous/wiener.hpp:64:12: error: no matching function for call to 'check_not_nan'

Any suggestions?

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

--

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+unsubscribe@googlegroups.com.

[Mike Lawrence]

On Joachim's suggestion, I was able to figure out from the policy refactoring commit back in August (http://goo.gl/yG9ug6) what needed to be changed in wiener.hpp. If anyone is interested, here's the diff:

https://gist.github.com/mike-lawrence/7866686/revisions

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[Bob Carpenter]

Thanks for sorting it out and contributing this back.

Do you have any interest in wrapping this up so we can add it
directly into Stan? Or is it something too specific for a general
tool?

I'm afraid we don't have a plug-in architecture for Stan yet
that makes it easy to just add in new code without merging it
into the Stan code base, which is a pain to maintain (even in
addition to Stan's interfaces changing out from under you).

- Bob

> requires three files, all available here: http://www.cidlab.com/supp.__php?o=wienerstan

> <http://www.cidlab.com/supp.php?o=wienerstan>
> Put the files in the right location, overwrite where needed, remake Stan, now you have ~wiener(a,t,b,v).
>
> If anyone needs help getting this to run, let me know.
>
> I certainly wouldn't mind if one of the Stan developers would care to have a look at the wiener.hpp file I
> wrote.
> Suggestions to make it more style compliant are welcome.
>
>
> I'll absolutely take a look at it. Right now, I'm
> having to prioritize finishing Stan 2.0, but we should
> be able to get this incorporated for Stan 2.1, assuming
> we can figure out the test cases for it.
>
> I added an issue to remind us:
>

> https://github.com/stan-dev/__stan/issues/245 <https://github.com/stan-dev/stan/issues/245>

>
> - 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+unsubscribe@__googlegroups.com <mailto:stan-users%2Bunsu...@googlegroups.com>.
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

>
>
>
> --
> 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.

[Joachim Vandekerckhove]

Thread necromancy!

I'm trying to plug this wiener distribution into Stan again, but I'm missing something.  I have the templated C++ file that implements wiener_log() (http://www.cidlab.com/p/wiener.hpp) but I can't find all the places I need to register it.

These are my steps:

(1) plunk wiener.hpp in ~/stan/cmdstan/stan/src/stan/prob/distributions/univariate/continuous/

(2) Add 

    add("wiener_log",DOUBLE_T,DOUBLE_T,DOUBLE_T,DOUBLE_T,DOUBLE_T,DOUBLE_T);

 to ~/stan/cmdstan/stan/src/stan/gm/function_signatures.h

(3) Add 

    #include <stan/prob/distributions/univariate/continuous/wiener.hpp>

~/stan/cmdstan/stan/src/stan/prob/distributions/univariate/continuous.hpp

(4) make bin/libstan.a && make -j2 bin/stanc && make bin/print

Step (4) runs without errors, but Stan still doesn't recognize wiener_log() or ~wiener().

Ideas welcome!

Joachim

[Ben Goodrich]

On Sunday, January 11, 2015 at 5:14:12 AM UTC-5, Joachim Vandekerckhove wrote:

(2) Add 

    add("wiener_log",DOUBLE_T,DOUBLE_T,DOUBLE_T,DOUBLE_T,DOUBLE_T,DOUBLE_T);

 to ~/stan/cmdstan/stan/src/stan/gm/function_signatures.h

My guess is that you need to add more function signatures. This one above is only for scalar arguments. For example, the normal has all these vector_types:

for (size_t i = 0; i < vector_types.size(); ++i)
  for (size_t j = 0; j < vector_types.size(); ++j)
    for (size_t k = 0; k < vector_types.size(); ++k) {
      add("normal_ccdf_log",
          DOUBLE_T, // result
          vector_types[i], vector_types[j], vector_types[k]); // args
      add("normal_cdf",
          DOUBLE_T, // result
          vector_types[i], vector_types[j], vector_types[k]); // args
      add("normal_cdf_log",
          DOUBLE_T, // result
          vector_types[i], vector_types[j], vector_types[k]); // args
      add("normal_log",
          DOUBLE_T, // result
          vector_types[i], vector_types[j], vector_types[k]); // args
    }
add_binary("normal_rng");

Ben

[Marco Inacio]

You can try using the (already patched) version from the branch feature/issue-765-wiener-distribution: (download https://github.com/stan-dev/stan/archive/feature/issue-765-wiener-distribution.zip, note that it's probably a bad idea since it's under development, and might not produce correct results) and has some documentation included.

I'll try to finish this soon, since there is considerable demand for that distribution. So maybe it'll be in Stan 2.6.

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.

[Joachim Vandekerckhove]

Oh, awesome, I didn't know it was there.  I had a quick look at the wiener.hpp and it needs some fixing, I will try to do that in the next couple of days.  This is extremely helpful, thanks!

--
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/29IcLbQ__uQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Kyle Payne]

Hi all,

Complete Newbie here.

I am currently trying to find a way to use the wiener model in rstan. Is there a way to use this fork of the distro, and compile it so that it works with R?

Thanks,

Kyle

>         stan-users+unsubscribe@__googlegroups.com <mailto:stan-users%2Bunsubscrib...@googlegroups.com>.

>         For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.
>
>
>
> --
> 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/groups/opt_out.

--
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.

[Jona Sassenhagen]

Seconding this request; DDMs are becoming increasingly popular in cognitive neuroscience and it would be amazing to have them readily available in a language as powerful as Stan.

[Bob Carpenter]

You can just check out the branch and then follow the R
build instructions.

It's passed unit testing and may wind up in our next
release (2.7) if it passes code review, which I'm about
to do now.

- Bob

>> > stan-users+unsubscribe@__googlegroups.com <mailto:stan-users%2Bunsu...@googlegroups.com>.

[Jona Sassenhagen]

That's great to hear.