library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores=5)

vs.data <- list(I=I, T=T, Tmax=Tmax,
                Obs=Obs, Time=Time, Dose=Dose,
                T0=T0, slope_mu_mu=slope_mu_mu,slope_mu_std=slope_mu_std,
                obs_rel=obs_rel, obs_mean_t1=obs_mean_t1,
                obs_std_t1=obs_std_t1)

vs.fit <- stan(file="varSlopes.stan", data=vs.data, iter=2000, chains=5)

print(vs.fit,pars=c("var_innov","treat_eff","T0_err_std","slope_mu",
                    "slope_std","lp__"))
plot(vs.fit,pars=c("var_innov","treat_eff","T0_err_std","slope_mu",
                    "slope_std","lp__"))
pairs(vs.fit,pars=c("var_innov","treat_eff","T0_err_std","slope_mu",
                    "slope_std","lp__"))

monitor(extract(vs.fit,pars="slope",permuted=FALSE))
innov_est <- extract(vs.fit,pars="innov")$innov
innov_mean <- apply(innov_est,c(2,3),mean)
pairs(innov_mean)
cor(innov_mean)
summary(as.vector(cor(innov_mean)))
hist(as.vector(cor(innov_mean)))