Thankfully @aammd had a reproducible example. I ran it (with the seeds and number of iterations) changed and with 4 chains and the results are alarming.
CmdStanR gave
> fit$summary()
# A tibble: 603 x 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 a 0.523 0.570 0.203 0.223 0.213 0.760 3.84 4.37 14.9
2 theta 1.71 1.11 1.59 1.32 0.130 4.47 3.44 4.42 17.3
It aslo has low EBFMI and 90% of the iterations hitting a maxtreedepth warning.
RStan gave
> print(s,pars = c("a","theta"))
Inference for Stan model: betabin.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a 0.31 0.00 0.01 0.29 0.30 0.31 0.31 0.32 3024 1
theta 19.37 0.03 1.56 16.36 18.33 19.37 20.37 22.57 2730 1
(True paramters are a=0.3, theta=40)
This looks a great deal like a bug or a performance reversion. Rstan also sampled WAAAAAY faster than CmdStanR, which is not a good sign!
I obviously have no idea what’s wrong so I’m just gonna tag people. @stevebronder was the last person to touch the beta_binomial C++ code so maybe he has a clue. It also might be an algorithm change between. It might be a change in the algorithm since 2.19 (@betanalpha). It might be a weird cmdstan thing (@mitzimorris @jgabry) It might also be a stanc3 thing (@seantalts).
Thanks for the report Andrew!