Sampling does not converge with default NUTS sampler

[jCalderTravis]

Thank you to the developers for a very exciting new code package, and for all the associated documentation.

We have experienced similar convergence issues to those that have been discussed elsewhere (#306, #323, #413, #471). If I understand correctly, these are likely caused by issues with the "t" (non-decision time) parameter, when the value of the t parameter is close to the value of a response time (RT) observed in the dataset.

In the linked discussions the following solutions have been proposed:

1. Narrowing the priors, in particular those associated with t (Large number of divergences in simple model #306, Tips for converging variance estimates in hierarchical models #413, hierarchical DDM does not seem to converge at all on 0.2.2 #471)
2. Set a very small initial value for the t chains (Large number of divergences in simple model #306, non-decision time does not converge #323, hierarchical DDM does not seem to converge at all on 0.2.2 #471)
3. Switching to a slice sampler from the default gradient-based sampler (non-decision time does not converge #323, Tips for converging variance estimates in hierarchical models #413)
4. Switch to the approximate likelihoods (Trouble Getting Hierarchical Model to Work with All Parameters Hierarchical #411)
5. Adding a constant to all RTs observed in the dataset (hierarchical DDM does not seem to converge at all on 0.2.2 #471)

My collaborator has found that the slice sampler (option 3) works very well, although unfortunately the slice sampler is prohibitively slow for hierarchical models (an open issue; #388).

Regarding setting initial values for the t parameter to very small values, I am unsure about this approach for two reasons. First, I wondered whether this could bias the Gelman-Rubin statistic: If all chains start at the same place then even if sampling fails, between-chain variance will be very small making it appear as if convergence has been achieved. Second, real non-decision time values can be quite high, depending on the task. For example, van den Berg et al. (2016) report non-decision times of about 400ms. Certainly real datasets could contain lapse trials with RTs shorter than this. If the sampling will encounter issues as soon as the non-decision time parameter approaches the shortest RTs, and the chain is initialised with a very low value for t, then the true value for t will likely never be achieved.

I wondered whether these concerns with setting the initial value for t are legitimate, or whether I might have misunderstood something. I would be grateful for any thoughts.

[AlexanderFengler]

Hi @jCalderTravis ,

the issues we are experiencing with the t parameter are somewhat specific to the analytical likelihood and, to the best of our current understanding, they stem from gradient instabilities when the min_rt ~ t.

Our current attempt to mitigate this is by starting at very small parameter values t and "hoping" that the sampler will not go traverse these "dangerous" reasons because they should a priori be relatively unlikely.

Simultaneously, we are trying to figure out how to improve the likelihood formulation to have more robust gradients in those regions.

The suggestion to use the slice sampler, derives from the insight that it is really the gradient, not the basic likelihood evaluation, that is a problem here. Slice samplers don't use gradients (making the slice sampler worse in certain respects ofc including speed, otherwise we wouldn't care about gradients), so shouldn't suffer the same fate.

Your concerns about equal starting points have some merit, let me share a few thoughts below.

1. Equal starting points is just one of two level of stochasticity that is applied through the sampler, so the chains are otherwise not the same.
2. Sampling across starting points is mostly valuable (it is not uniquely able to identify it) if you suspect the posterior to be multimodal, which is, in my experience, somewhat rate, especially for the t parameter, but can be a concern. Even if you start with the same starting points, you are in principle able to find multiple modes of course, but there can certainly be cases were it is just to unlikely to reach a particular mode from a given starting point, even though that mode has significant posterior probability. I want to mention that you can de facto set your custom starting points in multiple places, and override those defaults, so it's not a hard limit we are imposing, just our current best attempt at defaults. First, you can do so when you define your own priors for any given parameter, second you can pass the initvals argument to the .sample() function.
3. Concerning your point regarding early lapse trials, I agree with you that this can be an issue. As I said we are trying to make the likelihood formulation more robust to this. Also let me mention again that the approx_differentiable likelihoods tend to be a bit more robust to this issue already.

Best,
Alex

[jCalderTravis]

Hi @AlexanderFengler,

Thank you for your reply and your thoughts. I would like to stick to non-approximate likelihoods, so I appreciate that you are working on improving the formulation of the analytical likelihood function, and I will make sure to follow HSSM developments. For now we will think carefully about the priors we use and the option of customising the initial starting points, as you suggest.

Best,
Josh