Beta-binomial model


#1

Please share your Stan program and accompanying data if possible.


I’ve been away from stan for a bit. I’m trying to write what should be a straighforward, simple program to sample a binomial where the parameter for probability of success R follows a beta distribution with known parameters, alpha and beta. I use Rstan. Apparently I cannot figure out how to declare the random outcome integer variable M correctly. Here is the .R file and the .stan file:

#Beta-Binimial models preparatory for Jolly-Seber-Lincoln-Petersen

#N0 = 1000; #number marked and released
#A =399; #Beta alpha param, # prior successes from genetic assignment
#B = 97; #Beta beta param, # prior failures from genetic assignment

stanDat <- list(N0 = 1000, A = 399, B = 97)

#sample posterior:
BetaBin1fit <- stan(file=“BetaBin1.stan”, data = stanDat,iter = 100, chains = 1)

print(BetaBin1fit)

//BetaBin1.stan

data{
int<lower = 0> N0;

real<lower = 0> A;
real<lower = 0> B;
}

parameters{
vector[] M;
real<lower = 0.70, upper = 0.95> R;
}

model {
R ~ beta(A,B);
M~binomial(N0, R);
}

where should I put the declaration for M and how should it be written. I have tried several approaches but none work.

Mick


#2

The reason you’re running into trouble is that Stan does not admit integer parameters (has to do with differentiability of the target with respect to all parameters in the model).

To get samples from R \sim \text{beta}(A, B); M \mid R \sim \text{binomial}(N_0, R) , this will work:

data{
  int<lower = 0> N0;
  real<lower = 0> A;
  real<lower = 0> B;
}
parameters{
  real<lower = 0.70, upper = 0.95> R;
}
model {
  R ~ beta(A, B);
}
generated quantities{
  int M = binomial_rng(N0, R);
}

image


#3

Much thanks for this solution!