add R/bugs/jags model

.gitignore

@@ -3,3 +3,5 @@
 *.pyc
 *.pickle
 *.png
+*.pdf
+radon.coefs.from.bugs.csv

radon.wishart.R

@@ -0,0 +1,42 @@
+library (MCMCpack)
+library(R2jags)
+
+srrs <- read.csv("srrs.csv",header=FALSE,as.is=TRUE)
+srrs <- srrs[order(srrs[,1]),]
+n <- nrow(srrs)
+county <- as.integer(as.factor(srrs[,1]))
+J <- length(unique(county))
+X <- cbind(rep(1,nrow(srrs)),srrs[,3])
+radon <- srrs[,2]
+y <- log(ifelse (radon==0, .1, radon))
+K <- ncol(X)
+W <- diag(K)
+
+
+model.names <- list ("n", "J","K","y","county","X","W")
+parameters.to.save <- c("B","mu","sigma.y","sigma.B","rho.B","e.y")
+model.inits <- function () {
+    list(B.raw=array(rnorm(J*K), c(J,K)),
+         mu.raw=rnorm(K),
+         sigma.y=runif(1),
+         Tau.B.raw=rwish(K+1,diag(K)),
+         xi=runif(K))
+}
+
+wd <- getwd()
+jags.time <- system.time(radon.model <- jags(data=model.names, inits=model.inits, parameters.to.save=parameters.to.save,
+                                             model.file=paste(wd,"radon.wishart.bug",sep="/"),DIC=FALSE,
+                                             n.chains=3, n.iter=10000, n.burnin=3000))
+print(jags.time)
+
+pdf("radon.wishart.bugs.pdf")
+plot(radon.model)
+dev.off()
+
+radon.summary <- radon.model[["BUGSoutput"]][["summary"]]
+radon.means <- radon.summary[,"mean",drop=F]
+radon.coefs <- cbind(radon.means[1:J],radon.means[1:J +J])
+e.y <- radon.model[["BUGSoutput"]][["sims.matrix"]][,grep('e.y', colnames(radon.model[["BUGSoutput"]][["sims.matrix"]]))]
+radon.rsq <- 1 - mean(apply(e.y,1,var)) / var(y)
+
+write.csv(radon.coefs,"radon.coefs.from.bugs.csv")

radon.wishart.bug

@@ -0,0 +1,30 @@
+model {
+  for (i in 1:n){
+    y[i] ~ dnorm(y.hat[i], tau.y)
+    y.hat[i] <- inprod(B[county[i],],X[i,])
+    e.y[i] <- y[i] - y.hat[i]
+  }
+  tau.y <- pow(sigma.y, -2)
+  sigma.y ~ dunif(0, 100)
+
+  for (j in 1:J) {
+      for(k in 1:K) {
+        B[j,k] <- xi[k] * B.raw[j,k]
+      }
+      B.raw[j,1:K] ~ dmnorm(mu.raw[], Tau.B.raw[,])
+  }
+  for(k in 1:K) {
+    mu[k] <- xi[k] * mu.raw[k]
+    mu.raw[k] ~ dnorm(0,.0001)
+    xi[k] ~ dunif(0,100)
+  }
+  Tau.B.raw[1:K,1:K] ~ dwish(W[,],df)
+  df <- K + 1
+  Sigma.B.raw[1:K,1:K] <- inverse(Tau.B.raw[,])
+  for(k in 1:K) {
+    for(k.prime in 1:K) {
+      rho.B[k,k.prime] <- Sigma.B.raw[k,k.prime] / sqrt(Sigma.B.raw[k,k] * Sigma.B.raw[k.prime,k.prime])
+    }
+    sigma.B[k] <- abs(xi[k]) * sqrt(Sigma.B.raw[k,k])
+  }
+}