Merge pull request #44 from bsvars/43-cran-10-additional-issues

43 cran 10 additional issues
bsvars · Aug 10, 2024 · efa0beb · efa0beb
2 parents a2c55e1 + 55a62a5
commit efa0beb
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 18 deletions.
diff --git a/src/restrictions_zero.cpp b/src/restrictions_zero.cpp
@@ -34,13 +34,15 @@ arma::colvec zero_restrictions(
     const arma::colvec            vec_structural
 ) {
   int N  = Z(0).n_cols;
-  mat A0 = reshape(vec_structural.rows(0, N*N-1), N, N);
+
+  mat A0 = reshape(vec_structural.rows(0, N * N - 1), N, N);
 
   arma::field<arma::mat> ZF = ZIRF(Z, inv(A0.t()));
 
-  colvec z;
+  vec z;
   for (int j=0; j<ZF.n_elem; j++) {
-    z = join_vert(z, ZF(j).col(j));
+    mat ZF_j = ZF(j);
+    if (ZF_j.n_rows > 0) z = join_vert(z, ZF_j.col(j));
   }
 
   return z;

diff --git a/src/sample_hyper.cpp b/src/sample_hyper.cpp
@@ -1,4 +1,5 @@
 
+#define ARMA_WARN_LEVEL 1
 #include <RcppArmadillo.h>
 
 #include "utils.h"
@@ -108,25 +109,26 @@ double log_ml(
   int N = Y.n_cols;
 
   double log_ml    = 0;
+
   mat    inv_Omega = diagmat(1 / Omega.diag());
+  mat    XX        = X.t() * X + inv_Omega;
 
-  try {
-    mat Bhat = inv_sympd(X.t() * X + inv_Omega) * (X.t() * Y + inv_Omega * b);
-    mat ehat = Y - X * Bhat;
-
-    log_ml += - N * T / 2.0 * log(M_PI);
-    log_ml += log_mvgamma(N, (T + d) / 2.0);
-    log_ml += -log_mvgamma(N, d / 2.0);
-    log_ml += - N / 2.0 * log_det_sympd(Omega);
-    log_ml += d / 2.0 * log_det_sympd(Psi);
-    log_ml += - N / 2.0 * log_det_sympd(X.t() * X + inv_Omega);
-    mat A   = Psi + ehat.t() * ehat + (Bhat - b).t() * inv_Omega * (Bhat - b);
-    log_ml += - (T + d) / 2.0 * log_det_sympd(A);
-
-  } catch(...) {
-    log_ml = -1e+10;
+  if (!Omega.is_sympd() or !Psi.is_sympd() or !XX.is_sympd()) {
+    return -1e10;
   }
 
+  mat Bhat = solve(XX, X.t() * Y + inv_Omega * b, solve_opts::likely_sympd);
+  mat ehat = Y - X * Bhat;
+
+  log_ml  += - N * T / 2.0 * log(M_PI);
+  log_ml  += log_mvgamma(N, (T + d) / 2.0);
+  log_ml  += -log_mvgamma(N, d / 2.0);
+  log_ml  += - N / 2.0 * log_det_sympd(Omega);
+  log_ml  += d / 2.0 * log_det_sympd(Psi);
+  log_ml  += - N / 2.0 * log_det_sympd(XX);
+  mat A    = Psi + ehat.t() * ehat + (Bhat - b).t() * inv_Omega * (Bhat - b);
+  log_ml  += - (T + d) / 2.0 * log_det_sympd(A);
+
   return log_ml;
 }