Rstan works but cmdstan does not for a beta binomial model

Interfaces CmdStan
1

I am trying to fit a betabinomial model in cmdstan and am finding problems relative to rstan, which works well.

Here is the model, which i store in a file called simple_betabin.stan:

data{
    int W;
    int F[W];
    int R[W];
}
parameters{
    real<lower=0,upper=1> a;
    real<lower=0> theta;
}
model{
    vector[W] pbar;
    theta ~ exponential( 3 );
    a ~ beta( 3 , 7 );
    R ~ beta_binomial(F ,  a * theta, (1 - a) * theta );
}
generated quantities{
    vector[W] log_lik;
    vector[W] y_hat;
    for ( i in 1:W ) {
        log_lik[i] = beta_binomial_lpmf( R[i] | F[i] , a * theta, (1 - a) * theta  );
        y_hat[i] = beta_binomial_rng(F[i] , a * theta, (1 - a) * theta );
    }
}

If I fit this model in cmdstan, it fits very slowly and also gives a very poor result: Rhat of 2.2, ESS of 2.5
On the other hand, fitting in Rstan is faster and works better, with a perfect Rhat and plenty of effective samples.

library(cmdstanr)
mod <- cmdstan_model(stan_file = "code/simple_betabin.stan")

## generate data
probs <- rbeta(300, 0.3 * 40, 0.7 * 40)
trials <- rep(450, 300)
res <- rbinom(300, trials, probs)

d <- list(W = 300, F = trials, R = res)

fit <- mod$sample(
  data = d, 
  seed = 123, 
  num_chains = 2, 
  num_cores = 2
)

fit$summary()

## compare rstan

m <- rstan::stan_model(file = "code/simple_betabin.stan")

s <- rstan::sampling(m, data = d, iter = 1000, chains = 2, cores = 2)
  • Operating System: MacOS X Mojave 10.14.6
  • CmdStan Version: 0.0.0.9
  • Compiler/Toolkit:
  • rstan version: 2.19.2
2 Likes
2

I noticed that CmdStanR and RStan are not being used with quite the same settings:

  • The seed is set with CmdStanR (to 123), but isn’t explicitly set to anything with RStan. Sometimes models will seem well-behaved with some seeds but show pathologies with others.

  • Judging from Getting started with CmdStanR, the default total number of iterations is 2000, with half of those being for warmup. Your call to rstan::sampling sets the total number of iterations (again, including warmup!) to only 1000.

Both these issues make it difficult to make an apples-to-apples comparison.

1 Like
3

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!

2 Likes
4

Yikes! I just put up a revert to the changes I made in beta-binomial.

I ran the example on the revert and it gives much better answers

                    Mean     MCSE   StdDev        5%       50%       95%    N_Eff  N_Eff/s    R_hat
a                3.0e-01  1.5e-04  5.7e-03   2.9e-01   3.0e-01   3.1e-01  1.6e+03  1.7e+03  1.0e+00
theta            2.3e+01  4.7e-02  1.9e+00   2.0e+01   2.3e+01   2.6e+01  1.6e+03  1.7e+03  1.0e+00
lp__            -8.1e+04  3.2e-02  1.0e+00  -8.1e+04  -8.1e+04  -8.1e+04  9.7e+02  1.1e+03  1.0e+00
7 Likes
5

Wooooooo!! Thanks Steve! And thanks again @aammd for finding it!!

1 Like
6

I’m gonna guess it’s in stanc3 based mostly on the severity - I can take a look tomorrow. Definitely need to get to the bottom of this before the next release mid Jan. Thanks for replicating and tagging me :)

1 Like
7

Turns out it was a change to the function between 2.19 and now. Stanc3 continues to be fab!

2 Likes
8

Thank you everyone! Both for the quick resolution and for creating this tool in the first place!

Could someone give me a hint how to have the version including @stevebronder 's change? should I wait for another release?

One more question since I have all of your ears: why is theta estimated to be half its original value?

9

This bugfix will be included in the next cmdstan release, which will probably be available on the 18th of January.

You can get it now (well technically in a few hours) by cloning cmdstan directly from github by running:

git clone --recursive https://github.com/stan-dev/cmdstan.git

