#Load packages
library(rstan)
## Loading required package: StanHeaders
##
## rstan version 2.32.6 (Stan version 2.32.2)
## 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)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
library(tidyverse)
## āā Attaching core tidyverse packages āāāāāāāāāāāāāāāāāāāāāāāā tidyverse 2.0.0 āā
## ā dplyr 1.1.4 ā readr 2.1.5
## ā forcats 1.0.0 ā stringr 1.5.1
## ā ggplot2 3.5.1 ā tibble 3.2.1
## ā lubridate 1.9.3 ā tidyr 1.3.1
## ā purrr 1.0.2
## āā Conflicts āāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāā tidyverse_conflicts() āā
## ā tidyr::extract() masks rstan::extract()
## ā dplyr::filter() masks stats::filter()
## ā dplyr::lag() masks stats::lag()
## ā¹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(bayesplot)
## This is bayesplot version 1.11.1
## - 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
library(ggplot2)
library(patchwork)
library(modeest)
#a)Simulate the posterior model of Ī¼
# STEP 1: DEFINE the model
nn_model <- "
data {
real Y[5];
}
parameters {
real mu;
}
model {
Y ~ normal(mu,8);
mu ~ normal(-14, 2);
}
"
# STEP 2: SIMULATE the posterior
nn_sim <- stan(model_code = nn_model, data = list(Y = c(-10.1,5.5,0.1,-1.4,11.5)),
chains = 4, iter = 5000*2, seed = 84131)
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.028 seconds (Warm-up)
## Chain 1: 0.025 seconds (Sampling)
## Chain 1: 0.053 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 2e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.026 seconds (Warm-up)
## Chain 2: 0.026 seconds (Sampling)
## Chain 2: 0.052 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.01 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.027 seconds (Warm-up)
## Chain 3: 0.026 seconds (Sampling)
## Chain 3: 0.053 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.01 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.027 seconds (Warm-up)
## Chain 4: 0.03 seconds (Sampling)
## Chain 4: 0.057 seconds (Total)
## Chain 4:
as.array(nn_sim, pars = "mu") %>%
head()
## , , parameters = mu
##
## chains
## iterations chain:1 chain:2 chain:3 chain:4
## [1,] -9.716761 -11.46952 -10.267167 -10.489109
## [2,] -11.051514 -11.00563 -10.338113 -10.489109
## [3,] -12.333823 -11.08624 -10.507266 -9.583469
## [4,] -12.118245 -12.60533 -10.287138 -9.581286
## [5,] -11.645466 -12.60533 -8.974546 -9.365798
## [6,] -9.467670 -10.00483 -9.196780 -11.237714
#b)trace and density plots
# Trace plots of the 4 Markov chains
mcmc_trace(nn_sim, pars = "mu", size = 0.1) +
labs (x='Iteration',y='mu value',
title='Trace plot - mu')
# Histogram and density of the Markov chain values
h<-mcmc_hist(nn_sim, pars = "mu") +
labs (x='mu',y='count',
title='Histogram of the Markov chain values')
d<-mcmc_dens(nn_sim, pars = "mu") +
labs (x='mu',y='density',
title='Density: Markov chain values')
h+d
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#c)the most posterior plausible value of Ī¼
summary(nn_sim)
## $summary
## mean se_mean sd 2.5% 25% 50% 75%
## mu -10.41199 0.020240937 1.7515587 -13.81225 -11.594990 -10.399089 -9.219436
## lp__ -9.34021 0.007534389 0.7114367 -11.31873 -9.509052 -9.068771 -8.889390
## 97.5% n_eff Rhat
## mu -6.990454 7488.385 1.0001004
## lp__ -8.837261 8916.132 0.9999944
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## mu -10.453798 1.7622917 -13.89320 -11.656695 -10.427440 -9.240806
## lp__ -9.346772 0.7131484 -11.33699 -9.526227 -9.077906 -8.891511
## stats
## parameter 97.5%
## mu -7.095868
## lp__ -8.837275
##
## , , chains = chain:2
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## mu -10.391109 1.7495596 -13.88998 -11.555318 -10.365428 -9.190345
## lp__ -9.338976 0.7155019 -11.27400 -9.501367 -9.067792 -8.887639
## stats
## parameter 97.5%
## mu -7.008598
## lp__ -8.837162
##
## , , chains = chain:3
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## mu -10.391333 1.7677099 -13.79405 -11.582347 -10.389337 -9.199179
## lp__ -9.349447 0.7378461 -11.36575 -9.514434 -9.069344 -8.887941
## stats
## parameter 97.5%
## mu -6.901348
## lp__ -8.837318
##
## , , chains = chain:4
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## mu -10.411716 1.7261612 -13.75336 -11.577911 -10.418617 -9.247527
## lp__ -9.325646 0.6779308 -11.23872 -9.495355 -9.059019 -8.889619
## stats
## parameter 97.5%
## mu -7.009440
## lp__ -8.837272
#d)Specify the posterior model of Ī¼. How does your MCMC approximation compare?
# Histogram and density of the Markov chain values, mean and mode (posterior model)
#reorder data
chains<- data.frame(as.array(nn_sim, pars = "mu"))
vchains<-c(chains[,1],chains[,2],chains[,3],chains[,4])
dchains<-data.frame(vchains)
# Histogram of the Markov chain values
hist(dchains$vchains,
prob = TRUE,
xlim = c(-16, -3),
col = "royalblue3",
border = "blue",
main = "Histogram and Density of MCMC values - Mode and Mean of Posterior",
xlab = "mu",
las = 1)
# Density of the Markov chain values
lines(density(dchains$vchains), # density plot
lwd = 2, # thickness of line
col = "black")
# Mean (posteriori model)
segments(x0 = -10.4, y0 = 0, x1 = -10.4,
y1 = dnorm(-10.4,-10.4,1.7457),
col = "red", lwd = 2)
# Legend
legend(x ="topright", # UbicaciĆ³n de la leyenda
c("Density of MCMC values","Mean of Posterior"),
col = c("black", "red"),
lwd = c(1, 1, 1),
bty = "n")
#Density of the Markov chain values and density Posterior, mean and mode (posterior model)
# Density posterior moodel Normal(8.647,0.5715^2)
x = -16:-3
curve(dnorm(x, mean=-10.4, sd=1.7457),-16,-3, # density plot
xlab = "mu",
ylab = "density",
main = "Density of MCMC values and Posterior - Mode and Mean of Posterior",
lwd = 5, # thickness of line
col = "orange")
# Density of the Markov chain values
lines(density(dchains$vchains), # density plot
lwd = 2, # thickness of line
col = "black")
# Mean (posteriori model)
segments(x0 = -10.4, y0 = 0, x1 = -10.4,
y1 = dnorm(-10.4,-10.4,1.7457),
col = "red", lwd = 2)
# Legend
legend(x ="topright", # UbicaciĆ³n de la leyenda
c("Density of Posterior","Density of MCMC values", "Mean of Posterior"),
col = c("orange","black", "red"),
lwd = c(1,1, 1, 1),
bty = "n")
#Density plot of the Markov chain values and posterior density with mcmc_dens and stat_function
mcmc_dens(nn_sim, pars = "mu") +
labs (x='mu',y='density',
title='Density of MCMC values and Posterior') +
stat_function (
fun=dnorm,
args=list(mean=-10.4,sd=1.7457),
lwd = 1,
color="black"
)
