example-runModel7.R

Revision 1 - 6/9/09 at 9:03 am by eric.ward

Next rev
Back to revision history for example-runModel7.R
This file is part of the project Individual Diet Modeling
library(runjags)
library(R2WinBUGS)
library(R2jags)
library(BRugs)
library(gtools)
library(gdata)

num.iso = 2   # number of isotopes
num.prey = 3  # number of prey species
num.pop = 3   # number of consumer populations (regions)
# prey * isotopes
u = array(0,dim=c(num.prey,num.iso,num.pop))
sigma2 = array(0,dim=c(num.prey,num.iso,num.pop))


   X = as.matrix(read.table("test_fake_data.txt",header=F))
   Region = sort(rep(seq(1,3),30))
   N = dim(X)[1]
   MU = read.table("test_fake_means.txt",header=F)
   SIG = read.table("test_fake_sds.txt",header=F)
   # these next 5 lines created so we won't have to change any of the
   # args passed into the code
   u = array(0,dim=c(num.prey,num.iso,num.pop))
   sigma2 = array(0,dim=c(num.prey,num.iso,num.pop))
   for(i in 1:3) {
      for(j in 1:2) {
         u[,j,i] = MU[,j]
         sigma2[,j,i] = SIG[,j]^2
      }
   }
   
# These are the parameters that need to be set by the user
  mcmc.chainLength <- as.integer(10000)  # post-burn
  mcmc.burn <- as.integer(5000)
  mcmc.thin = 1;
  mcmc.adapt = 5000
  mcmc.chains = 3       # needs to be at least 2 for DIC
  dic.samples = 5000
  alpha=rep(1,num.prey)
  alphaMix = c(1,1)
  R = diag(num.prey)
  dic = 0
  calc.DIC = TRUE   
        

        
  jags.data = list("u", "sigma2", "N", "num.prey", "num.iso", "num.pop","X","alpha","Region")  
  jags.params = c("p")
  jags.inits <- function(){
     list("mu"=runif(3), p.transform=matrix(runif(num.pop*num.prey),num.pop,num.prey), region.sig=0.5, ind.sig=0.5, p.ind = matrix(runif(N*num.prey), N, num.prey))
  }
  model.loc = "Model 7-clr.txt"
  jags.1 = jags(jags.data, inits = jags.inits, parameters.to.save= jags.params, model.file=model.loc, n.chains = mcmc.chains, n.burnin = mcmc.burn, n.thin = mcmc.thin, n.iter = mcmc.chainLength, DIC = TRUE, n.adapt = mcmc.adapt) 
Sculpin 0.2 | xhtml | problems or comments? | report bugs