Sampling time using beta_binomial of rstan version 2.19.2 and 2.19.3

Developers
1

Hi,

I’m having an issue using the beta_binomial distribution.

I wrote a stan model where the outcome is modelled by a beta_binomial.
With rstan 2.19.2 it runs without any problems.
Now, I updated rstan to 2.19.3 and I’m getting on same data and model a warning about low Tail ESS.

Therefore I minimized my model to check the beta_binomial and compared the result on both versions.

data {
    int<lower=0> N; 
    int y[N];
    int group[N];
}
parameters {
	real<lower=0> a[2];
	real<lower=0> b[2];
}
model {
	y ~ beta_binomial(60, a[group], b[group]);
	
	a ~ exponential(0.5);
	b ~ exponential(0.5);
}
generated quantities {
	vector[N] ppc_y;
	vector[N] log_lik;
	
	for(i in 1:N){
		log_lik[i] = beta_binomial_lpmf(y[i] | 60, a[group[i]], b[group[i]]);
		ppc_y[i] = beta_binomial_rng(60, a[group[i]], b[group[i]]);
	}
}

I’m using fake data.csv (16.4 KB) including 500 data points for two groups and sampling parameters are iter=2000 and chains=4.
The fit with version *.2 took ~10 sec per chain.
With version *.3 ~280 sec per chain and there are 4000 transition divergents, high R-hat values and low Bulk/Tail ESS.

Using the target += notation:

data {
    int<lower=0> N; 
    int y[N]; 
	int group[N];
}
parameters {
	real<lower=0> a[2];
	real<lower=0> b[2];
}
model {
	target += beta_binomial_lpmf(y |60, a[group], b[group]);
	a ~ exponential(0.5);
	b ~ exponential(0.5);
}
generated quantities {
	vector[N] ppc_y;
	vector[N] log_lik;
	
	for(i in 1:N){
		log_lik[i] = beta_binomial_lpmf(y[i] | 60, a[group[i]], b[group[i]]);
		ppc_y[i] = beta_binomial_rng(60, a[group[i]], b[group[i]]);
	}
}

The result is:
2.19.2 ~12 sec
2.19.3 ~ 16 sec

Have I done anything wrong or is there an issue with the sampling statement?
Is it right that the sampling using *.3 in comparison with *.2 could take longer even with the target notation?
And regarding my first issue, is it also possible to get low ESS using version 2.19.3 where everything is fine with 2.19.2?

R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17134)

2

My results for the original model with sampling statements is:

2.19.2: 8s/chain running 4 chains in parallel
2.19.3: 8s/chain running 4 chains in parallel

So no change. That’s on my new Macbook Pro.

Are you sure that Makevars didn’t change in the interim and you’re still running with -O3 optimization and -march=native? My ~/.R/Makevars file has this:

CXXFLAGS= -O3 -mtune=native $(LTO) -w

for flags.

If that’s not the issue, it sounds like some other component might not be installed properly, at which point I’d recommend posting here under the RStan tag or as an issue in the RStan repo (stan-dev/rstan).

P.S. We’ve usually found a parameterization of beta or beta-binomial in terms of mean plus count to be the easiest to manage, so that you’d have something like:

parameters {
  vector<lower = 0, upper = 1>[2] theta;
  vector<lower = 0>[2] kappa;
  ...
transformed parameters {
  vector[2] a = theta .* kappa;
  vector[2] b = (1 - theta) .* kappa;
  ...
model {
  kappa ~ exponential(1);
  theta ~ beta(2, 2);  // or even 1, 1 if you'd like to keep it uniform

It let’s the count (kappa) vary independently of the mean (theta) in the prior, whereas the a, b parameterization tightly couples them to define the mean (which is usually way better characterized than the variance in a beta).