Would appreciate your clarification of the following issues that seem confusing.
In the documentation (eg, p.29 of cmdstan-guide-2.19.1.pdf), the reported divergent__ row has the NEff equal to the 1000 iternations saved (but in my case below, I have -nan) and also Rhat = nan (my case: 9.99e-01). Are these NEff and Rhat, in either case, meaningful at all?
The Mean in my case is non-zero but very small (5.00e-03). Is this merely due to trying to report an integer value in floating point? (In other words, that reported value should be taken as 0.)
Don’t worry about Rhat for the divergent__ parameter, that’s for the parameters mostly.
5e-3 is the fraction of post-warmup trajectories that were divergent in all chains, so there were 5e-3 * 1600 = 8 divergences across all four chains.
This model looks like it took way longer than the previous one to run? Nearly 12 hours four 1600 draws vs. 7 hours for 4000 draws? And the treedepth is way larger and stepsize way smaller. Was this with more data or something?
Was the adapt_delta higher in this run? I guess if it has to be higher to keep divergences down then there hasn’t been any real advantage to switching to the dense adaptation.
The warning divergences give is that you might not be sampling the distribution you think you are. I guess since you have results from the diagonal adaptation which isn’t producing divergences, you could make predictions with each of them with example data and see if there are differences.
If you can check your results against the model run with diagonal adaptation. Rhat occasionally ends up less than 1.0. I don’t think what you’re seeing in your results here it’s a thing to worry about.
I have been assuming that the 1.01 threshold for Rhat you mentioned before is equivalent to a 0.99 threshold. So only the deviation from 1 matters, right?
Thanks for the clarification and additional information. I need to look into the issues involved more carefully and systematically. Let me simply reply to your question about the slow run time: I suspect that my changing the p[3]'s true parameter value resulted in divergences, which typically would lead to a much slower run time.
The p[3] of my model has been hard to estimate reliably – usually overestimating the true value (0.35) to 0.43 (but other than that, the run time then seemed all right and I didn’t have divergences). I wanted to understand to what extent the hard-to-estimate p[3] is driven by my prior (with a mean at 0.5). So in the recent runs, I set the true value of p[3] to 0.68 (above 0.5) to see whether the estimate would fall somehwere in between, which I believe would suggest the p[3] cannot be well identified by the data and thus the estimate is always pulled toward the prior’s mean (0.5). p[1:3] are the parameters of a sigmoid function modeled with inv_logit(). I guess the nonlinearity might make it hard to be accurately estimated.
I think Rhat is supposed to be greater than or equal to 1.0. I think it being less than one just happens randomly, but isn’t a big deal. So less than or greater than one are interpreted differently.
Only suggestion I have is maybe generate a dataset that doesn’t take 12 hours to run haha. It’s gonna be hard to make much headway with calculations that slow.
Thanks for your clarification and suggestion. Make a lot of sense. Owing to some other considerations, for now I still use 10,000 simulated data points.
I am returning to more conservative choices of warmup (1000) and delta (0.99) to first ensure no divergence. Just finished a comparison of the dense vs diag runs on the slower cluster I have access to. I believe the results confirm that it’s fine (ie, similar results) and benefitial (ie, shorter run time) to use the dense metric.
Would appreciate your clarification of one thing though. In the diag results below, the rows from stepsize to n_leapfrog have MCSE and N_Eff reported as -nan, which is different from those reported in the dense results. Is this difference something to worry about by me as a model builder?
Nah, don’t worry about the stats for stepsize or n_leapfrog. It’s just the distribution of stepsizes and number of leapfrog steps used in each trajectory. NaNs just indicate every HMC trajectory used the same stepsize and n_leapfrog in the diagonal case.
I dunno why the MCSE for dense isn’t NaNs too for stepsize, cause that should be constant. I’d assume if you did the stansummary over 1 chain it still would be.
Warmup took (21000, 20885, 20857, 22902) seconds, 23.8 hours total
Sampling took (7961, 4248, 7934, 8000) seconds, 7.82 hours total
Ouch these warmup times are rough. You might benefit from doing a calculation with save_warmup=1 and examining what’s happening at warmup. I’m not sure what you might find but maybe there’s something there.
Thanks for the suggestions in your reply to my another post. I am working on the plotting as you suggest.
Meanwhile, I have a look at the effect of changing the max_depth and realize that setting it to the 9 or 10 that I had been using leads to 600 or 121 transitions hitting the limit, resp. (when running 3 x 200 = 600 transitions). Increasing to max_depth=11 or above avoids hitting the limit completely.
The max_depth=9 case also has the issue of split R-hat greater than 1.1: gw_sd[3], gw_sd[5]. But all the parameters of interest (reported below) have R-hat ~ 1.0 and the run time is significantly shorter than in the max_depth=10 case (where no R-hat > 1.1 at all).
Should R-hat > 1.1 be avoided at all cost even when(i) the problem does not affect the good Rhat of the parameters of interest; (ii) their mean estimates remain very similar; (iii) ignoring the issue gives a much shorter run time ?
I have looked at max_depth=11 to 13 as well. The warmup time increases quite substantially even though the final sampling time is very similar and so are the mean estimates, N_Eff, and Rhat of the parameters of interest.
Is it still worth trying much bigger max_depth, like 20, 30, …?
Had a look at the CmdStan manual concerning init_buffer, term_buffer, and window, and the adaptation process. I guess I have a basic understanding. Tried to look around (like the Stan user guide, forum, …) and found only very limited related info (such as this post). Another discussion by @aaronjg sounds very relevant but most of the responses in that post are at a level too technical for me to comprehend.
My take from that post is that I can reduce the overall warmup time by forcing more iterations into the fast init_buffer and term_buffer intervals (eg, 200, instead of the default 75). Am I getting this right? Is it useful to adjust the stepsize at the same time and if so, how?
Would appreciate very much your feedback. Many thanks.
Yes, unfortunately. You might look at what the different chains are doing, or what your estimates of these parameters are. You might be able to get around this with tighter priors if you can justify that in your modeling.
No, if your Rhats are good at 11 and you’re not getting treedepth warnings, no need to make the max treedepth deeper. Also this is a log scale sorta thing. Max treedepth = 10 means that NUTS will use at most 2^10 leapfrog steps to generate a sample. 11 means NUTS will use at most 2^11 leapfrog steps, etc.
It’s only a guess. It may or may not get you anywhere, but if adaptation starts before you’re near the solution then the adaptation itself might screw things up. I would just leave the stepsize adaptation at default and first try making the init_buffer larger.
You can get the stepsize from the output csv files. If you make a plot of that you can search to see if anything weird happens (for instance, the stepsize plummets to 1e-10 in part of warmup). It looks like for the max_treedepth = 11 chain, the final stepsize you take is 1e-2, so that’s what to compare things against. I haven’t really looked into this much, but I’d guess you should stay within a couple orders of magnitude of this through the rest of adaptation. I could easily be wrong there though – just a guess.
Would appreciate feedback on whether the following weird experience is a bug or my mistake.
I cannot get the warmups into the csv files (eg, samples1.csv (1.4 MB) ) even when I used the save_warmup=1 option and the stansummary report cmdStan0.5.2a_bkgd.sh.o6080738.txt (132.9 KB) suggests that this was successfully done:
3 chains: each with iter=(6,6,6); warmup=(8,8,8); thin=(1,1,1); 42 iterations saved.
...
...
Warning: non-fatal error reading adapation data
Warning: non-fatal error reading adapation data
Warning: non-fatal error reading adapation data
2 of 42 (4.8%) transitions ended with a divergence.
In the shell script (cmdStan0.5.2a_bkgd.sh.txt (3.2 KB) ), I previously used an environment variable but now hard-coded the save_warmup=1 option. This still cannot put the warmups into the csv files:
for i in {1..3}
do
nohup ../SSM0.5.2a sample algorithm=hmc engine=nuts max_depth=$MAXDEPTH metric=$METRIC num_samples=$N_ITER num_warmup=$N_WARMUP save_warmup=1 adapt delta=$DELTA init_buffer=$INIT_BUFFER term_buffer=$TERM_BUFFER window=$WINDOW random seed=$SEED id=$i data file=../SSM_data_J90d.dat output file=/scratch/city/${NOW}_${JOB_ID}/samples$i.csv > $HOME/$SHELLSCRIPT.${NOW}_${JOB_ID}.$i.out 2> $HOME/$SHELLSCRIPT.${NOW}_${JOB_ID}.$i.err &
done
jobs
wait
The above was for cmdstan 2.19.1. For rstan 2.18.2, I also tried save_warmup = TRUE (or omitting save_warmup, whose default should be TRUE) to obtain the model fit. The attempt to plot with mcmc_trace() also shows no warmup at all. (Although here, I don’t have the csv files and cannot look into them to confirm that the option did not really save the warmups.)
Lines 40-47 in samples.csv are warmup samples. Everything before the “# Adaptation terminated”. You can parse this input using read_stan_csv in R or do it manually with read_csv or whatnot.
You have an incredibly large output. If you have things in the transformed parameters block you might move them to the model block so they don’t print.
After comparing many combinations of the tuning parameters (warmup, iter, window , term_buffer , max_treedepth , metric=dense_e vs. diag_e , adapt_delta), the following combinations seem to give the shortest run time while keeping the sampling results at an acceptable level:
For a given amount of simulated data points (9K or 5.4K), diag_e runs a bit faster than dense_e mainly because the former can take a lower level of adapt_delta (0.82, instead of 0.982) without leading to divergences.
But as far as the parameters of interest are concerned, the N_Effs from diag_e seem to be more diverse, with the lowest being acceptable but substantially lower than the lowest of those from dense_e.
The default choices of window, init_buffer, and term_buffer (=25, 75, 50) appear to be quite good: alternative choices do not lead to dramatic improvements, especially when the sample size is smaller (ie, 5.4K, rather than 90K).
(window, term_buffer) = (35, 150) performs slightly better. Not sure of the reason but my guess is that a larger window (35, rather than 25) reduces the iterations spent on the last and slowiest stage of the adaptation process, where further search for improvement has little incremental gains. On the other hand, a larger term_buffer (150, rather than 50) allocates more iternations to optimally adjust the initial typical set of estimates according to the findings from the adaptation process.
Despite the rationalization given above, increasing (window, term_buffer) beyond (35, 150) would not necessarily improve further owing to the complicated adaptation process. There is simply no obvious relation between the run time and these tuning parameters. Apparently, they interact in a complicated way with the warmup and iter to determine the run time.
During the course of my comparison, I have the impression that choosing an unnecessarily high level of adapt_delta often results in hitting the max_treedepth (many such instances can lead to a rather long run time). So to avoid divergences, it’s about finding the suitable level of adapt_delta rather than the higher the better.
This makes me think dense would still be worth doing. But it sounds like it’s been pretty annoying to get the dense metric to run without divergences, which is no good.
Thanks for taking the time to report this stuff back!
Divergences imply that your model features strong, possibly spatially-varying, curvature that compromises the symplectic integration. Increasing adapt_deltamight help, but at the cost of more expensive transitions. Really the ideal is to identify exactly which features of your model are inducing the strong curvature (often strong correlations from weak or non identifiabilities between parameters in the likelihood function) and remedy them with modification to your model implementation (reparameterization, rescaling of data, stronger priors, additional data, etc.
In other words increasing adapt_delta is a bandaid – it might help for paper cuts but slapping on more bandaid when you’re hemorrhaging isn’t going to be of much use to anyone.