jhoupt / dbda2estan Goto Github PK

First, I'd put the Stan programs in their own files. It makes it easier to line up error messages and lets you reuse the model in other scripts.

OneOddGroupModelComp2E_Stan.R

  real<lower=0> kappa[nCond];
  for ( j in 1:nCond ) {
    kappa[j] <- kappaMinusTwo[j] + 2;
  }

can be

  vector<lower=0>[nCond] kappa = kappaMinusTwo + 2;

if you also declare kappaMinusTwo as a vector.

Same thing for a and b local variables in the model block, though you're much better off for efficiency defining true constants as transformed data (then they don't waste time during derivative calculations that go unused).

Pretty much everything you're doing can be vectorized, so this

for ( s in 1:nSubj ) {
    nCorrOfSubj[s] ~ binomial( nTrlOfSubj[s],  theta[s] );

is just

nCorrOfSubj ~ binomial(nTrlOfSubj, theta);

beta_log is deprecated in favor of beta_lpdf.

<- is deprecated in favor of =. The deprecations matter because they'll go away in a future release.

The bigger problem is that you almost never want to work on the exponentiated scale, so this is a big red flag:

prob1 <- modelProb1*exp(beta_log( theta[s], aBeta, bBeta ));

You have to stay on the log scale to prevent underflow and go back to that log_sum_exp form, which is right, though log1m(u) is more stable than log(1 - u) and there's a log_mix function in Stan that turns this all into a one-liner.

Vectorize

kappaMinusTwo ~ gamma( 2.618 , 0.0809 );

but that is a suspicious prior both in its decimal points of precision and where its mean and tail go.

omega ~ beta(aP, bP);

You also don't need all the spaces. There's a style guide in the Stan manual that basically follows C++ programming standards.

modelProb1 ~ beta(1,1);

is redundant because beta(1, 1) is uniform. It won't be too bad as Stan will drop it from everything other than error checking (1 > 0, 1 > 0, modelProb in (0, 1)).