Likelihood with infinite sum - RAM goes through the roof

Modeling
1

Hi everyone,

TL;DR

  • I have likelihood function that has an infinite sum;

  • This infinite sum needs to be truncated;

  • I use a while loop and commonly used tolerance rule to do this;

  • Can’t sample from the model and RAM goes through the roof.

I am trying to write a function to compute the density of the first passage time of a Wiener process over one of two absorbing boundaries.

There is already a function for this in Stan (wiener_lpdf), however, I am trying an alternative that does not involve too many decisions between different forms to compute the density of the first passage time.

My goal is to understand how well I can fit this model to response times and binary choice data from psychological experiments. Also to see how it compares to the already implemented version in Stan, how well (or not) it scales when I have data from more people, trials, etc.

To achieve this I am taking the so called large time representation from Navarro and Fuss (2009, equation 1):

WFPT( t | \delta, \alpha, \beta)=\frac{\pi}{\alpha^{2}} \exp \left(-\delta \alpha \beta-\frac{\delta^{2} t}{2}\right) \sum_{k=1}^{\infty} k \exp \left(-\frac{k^{2} \pi^{2} t}{2 \alpha^{2}}\right) \sin (k \pi \beta)

Where t is a shifted response time (t = RT - \tau, where \tau is a shift parameter), \alpha the distance between lower and upper boundaries, \beta relative starting point of evidence accumulation and \delta the drift rate of the process.

This function computes the density for when the process has terminated in the lower boundary. For the upper boundary one has to replace \beta and \delta by 1 - \beta and -\delta respectively.

One of the challenges to compute the likelihood from this function is that we have an infinite sum that needs to be truncated.

One way this can be done is by checking at each iteration of the sum (from \kappa = 1, ... , \infty), if the current term and previous term are both less than a prespecified tolerance. (for example, Ratcliff and Tuerlinckx (2002, appendix B) tolerance is 10^{-29} \times the sum at step \kappa).

I have tried to implement this rule using a while loop in Stan, however, for a very simple model and data, sampling just requires too much RAM and does not do any sampling. I tried both rstan and CmdStan and the issue remains.

Bellow you’ll find my Stan code and you can reproduce the problem in R using https://bitbucket.org/snippets/tcabaco/qnBngo .

Thank you for your help!

functions {
  real wfpt_lpdf (real t, real alpha, real beta, real delta) {
    real log_lik = 0;
    real wiener_error = 10 ^ -29;
    // make sure the parameterization match with Stan's. 
    real T_beta = 1 - beta;
    real T_delta = - delta;
    
    real temp1 = log(pi()) - 2 * log(alpha);
    real temp2 = -(T_delta * alpha * T_beta) - 0.5 * (square(T_delta) * t);
    
/*
* Truncating infinite sum using a tolerance rule
* From Ratcliff and Tuerlickx (2002)
*
* given k steps in a sum from k = 1 to infinity
*
* we stop the sum when the sum at steps k and k-1 are both 
* bellow 10^-29 times the current sum at step k
*
* I should correct this asap
*/
    real temp3 = 0; // cumulative sum

    {
      real mem1 = 1; // sum at step k
      real mem2 = 1; // sum at step k - 1
      real temp_sum; // temporary sum at step k
      int k = 1;
      
      while(mem1 >= temp3 * wiener_error && mem2 >= temp3 * wiener_error) {
        temp_sum = k * 
                   exp( 
                        - (square(k) * square(pi()) * t) / 
                        (2 * square(alpha))
                      ) * 
                   sin (k * pi() * T_beta);
        temp3 += temp_sum;
        mem2 = mem1;
        mem1 = temp_sum;
        k += 1;
      }
    }
    
    log_lik += temp1 + temp2 + log(temp3);
    
    return log_lik;
  }  
}
3 Likes
2

My first guess is that your approximation doesn’t converge, so the termination condition of your while loop never triggers.Have you tried debugging this function outside of Stan, say in R? It would be easier to double check that your coding does what’s expected, such as monitoring convergence.

If you are sticking to Stan, I’d suggest setting a maximum k and changing the while loop to:

while(mem1 >= temp3 * wiener_error && mem2 >= temp3 * wiener_error && k < max_k)

You could then try to make it print your variables if max_k is reached.

6 Likes
3

Thank you for the response.

Yes, after adding a print statement to see the number of k terms I noticed that it was going over to the hundreds of thousands.

I then did some testing and I got the model to converge by setting good initial values.

“Converge” here is a strong statement, but it is sampling at least and now I can start evaluating the likelihood implementation.

Thanks again!

1 Like