top of page

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida

you are kindly asked to include the complete citation if you used this material in a publication

â€‹

â€‹

â€‹

Code 10.22 Normal autoregressive model AR(1) for accessing the evolution of the number of sunspots through the years.

===================================================================================

require(R2jags)
require(jagstools)

â€‹

# Data

â€‹

# Prepare data to JAGS
y <- round(sunspot[,2])
t <- seq(1700,2015,1)
N <- length(y)

sun_data <- list(Y = y,                                           # Response variable
N = N)                                         # Sample size

â€‹

# Fit
AR1_NORM<-"model{
# Priors
sd2 ~ dgamma(1e-3,1e-3)
tau <- 1/sd2
sd <- sqrt(sd2)

â€‹

for(i in 1:2){
phi[i] ~dnorm(0,1e-2)
}

â€‹

mu[1] <- Y[1]

â€‹

# Likelihood function
for (t in 2:N) {
Y[t] ~ dnorm(mu[t],tau)
mu[t] <- phi[1] + phi[2] * Y[t-1]
}

â€‹

# Prediction
for (t in 1:N){
Yx[t]~dnorm(mu[t],tau)
}
}"

â€‹

# Generate initial values for mcmc
inits <- function () {
list(phi = rnorm(2, 0, 0.1))
}

â€‹

# Identify parameters
# Include "Yx" in the list bellow if you intend to generate plots

params <- c("sd", "phi", "Yx")

â€‹

# Run mcmc
jagsfit <- jags(data = sun_data,
inits = inits,
parameters = params,
model = textConnection(AR1_NORM),
n.thin = 1,
n.chains = 3,
n.burnin = 3000,
n.iter = 5000)

â€‹

# Output
print(jagsfit,intervals = c(0.025, 0.975),justify = "left", digits=2)

===================================================================================

â€‹

â€‹

â€‹

â€‹

Code 10.23 Model shown in Fiugre 10.21

===============================
require(ggplot2)

require(jagstools)

# Format data for ggplot
sun <- data.frame(x=t,y=y)                                                                          # original data
yx <- jagsresults(x=jagsfit, params=c("Yx"))                                              # fitted values
gdata <- data.frame(x = t, mean = yx[,"mean"],lwr1=yx[,"25%"],
lwr2=yx[,"2.5%"],upr1=yx[,"75%"],upr2=yx[,"97.5%"])

# Plot
ggplot(sun,aes(x=t,y=y))+
geom_ribbon(data=gdata,aes(x=t,ymin=lwr1, ymax=upr1,y=NULL),
alpha=0.95, fill=c("#969696"),show.legend=FALSE)+
geom_ribbon(data=gdata,aes(x=t,ymin=lwr2, ymax=upr2,y=NULL),
alpha=0.75, fill=c("#d9d9d9"),show.legend=FALSE)+
geom_line(data=gdata,aes(x=t,y=mean),colour="black",
linetype="solid",size=0.5,show.legend=FALSE)+
geom_point(colour="orange2",size=1.5,alpha=0.75)+
theme_bw() + xlab("Year") + ylab("Sunspots") +
theme(legend.background = element_rect(fill="white"),
legend.key = element_rect(fill = "white",color = "white"),
plot.background = element_rect(fill = "white"),
legend.position="top",
axis.title.y = element_text(vjust = 0.1,margin=margin(0,10,0,0)),
axis.title.x = element_text(vjust = -0.25),
text = element_text(size = 25,family="serif"))+
geom_hline(aes(yintercept=0),linetype="dashed",colour="gray45",size=1.25)

===================================================================================

Anchor 1

Output on screen:

â€‹

Inference for Bugs model at "4", fit using jags,

3 chains, each with 5000 iterations (first 3000 discarded)

n.sims = 6000 iterations saved

mu.vect    sd.vect         2.5%       97.5%  Rhat     n.eff

phi[1]               13.42        3.12          7.29         19.54        1     3900

phi[2]                 0.83        0.03          0.76           0.89        1     6000

sd                     35.81        1.41         33.22         38.67       1     6000

deviance      3150.05        2.46     3147.19     3156.26        1     6000

For each parameter, n.eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)

pD = 3.0 and DIC = 3153.1

DIC is an estimate of expected predictive error (lower deviance is better).

bottom of page