########## R script: rstanTest ##########

# For testing whether the rstan package is
# successfully installed.

# Last changed: 02 APR 2019

# Suppress warning messages:

options(warn=-1)

# Load required package:

library(rstan) 

# Set flag for code compilation (needed if 
# running script first time in current session) :

compileCode <- TRUE

# Generate some Bernoulli data:

pTrue <- 0.37 ; set.seed(1) 
x <- rbinom(100,1,pTrue) 

# Specify model in Stan:

BernoulliModel <- "data
                   {
                      int x[100];
                   } 
                   parameters
                   {
                      real<lower=0,upper=1> p;
                   }
                   model
                   {
                      x ~ bernoulli(p); 
                      p ~ uniform(0,1);
                   }"

# Compile code for model if required:

if (compileCode)
   stanCompilObj <- stan(model_code = BernoulliModel,
                      data = list(x = x),
                      iter = 1,chains = 1)

stanObj <- stan(model_code = BernoulliModel,
                data = list(x = x),warmup = 1000,
                iter = 2000,chains = 1,
                refresh = 100,
                fit = stanCompilObj)
pMCMC <- extract(stanObj,"p",permuted = FALSE)

plot(pMCMC,type="n",bty="l",xlab="index",ylab="MCMC draws from posterior of p",
     cex.lab=1.5) 
lines(pMCMC,col="darkgreen")

########## End of rstanTest ##########