Merge pull request #2984 from stan-dev/fix-inv-wishart-cholesky-rng

Fix inv wishart cholesky rng

Author: andrjohns

stan/math/prim/prob/inv_wishart_cholesky_rng.hpp

@@ -3,9 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
-#include <stan/math/prim/fun/mdivide_left_tri.hpp>
-#include <stan/math/prim/prob/wishart_cholesky_rng.hpp>
-#include <stan/math/prim/prob/wishart_rng.hpp>
+#include <stan/math/prim/fun/mdivide_left_tri_low.hpp>
 
 namespace stan {
 namespace math {
@@ -16,6 +14,9 @@ namespace math {
  * from the inverse Wishart distribution with the specified degrees of freedom
  * using the specified random number generator.
  *
+ * Axen, Seth D. "Efficiently generating inverse-Wishart matrices and their
+ * Cholesky factors." arXiv preprint arXiv:2310.15884 (2023).
+ *
  * @tparam RNG Random number generator type
  * @param[in] nu scalar degrees of freedom
  * @param[in] L_S lower Cholesky factor of the scale matrix
@@ -38,8 +39,15 @@ inline Eigen::MatrixXd inv_wishart_cholesky_rng(double nu,
   check_positive(function, "Cholesky Scale matrix", L_S.diagonal());
   check_positive(function, "columns of Cholesky Scale matrix", L_S.cols());
 
-  MatrixXd L_Sinv = mdivide_left_tri<Eigen::Lower>(L_S);
-  return mdivide_left_tri<Eigen::Lower>(wishart_cholesky_rng(nu, L_Sinv, rng));
+  MatrixXd B = MatrixXd::Zero(k, k);
+  for (int j = 0; j < k; ++j) {
+    for (int i = 0; i < j; ++i) {
+      B(j, i) = normal_rng(0, 1, rng);
+    }
+    B(j, j) = std::sqrt(chi_square_rng(nu - k + j + 1, rng));
+  }
+
+  return mdivide_left_tri_low(B, L_S);
 }
 
 }  // namespace math

test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp

@@ -91,14 +91,15 @@ TEST(ProbDistributionsInvWishartCholesky, SpecialRNGTest) {
   using stan::math::inv_wishart_cholesky_rng;
   using stan::math::multiply_lower_tri_self_transpose;
 
-  boost::random::mt19937 rng(1234U);
+  boost::random::mt19937 rng(92343U);
   int N = 1e5;
   double tol = 0.1;
   for (int k = 1; k < 5; k++) {
-    MatrixXd sigma = MatrixXd::Identity(k, k);
+    MatrixXd L = MatrixXd::Identity(k, k);
     MatrixXd Z = MatrixXd::Zero(k, k);
     for (int i = 0; i < N; i++) {
-      Z += stan::math::crossprod(inv_wishart_cholesky_rng(k + 2, sigma, rng));
+      Z += multiply_lower_tri_self_transpose(
+          inv_wishart_cholesky_rng(k + 2, L, rng));
     }
     Z /= N;
     for (int j = 0; j < k; j++) {
@@ -111,3 +112,35 @@ TEST(ProbDistributionsInvWishartCholesky, SpecialRNGTest) {
     }
   }
 }
+
+TEST(ProbDistributionsInvWishartCholesky, compareToInvWishart) {
+  // Compare the marginal mean
+
+  using Eigen::MatrixXd;
+  using Eigen::VectorXd;
+  using stan::math::inv_wishart_cholesky_rng;
+  using stan::math::inv_wishart_rng;
+  using stan::math::multiply_lower_tri_self_transpose;
+  using stan::math::qr_thin_Q;
+
+  boost::random::mt19937 rng(92343U);
+  int N = 1e5;
+  double tol = 0.05;
+  for (int k = 1; k < 5; k++) {
+    MatrixXd L = qr_thin_Q(MatrixXd::Random(k, k)).transpose();
+    L.diagonal() = stan::math::abs(L.diagonal());
+    MatrixXd sigma = multiply_lower_tri_self_transpose(L);
+    MatrixXd Z_mean = sigma / (k + 3);
+    MatrixXd Z_est = MatrixXd::Zero(k, k);
+    for (int i = 0; i < N; i++) {
+      Z_est += multiply_lower_tri_self_transpose(
+          inv_wishart_cholesky_rng(k + 4, L, rng));
+    }
+    Z_est /= N;
+    for (int j = 0; j < k; j++) {
+      for (int i = 0; i < j; i++) {
+        EXPECT_NEAR(Z_est(i, j), Z_mean(i, j), tol);
+      }
+    }
+  }
+}