###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 <- 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, 1, 1, 1) VbetaInv <- diag(1/c(tsigmaP0, tsigmaP1)) ##Variance of betas for proxies VbetaInv.t <- diag(1/rep(tsigmaT, 4)) ##Variance of betas in process stage qP <- 3 ##For inverse gamma rP <- 1 qT <- 3 rT <- 1 sigma2P6 <- 1/64 xi1r <- 0.0001 ##Tree|T xi2r <- 0.0001 ##Tree|T xi1p <- 0.0001 ##poll|T xi2p <- 0.0001 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/new9/commonv.R') ###################################################################### # Storage Setting # ###################################################################### nsave <- ngibbs-nburn Tsave <- matrix(0, nsave, Ty1) #Psave <- matrix(0, nsave, np) #Bsave <- matrix(0, nsave, nz) Vsave <- matrix(0, nsave, Ty) betasave <- matrix(0, ngibbs, 4) ### (beta0,beta1,beta2,beta3) beta0save <- matrix(0, ngibbs, nTr) ### (beta01, beta02, beta03) beta1save <- matrix(0, ngibbs, nTr) ### (beta11, beta12, beta13) phisave <- matrix(0, ngibbs, 4) ## (phi1r, phi2r, phi1p, phi2p, phi1t, phi2t) sigma2Psave <- matrix(0, ngibbs, 1) ## (sigma2P3--sigma2P5) sigma2Tsave <- rep(0, ngibbs) ###################################################################### # MCMC sampling # ###################################################################### flag1 <- flag2 <- flag3 <- flag4 <- 0 #how often does M-H accept the new value saveind <- 1 for (j in 1:ngibbs){ ################################################################### ### sample T1 tmp <- cbind(1, Sol, V, Co2)%*%betat tmp1 <- t(M.tree) %*% VPinv.r D1.e <- tmp1[1:Ty1,] D1.t <- VTinv[1:Ty1,] tmp1 <- tmp1%*%M.tree V11.e <- tmp1[1:Ty1,1:Ty1] V11.t <- VTinv[1:Ty1, 1:Ty1] V12.e <- tmp1[1:Ty1, (Ty1+1):Ty] V12.t <- VTinv[1:Ty1, (Ty1+1):Ty] Tcovinv <- sum(beta11^2)*V11.e + V11.t tmp11 <- D1.e %*% ((p.tree-rep(1,Ty)%*%t(beta01))%*%beta11)-sum(beta11^2)*V12.e%*%tempt2 Tb <- tmp11 + 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)) #T0 <- T #T1 <- T[1:Ty1] T <- c(T1,tempt2) ################################################################### ### sample Volcanism Vcovinv <- as.spam(diag(1/(sigma2P6*as.vector(V^2)))+betat[3]^2*VTinv) Vb <- Vol/(sigma2P6*V^2) + betat[3]*VTinv%*%(T-betat[1]-betat[2]*Sol-betat[4]*Co2) Vmn <- solve(Vcovinv, Vb) V <- Vmn + backsolve(chol(Vcovinv), rnorm(Ty)) ####################################################### # 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 beta0, beta1, beta2, beta3 Ftild <- cbind(1, Sol, V, Co2) 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(4)) ####################################################### # 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 phi1t, phi2t resid.t <- T - Ftild %*% betat 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 sigma2P6 ################################################################### ### 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) beta1save[j,] <- c(beta11) phisave[j,] <- c(phi1r, phi2r, phi1t, phi2t) sigma2Psave[j,] <- sigma2P3 sigma2Tsave[j] <- sigma2T if (j > nburn){ Tsave[saveind,] <- T1 #Psave[saveind,] <- P[,1] #Bsave[saveind,] <- B[,1] Vsave[saveind,] <- V 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/new9/main2v.RData')