# continue with cyc6, but Gibbs sample V, W1, W2, W3, and the states. mod <- mod.cyc6 data <- A535$Region.1 n.data <- length(data) # number of time points MC <- 100 # number of MCMC iterations keep <- 1:MC n.keep <- length(keep) # number of iterations to save # PARAMETERS # Theta: state vector # Tau: observation precision # Winvpoly: locally linear precision # Winvheart: heart precision # Winvresp: respiration precision # Winvweird: weird precision # PRIORS tau.a <- 1; tau.b <- 1 winvpoly.a <- 1; winvpoly.b <- 1 winvheart.a <- 1; winvheart.b <- 1 winvresp.a <- 1; winvresp.b <- 1 winvweird.a <- 1; winvweird.b <- 1 # STORAGE FOR MCMC OUTPUT Theta.Gibbs <- array ( NA, dim = c ( n.keep, n.data+1, length(mod$m0) ) ) Tau.Gibbs <- rep ( NA, n.keep ) Winvpoly.Gibbs <- rep ( NA, n.keep ) Winvheart.Gibbs <- rep ( NA, n.keep ) Winvresp.Gibbs <- rep ( NA, n.keep ) Winvweird.Gibbs <- rep ( NA, n.keep ) # THE SAMPLER for ( i in 1:MC ) { %print ( paste ( "beginning loop", i ) ) # sample theta filt <- dlmFilter ( data, mod, simplify=TRUE ) theta <- dlmBSample ( filt ) theta.pred <- theta[-n.data+1,] %*% t(mod$GG) theta.res <- theta[-1,] - theta.pred # sample tau fit <- theta[-1,] %*% t(mod$FF) rss <- sum ( (data-fit)^2 ) tau <- rgamma ( 1, shape=tau.a+n.data/2, rate=tau.b+rss/2 ) mod$V <- 1/tau # sample winvpoly rss <- sum ( theta.res[,2]^2 ) winvpoly <- rgamma ( 1, shape = winvpoly.a + n.data/2, rate = winvpoly.b + rss/2 ) mod$W[2,2] <- 1/winvpoly # sample winvheart rss <- sum ( theta.res[,3:4]^2 ) winvheart <- rgamma ( 1, shape = winvheart.a + n.data, rate = winvheart.b + rss/2 ) diag(mod$W)[3:4] <- 1/winvheart # sample winvresp rss <- sum ( theta.res[,5:8]^2 ) winvresp <- rgamma ( 1, shape = winvresp.a + 2*n.data, rate = winvresp.b + rss/2 ) diag(mod$W)[5:8] <- 1/winvresp # sample winvweird rss <- sum ( theta.res[,9:12]^2 ) winvweird <- rgamma ( 1, shape = winvweird.a + 2*n.data, rate = winvweird.b + rss/2 ) diag(mod$W)[9:12] <- 1/winvweird # keep this iteration? if ( !is.na ( j <- match ( i, keep ) ) ) { Theta.Gibbs[j,,] <- theta Tau.Gibbs[j] <- tau Winvpoly.Gibbs[j] <- winvpoly Winvheart.Gibbs[j] <- winvheart Winvresp.Gibbs[j] <- winvresp Winvweird.Gibbs[j] <- winvweird } %print ( paste ( "ending loop", i ) ) }