How much to log in a custom made pdf (for the linear ballistic accumulator)

Modeling
,
1

Hi,
I’m trying to implement a custom made for the linear ballistic accumulator with a drift from a normal+, to fit responses and reaction times. I’m taking the pdfs and cdf of each accumulator from here single-lba.r (dlba_norm_core and plba_norm_core), and then basically the PDF of the model is (for two accumulator) pdf_r * (1-cdf_n), where r is the accumulator that finished first (and “gave” the answer.)

So this is (a simplification) of the code for the pdf of one accumulator in R :

eta_i is real, b>0, A >b, s_i>0 , t>0

pdf <- function(t=.3, A=1, b=1.4, eta_i=3.5, s_i=1){
  Phi <- function(x) pnorm(x,0,1) 
  npdf <- function(x) dnorm(x,0,1)
  denom <- Phi(eta_i/s_i)
  zs <- t*s_i
  zu <- t*eta_i
  chiminuszu <- b-zu
   chizu <- chiminuszu/zs
   chizumax <- (chiminuszu-A)/zs
   pdf <- (eta_i*(Phi(chizu)-Phi(chizumax)) + s_i*(npdf(chizumax)-npdf(chizu)))/(A*denom)
   return(pdf)
}

I need to get the log(pdf), so my question is how “inside” should I have the logs. There are many things that can under/overflow very quickly, if eta_i is very large positive or negative, if s_i is very small or of b is very large. My first instinct was too transform everything into a sum of terms minus another sum of terms, and then
the lpdf becomes log_diff_exp(log_sum_exp(logterms1),log_sum_exp(logterms2)) - log(A) - normal_lcdf(eta_i/s_i|0,1)

But then the log of the Phi gets to 0 or to -Inf and the eta_i can be negative during the sampling, and then I have a lot of if clauses to catch all these cases and avoid NaNs.
My lpdf of the accumulator ends up being horribly long, and even worst for the lccdf. Am I on the right track? My custom made pdfs and cdfs seem to give the same answers as the package single-lba.r, except for extreme values.

But I was wondering if what I did makes sense and if there is a way to optimize my code. (I want to make this hierarchical as a next step).

A second question is when (or if) I should use log(Phi_approx(.)) instead of the normal_lcdf (that are inside the lpdf and lccdf of the acummulator), (when should I use this function if at all?)

I’m attaching the code just in case. tLBA.stan (6.8 KB)

And this is to test the model and to test the lpdf and lccdf: test_tnormLBA.R (4.6 KB)

Bests,
Bruno

1 Like
2

Less ugly (and working version) of the code here:tLBA.stan (6.2 KB)

But all my questions are still relevant

1 Like
3

You want to do as much as possible on the log scale, so I think your first instict is right.

We don’t have good examples in the code as we’re using standard CDF and CCDF algorithms and then just naively taking logs.

It uses a logistic approximation that’s much more efficient to calculate and pretty close. So I’d try it both ways and see if you get the same answer and then use Phi_approx where possible. That’d give you more opportunites to roll in some log terms.

4

Thanks for the answers.

Why would Phi_approx allow me to roll in more log terms? normal_lcdf is already logged, isn’t it?
And related to this, would Phi_approx allow me to go more negative than -37.5 and more positive than 8.25 without under/overflowing? The manual is not very clear in this, and it seems not…

And since I have you here, is there a way to approximate log(Phi(b) - Phi(b-a)) to avoid over/underflow, this is the most problematic aspect of my pdf.

5

The implementation of normal_lcdf just applies log() to normal_cdf, so it doesn’t operate on the log scale. The logistic form is simple enough that you could easily implement it stably on the log scale.

6

Ok, I understand, and any insight about log(Phi(b) - Phi(b-a)) ? Or I can’t do more than log(Phi_approx(b) - Phi_approx(b-a))? (Or this is probably more for stackexchange than for here).

7

I’d try substituting the definition for Phi_approx into

log(Phi_approx(b) - Phi_approx(b-a))

and see if you can reduce it at all to keep everything on the log scale. We have a log_diff_exp function that can help, as this will be:

log_diff_exp(log(Phi_approx(b)), log(Phi_approx(b - a)))

so then it’s just a matter of figuring out if log(Phi_approx(u)) can be reduced to something more stable.