###number of MCMC runs and number of burn in time ngibbs <- 3000 nburn <- 2000 ###The threshold of frequency for tree/poll freq.tree <- 10 freq.poll <- c(10, 100) ###How often is poll collected poll.ind <- seq(30,1150,30) #poll.ind <- seq(10,1150-10,10) #poll.ind <- 1:1150 ###################################################################### # Hyperparameters # ###################################################################### #mu <- 0 tsigmaT <- 1 ##variance for beta's in the process model tsigmaP0 <- 1 ## variance for beta's in the data model tsigmaP1 <- 1 mubetap <- c(0, 1) mubetat <- c(0) VbetaInv <- diag(1/c(tsigmaP0, tsigmaP1)) ##Variance of betas for proxies VbetaInv.t <- diag(1/tsigmaT) ##Variance of betas in process stage qP <- 3 ##For inverse gamma rP <- 1 qT <- 3 rT <- 1 xi1r <- 0.0001 ##Tree|T xi2r <- 0.0001 ##Tree|T xi1p <- 0.01 ##poll|T xi2p <- 0.01 xi1b <- 0.005 ##bore|T xi2b <- 0.005 xi1t <- 0.001 ##T process xi2t <- 0.001 set.seed(123) source('/fs/image/home/boli/Hockey/HBM/combination/new6/common.R') betat <- beta0 ###################################################################### # Storage Setting # ###################################################################### nsave <- ngibbs-nburn Tsave <- matrix(0, nsave, Ty1) #Psave <- matrix(0, nsave, np) #Bsave <- matrix(0, nsave, nz) #Ssave <- matrix(0, nsave, Ty) betasave <- matrix(0, ngibbs, 1) ### (beta0,beta1,beta2,beta3) beta0save <- matrix(0, ngibbs, nTr+nPl) ### (beta01, beta02, beta03) beta1save <- matrix(0, ngibbs, nTr+nPl) ### (beta11, beta12, beta13) phisave <- matrix(0, ngibbs, 6) ## (phi1r, phi2r, phi1p, phi2p, phi1t, phi2t) sigma2Psave <- matrix(0, ngibbs, 2) ## (sigma2P3--sigma2P5) sigma2Tsave <- rep(0, ngibbs) ###################################################################### # MCMC sampling # ###################################################################### flag1 <- flag2 <- flag3 <- flag4 <- 0 #how often does M-H accept the new value saveind <- 1 M.poll <- as.matrix(M.poll) for (j in 1:ngibbs){ ################################################################### ### sample T1 tmp <- rep(betat, nrow(VTinv)) tmp1 <- t(M.tree) %*% VPinv.r tmp2 <- t(M.poll) %*% VPinv.p D1.e <- tmp1[1:Ty1,] D1.p <- tmp2[1:Ty1,] D1.t <- VTinv[1:Ty1,] tmp1 <- tmp1%*%M.tree tmp2 <- tmp2%*%M.poll V11.e <- tmp1[1:Ty1,1:Ty1] V11.p <- tmp2[1:Ty1,1:Ty1] V11.t <- VTinv[1:Ty1, 1:Ty1] V12.e <- tmp1[1:Ty1, (Ty1+1):Ty] V12.p <- tmp2[1:Ty1, (Ty1+1):Ty] V12.t <- VTinv[1:Ty1, (Ty1+1):Ty] Tcovinv <- sum(beta11^2)*V11.e + sum(beta12^2)*V11.p + V11.t tmp11 <- D1.e %*% ((p.tree-rep(1,Ty)%*%t(beta01))%*%beta11)-sum(beta11^2)*V12.e%*%tempt2 tmp21 <- D1.p %*% ((p.poll-rep(1,np)%*%t(beta02))%*%beta12)-sum(beta12^2)*V12.p%*%tempt2 Tb <- tmp11 + tmp21 + D1.t%*%tmp - V12.t%*%tempt2 Tmn <- solve(0.5*(Tcovinv+t(Tcovinv)), Tb) T1 <-Tmn + backsolve(chol(0.5*(Tcovinv+t(Tcovinv))), rnorm(Ty1)) T <- c(T1,tempt2) ####################################################### # Sample coefficients # ####################################################### ################################################################### ### sample beta01, beta11 Ttild.r <- cbind(1, M.tree%*%T) tmp <- t(Ttild.r) %*% VPinv.r betacov <- solve(tmp %*% Ttild.r + VbetaInv) betamn <- betacov%*%(tmp%*%p.tree + as.vector(VbetaInv%*%mubetap)) beta <- betamn+crossprod(chol(0.5*(betacov+t(betacov))), matrix(rnorm(2*nTr),2,nTr)) ##no change beta01 <- beta[1,] beta11 <- beta[2,] ################################################################### ### sample beta02, beta12 Ttild.p <- cbind(1, M.poll%*%T) tmp <- t(Ttild.p) %*% VPinv.p betacov <- solve(tmp %*% Ttild.p + VbetaInv) ##no change betamn <- betacov%*%(tmp%*%P + as.vector(VbetaInv%*%mubetap)) beta <- betamn+crossprod(chol(0.5*(betacov+t(betacov))), matrix(rnorm(2*nPl),2,nPl)) ##no change beta02 <- beta[1,] beta12 <- beta[2,] ################################################################### ### sample beta0, beta1, beta2, beta3 Ftild <- rep(1, nrow(VTinv)) tmp <- t(Ftild) %*% VTinv betacov <- solve(tmp %*% Ftild + VbetaInv.t) betamn <- betacov%*%(tmp%*%T + VbetaInv.t%*%mubetat) betat <- betamn+crossprod(chol(0.5*(betacov+t(betacov))),rnorm(1)) ####################################################### # Sample time correlation coefficients # ####################################################### ################################################################### ### sample phi1r, phi2r resid.r <- p.tree - Ttild.r %*% rbind(beta01, beta11) tmp.r <- phi4ar2(0, resid.r, phi1r, phi2r, xi1r, xi2r, sigma2P3, VPinv.r, flag1, Ty,nTr) phi1r <- tmp.r$phi1 phi2r <- tmp.r$phi2 VPinv.r <- tmp.r$VPinv flag1 <- tmp.r$flag ### sample phi1p, phi2p resid.p <- P - Ttild.p %*% rbind(beta02, beta12) tmp.p <- phi4ar2(0, resid.p, phi1p, phi2p, xi1p, xi2p, sigma2P4, VPinv.p, flag2, np,nPl) phi1p <- tmp.p$phi1 phi2p <- tmp.p$phi2 VPinv.p <- tmp.p$VPinv flag2 <- tmp.p$flag ### sample phi1t, phi2t resid.t <- T - betat*Ftild tmp.t <- phi4ar2(0, resid.t, phi1t, phi2t, xi1t, xi2t, sigma2T, VTinv, flag4, Ty,1) ##for process the last parameter is 1 phi1t <- tmp.t$phi1 phi2t <- tmp.t$phi2 VTinv <- tmp.t$VPinv flag4 <- tmp.t$flag ####################################################### # Sample variance paramters # ####################################################### ################################################################### ### sample sigma2P3 sigma2Pold <- sigma2P3 qq <- qP+nTr*Ty/2 rr <- 1/(1/rP+0.5*sigma2Pold*sum(diag(t(resid.r)%*%VPinv.r%*%resid.r))) sigma2P3 <- 1/rgamma(1, qq,, rr) VPinv.r <- VPinv.r*sigma2Pold/sigma2P3 ################################################################### ### sample sigma2P4 sigma2Pold <- sigma2P4 qq <- qP+nPl*np/2 rr <- 1/(1/rP+0.5*sigma2Pold*sum(diag(t(resid.p)%*%VPinv.p%*%resid.p))) sigma2P4 <- 1/rgamma(1, qq,, rr) VPinv.p <- VPinv.p*sigma2Pold/sigma2P4 ################################################################### ### sample sigma2T sigma2Told <- sigma2T qq <- qT+Ty/2 rr <- 1/(1/rT+0.5*sigma2Told*t(resid.t)%*%VTinv%*%resid.t) sigma2T <- 1/rgamma(1, qq,, rr) VTinv <- VTinv*sigma2Told/sigma2T ################################################################### ### save posteriors betasave[j,] <- betat beta0save[j,] <- c(beta01, beta02) beta1save[j,] <- c(beta11, beta12) phisave[j,] <- c(phi1r, phi2r, phi1p, phi2p, phi1t, phi2t) sigma2Psave[j,] <- c(sigma2P3,sigma2P4) sigma2Tsave[j] <- sigma2T if (j > nburn){ Tsave[saveind,] <- T1 #Psave[saveind,] <- P[,1] #Bsave[saveind,] <- B[,1] #Ssave[saveind,] <- S saveind <- saveind+1 } print(j) } #end of MC rm(Tcovinv, M.tree, erfc) #rm(VPnew, VPnewInv, VTnew, VTnewInv, Tcovinv, Scovinv) #rm(M.tree, coeffM) save.image('/fs/image/home/boli/Hockey/HBM/combination/new6/main8.RData')