From 24759b4e6e8f80557aef315ac35b46c0538737a8 Mon Sep 17 00:00:00 2001 From: Franzi2114 Date: Tue, 1 Aug 2023 10:39:22 +0200 Subject: [PATCH] changed test --- stan/math/prim/functor/hcubature.hpp | 25 ++++---- stan/math/prim/prob/wiener_full_lpdf.hpp | 3 +- .../unit/math/prim/functor/hcubature_test.cpp | 62 +++++++++---------- 3 files changed, 43 insertions(+), 47 deletions(-) diff --git a/stan/math/prim/functor/hcubature.hpp b/stan/math/prim/functor/hcubature.hpp index 4f0867d6ad3..4066d4999a2 100644 --- a/stan/math/prim/functor/hcubature.hpp +++ b/stan/math/prim/functor/hcubature.hpp @@ -192,10 +192,9 @@ template std::pair gauss_kronrod(const F& integrand, const double a, const double b, const ParsPairT& pars_pair) { - std::vector c(1, 0); - std::vector cp(1, 0); - std::vector cm(1, 0); - c[0] = 0.5 * (a + b); + Eigen::VectorXd c {{0.5 * (a + b)}}; + Eigen::VectorXd cp(1); + Eigen::VectorXd cm(1); double delta = 0.5 * (b - a); double f0 = math::apply( [&integrand, &c](auto&&... args) { return integrand(c, args...); }, @@ -308,10 +307,10 @@ std::tuple integrate_GenzMalik( double twelvef1 = 12 * f1; double maxdivdiff = 0.0; - std::vector divdiff(dim); - std::vector p2(dim); - std::vector p3(dim); - std::vector cc(dim, 0); + Eigen::VectorXd divdiff(dim); + Eigen::VectorXd p2(dim); + Eigen::VectorXd p3(dim); + Eigen::VectorXd cc(dim); for (std::size_t i = 0; i != dim; i++) { for (std::size_t j = 0; j != dim; j++) { @@ -352,7 +351,7 @@ std::tuple integrate_GenzMalik( f3 += f3i; divdiff[i] = fabs(f3i + twelvef1 - 7 * f2i); } - std::vector p4(dim); + Eigen::VectorXd p4(dim); double f4 = 0.0; for (std::size_t i = 0; i != points[2].cols(); i++) { for (std::size_t j = 0; j != dim; j++) { @@ -367,7 +366,7 @@ std::tuple integrate_GenzMalik( f4 += temp; } double f5 = 0.0; - std::vector p5(dim); + Eigen::VectorXd p5(dim); for (std::size_t i = 0; i != points[3].cols(); i++) { for (std::size_t j = 0; j != dim; j++) { p5[j] = deltac[j] * points[3](j, i); @@ -521,12 +520,10 @@ double hcubature(const F& integrand, const ParsTuple& pars, const int dim, auto&& box = ms[err_idx]; double w = (box.b_[box.kdiv_] - box.a_[box.kdiv_]) / 2; - Eigen::VectorXd ma - = Eigen::Map(box.a_.data(), box.a_.size()); + Eigen::VectorXd ma = Eigen::Map(box.a_.data(), box.a_.size()); ma[box.kdiv_] += w; - Eigen::VectorXd mb - = Eigen::Map(box.b_.data(), box.b_.size()); + Eigen::VectorXd mb = Eigen::Map(box.b_.data(), box.b_.size()); mb[box.kdiv_] -= w; double result_1; diff --git a/stan/math/prim/prob/wiener_full_lpdf.hpp b/stan/math/prim/prob/wiener_full_lpdf.hpp index b0431524cb9..88e8e5246db 100644 --- a/stan/math/prim/prob/wiener_full_lpdf.hpp +++ b/stan/math/prim/prob/wiener_full_lpdf.hpp @@ -288,7 +288,7 @@ inline ReturnT wiener_full_lpdf(const T_y& y, const T_a& a, const T_t0& t0, T_st0>::value) { return 0; } - + using T_y_ref = ref_type_t; using T_a_ref = ref_type_t; using T_v_ref = ref_type_t; @@ -392,6 +392,7 @@ inline ReturnT wiener_full_lpdf(const T_y& y, const T_a& a, const T_t0& t0, } } + const double log_error_density = log(1e-6); // precision for density const double error_bound = precision_derivatives; // precision for // derivatives (controllable by user) diff --git a/test/unit/math/prim/functor/hcubature_test.cpp b/test/unit/math/prim/functor/hcubature_test.cpp index 26e0b647a91..2b1df83761d 100644 --- a/test/unit/math/prim/functor/hcubature_test.cpp +++ b/test/unit/math/prim/functor/hcubature_test.cpp @@ -105,8 +105,8 @@ stan::return_type_t f7(const T_x& x, double a) { template void test_integration(const F& f, const ArgsTupleT& pars, int dim, - std::vector a, std::vector b, int maxEval, - double reqAbsError, std::vector reqRelError, + const Eigen::VectorXd& a, const Eigen::VectorXd& b, int maxEval, + double reqAbsError, const Eigen::VectorXd& reqRelError, double val) { using stan::math::hcubature; @@ -122,46 +122,44 @@ TEST(StanMath_hcubature_prim, test1) { // https://www.quantargo.com/help/r/latest/packages/cubature/2.0.4.1/hcubature int dim = 1; - std::vector a = {0.0}; - std::vector b = {1.0}; - std::vector reqRelError = {1e-4, 1e-6, 1e-7}; - test_integration(hcubature_test::f1>, std::make_tuple(), + const Eigen::VectorXd a {{0.0}}; + const Eigen::VectorXd b {{1.0}}; + const Eigen::VectorXd reqRelError {{1e-4, 1e-6, 1e-7}}; + test_integration(hcubature_test::f1, std::make_tuple(), dim, a, b, 6000, 0.0, reqRelError, 0.841471); dim = 2; - a = {0.0, 0.0}; - b = {1.0, 1.0}; - reqRelError = {1e-4, 1e-6, 1e-7}; - test_integration(hcubature_test::f2>, std::make_tuple(), - dim, a, b, 6000, 0.0, reqRelError, 0.7080734); + const Eigen::VectorXd a_2 {{0.0, 0.0}}; + const Eigen::VectorXd b_2 {{1.0, 1.0}}; + test_integration(hcubature_test::f2, std::make_tuple(), + dim, a_2, b_2, 6000, 0.0, reqRelError, 0.7080734); - reqRelError = {1e-4}; - test_integration(hcubature_test::f3>, - std::make_tuple(0.50124145262344534123412), dim, a, b, 10000, - 0.0, reqRelError, 0.1972807); + const Eigen::VectorXd reqRelError_2 {{1e-4}}; + test_integration(hcubature_test::f3, + std::make_tuple(0.50124145262344534123412), dim, a_2, b_2, 10000, + 0.0, reqRelError_2, 0.1972807); // (Gaussian centered at 1/2) - reqRelError = {1e-4, 1e-6, 1e-7}; - test_integration(hcubature_test::f4>, - std::make_tuple(0.1), dim, a, b, 6000, 0.0, reqRelError, 1); + test_integration(hcubature_test::f4, + std::make_tuple(0.1), dim, a_2, b_2, 6000, 0.0, reqRelError, 1); dim = 3; - a = {0.0, 0.0, 0.0}; - b = {1.0, 1.0, 1.0}; - reqRelError = {1e-4, 1e-6}; - test_integration(hcubature_test::f5>, std::make_tuple(), - dim, a, b, 6000, 0.0, reqRelError, 1.00001); + const Eigen::VectorXd a_3 {{0.0, 0.0, 0.0}}; + const Eigen::VectorXd b_3 {{1.0, 1.0, 1.0}}; + const Eigen::VectorXd reqRelError_3 {{1e-4, 1e-6}}; + test_integration(hcubature_test::f5, std::make_tuple(), + dim, a_3, b_3, 6000, 0.0, reqRelError_3, 1.00001); - reqRelError = {1e-4, 1e-6, 1e-8}; - test_integration(hcubature_test::f6>, std::make_tuple(), - dim, a, b, 6000, 0.0, reqRelError, 1); + const Eigen::VectorXd reqRelError_4 {{1e-4, 1e-6, 1e-8}}; + test_integration(hcubature_test::f6, std::make_tuple(), + dim, a_3, b_3, 6000, 0.0, reqRelError_4, 1); // (Tsuda's example) dim = 4; - a = {0.0, 0.0, 0.0, 0.0}; - b = {1.0, 1.0, 1.0, 1.0}; - reqRelError = {1e-4, 1e-6}; - test_integration(hcubature_test::f7>, - std::make_tuple((1 + sqrt(10.0)) / 9.0), dim, a, b, 20000, - 0.0, reqRelError, 0.999998); + const Eigen::VectorXd a_4 {{0.0, 0.0, 0.0, 0.0}}; + const Eigen::VectorXd b_4 {{1.0, 1.0, 1.0, 1.0}}; + test_integration(hcubature_test::f7, + std::make_tuple((1 + sqrt(10.0)) / 9.0), dim, a_4, b_4, 20000, + 0.0, reqRelError_3, 0.999998); + }