fixes to gibbs step

NAMESPACE

@@ -14,6 +14,6 @@ importFrom("RcppParallel", "CxxFlags", "RcppParallelLibs")
 
 import(rstantools)
 
-export("bsam", "bsam_model")
+export("bsam", "stanmodels")
 
 useDynLib(bsam, .registration = TRUE)

R/bsam.R

@@ -1,15 +1,16 @@
 ## BSAM: Bayesian Structural After Measurement
 ## uses a new Stan file, then uses blavaan for everything else
 
-bsam <- function(...) {
+bsam <- function(..., ngibbs = 50L) {
 
   dotdotdot <- list(...)
 
   if (!("mcmcextra" %in% names(dotdotdot))) {
-    mcmcextra <- list(dosam = TRUE)
+    mcmcextra <- list(dosam = TRUE, data = list(ngibbs = ngibbs, fullpsi_c = 0))
   } else {
     mcmcextra <- dotdotdot$mcmcextra
     mcmcextra$dosam <- TRUE
+    mcmcextra$data <- c(mcmcextra$data, list(ngibbs = ngibbs, fullpsi_c = 0L))
   }
 
   mc <- match.call()

R/stanmodels.R

@@ -7,7 +7,7 @@ stanmodels <- c("stanmarg_bsam")
 Rcpp::loadModule("stan_fit4stanmarg_bsam_mod", what = TRUE)
 
 # instantiate each stanmodel object
-bsam_model <- sapply(stanmodels, function(model_name) {
+stanmodels <- sapply(stanmodels, function(model_name) {
   # create C++ code for stan model
   stan_file <- if(dir.exists("stan")) "stan" else file.path("inst", "stan")
   stan_file <- file.path(stan_file, paste0(model_name, ".stan"))

inst/stan/stanmarg_bsam.stan

@@ -1261,7 +1261,6 @@ generated quantities { // these matrices are saved in the output but do not figu
   vector[len_free[4]] b_primn;
   vector[len_free[14]] alpha_primn;
   array[Ng] int matdim;
-  array[Ng] vector[m] corrfac; // correction for std.lv
 
   array[Ntot] vector[m] eta;
 
@@ -1350,7 +1349,6 @@ generated quantities { // these matrices are saved in the output but do not figu
     b_prior_prec[g] = fill_matrix(pow(to_vector(b_sd), -2), B_skeleton[g], w4skel, g_start4[g,1], g_start4[g,2]);
     Psi_prior_shape[g] = fill_matrix(to_vector(psi_sd_shape), Psi_skeleton[g], w9skel, g_start9[g,1], g_start9[g,2]);
     Psi_prior_rate[g] = fill_matrix(to_vector(psi_sd_rate), Psi_skeleton[g], w9skel, g_start9[g,1], g_start9[g,2]);
-    corrfac[g] = rep_vector(1, m);
 
     if (g == Ng) {
       matdim[g] = (len_free[4] - g_start4[g, 1] + 1) + (len_free[14] - g_start14[g, 1] + 1);
@@ -1434,16 +1432,12 @@ generated quantities { // these matrices are saved in the output but do not figu
       Dchol = cholesky_decompose(D);
       d = to_vector(Psi0_inv * IBinv * Alpha[g]);
 
-      if (Ndum_x[mm] == 0) {
-        for (ridx in r1:r2) {
-	  eta[ridx] = multi_normal_cholesky_rng(D * (d + Lamt_Thet_inv * (YX[ridx] - to_vector(Nu[g]))), Dchol) ./ corrfac[g];
-	}
-      } else {
-	for (ridx in r1:r2) {
-	  eta[ridx] = multi_normal_cholesky_rng(D * (d + to_vector(Psi0_inv * IBinv * B[g, , dum_lv_x_idx[mm, 1:Ndum_x[mm]]] * YX[ridx, dum_ov_x_idx[mm, 1:Ndum_x[mm]]]) + Lamt_Thet_inv * (YX[ridx] - to_vector(Nu[g]))), Dchol) ./ corrfac[g];
-
+      for (ridx in r1:r2) {
+	eta[ridx] = multi_normal_cholesky_rng(D * (d + Lamt_Thet_inv * (YX[ridx] - to_vector(Nu[g]))), Dchol);
+
+	if (Ndum[mm] > 0) {
 	  eta[ridx, dum_lv_idx[mm, 1:Ndum[mm]]] = YX[ridx, dum_ov_idx[mm, 1:Ndum[mm]]];
-        }
+	}
       }
 
       // sample alpha, beta
@@ -1488,7 +1482,6 @@ generated quantities { // these matrices are saved in the output but do not figu
       params = multi_normal_rng(Dinv * FVz, Dinv);
 
       // now put parameters in free parameter vectors
-      // FIXME? fill_matrix works by columns instead of rows?
       pidx = 1;
       for (r in 1:m) {
 	real askel = Alpha_skeleton[g, r, 1];
@@ -1552,11 +1545,8 @@ generated quantities { // these matrices are saved in the output but do not figu
 	    Psi[gg, srow:erow, srow:erow] = inv_wishart_rng((r2 - r1 + 1) + Psi_prior_shape[gg, srow, srow], residcp[srow:erow, srow:erow] + Psi_prior_rate[gg, srow:erow, srow:erow]);
 	  } else {
 	    // line numbers are off by about 156 here
-	    if (is_inf(Psi_skeleton[gg, srow, srow])) {
-	      Psi[gg, srow, srow] = inv_gamma_rng(.5 * (r2 - r1 + 1) + Psi_prior_shape[gg, srow, srow], .5 * residcp[srow, srow] + Psi_prior_rate[gg, srow, srow]);
-	    } else if (Psi_skeleton[gg, srow, srow] == 1) {
-	      corrfac[gg, srow] = pow(inv_gamma_rng(.5 * (r2 - r1 + 1) + Psi_prior_shape[gg, srow, srow], .5 * residcp[srow, srow] + Psi_prior_rate[gg, srow, srow]), .5);
-	    }
+	    // even if Psi was originally fixed, we sample in this step
+	    Psi[gg, srow, srow] = inv_gamma_rng(.5 * (r2 - r1 + 1) + Psi_prior_shape[gg, srow, srow], .5 * residcp[srow, srow] + Psi_prior_rate[gg, srow, srow]);
 	  }
 	}
       }