/
code_10.5.R
77 lines (64 loc) · 2.24 KB
/
code_10.5.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
# 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.5 - Multivariate normal model in R using JAGS for accessing the relationship
# between period, luminosity, and color in early-type contact binaries
library(R2jags)
# Data
PLC <- read.csv("https://raw.githubusercontent.com/astrobayes/BMAD/master/data/Section_10p3/PLC.csv", header = T)
# Prepare data for JAGS
nobs = nrow(PLC) # number of data points
x1 <- PLC$logP # log period
x2 <- PLC$V_I # V-I color
y <- PLC$M_V # V magnitude
type <- as.numeric(PLC$type) # type NC/GC
X <- model.matrix(~ 1 + x1+x2) # covariate matrix
K <- ncol(X) # number of covariates per type
jags_data <- list(Y = y,
X = X,
K = K,
type = type,
N = nobs)
# Fit
NORM <-"model{
# Shared hyperprior
tau0 ~ dgamma(0.001,0.001)
mu0 ~ dnorm(0,1e-3)
# Diffuse normal priors for predictors
for(j in 1:2){
for (i in 1:K) {
beta[i,j] ~ dnorm(mu0, tau0)
}
}
# Uniform prior for standard deviation
for(i in 1:2) {
tau[i] <- pow(sigma[i], -2) #precision
sigma[i] ~ dgamma(1e-3, 1e-3) #standard deviation
}
# Likelihood function
for (i in 1:N){
Y[i]~dnorm(mu[i],tau[type[i]])
mu[i] <- eta[i]
eta[i] <- beta[1, type[i]] * X[i, 1] + beta[2, type[i]] * X[i, 2] +
beta[3, type[i]] * X[i, 3]
}
}"
# Determine initial values
inits <- function () {
list(beta = matrix(rnorm(6,0, 0.01),ncol=2))
}
# Identify parameters
params <- c("beta", "sigma")
# Fit
jagsfit <- jags(data = jags_data,
inits = inits,
parameters = params,
model = textConnection(NORM),
n.chains = 3,
n.iter = 5000,
n.thin = 1,
n.burnin = 2500)
## Output
print(jagsfit,justify = "left", digits=2)