Thank you @avehtari and @bbbales2 for helping to integrate DS(G,m) with STAN…your advice helped us implement the DS(G,m) model described in the “ModelDescription.pdf” attachment of the original post:

```
library(rstan)
### MC Options
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
########
# Rat Data from Gelman's BDA (Chptr 5)
J <- 71
y <- as.integer(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
1, 1, 1, 1, 1, 1, 1, 3, 2, 2, 2, 2, 2, 2, 2, 2,
2, 1, 5, 2, 5, 2, 7, 7, 3, 3, 2, 9, 10, 4, 4,
4, 4, 4, 4, 4, 10, 4, 4, 4, 5, 11, 12, 5, 5, 6,
5, 6, 6, 6, 6, 16, 15, 15, 9,4))
n <- as.integer(c(20, 20, 20, 20, 20, 20, 20, 19, 19, 19, 19, 18,
18, 17, 20, 20, 20, 20, 19, 19, 18, 18, 27, 25, 24,
23, 20, 20, 20, 20, 20, 20, 10, 49, 19, 46, 17, 49,
47, 20, 20, 13, 48, 50, 20, 20, 20, 20, 20, 20, 20,
48, 19, 19, 19, 22, 46, 49, 20, 20, 23, 19, 22, 20,
20, 20, 52, 46, 47, 24,14))
#####################
####### STAN
# The model specification
rat.model2.string <- "
functions{
real DSbeta_lpdf(real theta , real alpha, real beta, real c1,
real c2, real c3, real c4){
real lprobBeta;
real Goftheta;
real lprobDS;
Goftheta = exp(beta_lcdf(theta | alpha, beta));
lprobBeta = beta_lpdf(theta|alpha, beta);
lprobDS = log(1 + c1*sqrt(12)*(Goftheta-0.5) +
c2*sqrt(5)*(6*Goftheta^2 - 6*Goftheta +1) +
c3*sqrt(7)*(20*Goftheta^3 - 30*Goftheta^2 +12*Goftheta -1)+
c4*3*(70*Goftheta^4 -140*Goftheta^3+90*Goftheta^2-20*Goftheta+1));
return(lprobBeta + lprobDS);
}
}
data {
int<lower=0> J;
int<lower=0> n[J];
int<lower=0> y[J];
}
parameters {
real<lower=0> alpha;
real<lower=0> beta;
real<lower=0,upper=1> theta[J];
real c1;
real c2;
real c3;
real c4;
}
model {
theta[J] ~ DSbeta(alpha, beta, c1, c2, c3, c4);
y ~ binomial(n, theta);
alpha ~ normal(2.4, 5);
beta ~ normal(14,5);
c1 ~ normal(0,5);
c2 ~ normal(0,5);
c3 ~ normal(0,5);
c4 ~ normal(0,5);
}"
# Running the model
stan.mcmc.rat <- stan(model_code= rat.model2.string,
data= list(J = J, y = y, n = n),
pars=c("c1", "c2", "c3", "c4"),
chains=3, iter=10000, warmup=2000,
control = list(adapt_delta = 0.99))
```

Although the code is running, it is not producing desirable output for the parameters…@deepstan and I suspect that there is something we have missed in the construction of the model. We respectfully request your assistance in verifying whether we have properly coded the model in STAN.