Good morning all
I have a really simple question! How does one specify the following priors in jags?
tau1 ~ dgamma(a, b)
and now for some positive value 'k'
tau2 ~ dgamma(a, b) truncated on the interval tau1/k < tau2 < tau1
I've tried the following:
method 1
--------------
tau1 ~ dgamma(a, b)
tau2 ~ dgamma(a, b) T(tau1/K2, tau1)
This does not work since T() has to have deterministic values for the lower and upper bound entries of T()
Given the above comment. I them transform from tau2 to tau2t where tau2t = tau2/tau1. In this case tau2t ~G(a, b*tau1) truncated on the interval 1/k < tau2t < 1.
tau1~dgamma(a,b)
hi<- b*tau1
tau2t~dgamma(a,hi)T(1/K2, 1) #I've also tried tau2t~dgamma(a, b*tau1)T(0.02, 1)
tau2<-tau1*tau2t
This does not work either! I get the following error:
"Error in update.jags(object, n.iter, ...) : Error in node y[1]
Failure to calculate log density"
Does anyone have any ideas?
regards
The complete code is included below.
Its a really simple model:
Basically a linear model where some of the data points have inflated error variance!
#The data
hills<- DAAG::hills
X<-cbind(1, hills$dist, hills$climb)
y<-matrix(hills$time, ncol=1)
#the jags code
writeLines("
model{
#prior specifications
b0~dnorm(0, inv.v)
b1~dnorm(0, inv.v)
b2~dnorm(0, inv.v)
tau1~dgamma(a,b)
tau2t~dgamma(a,b*tau1)T(1/K2, 1)
tau2<-tau1*tau2t
#the likelihood of the data
for (i in 1:n1){
mu[i]<- b0 + b1*x1[i] + b2*x2[i]
y[i]~ dnorm(mu[i], tau1)
}
for (i in (n1+1):n){
mu[i]<- b0 + b1*x1[i] + b2*x2[i]
y[i]~ dnorm(mu[i], tau2)
}
}
", con = "model2.txt")
#---------------------------------------------------------------------
#inputs for the Bayesian method
niter <- 50000
nburn <- 20000
nthin <- 1
nchains <- 1
nkeep<-(niter-nburn)*nchains
# Input data
index<- c(1, 18, 33)
Data <- list(x1=c(X[-index,2],X[index,2]), x2=c(X[-index,3],X[index,3]),
y=c(y[-index],y[index]), n1=32, n=35 , a=0.001, b=0.001, inv.v=1/1000, K2=50)
# Initial values
inits <- function()
{
list("b0"=rnorm(1),"b1"=rnorm(1),"b2"=rnorm(1), "tau1"=10, "tau2t"=0.5)
}
inits1=inits(); #inits2=inits(); inits3=inits()
Inits=list(inits1)#,inits2,inits3)
# Paramaters to be monitored
parameters <- c("b0", "b1", "b2", "tau1","tau2")
modfile = "model2.txt"
rm(Runs2)
Runs2 <- jags(data=Data, inits=Inits, parameters.to.save=parameters,
model.file=modfile,
n.iter = niter, n.burnin = nburn, n.thin = nthin,
n.chains = nchains, DIC = FALSE, verbose=FALSE, parallel=TRUE)
Runs2