sampling tauSq before sampling gammas

Author: zhizuio

R/GPTCM.R

@@ -26,7 +26,7 @@
 #' common (\code{dirichlet = TRUE}) or alternative (\code{dirichlet = FALSE})
 #' parametrization of the Dirichlet regression model
 #' @param hyperpar a list of relevant hyperparameters
-#' @param BVS logical value for implementing Bayesian variable selection//' 
+#' @param BVS logical value for implementing Bayesian variable selection
 #' @param threads maximum threads used for parallelization. Default is 1
 #' @param kappaIGamma logical value for using inverse-gamma prior (\code{TRUE})
 #' or gamma prior (\code{FALSE}) for Weibull's shape parameter

src/arms_gibbs.cpp

@@ -240,7 +240,7 @@ void ARMS_Gibbs::arms_gibbs_beta(
     // mydata->xr = xr;
     // mydata->vA = vA;
     // mydata->vB = vB;
-    // mydata->tauSq = tauSq;
+    mydata->tau0Sq = tau0Sq;
     mydata->kappa = kappa;
     mydata->datTheta = datTheta.memptr();
     mydata->datMu = datMu.memptr();
@@ -257,10 +257,7 @@ void ARMS_Gibbs::arms_gibbs_beta(
         // Gibbs sampling
         // mydata->tauSq = sampleW(tauA, tauB, currentPars.col(l));
         // mydata->tau0Sq = sampleW0(tauA, tauB, currentPars(0,l));
-        
-        tau0Sq = sampleV(tau0A, tau0B, currentPars.row(0).t());
-        mydata->tau0Sq = tau0Sq;
-        tauSq[l] = sampleV(tauA, tauB, currentPars.submat(1,l,p,l));
+
         mydata->tauSq = tauSq[l];
         // tauSq_tmp[l] = sampleW(tauA, tauB, currentPars.col(l));
         // mydata->tauSq = tauSq_tmp[l];
@@ -556,7 +553,7 @@ void ARMS_Gibbs::arms_gibbs_zeta(
     mydata->p = p;
     mydata->L = L;
     mydata->N = N;
-    // mydata->w0Sq = w0Sq;
+    mydata->w0Sq = w0Sq;
     // mydata->wSq = wSq;
     //mydata->phi = phi;
     //mydata->dirichlet = dirichlet,
@@ -576,9 +573,6 @@ void ARMS_Gibbs::arms_gibbs_zeta(
     {
         // Gibbs sampling
         // mydata->w0Sq = sampleW0(w0A, w0B, currentPars(0,l));
-        w0Sq = sampleV(w0A, w0B, currentPars.row(0).t());
-        mydata->w0Sq = w0Sq;
-        wSq[l] = sampleV(wA, wB, currentPars.submat(1,l,p,l));
         mydata->wSq = wSq[l];
         // wSq_tmp[l] = sampleW(wA, wB, currentPars.submat(1,l,p,l));
         // mydata->wSq = wSq_tmp[l];

src/drive.cpp

@@ -481,6 +481,16 @@ Rcpp::List run_mcmc(
                         std::string(cTotalLength * (1. - (m + 1.) / nIter), '-') << // printing empty part
                         "] " << (int)((m + 1.) / nIter * 100.0) << "%\r";             // printing percentage
 
+
+        // update \beta's and \zeta's variacnes tau0Sq, tauSq, w0Sq, wSq
+        tau0Sq = sampleV(hyperpar->tau0A, hyperpar->tau0B, betas.row(0).t());
+        w0Sq = sampleV(hyperpar->w0A, hyperpar->w0B, zetas.row(0).t());
+        for (unsigned int l=0; l<L; ++l)
+        {
+            tauSq[l] = sampleV(hyperpar->tauA, hyperpar->tauB, betas.submat(1,l,p,l));
+            wSq[l] = sampleV(hyperpar->wA, hyperpar->wB, zetas.submat(1,l,p,l));
+        }
+
         // update \xi's variance vSq
         v0Sq = sampleV0(hyperpar->v0A, hyperpar->v0B, xi[0]);
         // hyperpar->v0Sq = v0Sq;