library(rstan)
## Loading required package: StanHeaders
## Loading required package: ggplot2
## rstan (Version 2.21.7, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
library(extraDistr)
library(loo)
## This is loo version 2.5.1
## - Online documentation and vignettes at mc-stan.org/loo
## - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
##
## Attaching package: 'loo'
## The following object is masked from 'package:rstan':
##
## loo
library(coda)
##
## Attaching package: 'coda'
## The following object is masked from 'package:rstan':
##
## traceplot
library(bayesplot)
## This is bayesplot version 1.9.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
Read in data
library(readxl)
BayesGain<-readxl::read_excel(path = "/Users/katarinaploch/Downloads/Bayes Stats Assignments /BayesGain Final.xlsx")
n <- NROW(BayesGain)
data.list <- with(BayesGain, list(ycat=AS, n = nrow(BayesGain)))
modelString ="
data {
int<lower=0> n;
int<lower=1, upper=9> ycat[n];
}
parameters {
simplex[9] theta; //simplex indicates a vector of real numbers that sum to 1
}
model {
ycat ~ categorical(theta);
theta ~ dirichlet([1,1,1,1,1,1,1,1,1]');
// The hyperparameter values are the concentration hyperparameters. Here all outcomes are equally likely to have been observed.
}
"
nChains = 4
nIter= 30000
thinSteps = 10
burnInSteps = floor(nIter/2)
MultiFit= stan(data=data.list,model_code=modelString,chains=nChains,
iter=nIter,warmup=burnInSteps,thin=thinSteps)
##
## SAMPLING FOR MODEL '0ed4ee8acce8b691d156bf50858d5fa1' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 4.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 30000 [ 0%] (Warmup)
## Chain 1: Iteration: 3000 / 30000 [ 10%] (Warmup)
## Chain 1: Iteration: 6000 / 30000 [ 20%] (Warmup)
## Chain 1: Iteration: 9000 / 30000 [ 30%] (Warmup)
## Chain 1: Iteration: 12000 / 30000 [ 40%] (Warmup)
## Chain 1: Iteration: 15000 / 30000 [ 50%] (Warmup)
## Chain 1: Iteration: 15001 / 30000 [ 50%] (Sampling)
## Chain 1: Iteration: 18000 / 30000 [ 60%] (Sampling)
## Chain 1: Iteration: 21000 / 30000 [ 70%] (Sampling)
## Chain 1: Iteration: 24000 / 30000 [ 80%] (Sampling)
## Chain 1: Iteration: 27000 / 30000 [ 90%] (Sampling)
## Chain 1: Iteration: 30000 / 30000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.763123 seconds (Warm-up)
## Chain 1: 0.961254 seconds (Sampling)
## Chain 1: 1.72438 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '0ed4ee8acce8b691d156bf50858d5fa1' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.8e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 30000 [ 0%] (Warmup)
## Chain 2: Iteration: 3000 / 30000 [ 10%] (Warmup)
## Chain 2: Iteration: 6000 / 30000 [ 20%] (Warmup)
## Chain 2: Iteration: 9000 / 30000 [ 30%] (Warmup)
## Chain 2: Iteration: 12000 / 30000 [ 40%] (Warmup)
## Chain 2: Iteration: 15000 / 30000 [ 50%] (Warmup)
## Chain 2: Iteration: 15001 / 30000 [ 50%] (Sampling)
## Chain 2: Iteration: 18000 / 30000 [ 60%] (Sampling)
## Chain 2: Iteration: 21000 / 30000 [ 70%] (Sampling)
## Chain 2: Iteration: 24000 / 30000 [ 80%] (Sampling)
## Chain 2: Iteration: 27000 / 30000 [ 90%] (Sampling)
## Chain 2: Iteration: 30000 / 30000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.759316 seconds (Warm-up)
## Chain 2: 0.837041 seconds (Sampling)
## Chain 2: 1.59636 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '0ed4ee8acce8b691d156bf50858d5fa1' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.5e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 30000 [ 0%] (Warmup)
## Chain 3: Iteration: 3000 / 30000 [ 10%] (Warmup)
## Chain 3: Iteration: 6000 / 30000 [ 20%] (Warmup)
## Chain 3: Iteration: 9000 / 30000 [ 30%] (Warmup)
## Chain 3: Iteration: 12000 / 30000 [ 40%] (Warmup)
## Chain 3: Iteration: 15000 / 30000 [ 50%] (Warmup)
## Chain 3: Iteration: 15001 / 30000 [ 50%] (Sampling)
## Chain 3: Iteration: 18000 / 30000 [ 60%] (Sampling)
## Chain 3: Iteration: 21000 / 30000 [ 70%] (Sampling)
## Chain 3: Iteration: 24000 / 30000 [ 80%] (Sampling)
## Chain 3: Iteration: 27000 / 30000 [ 90%] (Sampling)
## Chain 3: Iteration: 30000 / 30000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.753476 seconds (Warm-up)
## Chain 3: 0.994665 seconds (Sampling)
## Chain 3: 1.74814 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '0ed4ee8acce8b691d156bf50858d5fa1' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.5e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 30000 [ 0%] (Warmup)
## Chain 4: Iteration: 3000 / 30000 [ 10%] (Warmup)
## Chain 4: Iteration: 6000 / 30000 [ 20%] (Warmup)
## Chain 4: Iteration: 9000 / 30000 [ 30%] (Warmup)
## Chain 4: Iteration: 12000 / 30000 [ 40%] (Warmup)
## Chain 4: Iteration: 15000 / 30000 [ 50%] (Warmup)
## Chain 4: Iteration: 15001 / 30000 [ 50%] (Sampling)
## Chain 4: Iteration: 18000 / 30000 [ 60%] (Sampling)
## Chain 4: Iteration: 21000 / 30000 [ 70%] (Sampling)
## Chain 4: Iteration: 24000 / 30000 [ 80%] (Sampling)
## Chain 4: Iteration: 27000 / 30000 [ 90%] (Sampling)
## Chain 4: Iteration: 30000 / 30000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.797262 seconds (Warm-up)
## Chain 4: 1.04795 seconds (Sampling)
## Chain 4: 1.84521 seconds (Total)
## Chain 4:
stan_plot(MultiFit, pars=c("theta"))
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)
stan_trace(MultiFit,inc_warmup=T, pars=c("theta"))
stan_dens(MultiFit, pars=c("theta"))
stan_ac(MultiFit, pars=c("theta"))
print(MultiFit,pars=c("theta"))
## Inference for Stan model: 0ed4ee8acce8b691d156bf50858d5fa1.
## 4 chains, each with iter=30000; warmup=15000; thin=10;
## post-warmup draws per chain=1500, total post-warmup draws=6000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## theta[1] 0.17 0 0.02 0.13 0.15 0.17 0.18 0.21 5821 1
## theta[2] 0.15 0 0.02 0.12 0.14 0.15 0.17 0.19 5754 1
## theta[3] 0.16 0 0.02 0.13 0.15 0.16 0.17 0.20 5742 1
## theta[4] 0.13 0 0.02 0.10 0.12 0.13 0.14 0.17 5763 1
## theta[5] 0.14 0 0.02 0.10 0.12 0.13 0.15 0.17 6150 1
## theta[6] 0.09 0 0.01 0.06 0.08 0.09 0.09 0.11 5887 1
## theta[7] 0.09 0 0.01 0.06 0.08 0.09 0.10 0.12 6149 1
## theta[8] 0.06 0 0.01 0.04 0.05 0.06 0.07 0.08 5951 1
## theta[9] 0.01 0 0.01 0.01 0.01 0.01 0.02 0.03 5774 1
##
## Samples were drawn using NUTS(diag_e) at Tue Oct 11 20:22:23 2022.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
library(coda)
newMultiFit<- As.mcmc.list(MultiFit)
HPDinterval(newMultiFit, prob=0.95)
## [[1]]
## lower upper
## theta[1] 1.327425e-01 0.20539802
## theta[2] 1.207722e-01 0.18975653
## theta[3] 1.291093e-01 0.20002635
## theta[4] 9.970627e-02 0.16295110
## theta[5] 1.058405e-01 0.17114519
## theta[6] 6.218472e-02 0.11400088
## theta[7] 6.250095e-02 0.11630344
## theta[8] 3.515281e-02 0.08075105
## theta[9] 5.055868e-03 0.02682651
## lp__ -8.511787e+02 -843.56317366
## attr(,"Probability")
## [1] 0.95
##
## [[2]]
## lower upper
## theta[1] 1.337481e-01 0.20499252
## theta[2] 1.192928e-01 0.18921009
## theta[3] 1.301460e-01 0.19995011
## theta[4] 9.960522e-02 0.16581580
## theta[5] 1.009104e-01 0.16799242
## theta[6] 5.781139e-02 0.11018548
## theta[7] 6.190776e-02 0.11453766
## theta[8] 3.674444e-02 0.08208742
## theta[9] 5.906006e-03 0.02813354
## lp__ -8.508687e+02 -843.58831723
## attr(,"Probability")
## [1] 0.95
##
## [[3]]
## lower upper
## theta[1] 1.312393e-01 0.20359029
## theta[2] 1.239632e-01 0.19113168
## theta[3] 1.278433e-01 0.19631911
## theta[4] 1.005959e-01 0.16366459
## theta[5] 1.033640e-01 0.16745725
## theta[6] 6.050637e-02 0.11287980
## theta[7] 6.308925e-02 0.11713055
## theta[8] 3.845971e-02 0.08304086
## theta[9] 4.298999e-03 0.02638565
## lp__ -8.505559e+02 -843.38810450
## attr(,"Probability")
## [1] 0.95
##
## [[4]]
## lower upper
## theta[1] 1.339860e-01 0.20479230
## theta[2] 1.212594e-01 0.18935979
## theta[3] 1.259914e-01 0.19894415
## theta[4] 9.793910e-02 0.16232645
## theta[5] 1.015954e-01 0.16581031
## theta[6] 6.064911e-02 0.11107420
## theta[7] 6.210433e-02 0.11632125
## theta[8] 3.798438e-02 0.08190471
## theta[9] 3.993838e-03 0.02650840
## lp__ -8.505774e+02 -843.64710506
## attr(,"Probability")
## [1] 0.95