#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
gp_model <- "
data {
int<lower = 0> Y[3];
}
parameters {
real<lower = 0> lambda;
}
model {
Y ~ poisson(lambda);
lambda ~ gamma(20, 5);
}
"
# STEP 2: SIMULATE the posterior
gp_sim <- stan(model_code = gp_model, data = list(Y = c(0,1,0)),
chains = 4, iter = 5000*2, seed = 84131)
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.16 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.044 seconds (Warm-up)
## Chain 1: 0.053 seconds (Sampling)
## Chain 1: 0.097 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.036 seconds (Warm-up)
## Chain 2: 0.038 seconds (Sampling)
## Chain 2: 0.074 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 2e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.02 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.036 seconds (Warm-up)
## Chain 3: 0.036 seconds (Sampling)
## Chain 3: 0.072 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 2e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.02 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.035 seconds (Warm-up)
## Chain 4: 0.043 seconds (Sampling)
## Chain 4: 0.078 seconds (Total)
## Chain 4:
as.array(gp_sim, pars = "lambda") %>%
head()
## , , parameters = lambda
##
## chains
## iterations chain:1 chain:2 chain:3 chain:4
## [1,] 3.727827 1.929448 2.122099 3.702171
## [2,] 4.033571 2.577971 2.227208 3.229174
## [3,] 2.473422 2.665157 2.373778 2.805807
## [4,] 2.308675 2.665157 2.949798 2.645933
## [5,] 2.345855 3.566446 3.509082 2.611989
## [6,] 2.388877 2.165358 3.509082 2.692876
#b)trace and density plots
# Trace plots of the 4 Markov chains
mcmc_trace(gp_sim, pars = "lambda", size = 0.1) +
labs (x='Iteration',y='lambda value',
title='Trace plot - lambda')
# Histogram and density of the Markov chain values
h<-mcmc_hist(gp_sim, pars = "lambda") +
labs (x='lambda',y='count',
title='Histogram of the Markov chain values')
d<-mcmc_dens(gp_sim, pars = "lambda") +
labs (x='lambda',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(gp_sim)
## $summary
## mean se_mean sd 2.5% 25% 50% 75%
## lambda 2.617061 0.006586424 0.5739810 1.630153 2.209710 2.574702 2.9742381
## lp__ -1.240702 0.008494703 0.7118391 -3.237515 -1.413263 -0.962999 -0.7838307
## 97.5% n_eff Rhat
## lambda 3.865794 7594.437 1.000367
## lp__ -0.733799 7022.104 1.000351
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## lambda 2.619808 0.5741498 1.639183 2.207181 2.574377 2.9759515
## lp__ -1.238223 0.6883728 -3.161966 -1.432394 -0.970706 -0.7841654
## stats
## parameter 97.5%
## lambda 3.8325850
## lp__ -0.7338434
##
## , , chains = chain:2
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## lambda 2.617336 0.5660456 1.640373 2.221380 2.5802976 2.9699064
## lp__ -1.228558 0.7139942 -3.239109 -1.382237 -0.9529592 -0.7828873
## stats
## parameter 97.5%
## lambda 3.8406321
## lp__ -0.7336622
##
## , , chains = chain:3
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## lambda 2.600924 0.5792129 1.597723 2.189538 2.5595662 2.9640474
## lp__ -1.256558 0.7332267 -3.398458 -1.425861 -0.9762372 -0.7869109
## stats
## parameter 97.5%
## lambda 3.8678122
## lp__ -0.7339998
##
## , , chains = chain:4
##
## stats
## parameter mean sd 2.5% 25% 50% 75%
## lambda 2.630173 0.5762214 1.634297 2.227473 2.5883325 2.9869224
## lp__ -1.239471 0.7109786 -3.227732 -1.413204 -0.9551701 -0.7802573
## stats
## parameter 97.5%
## lambda 3.9070489
## lp__ -0.7337052
#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(gp_sim, pars = "lambda"))
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(1, 5),
col = "royalblue3",
border = "blue",
main = "Histogram and Density of MCMC values - Mode and Mean of Posterior",
xlab = "lambda",
las = 1)
# Density of the Markov chain values
lines(density(dchains$vchains), # density plot
lwd = 2, # thickness of line
col = "black")
# Mode (posteriori model)
mode <- gammaMode(21, 8)
segments(x0 = mode, y0 = 0, x1 = mode,
y1 = dgamma(mode, 21, 8),
col = "green", lwd = 2)
# Mean (posteriori model)
mean <- 21/8
segments(x0 = mean, y0 = 0, x1 = mean,
y1 = dgamma(mean, 21, 8),
col = "red", lwd = 2)
# Legend
legend(x ="topright", # UbicaciĆ³n de la leyenda
c("Density of MCMC values","Mode of Posterior" , "Mean of Posterior"),
col = c("black", "green", "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 Gamma(21,8)
x = 1:5
curve(dgamma(x, 21, 8),1,5, # density plot
xlab = "lambda",
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")
# Mode (posteriori model)
mode <- gammaMode(21, 8)
segments(x0 = mode, y0 = 0, x1 = mode,
y1 = dgamma(mode, 21, 8),
col = "green", lwd = 2)
# Mean (posteriori model)
mean <- 21/8
segments(x0 = mean, y0 = 0, x1 = mean,
y1 = dgamma(mean, 21, 8),
col = "red", lwd = 2)
# Legend
legend(x ="topright", # UbicaciĆ³n de la leyenda
c("Density of Posterior","Density of MCMC values","Mode of Posterior" , "Mean of Posterior"),
col = c("orange","black", "green", "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(gp_sim, pars = "lambda") +
labs (x='lambda',y='density',
title='Density of MCMC values and Posterior') +
stat_function (
fun=dgamma,
args=list(shape=21,rate=8),
lwd = 1.5,
color="black"
)
