From 80cfead22127b526151d3ed9cbdd4c9be0dce202 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Thu, 7 Dec 2023 06:48:12 -0500 Subject: [PATCH 1/4] add new test --- .../prim/prob/inv_wishart_cholesky_rng.hpp | 18 ++++++--- .../prob/inv_wishart_cholesky_rng_test.cpp | 39 +++++++++++++++++-- 2 files changed, 49 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/prob/inv_wishart_cholesky_rng.hpp b/stan/math/prim/prob/inv_wishart_cholesky_rng.hpp index d8078e8907d..f1fed4dff62 100644 --- a/stan/math/prim/prob/inv_wishart_cholesky_rng.hpp +++ b/stan/math/prim/prob/inv_wishart_cholesky_rng.hpp @@ -3,9 +3,7 @@ #include #include -#include -#include -#include +#include 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(L_S); - return mdivide_left_tri(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 diff --git a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp index 65464672a07..f9be95c8fac 100644 --- a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp +++ b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp @@ -91,14 +91,14 @@ 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 +111,36 @@ 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::multi_normal_rng; + using stan::math::inv_wishart_rng; + using stan::math::lkj_corr_cholesky_rng; + using stan::math::qr_thin_R; + using stan::math::multiply_lower_tri_self_transpose; + + boost::random::mt19937 rng(92343U); + int N = 1e5; + double tol = 0.1; + for (int k = 1; k < 5; k++) { + MatrixXd sigma = inv_wishart_rng(15, MatrixXd::Identity(k, k), rng); + MatrixXd L = stan::math::cholesky_decompose(sigma); + MatrixXd Z_mean = sigma / (k + 3); + MatrixXd Z_est = MatrixXd::Zero(k, k); + for (int i = 0; i < N; i++) { + Z_est += inv_wishart_cholesky_rng(k + 4, L, rng); + } + Z_est /= N; + Z_est = multiply_lower_tri_self_transpose(Z_est); + 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); + } + } + } +} \ No newline at end of file From b1aedefada8ce27383b04364e60daefd8a873417 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Thu, 7 Dec 2023 07:02:23 -0500 Subject: [PATCH 2/4] Update inv_wishart_cholesky_rng_test.cpp --- .../prim/prob/inv_wishart_cholesky_rng_test.cpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp index f9be95c8fac..521249a926c 100644 --- a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp +++ b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp @@ -118,25 +118,23 @@ TEST(ProbDistributionsInvWishartCholesky, compareToInvWishart) { using Eigen::MatrixXd; using Eigen::VectorXd; using stan::math::inv_wishart_cholesky_rng; - using stan::math::multi_normal_rng; using stan::math::inv_wishart_rng; - using stan::math::lkj_corr_cholesky_rng; - using stan::math::qr_thin_R; 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.1; + double tol = 0.05; for (int k = 1; k < 5; k++) { - MatrixXd sigma = inv_wishart_rng(15, MatrixXd::Identity(k, k), rng); - MatrixXd L = stan::math::cholesky_decompose(sigma); + 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 += inv_wishart_cholesky_rng(k + 4, L, rng); + Z_est += multiply_lower_tri_self_transpose(inv_wishart_cholesky_rng(k + 4, L, rng)); } Z_est /= N; - Z_est = multiply_lower_tri_self_transpose(Z_est); 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); From 80168bb90c39c6350082a9d1108dc0b6f14cae75 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 7 Dec 2023 07:04:18 -0500 Subject: [PATCH 3/4] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/prob/inv_wishart_cholesky_rng.hpp | 6 +++--- .../math/prim/prob/inv_wishart_cholesky_rng_test.cpp | 10 ++++++---- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/stan/math/prim/prob/inv_wishart_cholesky_rng.hpp b/stan/math/prim/prob/inv_wishart_cholesky_rng.hpp index f1fed4dff62..a3376db5326 100644 --- a/stan/math/prim/prob/inv_wishart_cholesky_rng.hpp +++ b/stan/math/prim/prob/inv_wishart_cholesky_rng.hpp @@ -14,9 +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). - * + * 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 diff --git a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp index 521249a926c..536911e174d 100644 --- a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp +++ b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp @@ -98,7 +98,8 @@ TEST(ProbDistributionsInvWishartCholesky, SpecialRNGTest) { MatrixXd L = MatrixXd::Identity(k, k); MatrixXd Z = MatrixXd::Zero(k, k); for (int i = 0; i < N; i++) { - Z += multiply_lower_tri_self_transpose(inv_wishart_cholesky_rng(k + 2, L, rng)); + Z += multiply_lower_tri_self_transpose( + inv_wishart_cholesky_rng(k + 2, L, rng)); } Z /= N; for (int j = 0; j < k; j++) { @@ -113,7 +114,7 @@ TEST(ProbDistributionsInvWishartCholesky, SpecialRNGTest) { } TEST(ProbDistributionsInvWishartCholesky, compareToInvWishart) { - // Compare the marginal mean + // Compare the marginal mean using Eigen::MatrixXd; using Eigen::VectorXd; @@ -132,12 +133,13 @@ TEST(ProbDistributionsInvWishartCholesky, compareToInvWishart) { 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 += 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); + EXPECT_NEAR(Z_est(i, j), Z_mean(i, j), tol); } } } From 428bd1569207fd819af0483e2e5295cb723b7902 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Thu, 7 Dec 2023 07:16:32 -0500 Subject: [PATCH 4/4] add new line at eof --- test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp index 521249a926c..ca1b847ec04 100644 --- a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp +++ b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp @@ -141,4 +141,4 @@ TEST(ProbDistributionsInvWishartCholesky, compareToInvWishart) { } } } -} \ No newline at end of file +}