If you are using cmdstanr you can also install it using install_cmdstan(repo_clone = TRUE).

The fix was merged to the math library but it still needs to pass some tests before it is available in Cmdstan (hence the need to wait a few more hours).

10

That’s because of your exponential(3) prior on theta, which has a mean of 1/3, pulling your estimate to the left! Hopefully a prior-posterior plot will show it’s bad!

1 Like
11

Thanks once again @anon75146577! I changed the model block to

model{
    vector[W] pbar;
    theta ~ normal( 3,0.5 );
    a ~ beta( 3 , 7 );
    R ~ beta_binomial(F ,  a * exp(theta), (1 - a) * exp(theta) );
}

And it works beautifully!
(EDIT that is, it estimates it as between 37 and 42)

1 Like
12

Hi Stan team,

Having some trouble with some code that I swear used to work fine and ran across this thread. I tried the test code provided above (copied at bottom of this post for convenience.)

I recently got a new computer, so unfortunately I don’t know what conditions I was running under before. Code is running /slow/ and constantly exceeds maximum treedepth. I just installed rstan from CRAN (version 2.13.3).

It sounds like it could be the same bug?

library(rstan)
bb.mod <- "data{
  int W;
  int F[W];
  int R[W];
}
parameters{
  real<lower=0,upper=1> a;
  real<lower=0> theta;
}
model{
    vector[W] pbar;
    theta ~ normal( 3,0.5 );
    a ~ beta( 3 , 7 );
    R ~ beta_binomial(F ,  a * exp(theta), (1 - a) * exp(theta) );
}"

## generate data
probs <- rbeta(300, 0.3 * 40, 0.7 * 40)
trials <- rep(450, 300)
res <- rbinom(300, trials, probs)

d <- list(W = 300, F = trials, R = res)



## compare rstan

m <- stan_model(model_code=bb.mod)

s <- sampling(m, data = d, iter = 1000, chains = 2, cores = 2)

Warning messages:
1: There were 941 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded 
2: Examine the pairs() plot to diagnose sampling problems
 
3: The largest R-hat is 2.34, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat 
4: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess 
5: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess 
> packageVersion("rstan")
[1] ‘2.19.3’
13

Do you mean 2.19.3?

15

Yes. That was a typo, sorry.
> packageVersion(“rstan”)
[1] ‘2.19.3’

16

Can you see if you have the same problem using cmdstanr?

17

I was finally able to get cmdstan/cmdstanr up and running on my Windows machine. With cmdstan it runs fast (~1 sec vs ~1 min with rstan) and diagnostics show no problem with convergence.)

cmdstan output:

> fit$summary()
# A tibble: 3 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 lp__     -81249.    -81248.    1.08    0.741   -81251     -81248.     1.00     905.      NA 
2 a             0.302      0.302 0.00462 0.00471      0.295      0.309  1.00    1634.    1296.
3 theta         3.59       3.59  0.0914  0.0944       3.43       3.73   1.00    1657.    1165.

rstan output:

> summary(s)$summary
             mean   se_mean        sd        2.5%         25%         50%        75%      97.5%
a       0.6329534 0.2448504 0.2449837   0.3814097   0.3886735   0.6351125  0.8778146  0.8801354
theta   1.6253428 0.7332779 0.7337782   0.8736830   0.8923865   1.6099964  2.3577879  2.3985647
lp__  -15.1753221 8.9313636 8.9359453 -24.1864680 -24.1042734 -15.1780398 -6.2443677 -6.1329265
         n_eff      Rhat
a     1.001089 120.85057
theta 1.001365  73.14958
lp__  1.001026 224.31933
18

Why are the results so radically different?

19

@maxbiostat: if you mean relative to the OP’s results, I used the change mid-thread where theta -> log(theta). The cmdstan results for a and lp__ are very similar to what the OP reported as good results. The point of my question is what is going on with rstan for this model; certainly the rstan results are not comparable to what I got using cmdstan.

20

I meant why is this

so different from this

21

That’s my question. Why doesnt rstan return reasonable results for this model.

The scenario is almost the same as the OPs, except I am having problems with rstan rather than cmdstan.