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

â€‹

â€‹

â€‹

# Data from  Code 10.22

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

â€‹

â€‹

â€‹

â€‹

â€‹

Code 10.24 Negative binomial model for assessing the evolution of the number of sunspots through the years

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

# Fit
AR1_NB<-"model{
for(i in 1:2){
phi[i] ~ dnorm(0,1e-2)
}

â€‹

theta ~ dgamma(0.001,0.001)
mu[1] <- Y[1]

â€‹

# Likelihood function
for (t in 2:N) {
Y[t] ~ dnegbin(p[t],theta)
p[t] <- theta/(theta+mu[t])
log(mu[t]) <- phi[1] + phi[2]*Y[t-1]
}

â€‹

for (t in 1:N){
Yx[t] ~ dnegbin(px[t],theta)
px[t] <- theta/(theta+mu[t])
}
}"

â€‹

# Identify parameters
# Include Yx in the list bellow only if interested in prediction

params <- c("phi","theta")

â€‹

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

â€‹

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

â€‹

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

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

Anchor 1

Output on screen:

â€‹

â€‹

bottom of page