Hi all,

I’ve specified and validated a hierarchical GPR model. The validation involved running the model on simulated data that is similar in nature to experimental data I plan to analyze in the future.

My issue is with computation time and sampling efficiency. I cannot efficiently sample from the posterior without setting max_treedepth and adapt_delta to something like 15 and 0.99 respectively. Of course, these control settings drastically increase the amount of time it takes the to sample the posterior.

I am looking for some guidance regarding prior specifications on the GP hyperparameters that might result in sampling improvements. Or, even more generally, some procedural advice about how to diagnose issues what might be the issue would be greatly appreciated.

Currently, I specify a Weibull on the amplitude (variance) hyperparameter and a Cauchy on the volatility (the inverse of the lengthscale) hyperparameter.

I’ve played around with using a normal instead of a cauchy and with making the prior more informative but that does not appear to improve the sampling efficiency.

I’ve attached my stan code. Also, I’ve attached a Kruschke style diagram of my model for what it’s worth.

Also, I’ve attached pairs plots for a case in which some parameters are sample well and for a case in which some parameters are not sampled well (exceeding maximum treedepth).

I have a bundle of scripts and data that can reproduce my issues that I can link to (Google Drive of something) if necessary.

Thanks,

Please find attached gp_regression.stan (7.1 KB)