# edps 590BAY # Fall 2019 # C.J.Anderson # # Estimate the mean and variance of normaly distributed variable # # Illustration of jags input and output for # rjags # library(coda) library(rjags) ano <- read.table("anorexia_data.txt",header=T) ano$change <- ano$weight2 - ano$weight1 ############################################################################################ set.seed(234590) dataList <- list(y=ano$change, Ntotal=length(ano$change), meanY = mean(ano$change), sdY = sd(ano$change) ) Model1 = "model { for (i in 1:Ntotal){ y[i] ~ dnorm( mu, tau ) } mu ~ dnorm( meanY , 1/(100*sdY^2) ) # mean (y) and precision N(2.763889,1.568927e-06) tau <- pow(sigma, -2) # tau is presicion = 1/var(y) sigma ~ dunif( sdY/1000, sdY*1000 ) # variance (y) and uniform(.00798,7983.598) } " writeLines(Model1, con="Model1.txt") thetaInit = mean(ano$change) sigmaInit = sd(ano$change) initsList = list(mu=thetaInit,sigma=sigmaInit) jagsModel1 <- jags.model(file="Model1.txt", # compiles and intializes model data=dataList, inits=initsList, n.chains=4, n.adapt=500) # adaptive phase for 500 iterations update (jagsModel1, n.iter=500) # gets samples from posterior with a # burn in of 500 iterations # contains samples from all chains with 500 iterations Samples <- coda.samples(jagsModel1, variable.names=c("mu","tau","sigma"), n.iter=4000) # output summary information summary(Samples) plot(Samples)