diff --git a/stan/math/prim/functor/apply.hpp b/stan/math/prim/functor/apply.hpp index 48673e5b98a..e0c1c7c73bf 100644 --- a/stan/math/prim/functor/apply.hpp +++ b/stan/math/prim/functor/apply.hpp @@ -10,20 +10,26 @@ namespace math { namespace internal { /* * Invoke the functor f with arguments given in t and indexed in the index - * sequence I + * sequence I with other arguments possibly before or after * * @tparam F Type of functor * @tparam Tuple Type of tuple containing arguments + * @tparam PreArgs Parameter pack of arguments before the tuple * @tparam I Parameter pack of an index sequence going from 0 to * std::tuple_size::value - 1 inclusive * @param f functor callable * @param t tuple of arguments * @param i placeholder variable for index sequence + * @param pre_args parameter pack of arguments to place before elements in + * tuple. */ -template +template constexpr decltype(auto) apply_impl(F&& f, Tuple&& t, - std::index_sequence i) { - return f(std::forward(t))>(std::get(t))...); + std::index_sequence i, + PreArgs&&... pre_args) { + return std::forward(f)( + std::forward(pre_args)..., + std::forward(t))>(std::get(t))...); } } // namespace internal @@ -36,15 +42,19 @@ constexpr decltype(auto) apply_impl(F&& f, Tuple&& t, * * @tparam F Type of functor * @tparam Tuple Type of tuple containing arguments + * @tparam PreArgs Parameter pack of arguments before the tuple * @param f functor callable * @param t tuple of arguments + * @param pre_args parameter pack of arguments to place before elements in + * tuple. */ -template -constexpr decltype(auto) apply(F&& f, Tuple&& t) { +template +constexpr decltype(auto) apply(F&& f, Tuple&& t, PreArgs&&... pre_args) { return internal::apply_impl( std::forward(f), std::forward(t), std::make_index_sequence< - std::tuple_size>{}>{}); + std::tuple_size>{}>{}, + std::forward(pre_args)...); } } // namespace math diff --git a/stan/math/prim/functor/hcubature.hpp b/stan/math/prim/functor/hcubature.hpp index b342dc42e56..48b77c2a0fb 100644 --- a/stan/math/prim/functor/hcubature.hpp +++ b/stan/math/prim/functor/hcubature.hpp @@ -23,6 +23,8 @@ #include #include #include +#include +#include #include #include @@ -31,28 +33,29 @@ namespace math { namespace internal { -static constexpr double xd7[8] = {-9.9145537112081263920685469752598e-01, - -9.4910791234275852452618968404809e-01, - -8.6486442335976907278971278864098e-01, - -7.415311855993944398638647732811e-01, - -5.8608723546769113029414483825842e-01, - -4.0584515137739716690660641207707e-01, - -2.0778495500789846760068940377309e-01, - 0.0}; - -static constexpr double wd7[8] = {2.2935322010529224963732008059913e-02, - 6.3092092629978553290700663189093e-02, - 1.0479001032225018383987632254189e-01, - 1.4065325971552591874518959051021e-01, - 1.6900472663926790282658342659795e-01, - 1.9035057806478540991325640242055e-01, - 2.0443294007529889241416199923466e-01, - 2.0948214108472782801299917489173e-01}; - -static constexpr double gwd7[4] = {1.2948496616886969327061143267787e-01, - 2.797053914892766679014677714229e-01, - 3.8183005050511894495036977548818e-01, - 4.1795918367346938775510204081658e-01}; +static constexpr std::array xd7{ + -9.9145537112081263920685469752598e-01, + -9.4910791234275852452618968404809e-01, + -8.6486442335976907278971278864098e-01, + -7.415311855993944398638647732811e-01, + -5.8608723546769113029414483825842e-01, + -4.0584515137739716690660641207707e-01, + -2.0778495500789846760068940377309e-01}; + +static constexpr std::array wd7{ + 2.2935322010529224963732008059913e-02, + 6.3092092629978553290700663189093e-02, + 1.0479001032225018383987632254189e-01, + 1.4065325971552591874518959051021e-01, + 1.6900472663926790282658342659795e-01, + 1.9035057806478540991325640242055e-01, + 2.0443294007529889241416199923466e-01, + 2.0948214108472782801299917489173e-01}; + +static constexpr std::array gwd7{ + 1.2948496616886969327061143267787e-01, 2.797053914892766679014677714229e-01, + 3.8183005050511894495036977548818e-01, + 4.1795918367346938775510204081658e-01}; /** * Get the [x]-th lexicographically ordered set of [p] elements in [dim] @@ -62,119 +65,139 @@ static constexpr double gwd7[4] = {1.2948496616886969327061143267787e-01, * Vol. 3, No. 2, June 1977. User lucaroni from * https://stackoverflow.com/questions/561/how-to-use-combinations-of-sets-as-test-data#794 * - * @param c output vector + * @param[out] c output vector * @param dim dimension * @param p number of elements * @param x x-th lexicographically ordered set */ -inline void combination(std::vector& c, const int& dim, const int& p, - const int& x) { - size_t r, k = 0; - for (std::size_t i = 0; i < p - 1; i++) { - c[i] = (i != 0) ? c[i - 1] : 0; - do { +inline void combination(Eigen::Matrix& c, const int dim, + const int p, const int x) { + int r = 0; + int k = 0; + c[0] = 0; + for (; k < x; r = choose(dim - c[0], p - 1), k = k + r) { + c[0]++; + } + k = k - r; + for (int i = 1; i < p - 1; i++) { + c[i] = c[i - 1]; + for (; k < x; r = choose(dim - c[i], p - (i + 1)), k = k + r) { c[i]++; - r = choose(dim - c[i], p - (i + 1)); - k = k + r; - } while (k < x); + } k = k - r; } - if (p > 1) { - c[p - 1] = c[p - 2] + x - k; - } else { - c[0] = x; - } + c[p - 1] = (p > 1) ? c[p - 2] + x - k : x; } /** - * Compute a vector [p] of all [dim]-component vectors + * Compute a matrix [p] of all [dim]-component vectors * with [k] components equal to [lambda] and other components equal to zero. * + * @param[in,out] p matrix * @param k number of components equal to lambda * @param lambda scalar * @param dim dimension - * @param p vector of vectors */ -inline void combos(const int& k, const double& lambda, const int& dim, - std::vector>& p) { - std::vector c(k); +template +inline Eigen::Matrix combos( + const int k, const Scalar lambda, const int dim) { + Eigen::Matrix c(k); const auto choose_dimk = choose(dim, k); - for (std::size_t i = 1; i != choose_dimk + 1; i++) { - std::vector temp(dim, 0.0); - combination(c, dim, k, i); - for (std::size_t j = 0; j != k; j++) { - temp[c[j] - 1] = lambda; + Eigen::Matrix p + = Eigen::MatrixXd::Zero(dim, choose_dimk); + for (size_t i = 0; i < choose_dimk; i++) { + combination(c, dim, k, i + 1); + for (size_t j = 0; j < k; j++) { + p.coeffRef(c.coeff(j) - 1, i) = lambda; } - p.push_back(temp); } + return p; } /** - * Helper function for signcombos + * Compute a matrix [p] of all [dim]-component vectors + * with [k] components equal to [±lambda] and other components equal to zero + * (with all possible signs). * - * @param index helper vector + * @param[in,out] p matrix * @param k number of components equal to lambda * @param lambda scalar - * @param c ordered vector - * @param temp helper vector + * @param dim dimension */ -inline void increment(std::vector& index, const int& k, - const double& lambda, const std::vector& c, - std::vector& temp) { - if (index.size() == 0) { +template +inline Eigen::Matrix signcombos( + const int k, const Scalar lambda, const int dim) { + Eigen::Matrix c(k); + const auto choose_dimk = choose(dim, k); + Eigen::Matrix p + = Eigen::MatrixXd::Zero(dim, choose_dimk * std::pow(2, k)); + int current_col = 0; + const auto inner_iter_len = std::pow(2, k); + std::vector index; + index.reserve(inner_iter_len * (choose_dimk + 1)); + for (int i = 1; i != choose_dimk + 1; i++) { + combination(c, dim, k, i); index.push_back(false); - for (std::size_t j = 0; j != k; j++) { - temp[c[j] - 1] = lambda; - } - return; - } - int first_zero = 0; - while ((first_zero < index.size()) && index[first_zero]) { - first_zero++; - } - if (first_zero == index.size()) { - index.flip(); - for (std::size_t j = 0; j != index.size(); j++) { - temp[c[j] - 1] *= -1; - } - index.push_back(true); - temp[c[index.size() - 1] - 1] = -lambda; - } else { - for (std::size_t i = 0; i != first_zero + 1; i++) { - if (index[i]) { - index[i] = 0; + p(c.array() - 1.0, current_col).setConstant(lambda); + current_col += 1; + for (int j = 1; j != inner_iter_len; j++, current_col++) { + p.col(current_col) = p.col(current_col - 1); + int first_zero + = std::distance(std::begin(index), + std::find(std::begin(index), std::end(index), false)); + const std::size_t index_size = index.size(); + if (first_zero == index_size) { + index.flip(); + p(c.segment(0, index.size()).array() - 1, current_col).array() *= -1; + index.push_back(true); + p(c[index.size() - 1] - 1, current_col) = -lambda; } else { - index[i] = 1; + for (int h = 0; h != first_zero + 1; h++) { + index[h].flip(); + } + p(c.segment(0, first_zero + 1).array() - 1, current_col).array() *= -1; } - temp[c[i] - 1] *= -1; } + index.clear(); } + return p; } /** - * Compute a vector [p] of all [dim]-component vectors - * with [k] components equal to [±lambda] and other components equal to zero - * (with all possible signs). + * Compute the points and weights corresponding to a [dim]-dimensional + * Genz-Malik cubature rule * - * @param k number of components equal to lambda - * @param lambda scalar + * @param[in,out] points points for the last 4 GenzMalik weights + * @param[in,out] weights weights for the 5 terms in the GenzMalik rule + * @param[in,out] weights_low_deg weights for the embedded lower-degree rule * @param dim dimension - * @param p vector of vectors */ -inline void signcombos(const int& k, const double& lambda, const int& dim, - std::vector>& p) { - std::vector c(k); - const auto choose_dimk = choose(dim, k); - for (std::size_t i = 1; i != choose_dimk + 1; i++) { - std::vector temp(dim, 0.0); - combination(c, dim, k, i); - std::vector index; - index.clear(); - for (std::size_t j = 0; j != std::pow(2, k); j++) { - increment(index, k, lambda, c, temp); - p.push_back(temp); - } - } + +inline std::tuple< + std::array, 4>, + Eigen::Matrix, Eigen::Matrix> +make_GenzMalik(const int dim) { + std::array, 4> points; + Eigen::Matrix weights; + Eigen::Matrix weights_low_deg; + double twopn = std::pow(2, dim); + weights[0] + = twopn * ((12824.0 - 9120.0 * dim + 400.0 * dim * dim) * 1.0 / 19683.0); + weights[1] = twopn * (980.0 / 6561.0); + weights[2] = twopn * ((1820.0 - 400.0 * dim) * 1.0 / 19683.0); + weights[3] = twopn * (200.0 / 19683.0); + weights[4] = 6859.0 / 19683.0; + weights_low_deg[0] = twopn * ((729 - 950 * dim + 50 * dim * dim) * 1.0 / 729); + weights_low_deg[1] = twopn * (245.0 / 486); + weights_low_deg[2] = twopn * ((265.0 - 100.0 * dim) * 1.0 / 1458.0); + weights_low_deg[3] = twopn * (25.0 / 729.0); + points[0] = combos(1, std::sqrt(9.0 * 1.0 / 70.0), dim); + double l3 = std::sqrt(9.0 * 1.0 / 10.0); + points[1] = combos(1, l3, dim); + points[2] = signcombos(2, l3, dim); + points[3] = signcombos(dim, std::sqrt(9.0 * 1.0 / 19.0), dim); + return std::make_tuple(std::move(points), std::move(weights), + std::move(weights_low_deg)); } /** @@ -182,75 +205,46 @@ inline void signcombos(const int& k, const double& lambda, const int& dim, * for one dimension. * * @tparam F type of the integrand - * @tparam T_pars type of the parameters for the integrand + * @tparam T_a type of lower limit of integration + * @tparam T_b type of upper limit of integration + * @tparam ParsPairT type of the pair of parameters for the integrand * @param integrand function to be integrated * @param a lower limit of integration * @param b upper limit of integration - * @param pars parameters for the integrand + * @param pars_pair Pair of parameters for the integrand * @return numeric integral of the integrand and error */ -template -std::tuple gauss_kronrod(const F& integrand, const double& a, - const double& b, T_pars& pars) { - std::vector c(1, 0); - std::vector cp(1, 0); - std::vector cm(1, 0); - c[0] = 0.5 * (a + b); - double delta = 0.5 * (b - a); - double f0 = integrand(c, pars); - double I = f0 * wd7[7]; - double Idash = f0 * gwd7[3]; - for (std::size_t i = 0; i != 7; i++) { - double deltax = delta * xd7[i]; - cp[0] = c[0] + deltax; - cm[0] = c[0] - deltax; - double fx = integrand(cp, pars); - double temp = integrand(cm, pars); - fx += temp; +template +inline auto gauss_kronrod(const F& integrand, const T_a a, const T_b b, + const ParsPairT& pars_pair) { + using delta_t = return_type_t; + const delta_t c = 0.5 * (a + b); + const delta_t delta = 0.5 * (b - a); + auto f0 = math::apply([](auto&& integrand, auto&& c, + auto&&... args) { return integrand(c, args...); }, + pars_pair, integrand, c); + + auto I = f0 * wd7[7]; + auto Idash = f0 * gwd7[3]; + std::array deltax; + for (int i = 0; i < 7; ++i) { + deltax[i] = delta * xd7[i]; + } + for (auto i = 0; i != 7; i++) { + auto fx = math::apply( + [](auto&& integrand, auto&& c, auto&& delta, auto&&... args) { + return integrand(c + delta, args...) + integrand(c - delta, args...); + }, + pars_pair, integrand, c, deltax[i]); I += fx * wd7[i]; if (i % 2 == 1) { Idash += fx * gwd7[i / 2]; } } - double V = fabs(delta); + delta_t V = fabs(delta); I *= V; Idash *= V; - return std::make_tuple(I, fabs(I - Idash)); -} - -/** - * Compute the points and weights corresponding to a [dim]-dimensional - * Genz-Malik cubature rule - * - * @param dim dimension - * @param p points for the last 4 GenzMalik weights - * @param w weights for the 5 terms in the GenzMalik rule - * @param wd weights for the embedded lower-degree rule - */ -inline void make_GenzMalik(const int& dim, - std::vector>>& p, - std::vector& w, std::vector& wd) { - double l4 = std::sqrt(9 * 1.0 / 10); - double l2 = std::sqrt(9 * 1.0 / 70); - double l3 = l4; - double l5 = std::sqrt(9 * 1.0 / 19); - - double twopn = std::pow(2, dim); - - w[0] = twopn * ((12824 - 9120 * dim + 400 * dim * dim) * 1.0 / 19683); - w[1] = twopn * (980.0 / 6561); - w[2] = twopn * ((1820 - 400 * dim) * 1.0 / 19683); - w[3] = twopn * (200.0 / 19683); - w[4] = 6859.0 / 19683; - wd[3] = twopn * (25.0 / 729); - wd[2] = twopn * ((265 - 100 * dim) * 1.0 / 1458); - wd[1] = twopn * (245.0 / 486); - wd[0] = twopn * ((729 - 950 * dim + 50 * dim * dim) * 1.0 / 729); - - combos(1, l2, dim, p[0]); - combos(1, l3, dim, p[1]); - signcombos(2, l4, dim, p[2]); - signcombos(dim, l5, dim, p[3]); + return std::make_pair(I, fabs(I - Idash)); } /** @@ -258,118 +252,91 @@ inline void make_GenzMalik(const int& dim, * for more than one dimensions. * * @tparam F type of the integrand - * @tparam T_pars type of the parameters for the integrand - * @param integrand function to be integrated - * @param p - * @param w - * @param wd + * @tparam GenzMalik type of ----- + * @tparam T_a type of lower limit of integration + * @tparam T_b type of upper limit of integration + * @tparam ParsTupleT type of the tuple of parameters for the integrand + * @param[out] integrand function to be integrated + * @param[in] genz_malik tuple of GenzMalik weights * @param dim dimension of the multidimensional integral * @param a lower limit of integration * @param b upper limit of integration - * @param pars parameters for the integrand + * @param pars_tuple Tuple of parameters for the integrand * @return numeric integral of the integrand, error, and suggested coordinate to * subdivide next */ -template -std::tuple integrate_GenzMalik( - const F& integrand, std::vector>>& p, - std::vector& w, std::vector& wd, const int& dim, - const std::vector& a, const std::vector& b, T_pars& pars) { - std::vector c(dim, 0); - std::vector deltac(dim); - - for (std::size_t i = 0; i != dim; i++) { +template +inline auto integrate_GenzMalik(const F& integrand, const GenzMalik& genz_malik, + const int dim, + const Eigen::Matrix& a, + const Eigen::Matrix& b, + const ParsTupleT& pars_tuple) { + auto&& points = std::get<0>(genz_malik); + auto&& weights = std::get<1>(genz_malik); + auto&& weights_low_deg = std::get<2>(genz_malik); + using delta_t = return_type_t; + Eigen::Matrix c(dim); + Eigen::Matrix deltac(dim); + using Scalar = return_type_t; + for (auto i = 0; i != dim; i++) { + if (a[i] == b[i]) { + return std::make_tuple(Scalar(0.0), Scalar(0.0), 0); + } c[i] = (a[i] + b[i]) / 2; + deltac[i] = (b[i] - a[i]) / 2.0; } - - for (std::size_t i = 0; i != dim; i++) { - deltac[i] = fabs(b[i] - a[i]) / 2; - } - double v = 1.0; + delta_t v = 1.0; for (std::size_t i = 0; i != dim; i++) { v *= deltac[i]; } - - if (v == 0.0) { - return std::make_tuple(0.0, 0.0, 0); + Eigen::Matrix f = Eigen::Matrix::Zero(); + f.coeffRef(0) + = math::apply([](auto&& integrand, auto&& c, + auto&&... args) { return integrand(c, args...); }, + pars_tuple, integrand, c); + Eigen::Matrix divdiff(dim); + for (auto i = 0; i != dim; i++) { + auto p2 = deltac.cwiseProduct(points[0].col(i)); + auto f2i = math::apply( + [](auto&& integrand, auto&& c, auto&& p2, auto&&... args) { + return integrand(c + p2, args...) + integrand(c - p2, args...); + }, + pars_tuple, integrand, c, p2); + auto p3 = deltac.cwiseProduct(points[1].col(i)); + auto f3i = math::apply( + [](auto&& integrand, auto&& c, auto&& p3, auto&&... args) { + return integrand(c + p3, args...) + integrand(c - p3, args...); + }, + pars_tuple, integrand, c, p3); + f.coeffRef(1) += f2i; + f.coeffRef(2) += f3i; + divdiff[i] = fabs(f3i + 12.0 * f.coeff(0) - 7.0 * f2i); } - - double f1 = integrand(c, pars); - double f2 = 0.0; - double f3 = 0.0; - 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); - - for (std::size_t i = 0; i != dim; i++) { - for (std::size_t j = 0; j != dim; j++) { - p2[j] = deltac[j] * p[0][i][j]; - } - - for (std::size_t j = 0; j != dim; j++) { - cc[j] = c[j] + p2[j]; - } - double f2i = integrand(cc, pars); - for (std::size_t j = 0; j != dim; j++) { - cc[j] = c[j] - p2[j]; - } - double temp = integrand(cc, pars); - f2i += temp; - - for (std::size_t j = 0; j != dim; j++) { - p3[j] = deltac[j] * p[1][i][j]; - } - for (std::size_t j = 0; j != dim; j++) { - cc[j] = c[j] + p3[j]; - } - double f3i = integrand(cc, pars); - for (std::size_t j = 0; j != dim; j++) { - cc[j] = c[j] - p3[j]; - } - temp = integrand(cc, pars); - f3i += temp; - f2 += f2i; - f3 += f3i; - divdiff[i] = fabs(f3i + twelvef1 - 7 * f2i); + for (auto i = 0; i != points[2].cols(); i++) { + f.coeffRef(3) += math::apply( + [](auto&& integrand, auto&& cc, auto&&... args) { + return integrand(cc, args...); + }, + pars_tuple, integrand, c + deltac.cwiseProduct(points[2].col(i))); } - std::vector p4(dim); - double f4 = 0.0; - for (std::size_t i = 0; i != p[2].size(); i++) { - for (std::size_t j = 0; j != dim; j++) { - p4[j] = deltac[j] * p[2][i][j]; - } - for (std::size_t j = 0; j != dim; j++) { - cc[j] = c[j] + p4[j]; - } - double temp = integrand(cc, pars); - f4 += temp; - } - double f5 = 0.0; - std::vector p5(dim); - for (std::size_t i = 0; i != p[3].size(); i++) { - for (std::size_t j = 0; j != dim; j++) { - p5[j] = deltac[j] * p[3][i][j]; - } - - for (std::size_t j = 0; j != dim; j++) { - cc[j] = c[j] + p5[j]; - } - double temp = integrand(cc, pars); - f5 += temp; + for (auto i = 0; i != points[3].cols(); i++) { + f.coeffRef(4) += math::apply( + [](auto&& integrand, auto&& cc, auto&&... args) { + return integrand(cc, args...); + }, + pars_tuple, integrand, c + deltac.cwiseProduct(points[3].col(i))); } - double I = v * (w[0] * f1 + w[1] * f2 + w[2] * f3 + w[3] * f4 + w[4] * f5); - double Idash = v * (wd[0] * f1 + wd[1] * f2 + wd[2] * f3 + wd[3] * f4); - double E = fabs(I - Idash); + Scalar I = v * weights.dot(f); + Scalar Idash = v * weights_low_deg.dot(f.template head<4>()); + Scalar E = fabs(I - Idash); int kdivide = 0; - double deltaf = E / (std::pow(10, dim) * v); - for (std::size_t i = 0; i != dim; i++) { - double delta = divdiff[i] - maxdivdiff; + Scalar deltaf = E / (std::pow(10, dim) * v); + Scalar maxdivdiff = 0.0; + for (auto i = 0; i != dim; i++) { + Scalar delta = divdiff[i] - maxdivdiff; if (delta > deltaf) { kdivide = i; maxdivdiff = divdiff[i]; @@ -380,148 +347,168 @@ std::tuple integrate_GenzMalik( return std::make_tuple(I, E, kdivide); } +/** + * Compute the integral of the function to be integrated (integrand) from a to b + * for more than one dimensions. + * + * @tparam T_a Type of return_type_t 1 + * @tparam T_b Type of return_type_t 2 + * @param a lower bounds of the integral + * @param b upper bounds of the integral + * @param I value of the integral + * @param kdivide number of subdividing the integration volume + */ +template struct Box { - Box(const std::vector& a, const std::vector& b, double I, - double err, int kdivide) - : a(a), b(b), I(I), E(err), kdiv(kdivide) {} - bool operator<(const Box& box) const { return E < box.E; } - std::vector a; - std::vector b; - double I; - double E; - int kdiv; + template + Box(Vec1&& a, Vec2&& b, return_type_t I, int kdivide) + : a_(std::forward(a)), + b_(std::forward(b)), + I_(I), + kdiv_(kdivide) {} + Eigen::Matrix a_; + Eigen::Matrix b_; + return_type_t I_; + int kdiv_; }; } // namespace internal /** - * Compute the dim-dimensional integral of the function \f$f\f$ from \f$a\f$ to - \f$b\f$ within + * Compute the [dim]-dimensional integral of the function \f$f\f$ from \f$a\f$ + to \f$b\f$ within * specified relative and absolute tolerances or maximum number of evaluations. * \f$a\f$ and \f$b\f$ can be finite or infinite and should be given as vectors. * - * The prototype for \f$f\f$ is: - \verbatim - template - stan::return_type_t f(const T_x& x, const T_p& p) { - using T_x_ref = stan::ref_type_t; - T_x_ref x_ref = x; - stan::scalar_seq_view x_vec(x_ref); - my_params* pars = static_cast(p); - type_1 var_1 = (pars->par_1); - return ; - } - \endverbatim - * - * The parameters can be handed over to f as a struct: - \verbatim - struct my_params { - type_1 par_1; - type_2 par_2; - }; - \endverbatim - * * @tparam F Type of f - * @tparam T_pars Type of paramete-struct + * @tparam T_a Type of lower limit of integration + * @tparam T_b Type of upper limit of integration + * @tparam ParsTuple Type of parameter-tuple + * @tparam TAbsErr Type of absolute error + * @tparam TRelErr Type of relative error * * @param integrand a functor with signature given above - * @param pars parameters to be passed to f as a struct + * @param pars parameters to be passed to f as a tuple * @param dim dimension of the integral * @param a lower limit of integration as vector * @param b upper limit of integration as vector - * @param maxEval maximal number of evaluations + * @param max_eval maximal number of evaluations * @param reqAbsError absolute error * @param reqRelError relative error as vector - * @param val correct value of integral * - * @return The value of the dim-dimensional integral of \f$f\f$ from \f$a\f$ to - \f$b\f$. + * @return The value of the [dim]-dimensional integral of \f$f\f$ from \f$a\f$ + to \f$b\f$. * @throw std::domain_error no errors will be thrown. */ -template -double hcubature(const F& integrand, const T_pars& pars, const int& dim, - const std::vector& a, const std::vector& b, - int& maxEval, const double& reqAbsError, - const double& reqRelError) { - if (maxEval <= 0) { - maxEval = 1000000; - } - - double result, err; - int kdivide = 0; - - std::vector>> p(4); - std::vector w_five(5); - std::vector wd_four(4); +template +inline auto hcubature(const F& integrand, const ParsTuple& pars, const int dim, + const Eigen::Matrix& a, + const Eigen::Matrix& b, + const int max_eval, const TAbsErr reqAbsError, + const TRelErr reqRelError) { + using Scalar = return_type_t; + using eig_vec_a = Eigen::Matrix; + using eig_vec_b = Eigen::Matrix; + using namespace boost::math::quadrature; + + const int maxEval = max_eval <= 0 ? 1000000 : max_eval; + Scalar result; + Scalar err; + auto kdivide = 0; + std::tuple< + std::array, 4>, + Eigen::Matrix, Eigen::Matrix> + genz_malik; + + auto gk_lambda = [&integrand, &pars](auto&& c) { + return stan::math::apply( + [](auto&& integrand, auto&& c, auto&&... args) { + return integrand(c, args...); + }, + pars, integrand, c); + }; if (dim == 1) { std::tie(result, err) = internal::gauss_kronrod(integrand, a[0], b[0], pars); } else { - internal::make_GenzMalik(dim, p, w_five, wd_four); - std::tie(result, err, kdivide) = internal::integrate_GenzMalik( - integrand, p, w_five, wd_four, dim, a, b, pars); + genz_malik = internal::make_GenzMalik(dim); + std::tie(result, err, kdivide) + = internal::integrate_GenzMalik(integrand, genz_malik, dim, a, b, pars); } - int numevals + auto numevals = (dim == 1) ? 15 : 1 + 4 * dim + 2 * dim * (dim - 1) + std::pow(2, dim); - int evals_per_box = numevals; - double error = err; - double val = result; + + auto evals_per_box = numevals; + Scalar error = err; + Scalar val = result; if ((error <= fmax(reqRelError * fabs(val), reqAbsError)) || (numevals >= maxEval)) { - return val; + return result; } - std::priority_queue ms; - internal::Box box(a, b, result, err, kdivide); - ms.push(box); - numevals += 2 * evals_per_box; + using box_t = internal::Box; + std::vector ms; + ms.reserve(numevals); + ms.emplace_back(a, b, result, kdivide); + auto get_largest_box_idx = [](auto&& box_vec) { + auto max_it = std::max_element(box_vec.begin(), box_vec.end()); + return std::distance(box_vec.begin(), max_it); + }; + std::vector err_vec; + err_vec.reserve(numevals); + err_vec.push_back(err); while ((numevals < maxEval) - && (error > max(reqRelError * fabs(val), reqAbsError)) - && std::isfinite(val)) { - internal::Box box = ms.top(); - ms.pop(); - - double w = (box.b[box.kdiv] - box.a[box.kdiv]) / 2; - std::vector ma(box.a); - - ma[box.kdiv] += w; - std::vector mb(box.b); - mb[box.kdiv] -= w; - - double result_1, result_2, err_1, err_2; - double kdivide_1 = math::NOT_A_NUMBER; - double kdivide_2 = math::NOT_A_NUMBER; - + && (error > fmax(reqRelError * fabs(val), reqAbsError)) + && fabs(val) < INFTY) { + auto err_idx = get_largest_box_idx(err_vec); + auto&& box = ms[err_idx]; + auto w = (box.b_[box.kdiv_] - box.a_[box.kdiv_]) / 2; + eig_vec_a ma = Eigen::Map(box.a_.data(), box.a_.size()); + ma[box.kdiv_] += w; + eig_vec_b mb = Eigen::Map(box.b_.data(), box.b_.size()); + mb[box.kdiv_] -= w; + int kdivide_1{0}; + int kdivide_2{0}; + Scalar result_1; + Scalar result_2; + Scalar err_1; + Scalar err_2; if (dim == 1) { std::tie(result_1, err_1) - = internal::gauss_kronrod(integrand, ma[0], box.b[0], pars); + = internal::gauss_kronrod(integrand, ma[0], box.b_[0], pars); std::tie(result_2, err_2) - = internal::gauss_kronrod(integrand, box.a[0], mb[0], pars); + = internal::gauss_kronrod(integrand, box.a_[0], mb[0], pars); } else { std::tie(result_1, err_1, kdivide_1) = internal::integrate_GenzMalik( - integrand, p, w_five, wd_four, dim, ma, box.b, pars); + integrand, genz_malik, dim, ma, box.b_, pars); std::tie(result_2, err_2, kdivide_2) = internal::integrate_GenzMalik( - integrand, p, w_five, wd_four, dim, box.a, mb, pars); + integrand, genz_malik, dim, box.a_, mb, pars); } - internal::Box box1(ma, box.b, result_1, err_1, kdivide_1); - ms.push(box1); - internal::Box box2(box.a, mb, result_2, err_2, kdivide_2); - ms.push(box2); - val += box1.I + box2.I - box.I; - error += box1.E + box2.E - box.E; + box_t box1(std::move(ma), std::move(box.b_), result_1, kdivide_1); + box_t box2(std::move(box.a_), std::move(mb), result_2, kdivide_2); + result += result_1 + result_2 - box.I_; + err += err_1 + err_2 - err_vec[err_idx]; + ms[err_idx].I_ = 0; + err_vec[err_idx] = 0; + ms.push_back(std::move(box1)); + ms.push_back(std::move(box2)); + err_vec.push_back(err_1); + err_vec.push_back(err_2); numevals += 2 * evals_per_box; } - val = 0.0; - error = 0.0; + result = 0.0; + err = 0.0; - for (; !ms.empty(); ms.pop()) { - internal::Box box = ms.top(); - val += box.I; - error += box.E; + for (auto&& box : ms) { + result += box.I_; + } + for (auto err_i : err_vec) { + err += err_i; } - return val; + return result; } // hcubature } // namespace math diff --git a/stan/math/prim/meta/return_type.hpp b/stan/math/prim/meta/return_type.hpp index 79ef5816a9a..e4204343011 100644 --- a/stan/math/prim/meta/return_type.hpp +++ b/stan/math/prim/meta/return_type.hpp @@ -146,6 +146,26 @@ struct scalar_lub, std::complex> { template using scalar_lub_t = typename scalar_lub::type; +namespace internal { +template +struct return_type_impl { + using type = double; +}; + +template +struct return_type_impl { + using type + = scalar_lub_t, typename return_type_impl::type>; +}; + +template +struct return_type_impl, Ts...> { + using type = scalar_lub_t::type, + typename return_type_impl::type>; +}; + +} // namespace internal + /** * Template metaprogram to calculate the base scalar return type * resulting from promoting all the scalar types of the template @@ -184,15 +204,7 @@ using scalar_lub_t = typename scalar_lub::type; * @ingroup type_trait */ template -struct return_type { - using type = double; -}; - -template -struct return_type { - using type - = scalar_lub_t, typename return_type::type>; -}; +struct return_type : internal::return_type_impl...> {}; /** * Convenience type for the return type of the specified template diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 70aafe82422..598663b71e2 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -369,8 +369,10 @@ #include #include #include -#include +#include #include +#include +#include #include #include #include diff --git a/stan/math/prim/prob/wiener5_lpdf.hpp b/stan/math/prim/prob/wiener5_lpdf.hpp new file mode 100644 index 00000000000..075b54906b6 --- /dev/null +++ b/stan/math/prim/prob/wiener5_lpdf.hpp @@ -0,0 +1,846 @@ +#ifndef STAN_MATH_PRIM_PROB_WIENER5_LPDF_HPP +#define STAN_MATH_PRIM_PROB_WIENER5_LPDF_HPP + +#include + +namespace stan { +namespace math { +namespace internal { + +enum GradientCalc { OFF = 0, ON = 1 }; + +/** + * Calculate the 'error_term' term for a wiener5 density or gradient + * + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_v type of drift rate + * @tparam T_w type of relative starting point + * @tparam T_sv type of inter-trial variability in v + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param v_value The drift rate + * @param w_value The relative starting point + * @param sv The inter-trial variability of the drift rate + * @return 'error_term' term + */ +template +inline auto wiener5_compute_error_term(T_y&& y, T_a&& a, T_v&& v_value, + T_w&& w_value, T_sv&& sv) noexcept { + const auto w = 1.0 - w_value; + const auto v = -v_value; + const auto sv_sqr = square(sv); + const auto one_plus_svsqr_y = 1 + sv_sqr * y; + const auto two_avw = 2.0 * a * v * w; + const auto two_log_a = 2.0 * log(a); + return stan::math::eval((sv_sqr * square(a * w) - two_avw - square(v) * y) + / 2.0 / one_plus_svsqr_y + - two_log_a - 0.5 * log(one_plus_svsqr_y)); +} + +/** + * Calculate the 'n_terms_small_t' term for a wiener5 density or gradient + * + * @tparam Density Whether the calculation is for the density + * @tparam GradW Whether the calculation is for gradient w.r.t. 'w' + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_w_value type of relative starting point + * @tparam T_err type of error + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param w_value The relative starting point + * @param error The error tolerance + * @return 'n_terms_small_t' term + */ +template +inline auto wiener5_n_terms_small_t(T_y&& y, T_a&& a, T_w_value&& w_value, + T_err&& error) noexcept { + const auto two_error = 2.0 * error; + const auto y_asq = y / square(a); + const auto two_log_a = 2 * log(a); + const auto log_y_asq = log(y) - two_log_a; + const auto w = 1.0 - w_value; + + const auto n_1_factor = Density ? 2 : 3; + const auto n_1 = (sqrt(n_1_factor * y_asq) + w) / 2.0; + auto u_eps = (Density || GradW) + ? fmin(-1.0, LOG_TWO + LOG_PI + 2.0 * log_y_asq + two_error) + : fmin(-3.0, (log(8.0) - log(27.0) + LOG_PI + 4.0 * log_y_asq + + two_error)); + const auto arg_mult = (Density || GradW) ? 1 : 3; + const auto arg = -arg_mult * y_asq * (u_eps - sqrt(-2.0 * u_eps - 2.0)); + + const auto n_2 + = (arg > 0) ? GradW ? 0.5 * (sqrt(arg) + w) : 0.5 * (sqrt(arg) - w) : n_1; + + return ceil(fmax(n_1, n_2)); +} + +/** + * Calculate the 'n_terms_large_t' term for a wiener5 density + * + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_w_value type of relative starting point + * @tparam T_err type of error + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param w_value The relative starting point + * @param error The error tolerance + * @return 'n_terms_large_t' term + */ +template +inline auto wiener5_density_large_reaction_time_terms(T_y&& y, T_a&& a, + T_w_value&& w_value, + T_err&& error) noexcept { + const auto y_asq = y / square(a); + const auto log_y_asq = log(y) - 2 * log(a); + static constexpr double PI_SQUARED = pi() * pi(); + auto n_1 = 1.0 / (pi() * sqrt(y_asq)); + const auto two_log_piy = -2.0 * (LOG_PI + log_y_asq + error); + auto n_2 + = (two_log_piy >= 0) ? sqrt(two_log_piy / (PI_SQUARED * y_asq)) : 0.0; + return ceil(fmax(n_1, n_2)); +} + +/** + * Calculate the 'n_terms_large_t' term for a wiener5 gradient + * + * @tparam GradW Whether the calculation is for gradient w.r.t. 'w' + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_w_value type of relative starting point + * @tparam T_err type of error + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param w_value The relative starting point + * @param error The error tolerance + * @return 'n_terms_large_t' term + */ +template +inline auto wiener5_gradient_large_reaction_time_terms(T_y&& y, T_a&& a, + T_w_value&& w_value, + T_err&& error) noexcept { + const auto y_asq = y / square(a); + const auto log_y_asq = log(y) - 2 * log(a); + static constexpr double PI_SQUARED = pi() * pi(); + const auto n_1_factor = GradW ? 2 : 3; + auto n_1 = sqrt(n_1_factor / y_asq) / pi(); + const auto two_error = 2.0 * error; + const auto u_eps_arg + = GradW ? log(4.0) - log(9.0) + 2.0 * LOG_PI + 3.0 * log_y_asq + two_error + : log(3.0) - log(5.0) + LOG_PI + 2.0 * log_y_asq + error; + const auto u_eps = fmin(-1, u_eps_arg); + const auto arg_mult = GradW ? 1 : (2.0 / PI_SQUARED / y_asq); + const auto arg = -arg_mult * (u_eps - sqrt(-2.0 * u_eps - 2.0)); + auto n_2 = GradW ? ((arg > 0) ? sqrt(arg / y_asq) / pi() : n_1) + : ((arg > 0) ? sqrt(arg) : n_1); + return ceil(fmax(n_1, n_2)); +} + +/** + * Calculate the 'result' term and its sign for a wiener5 density or gradient + * + * @tparam Density Whether the calculation is for the density + * @tparam GradW Whether the calculation is for gradient w.r.t. 'w' + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_w type of relative starting point + * @tparam T_nsmall type of term number_small_terms + * @tparam T_nlarge type of term number_large_terms + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param w_value The relative starting point + * @param n_terms_small_t The n_terms_small_t term + * @param n_terms_large_t The n_terms_large_t term + * @return 'result' sum and its sign + */ +template +inline auto wiener5_log_sum_exp(T_y&& y, T_a&& a, T_w&& w_value, + T_nsmall&& n_terms_small_t, + T_nlarge&& n_terms_large_t) noexcept { + const auto y_asq = y / square(a); + const auto w = 1.0 - w_value; + const bool small_n_terms_small_t + = Density ? (2 * n_terms_small_t <= n_terms_large_t) + : (2 * n_terms_small_t < n_terms_large_t); + const auto scaling = small_n_terms_small_t ? inv(2.0 * y_asq) : y_asq / 2.0; + using ret_t = return_type_t; + ret_t fplus = NEGATIVE_INFTY; + ret_t fminus = NEGATIVE_INFTY; + int current_sign; + if (small_n_terms_small_t) { + constexpr double mult = Density ? 1.0 : 3.0; + if (GradW) { + for (auto k = n_terms_small_t; k >= 1; k--) { + const auto w_plus_2_k = w + 2.0 * k; + const auto w_minus_2_k = w - 2.0 * k; + const auto square_w_plus_2_k_minus_offset = square(w_plus_2_k) - y_asq; + if (square_w_plus_2_k_minus_offset > 0) { + const auto summand_plus = log(square_w_plus_2_k_minus_offset) + - square(w_plus_2_k) * scaling; + fplus = log_sum_exp(fplus, summand_plus); + } else if (square_w_plus_2_k_minus_offset < 0) { + const auto summand_plus = log(-square_w_plus_2_k_minus_offset) + - square(w_plus_2_k) * scaling; + fminus = log_sum_exp(fminus, summand_plus); + } + const auto square_w_minus_2_k_minus_offset + = square(w_minus_2_k) - y_asq; + if (square_w_minus_2_k_minus_offset > 0) { + const auto summand_minus = log(square_w_minus_2_k_minus_offset) + - square(w_minus_2_k) * scaling; + fplus = log_sum_exp(fplus, summand_minus); + } else if (square_w_minus_2_k_minus_offset < 0) { + const auto summand_minus = log(-square_w_minus_2_k_minus_offset) + - square(w_minus_2_k) * scaling; + fminus = log_sum_exp(fminus, summand_minus); + } + } + const auto square_w_minus_offset = square(w) - y_asq; + if (square_w_minus_offset > 0) { + const auto new_val = log(square_w_minus_offset) - square(w) * scaling; + fplus = log_sum_exp(fplus, new_val); + } else if (square_w_minus_offset < 0) { + const auto new_val = log(-square_w_minus_offset) - square(w) * scaling; + fminus = log_sum_exp(fminus, new_val); + } + } else { + for (auto k = n_terms_small_t; k >= 1; k--) { + const auto w_plus_2_k = w + 2.0 * k; + const auto w_minus_2_k = w - 2.0 * k; + const auto summand_plus + = mult * log(w_plus_2_k) - square(w_plus_2_k) * scaling; + fplus = log_sum_exp(fplus, summand_plus); + const auto summand_minus + = mult * log(-w_minus_2_k) - square(w_minus_2_k) * scaling; + if (fminus <= NEGATIVE_INFTY) { + fminus = summand_minus; + } else if (summand_minus <= NEGATIVE_INFTY) { + continue; + } else if (fminus > summand_minus) { + fminus = fminus + log1p_exp(summand_minus - fminus); + } else { + fminus = summand_minus + log1p_exp(fminus - summand_minus); + } + } + const auto new_val = mult * log(w) - square(w) * scaling; + fplus = log_sum_exp(fplus, new_val); + } + } else { // for large t + constexpr double mult = (Density ? 1 : (GradW ? 2 : 3)); + for (auto k = n_terms_large_t; k >= 1; k--) { + const auto pi_k = k * pi(); + const auto check = (GradW) ? cos(pi_k * w) : sin(pi_k * w); + if (check > 0) { + fplus = log_sum_exp( + fplus, mult * log(k) - square(pi_k) * scaling + log(check)); + } else if ((GradW && check < 0) || !GradW) { + fminus = log_sum_exp( + fminus, mult * log(k) - square(pi_k) * scaling + log(-check)); + } + } + } + current_sign = (fplus < fminus) ? -1 : 1; + if (fplus == NEGATIVE_INFTY) { + return std::make_pair(fminus, current_sign); + } else if (fminus == NEGATIVE_INFTY) { + return std::make_pair(fplus, current_sign); + } else if (fplus > fminus) { + return std::make_pair(log_diff_exp(fplus, fminus), current_sign); + } else if (fplus < fminus) { + return std::make_pair(log_diff_exp(fminus, fplus), current_sign); + } else { + return std::make_pair(ret_t(NEGATIVE_INFTY), current_sign); + } +} + +/** + * Calculate the wiener5 density + * + * @tparam NaturalScale Whether to return the density on natural (true) or + * log-scale (false) + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_w type of relative starting point + * @tparam T_v type of drift rate + * @tparam T_sv type of inter-trial variability in v + * @tparam T_err type of log error + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param v_value The drift rate + * @param w_value The relative starting point + * @param sv The inter-trial variability of the drift rate + * @param err The log error tolerance + * @return density + */ +template +inline auto wiener5_density(const T_y& y, const T_a& a, const T_v& v_value, + const T_w& w_value, const T_sv& sv, + T_err&& err = log(1e-12)) noexcept { + const auto error_term + = wiener5_compute_error_term(y, a, v_value, w_value, sv); + const auto error = (err - error_term); + const auto n_terms_small_t + = wiener5_n_terms_small_t( + y, a, w_value, error); + const auto n_terms_large_t + = wiener5_density_large_reaction_time_terms(y, a, w_value, error); + + auto res = wiener5_log_sum_exp( + y, a, w_value, n_terms_small_t, n_terms_large_t) + .first; + if (2 * n_terms_small_t <= n_terms_large_t) { + auto log_density = error_term - 0.5 * LOG_TWO - LOG_SQRT_PI + - 1.5 * (log(y) - 2 * log(a)) + res; + return NaturalScale ? exp(log_density) : log_density; + } else { + auto log_density = error_term + res + LOG_PI; + return NaturalScale ? exp(log_density) : log_density; + } +} + +/** + * Calculate the derivative of the wiener5 density w.r.t. 't' + * + * @tparam WrtLog Whether to return the derivative w.r.t. + * the natural (true) or log-scale (false) density + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_v type of drift rate + * @tparam T_w type of relative starting point + * @tparam T_sv type of inter-trial variability in v + * @tparam T_err type of log error + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param v_value The drift rate + * @param w_value The relative starting point + * @param sv The inter-trial variability of the drift rate + * @param err The log error tolerance + * @return Gradient w.r.t. t + */ +template +inline auto wiener5_grad_t(const T_y& y, const T_a& a, const T_v& v_value, + const T_w& w_value, const T_sv& sv, + T_err&& err = log(1e-12)) noexcept { + const auto two_log_a = 2 * log(a); + const auto log_y_asq = log(y) - two_log_a; + const auto error_term + = wiener5_compute_error_term(y, a, v_value, w_value, sv); + const auto w = 1.0 - w_value; + const auto v = -v_value; + const auto sv_sqr = square(sv); + const auto one_plus_svsqr_y = 1 + sv_sqr * y; + const auto density_part_one + = -0.5 + * (square(sv_sqr) * (y + square(a * w)) + + sv_sqr * (1 - (2.0 * a * v * w)) + square(v)) + / square(one_plus_svsqr_y); + const auto error = (err - error_term) + two_log_a; + const auto n_terms_small_t + = wiener5_n_terms_small_t( + y, a, w_value, error); + const auto n_terms_large_t + = wiener5_gradient_large_reaction_time_terms( + y, a, w_value, error); + auto wiener_res = wiener5_log_sum_exp( + y, a, w_value, n_terms_small_t, n_terms_large_t); + auto&& result = wiener_res.first; + auto&& newsign = wiener_res.second; + const auto error_log_density + = log(fmax(fabs(density_part_one - 1.5 / y), fabs(density_part_one))); + const auto log_density = wiener5_density( + y, a, v_value, w_value, sv, err - error_log_density); + if (2 * n_terms_small_t < n_terms_large_t) { + auto ans = density_part_one - 1.5 / y + + newsign + * exp(error_term - two_log_a - 1.5 * LOG_TWO - LOG_SQRT_PI + - 3.5 * log_y_asq + result - log_density); + return WrtLog ? ans * exp(log_density) : ans; + } else { + auto ans = density_part_one + - newsign + * exp(error_term - two_log_a + 3.0 * LOG_PI - LOG_TWO + + result - log_density); + return WrtLog ? ans * exp(log_density) : ans; + } +} + +/** + * Calculate the derivative of the wiener5 density w.r.t. 'a' + * + * @tparam WrtLog Whether to return the derivative w.r.t. + * the natural (true) or log-scale (false) density + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_v type of drift rate + * @tparam T_w type of relative starting point + * @tparam T_sv type of inter-trial variability in v + * @tparam T_err type of log error + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param v_value The drift rate + * @param w_value The relative starting point + * @param sv The inter-trial variability of the drift rate + * @param err The log error tolerance + * @return Gradient w.r.t. a + */ +template +inline auto wiener5_grad_a(const T_y& y, const T_a& a, const T_v& v_value, + const T_w& w_value, const T_sv& sv, + T_err&& err = log(1e-12)) noexcept { + const auto two_log_a = 2 * log(a); + const auto log_y_asq = log(y) - two_log_a; + const auto error_term + = wiener5_compute_error_term(y, a, v_value, w_value, sv); + const auto w = 1.0 - w_value; + const auto v = -v_value; + const auto sv_sqr = square(sv); + const auto one_plus_svsqr_y = 1 + sv_sqr * y; + const auto density_part_one + = (-v * w + sv_sqr * square(w) * a) / one_plus_svsqr_y; + const auto error = err - error_term + 3 * log(a) - log(y) - LOG_TWO; + + const auto n_terms_small_t + = wiener5_n_terms_small_t( + y, a, w_value, error); + const auto n_terms_large_t + = wiener5_gradient_large_reaction_time_terms( + y, a, w_value, error); + auto wiener_res = wiener5_log_sum_exp( + y, a, w_value, n_terms_small_t, n_terms_large_t); + auto&& result = wiener_res.first; + auto&& newsign = wiener_res.second; + const auto error_log_density = log( + fmax(fabs(density_part_one + 1.0 / a), fabs(density_part_one - 2.0 / a))); + const auto log_density = wiener5_density( + y, a, v_value, w_value, sv, err - error_log_density); + if (2 * n_terms_small_t < n_terms_large_t) { + auto ans + = density_part_one + 1.0 / a + - newsign + * exp(-0.5 * LOG_TWO - LOG_SQRT_PI - 2.5 * log(y) + + 2.0 * two_log_a + error_term + result - log_density); + return WrtLog ? ans * exp(log_density) : ans; + } else { + auto ans = density_part_one - 2.0 / a + + newsign + * exp(log(y) + error_term - 3 * (log(a) - LOG_PI) + result + - log_density); + return WrtLog ? ans * exp(log_density) : ans; + } +} + +/** + * Calculate the derivative of the wiener5 density w.r.t. 'v' + * + * @tparam WrtLog Whether to return the derivative w.r.t. + * the natural (true) or log-scale (false) density + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_v type of drift rate + * @tparam T_w type of relative starting point + * @tparam T_sv type of inter-trial variability in v + * @tparam T_err type of log error + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param v_value The drift rate + * @param w_value The relative starting point + * @param sv The inter-trial variability of the drift rate + * @param err The log error tolerance + * @return Gradient w.r.t. v + */ +template +inline auto wiener5_grad_v(const T_y& y, const T_a& a, const T_v& v_value, + const T_w& w_value, const T_sv& sv, + T_err&& err = log(1e-12)) noexcept { + auto ans = (a * (1 - w_value) - v_value * y) / (1.0 + square(sv) * y); + if (WrtLog) { + return ans * wiener5_density(y, a, v_value, w_value, sv, err); + } else { + return ans; + } +} + +/** + * Calculate the derivative of the wiener5 density w.r.t. 'w' + * + * @tparam WrtLog Whether to return the derivative w.r.t. + * the natural (true) or log-scale (false) density + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_v type of drift rate + * @tparam T_w type of relative starting point + * @tparam T_sv type of inter-trial variability in v + * @tparam T_err type of log error + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param v_value The drift rate + * @param w_value The relative starting point + * @param sv The inter-trial variability of the drift rate + * @param err The log error tolerance + * @return Gradient w.r.t. w + */ +template +inline auto wiener5_grad_w(const T_y& y, const T_a& a, const T_v& v_value, + const T_w& w_value, const T_sv& sv, + T_err&& err = log(1e-12)) noexcept { + const auto two_log_a = 2 * log(a); + const auto log_y_asq = log(y) - two_log_a; + const auto error_term + = wiener5_compute_error_term(y, a, v_value, w_value, sv); + const auto w = 1.0 - w_value; + const auto v = -v_value; + const auto sv_sqr = square(sv); + const auto one_plus_svsqr_y = 1 + sv_sqr * y; + const auto density_part_one + = (-v * a + sv_sqr * square(a) * w) / one_plus_svsqr_y; + const auto error = (err - error_term); + + const auto n_terms_small_t + = wiener5_n_terms_small_t( + y, a, w_value, error); + const auto n_terms_large_t + = wiener5_gradient_large_reaction_time_terms( + y, a, w_value, error); + auto wiener_res = wiener5_log_sum_exp( + y, a, w_value, n_terms_small_t, n_terms_large_t); + auto&& result = wiener_res.first; + auto&& newsign = wiener_res.second; + const auto log_density = wiener5_density( + y, a, v_value, w_value, sv, err - log(fabs(density_part_one))); + if (2 * n_terms_small_t < n_terms_large_t) { + auto ans = -(density_part_one + - newsign + * exp(result - (log_density - error_term) + - 2.5 * log_y_asq - 0.5 * LOG_TWO - 0.5 * LOG_PI)); + return WrtLog ? ans * exp(log_density) : ans; + } else { + auto ans + = -(density_part_one + + newsign * exp(result - (log_density - error_term) + 2 * LOG_PI)); + return WrtLog ? ans * exp(log_density) : ans; + } +} + +/** + * Calculate the derivative of the wiener5 density w.r.t. 'sv' + * + * @tparam WrtLog Whether to return the derivative w.r.t. + * the natural (true) or log-scale (false) density + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_v type of drift rate + * @tparam T_w type of relative starting point + * @tparam T_sv type of inter-trial variability in v + * @tparam T_err type of log error + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param v_value The drift rate + * @param w_value The relative starting point + * @param sv The inter-trial variability of the drift rate + * @param err The log error tolerance + * @return Gradient w.r.t. sv + */ +template +inline auto wiener5_grad_sv(const T_y& y, const T_a& a, const T_v& v_value, + const T_w& w_value, const T_sv& sv, + T_err&& err = log(1e-12)) noexcept { + const auto one_plus_svsqr_y = 1 + square(sv) * y; + const auto w = 1.0 - w_value; + const auto v = -v_value; + const auto t1 = -y / one_plus_svsqr_y; + const auto t2 = (square(a * w) + 2 * a * v * w * y + square(v * y)) + / square(one_plus_svsqr_y); + const auto ans = sv * (t1 + t2); + return WrtLog ? ans * wiener5_density(y, a, v_value, w_value, sv, err) + : ans; +} + +/** + * Utility function for replacing a value with a specified error value + * + * @tparam NestedIndex index of error position in tuple + * @tparam Scalar1 type of argument to be replaced + * @tparam Scalar2 type of error to replace + * + * @param arg argument to be replaced + * @param err argument to replace + */ +template +inline void assign_err(Scalar1 arg, Scalar2 err) { + arg = err; +} + +/** + * Utility function for replacing a value with a specified error value, + * overload for use when the value is stored within a tuple. + * + * @tparam NestedIndex index of element in tuple to be replaced + * @tparam Scalar type of error to replace + * @tparam TArgs type of arguments to be replaced + * + * @param args_tuple argument tuple to be replaced + * @param err argument to replace + */ +template +inline void assign_err(std::tuple& args_tuple, Scalar err) { + std::get(args_tuple) = err; +} + +/** + * Utility function for estimating a function with a given set of arguments, + * checking the result against a provided error tolerance, and re-estimating + * the function with the increased error if it fails. + * + * @tparam ErrIndex Position of the error argument in the provided arguments + * @tparam GradW7 Whether the gradient of w is computed + * @tparam NestedIndex Nested position if the error argument is within a tuple + * @tparam LogResult Whether result is on log- or on natural-scale + * @tparam F Type of functor + * @tparam T_err type of error + * @tparam ArgsTupleT Type of tuple of arguments for functor + * + * @param functor Function to apply + * @param err Error value to check against + * @param args_tuple Tuple of arguments to pass to functor + */ +template +inline auto estimate_with_err_check(F&& functor, T_err&& err, + ArgsTupleT&&... args_tuple) { + auto result = functor(args_tuple...); + auto log_fabs_result = LogResult ? log(fabs(result)) : fabs(result); + if (log_fabs_result < err) { + log_fabs_result = is_inf(log_fabs_result) ? 0 : log_fabs_result; + auto err_args_tuple = std::make_tuple(args_tuple...); + const auto new_error + = GradW7 ? err + log_fabs_result + LOG_TWO : err + log_fabs_result; + assign_err(std::get(err_args_tuple), new_error); + result + = math::apply([](auto&& func, auto&&... args) { return func(args...); }, + err_args_tuple, functor); + } + return result; +} +} // namespace internal + +/** + * Log-density function for the 5-parameter Wiener density. + * See 'wiener_lpdf' for more comprehensive documentation + * + * @tparam T_y type of scalar + * @tparam T_a type of boundary + * @tparam T_t0 type of non-decision time + * @tparam T_w type of relative starting point + * @tparam T_v type of drift rate + * @tparam T_sv type of inter-trial variability of drift rate + * @tparam T_precision type of precision + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param t0 The non-decision time + * @param w The relative starting point + * @param v The drift rate + * @param sv The inter-trial variability of the drift rate + * @param precision_derivatives Level of precision in estimation + * @return The log of the Wiener first passage time density with + * the specified arguments for upper boundary responses + */ +template +inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, + const T_w& w, const T_v& v, const T_sv& sv, + const double& precision_derivatives = 1e-4) { + using T_partials_return = partials_return_t; + using ret_t = return_type_t; + if (!include_summand::value) { + return ret_t(0.0); + } + using T_y_ref = ref_type_if_t::value, T_y>; + using T_a_ref = ref_type_if_t::value, T_a>; + using T_t0_ref = ref_type_if_t::value, T_t0>; + using T_w_ref = ref_type_if_t::value, T_w>; + using T_v_ref = ref_type_if_t::value, T_v>; + using T_sv_ref = ref_type_if_t::value, T_sv>; + + static constexpr const char* function_name = "wiener5_lpdf"; + + check_consistent_sizes(function_name, "Random variable", y, + "Boundary separation", a, "Drift rate", v, + "A-priori bias", w, "Nondecision time", t0, + "Inter-trial variability in drift rate", sv); + + T_y_ref y_ref = y; + T_a_ref a_ref = a; + T_t0_ref t0_ref = t0; + T_w_ref w_ref = w; + T_v_ref v_ref = v; + T_sv_ref sv_ref = sv; + + decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + decltype(auto) a_val = to_ref(as_value_column_array_or_scalar(a_ref)); + decltype(auto) v_val = to_ref(as_value_column_array_or_scalar(v_ref)); + decltype(auto) w_val = to_ref(as_value_column_array_or_scalar(w_ref)); + decltype(auto) t0_val = to_ref(as_value_column_array_or_scalar(t0_ref)); + decltype(auto) sv_val = to_ref(as_value_column_array_or_scalar(sv_ref)); + check_positive_finite(function_name, "Random variable", y_val); + check_positive_finite(function_name, "Boundary separation", a_val); + check_finite(function_name, "Drift rate", v_val); + check_less(function_name, "A-priori bias", w_val, 1); + check_greater(function_name, "A-priori bias", w_val, 0); + check_nonnegative(function_name, "Nondecision time", t0_val); + check_finite(function_name, "Nondecision time", t0_val); + check_nonnegative(function_name, "Inter-trial variability in drift rate", + sv_val); + check_finite(function_name, "Inter-trial variability in drift rate", sv_val); + + if (size_zero(y, a, t0, w, v, sv)) { + return ret_t(0.0); + } + const size_t N = max_size(y, a, t0, w, v, sv); + if (!N) { + return ret_t(0.0); + } + + scalar_seq_view y_vec(y_ref); + scalar_seq_view a_vec(a_ref); + scalar_seq_view t0_vec(t0_ref); + scalar_seq_view w_vec(w_ref); + scalar_seq_view v_vec(v_ref); + scalar_seq_view sv_vec(sv_ref); + const size_t N_y_t0 = max_size(y, t0); + + for (size_t i = 0; i < N_y_t0; ++i) { + if (y_vec[i] <= t0_vec[i]) { + std::stringstream msg; + msg << ", but must be greater than nondecision time = " << t0_vec[i]; + std::string msg_str(msg.str()); + throw_domain_error(function_name, "Random variable", y_vec[i], " = ", + msg_str.c_str()); + } + } + + const auto log_error_density = log(1e-6); + const auto log_error_derivative = log(precision_derivatives); + const double log_error_absolute_val = log(1e-12); + const T_partials_return log_error_absolute = log_error_absolute_val; + T_partials_return log_density = 0.0; + auto ops_partials + = make_partials_propagator(y_ref, a_ref, t0_ref, w_ref, v_ref, sv_ref); + + static constexpr double LOG_FOUR = LOG_TWO + LOG_TWO; + + // calculate density and partials + for (size_t i = 0; i < N; i++) { + // Calculate 4-parameter model without inter-trial variabilities (if + // sv_vec[i] == 0) or 5-parameter model with inter-trial variability in + // drift rate (if sv_vec[i] != 0) + + const auto y_value = y_vec.val(i); + const auto a_value = a_vec.val(i); + const auto t0_value = t0_vec.val(i); + const auto w_value = w_vec.val(i); + const auto v_value = v_vec.val(i); + const auto sv_value = sv_vec.val(i); + using internal::GradientCalc; + auto l_density = internal::estimate_with_err_check<5, 0, GradientCalc::OFF, + GradientCalc::OFF>( + [](auto&&... args) { + return internal::wiener5_density(args...); + }, + log_error_density - LOG_TWO, y_value - t0_value, a_value, v_value, + w_value, sv_value, log_error_absolute); + + log_density += l_density; + + const auto new_est_err = l_density + log_error_derivative - LOG_FOUR; + + // computation of derivative for t and precision check in order to give + // the value as deriv_y to edge1 and as -deriv_y to edge5 + const auto deriv_y + = internal::estimate_with_err_check<5, 0, GradientCalc::OFF, + GradientCalc::ON>( + [](auto&&... args) { + return internal::wiener5_grad_t(args...); + }, + new_est_err, y_value - t0_value, a_value, v_value, w_value, + sv_value, log_error_absolute); + + // computation of derivatives and precision checks + if (!is_constant_all::value) { + partials<0>(ops_partials)[i] = deriv_y; + } + if (!is_constant_all::value) { + partials<1>(ops_partials)[i] + = internal::estimate_with_err_check<5, 0, GradientCalc::OFF, + GradientCalc::ON>( + [](auto&&... args) { + return internal::wiener5_grad_a(args...); + }, + new_est_err, y_value - t0_value, a_value, v_value, w_value, + sv_value, log_error_absolute); + } + if (!is_constant_all::value) { + partials<2>(ops_partials)[i] = -deriv_y; + } + if (!is_constant_all::value) { + partials<3>(ops_partials)[i] + = internal::estimate_with_err_check<5, 0, GradientCalc::OFF, + GradientCalc::ON>( + [](auto&&... args) { + return internal::wiener5_grad_w(args...); + }, + new_est_err, y_value - t0_value, a_value, v_value, w_value, + sv_value, log_error_absolute); + } + if (!is_constant_all::value) { + partials<4>(ops_partials)[i] + = internal::wiener5_grad_v( + y_value - t0_value, a_value, v_value, w_value, sv_value, + log_error_absolute_val); + } + if (!is_constant_all::value) { + partials<5>(ops_partials)[i] + = internal::wiener5_grad_sv( + y_value - t0_value, a_value, v_value, w_value, sv_value, + log_error_absolute_val); + } + } // end for loop + return ops_partials.build(log_density); +} // end wiener_lpdf + +// ToDo: delete old wiener_lpdf implementation to use this one +// template +// inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, +// const T_w& w, const T_v& v, +// const double& precision_derivatives = 1e-4) { +// return wiener_lpdf(y, a, t0, w, v, 0, precision_derivatives); +//} // end wiener_lpdf + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/wiener_full_lpdf.hpp b/stan/math/prim/prob/wiener_full_lpdf.hpp new file mode 100644 index 00000000000..be4496c7484 --- /dev/null +++ b/stan/math/prim/prob/wiener_full_lpdf.hpp @@ -0,0 +1,616 @@ +#ifndef STAN_MATH_PRIM_PROB_WIENER_FULL_LPDF_HPP +#define STAN_MATH_PRIM_PROB_WIENER_FULL_LPDF_HPP + +#include +#include +#include + +namespace stan { +namespace math { +namespace internal { + +/** + * Calculate the derivative of the wiener7 density w.r.t. 'sw' + * + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_v type of drift rate + * @tparam T_w type of relative starting point + * @tparam T_sv type of inter-trial variability in v + * @tparam T_sw type of inter-trial variability in w + * @tparam T_err type of log error tolerance + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param v The drift rate + * @param w The relative starting point + * @param sv The inter-trial variability of the drift rate + * @param sw The inter-trial variability of the relative starting point + * @param log_error The log error tolerance + * @return Gradient w.r.t. sw + */ +template +inline auto wiener7_grad_sw(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, const T_sv& sv, const T_sw& sw, + T_err log_error) { + auto low = w - sw / 2.0; + const auto lower_value + = wiener5_density(y, a, v, low, sv, log_error); + auto high = w + sw / 2.0; + const auto upper_value + = wiener5_density(y, a, v, high, sv, log_error); + return 0.5 * (lower_value + upper_value) / sw; +} + +/** + * Helper function for agnostically calling wiener5 functions + * (to be integrated over) or directly calling wiener7 functions, + * accounting for the different number of arguments. + * + * Specialisation for wiener5 functions + * + * @tparam GradSW Whether the gradient of sw is computed + * @tparam F Type of Gradient/density functor + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_v type of drift rate + * @tparam T_w type of relative starting point + * @tparam T_sv type of inter-trial variability in v + * @tparam T_sw type of inter-trial variability in w + * @tparam T_err type of log error tolerance + * + * @param functor Gradient/density functor to apply + * @param y_diff A scalar variable; the reaction time in seconds without + * non-decision time + * @param a The boundary separation + * @param v The drift rate + * @param w The relative starting point + * @param sv The inter-trial variability of the drift rate + * @param sw The inter-trial variability of the relative starting point + * @param log_error The log error tolerance + * @return Functor applied to arguments + */ +template * = nullptr> +inline auto conditionally_grad_sw(F&& functor, T_y&& y_diff, T_a&& a, T_v&& v, + T_w&& w, T_sv&& sv, T_sw&& sw, + T_err&& log_error) { + return functor(y_diff, a, v, w, sv, log_error); +} + +/** + * Helper function for agnostically calling wiener5 functions + * (to be integrated over) or directly calling wiener7 functions, + * accounting for the different number of arguments. + * + * Specialisation for wiener7 functions + * + * @tparam GradSW Whether the gradient of sw is computed + * @tparam F Type of Gradient/density functor + * @tparam T_y type of scalar variable + * @tparam T_a type of boundary separation + * @tparam T_v type of drift rate + * @tparam T_w type of relative starting point + * @tparam T_sv type of inter-trial variability in v + * @tparam T_sw type of inter-trial variability in w + * @tparam T_err type of log error tolerance + * + * @param functor Gradient/density functor to apply + * @param y_diff A scalar variable; the reaction time in seconds without + * non-decision time + * @param a The boundary separation + * @param v The drift rate + * @param w The relative starting point + * @param sv The inter-trial variability of the drift rate + * @param sw The inter-trial variability of the relative starting point + * @param log_error The log error tolerance + * @return Functor applied to arguments + */ +template * = nullptr> +inline auto conditionally_grad_sw(F&& functor, T_y&& y_diff, T_a&& a, T_v&& v, + T_w&& w, T_sv&& sv, T_sw&& sw, + T_err&& log_error) { + return functor(y_diff, a, v, w, sv, sw, log_error); +} + +/** + * Implementation function for preparing arguments and functor to be passed + * to the hcubature() function for calculating wiener7 parameters via + * integration + * + * @tparam GradSW Whether the gradient of sw is computed + * @tparam GradW7 Whether the gradient of w is computed in the full model + * @tparam Wiener7FunctorT Type of functor + * @tparam T_err Type of error in hcubature + * @tparam Targs... Types of arguments in parameter pack + * + * @param wiener7_functor Gradient/density functor to apply + * @param hcubature_err Error tolerance for calculation + * @param args Additional arguments to be passed to the hcubature function + * @return Wiener7 density or gradient calculated by integration + */ +template +inline auto wiener7_integrate(const Wiener7FunctorT& wiener7_functor, + T_err&& hcubature_err, TArgs&&... args) { + const auto functor = [&wiener7_functor](auto&&... integration_args) { + return hcubature( + [&wiener7_functor](auto&& x, auto&& y, auto&& a, auto&& v, auto&& w, + auto&& t0, auto&& sv, auto&& sw, auto&& st0, + auto&& log_error) { + using ret_t = return_type_t; + scalar_seq_view x_vec(x); + const auto sw_val = GradSW ? 0 : sw; + const auto new_t0 = t0 + st0 * x_vec[(sw_val != 0) ? 1 : 0]; + if (y - new_t0 <= 0) { + return ret_t(0.0); + } else { + const auto new_w = w + sw_val * (x_vec[0] - 0.5); + return ret_t(conditionally_grad_sw( + wiener7_functor, y - new_t0, a, v, new_w, sv, sw, log_error)); + } + }, + integration_args...); + }; + return estimate_with_err_check<0, 8, GradW7, GradientCalc::ON>( + functor, hcubature_err, args...); +} +} // namespace internal + +/** \ingroup prob_dists + * The log of the first passage time density function for a (Wiener) + * drift diffusion model with up to 7 parameters, where + * \f$y\in \mathbb{R}_{+}\f$ is the reacion time, \f$a \in \mathbb{R}_{+}\f$ + * the boundary separation, \f$t_0 \in \mathbb{R}_{\geq 0}\f$ the non-decision + time, + * \f$w \in (0, 1)\f$ the relative starting point (aka a-priori bias), + * \f$v \in \mathbb{R}\f$ the drifte rate, \f$s_v \in + * \mathbb{R}_{\geq 0}\f$ the inter-trial variability of the drift rate, + * \f$s_w \in [0, 1)\f$ the inter-trial variability of the relative starting + * point, and \f$s_{t_0} \in \mathbb{R}_{\geq 0}\f$ the inter-trial variability + * of + * the non-decision time. + * + * + \f{eqnarray*}{ + y &\sim& \text{wiener_full}(a,t_0,w,v,s_v,s_w,s_{t_0}) \\ + \log(p(y|a,v,w,t_0,s_v,s_w,s_{t_0})) &=& \log(\frac{1}{s_{t_0}} + \int_{t_0}^{t_o + s_{t_0}} \frac{1}{s_{w}}\int_{w -0.5s_w}^{w + 0.5s_{w}} + \int_{-\infty}^{\infty} p_3(y-\tau_0|a,\nu,\omega) \\ + &&\times\frac{1}{\sqrt{2\pi s_\nu^2}} + \mathbb{e}^{-\frac{(\nu-v)^2}{2s_\nu^2}} \ d\nu \ d\omega \ d\tau_0) \\ + &=& \log(\frac{1}{s_{t_0}} + \int_{t_0}^{t_o + s_{t_0}} \frac{1}{s_{w}}\int_{w -0.5s_w}^{w + 0.5s_{w}} + M \times p_3(y-\tau_0|a,v,\omega) \ d\omega \ d\tau_0), + \f} + * where \f$M\f$ and \f$p_3()\f$ are defined, by using \f$t:=y-\tau_0\f$, as + \f{eqnarray*}{ + M &:=& \frac{1}{\sqrt{1+s^2_v t}} + \mathbb{e}^{av\omega+\frac{v^2t}{2}+\frac{s^2_v a^2 + \omega^2-2av\omega-v^2t}{2(1+s^2_vt)}} \text{ and} \\ p_3(t|a,v,w) &:=& + \frac{1}{a^2} \mathbb{e}^{-a v w -\frac{v^2t}{2}} f(\frac{t}{a^2}|0,1,w), \f} + * where \f$f(t^*=\frac{t}{a^2}|0,1,w)\f$ has two forms + \f{eqnarray*}{ + f_l(t^*|0,1,w) &=& \sum_{k=1}^{\infty} k\pi \mathbb{e}^{-\frac{k^2\pi^2t^*}{2}} + \sin{(k \pi w)}\text{ and} \\ + f_s(t^*|0,1,w) &=& \sum_{k=-\infty}^{\infty} \frac{1}{\sqrt{2\pi (t^*)^3}} + (w+2k) \mathbb{e}^{-\frac{(w+2k)^2}{2t^*}}, \f} + * which are selected depending on the number of components \f$k\f$, needed to + * guarantee a certain precision. + * + * Note that the parameterization for non-decision time and relative starting + point is as follows: + * \c t0 is the lower bound of the variability interval; + * \c w is the mean of the variability interval. + * + * See \b Details below for more details on how to use \c wiener_lpdf(). + * + * @tparam T_y type of scalar + * @tparam T_a type of boundary separation + * @tparam T_t0 type of non-decision time + * @tparam T_w type of relative starting point + * @tparam T_v type of drift rate + * @tparam T_sv type of inter-trial variability of drift rate + * @tparam T_sw type of inter-trial variability of relative starting point + * @tparam T_st0 type of inter-trial variability of non-decision time + * + * @param y A scalar variable; the reaction time in seconds + * @param a The boundary separation + * @param t0 The non-decision time + * @param w The relative starting point + * @param v The drift rate + * @param sv The inter-trial variability of the drift rate + * @param sw The inter-trial variability of the relative starting point + * @param st0 The inter-trial variability of the non-decision time + * @param precision_derivatives Level of precision in estimation of partial + derivatives + * @return The log of the Wiener first passage time density with + * the specified arguments for upper boundary responses + * @throw std::domain_error if non-decision time \c t0 is greater than reaction + time \c y. + * @throw std::domain_error if \c 1-sw/2 is smaller than or equal to \c w. + * @throw std::domain_error if \c sw/2 is larger than or equal to \c w. + * + * + * + * **Details** + * + * The function can be called by + @code + target += wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + @endcode + * or + @code + y ~ wiener_full(a, t0, w, v, sv, sw, st0); + @endcode + * + * By default \c wiener_lpdf() gives the log of the + * Wiener first-passage time probability density function for the \e upper + response + * boundary. To use the \e lower response boundary \c v and \c w must be changed + to + * \c -v and \c (1-w), respectively. + * + * \c sv, \c sw, \c st0, and \c t0 can be set to zero, indicating no inter-trial + * variability in \f$v\f$, no inter-trial variability in \f$w\f$, no + * inter-trial variability in \f$t_0\f$, and no non-decision time, respectively. + * If \c t0 is zero, \c st0 must be zero as well. For example, when no + inter-trial + * variability for the relative starting point is needed one can write something + like: + @code + target += wiener_lpdf(y, a, t0, w, v, sv, 0, st0) + @endcode + * If no inter-trial variability is needed at all one can write something like: + @code + target += wiener_lpdf(y, a, t0, w, v, 0, 0, 0) + @endcode + * If for some reason no non-decision time is assumed one can write something + like: + @code + target += wiener_lpdf(y, a, 0, w, v, sv, sw, 0) + @endcode + * If only inter-trial variability for the drift rate is needed can write + something like: + @code + target += wiener_lpdf(y, a, t0, w, v, sv) + @endcode + * + * To also control the precision in the estimation of the partial derivatives: + @code + target += wiener_lpdf(y, a, t0, w, v, sv, sw, st0, precision); + @endcode + * + * + * **References** + * - Blurton, S. P., Kesselmeier, M., & Gondan, M. (2017). The first-passage + * time distribution for the diffusion model with variable drift. + * *Journal of Mathematical Psychology, 76*, 7–12. + * https://doi.org/10.1016/j.jmp.2016.11.003 + * - Foster, K., & Singmann, H. (2021). Another Approximation of the + First-Passage + * Time Densities for the Ratcliff Diffusion Decision Model. + * *arXiv preprint arXiv:2104.01902* + * - Gondan, M., Blurton, S. P., & Kesselmeier, M. (2014). Even faster and even + * more accurate first-passage time densities and distributions for the Wiener + * diffusion model. *Journal of Mathematical Psychology, 60*, 20–22. + * https://doi.org/10.1016/j.jmp.2014.05.002 + * - Hartmann, R., & Klauer, K. C. (2021). Partial derivatives for the + * first-passage time distribution in Wiener diffusion models. + * *Journal of Mathematical Psychology, 103*, 102550. + * https://doi.org/10.1016/j.jmp.2021.102550 + * - Henrich, F., Hartmann, R., Pratz, V., Voss, A., & Klauer, K.C. (2023). + * The Seven-parameter Diffusion Model: An Implementation in Stan for Bayesian + * Analyses. *Behavior Research Methods*. + * - Navarro, D. J., & Fuss, I. G. (2009). Fast and accurate calculations for + * first-passage times in Wiener diffusion models. + * *Journal of Mathematical Psychology, 53*(4), 222–230. + * https://doi.org/10.1016/j.jmp.2009.02.003 + */ +template +inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, + const T_w& w, const T_v& v, const T_sv& sv, + const T_sw& sw, const T_st0& st0, + const double& precision_derivatives = 1e-4) { + using ret_t = return_type_t; + if (!include_summand::value) { + return ret_t(0); + } + + using T_y_ref = ref_type_if_t::value, T_y>; + using T_a_ref = ref_type_if_t::value, T_a>; + using T_v_ref = ref_type_if_t::value, T_v>; + using T_w_ref = ref_type_if_t::value, T_w>; + using T_t0_ref = ref_type_if_t::value, T_t0>; + using T_sv_ref = ref_type_if_t::value, T_sv>; + using T_sw_ref = ref_type_if_t::value, T_sw>; + using T_st0_ref = ref_type_if_t::value, T_st0>; + + using T_partials_return + = partials_return_t; + + static constexpr const char* function_name = "wiener_lpdf"; + check_consistent_sizes(function_name, "Random variable", y, + "Boundary separation", a, "Drift rate", v, + "A-priori bias", w, "Nondecision time", t0, + "Inter-trial variability in drift rate", sv, + "Inter-trial variability in A-priori bias", sw, + "Inter-trial variability in Nondecision time", st0); + + T_y_ref y_ref = y; + T_a_ref a_ref = a; + T_v_ref v_ref = v; + T_w_ref w_ref = w; + T_t0_ref t0_ref = t0; + T_sv_ref sv_ref = sv; + T_sw_ref sw_ref = sw; + T_st0_ref st0_ref = st0; + + decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + decltype(auto) a_val = to_ref(as_value_column_array_or_scalar(a_ref)); + decltype(auto) v_val = to_ref(as_value_column_array_or_scalar(v_ref)); + decltype(auto) w_val = to_ref(as_value_column_array_or_scalar(w_ref)); + decltype(auto) t0_val = to_ref(as_value_column_array_or_scalar(t0_ref)); + decltype(auto) sv_val = to_ref(as_value_column_array_or_scalar(sv_ref)); + decltype(auto) sw_val = to_ref(as_value_column_array_or_scalar(sw_ref)); + decltype(auto) st0_val = to_ref(as_value_column_array_or_scalar(st0_ref)); + check_positive_finite(function_name, "Random variable", y_val); + check_positive_finite(function_name, "Boundary separation", a_val); + check_finite(function_name, "Drift rate", v_val); + check_less(function_name, "A-priori bias", w_val, 1); + check_greater(function_name, "A-priori bias", w_val, 0); + check_nonnegative(function_name, "Nondecision time", t0_val); + check_finite(function_name, "Nondecision time", t0_val); + check_nonnegative(function_name, "Inter-trial variability in drift rate", + sv_val); + check_finite(function_name, "Inter-trial variability in drift rate", sv_val); + check_bounded(function_name, "Inter-trial variability in A-priori bias", + sw_val, 0, 1); + check_nonnegative(function_name, + "Inter-trial variability in Nondecision time", st0_val); + check_finite(function_name, "Inter-trial variability in Nondecision time", + st0_val); + + const size_t N = max_size(y, a, v, w, t0, sv, sw, st0); + if (N == 0) { + return ret_t(0); + } + scalar_seq_view y_vec(y_ref); + scalar_seq_view a_vec(a_ref); + scalar_seq_view v_vec(v_ref); + scalar_seq_view w_vec(w_ref); + scalar_seq_view t0_vec(t0_ref); + scalar_seq_view sv_vec(sv_ref); + scalar_seq_view sw_vec(sw_ref); + scalar_seq_view st0_vec(st0_ref); + const size_t N_y_t0 = max_size(y, t0, st0); + + for (size_t i = 0; i < N_y_t0; ++i) { + if (y_vec[i] <= t0_vec[i]) { + [&]() STAN_COLD_PATH { + std::stringstream msg; + msg << ", but must be greater than nondecision time = " << t0_vec[i]; + std::string msg_str(msg.str()); + throw_domain_error(function_name, "Random variable", y_vec[i], " = ", + msg_str.c_str()); + }(); + } + } + size_t N_beta_sw = max_size(w, sw); + for (size_t i = 0; i < N_beta_sw; ++i) { + if (unlikely(w_vec[i] - .5 * sw_vec[i] <= 0)) { + [&]() STAN_COLD_PATH { + std::stringstream msg; + msg << ", but must be smaller than 2*(A-priori bias) = " + << 2 * w_vec[i]; + std::string msg_str(msg.str()); + throw_domain_error(function_name, + "Inter-trial variability in A-priori bias", + sw_vec[i], " = ", msg_str.c_str()); + }(); + } + if (unlikely(w_vec[i] + .5 * sw_vec[i] >= 1)) { + [&]() STAN_COLD_PATH { + std::stringstream msg; + msg << ", but must be smaller than 2*(1-A-priori bias) = " + << 2 * (1 - w_vec[i]); + std::string msg_str(msg.str()); + throw_domain_error(function_name, + "Inter-trial variability in A-priori bias", + sw_vec[i], " = ", msg_str.c_str()); + }(); + } + } + // precision for density + const T_partials_return log_error_density = log(1e-6); + // precision for derivatives (controllable by user) + const auto error_bound = precision_derivatives; + const auto log_error_derivative = log(error_bound); + const T_partials_return absolute_error_hcubature = 0.0; + // eps_rel(Integration) + const T_partials_return relative_error_hcubature = .9 * error_bound; + const T_partials_return log_error_absolute = log(1e-12); + const int maximal_evaluations_hcubature = 6000; + T_partials_return log_density = 0.0; + auto ops_partials = make_partials_propagator(y_ref, a_ref, t0_ref, w_ref, + v_ref, sv_ref, sw_ref, st0_ref); + ret_t result = 0; + + // calculate density and partials + for (size_t i = 0; i < N; i++) { + if (sw_vec[i] == 0 && st0_vec[i] == 0) { + result += wiener_lpdf(y_vec[i], a_vec[i], t0_vec[i], w_vec[i], + v_vec[i], sv_vec[i], precision_derivatives); + continue; + } + const T_partials_return y_value = y_vec.val(i); + const T_partials_return a_value = a_vec.val(i); + const T_partials_return v_value = v_vec.val(i); + const T_partials_return w_value = w_vec.val(i); + const T_partials_return t0_value = t0_vec.val(i); + const T_partials_return sv_value = sv_vec.val(i); + const T_partials_return sw_value = sw_vec.val(i); + const T_partials_return st0_value = st0_vec.val(i); + const int dim = (sw_value != 0) + (st0_value != 0); + check_positive(function_name, + "(Inter-trial variability in A-priori bias) + " + "(Inter-trial variability in nondecision time)", + dim); + + Eigen::Matrix xmin = Eigen::VectorXd::Zero(dim); + Eigen::Matrix xmax = Eigen::VectorXd::Ones(dim); + if (st0_value != 0) { + xmax[dim - 1] = fmin(1.0, (y_value - t0_value) / st0_value); + } + + T_partials_return hcubature_err + = log_error_absolute - log_error_density + LOG_TWO + 1; + using internal::GradientCalc; + const auto params = std::make_tuple(y_value, a_value, v_value, w_value, + t0_value, sv_value, sw_value, st0_value, + log_error_absolute - LOG_TWO); + T_partials_return density + = internal::wiener7_integrate( + [](auto&&... args) { + return internal::wiener5_density(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2); + log_density += log(density); + hcubature_err = log_error_absolute - log_error_derivative + + log(fabs(density)) + LOG_TWO + 1; + + // computation of derivative for t and precision check in order to give + // the value as deriv_t to edge1 and as -deriv_t to edge5 + const T_partials_return deriv_t_7 + = internal::wiener7_integrate( + [](auto&&... args) { + return internal::wiener5_grad_t(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / density; + + // computation of derivatives and precision checks + T_partials_return derivative; + if (!is_constant_all::value) { + partials<0>(ops_partials)[i] = deriv_t_7; + } + if (!is_constant_all::value) { + partials<1>(ops_partials)[i] + = internal::wiener7_integrate( + [](auto&&... args) { + return internal::wiener5_grad_a(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / density; + } + if (!is_constant_all::value) { + partials<2>(ops_partials)[i] = -deriv_t_7; + } + if (!is_constant_all::value) { + partials<3>(ops_partials)[i] + = internal::wiener7_integrate( + [](auto&&... args) { + return internal::wiener5_grad_w(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / density; + } + if (!is_constant_all::value) { + partials<4>(ops_partials)[i] + = internal::wiener7_integrate( + [](auto&&... args) { + return internal::wiener5_grad_v(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / density; + } + if (!is_constant_all::value) { + partials<5>(ops_partials)[i] + = internal::wiener7_integrate( + [](auto&&... args) { + return internal::wiener5_grad_sv(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / density; + } + if (!is_constant_all::value) { + if (sw_value == 0) { + partials<6>(ops_partials)[i] = 0; + } else { + if (st0_value == 0) { + derivative = internal::estimate_with_err_check< + 6, 0, GradientCalc::OFF, GradientCalc::ON>( + [](auto&&... args) { return internal::wiener7_grad_sw(args...); }, + hcubature_err, y_value - t0_value, a_value, v_value, w_value, + sv_value, sw_value, log_error_absolute - LOG_TWO); + } else { + derivative = internal::wiener7_integrate( + [](auto&&... args) { return internal::wiener7_grad_sw(args...); }, + hcubature_err, params, 1, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2); + } + partials<6>(ops_partials)[i] = derivative / density - 1.0 / sw_value; + } + } + if (!is_constant_all::value) { + T_partials_return f; + if (st0_value == 0) { + partials<7>(ops_partials)[i] = 0; + } else if (y_value - (t0_value + st0_value) <= 0) { + partials<7>(ops_partials)[i] = -1 / st0_value; + } else { + const T_partials_return t0_st0 = t0_value + st0_value; + if (sw_value == 0) { + f = internal::estimate_with_err_check<5, 0, GradientCalc::OFF, + GradientCalc::ON>( + [](auto&&... args) { + return internal::wiener5_density(args...); + }, + log_error_derivative + log(st0_value), y_value - t0_st0, a_value, + v_value, w_value, sv_value, log_error_absolute - LOG_TWO); + } else { + const T_partials_return new_error = log_error_absolute - LOG_TWO; + auto params_st + = std::make_tuple(y_value, a_value, v_value, w_value, t0_st0, + sv_value, sw_value, 0.0, new_error); + f = internal::wiener7_integrate( + [](auto&&... args) { + return internal::wiener5_density(args...); + }, + hcubature_err, params_st, 1, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2.0); + } + partials<7>(ops_partials)[i] = -1 / st0_value + f / st0_value / density; + } + } + } + return result + ops_partials.build(log_density); +} +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/mix/prob/util.hpp b/test/unit/math/mix/prob/util.hpp new file mode 100644 index 00000000000..5b1dea8ad59 --- /dev/null +++ b/test/unit/math/mix/prob/util.hpp @@ -0,0 +1,48 @@ +#ifndef STAN_TEST_UNIT_MATH_MIX_PROB_UTIL_HPP +#define STAN_TEST_UNIT_MATH_MIX_PROB_UTIL_HPP + +#include + +class Wiener7MixArgs : public ::testing::Test { + public: + const char* function = "function"; + void SetUp() {} + + Eigen::VectorXd y{{1.0, 0.05, 3.0, 1.0, 1.0, 2.0, 1.0, 1e-14, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 100.0, 1e10}}; + Eigen::VectorXd a{{2.5, 2.5, 2.5, 0.2, 15.0, 1e-13, 2.5e10, 2.5, 2.5, 2.5, + 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5}}; + Eigen::VectorXd t0{{2.2, 0.01, 2.2, 0.2, 0.2, 1.99, 1e-13, 0.0, 0.2, 0.2, 0.2, + 0.2, 0.2, 0.2, 0.2, 0.2, 1e13, 0.2, 0.2}}; + Eigen::VectorXd w{{1e-13, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1e-13, 0.999999999999, + 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5}}; + Eigen::VectorXd v{{2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1e-12, 0.0, 2.0, 0.0, -5.0, + 8.0, 2.0, 2.0, 2.0, 2.0, 2.0, 4e10, -4e10}}; + Eigen::VectorXd sv{{0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 1e-14, 0.2, 0.2, + 0.2, 0.02, 9.0, 0.01, 0.2, 4, 0.2, 1e14}}; + Eigen::VectorXd sw{{1, 1e-13, 0.9999999999, 0.1, 0.4, 0.1, 0.1, 0.8, 0.1, 0.1, + 0.1, 0.1, 0.1, 0.1, 0.1, 0.02, 0.4, 0.1, 0.1}}; + Eigen::VectorXd st0{{0.0, 1e-13, 1e16, 0.3, 0.3, 0.3, 1, 10, 100, 0.3, 0.3, + 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.02, 5.0}}; +}; + +class Wiener5MixArgs : public ::testing::Test { + public: + const char* function = "function"; + void SetUp() {} + + Eigen::VectorXd y{{1.0, 0.05, 3.0, 1.0, 1.0, 2.0, 1.0, 1e-14, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 100.0, 1e10}}; + Eigen::VectorXd a{{2.5, 2.5, 2.5, 0.2, 15.0, 1e-13, 2.5e10, 2.5, 2.5, 2.5, + 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5}}; + Eigen::VectorXd t0{{2.2, 0.01, 2.2, 0.2, 0.2, 1.99, 1e-13, 0.0, 0.2, 0.2, 0.2, + 0.2, 0.2, 0.2, 0.2, 0.2, 1e13, 0.2, 0.2}}; + Eigen::VectorXd w{{0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1e-13, 0.999999999999, + 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5}}; + Eigen::VectorXd v{{2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1e-12, 0.0, 2.0, 0.0, -5.0, + 8.0, 2.0, 2.0, 2.0, 2.0, 2.0, 4e10, -4e10}}; + Eigen::VectorXd sv{{0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 1e-14, 0.2, 0.2, + 0.2, 0.02, 9.0, 0.01, 0.2, 4, 0.2, 1e14}}; +}; + +#endif diff --git a/test/unit/math/mix/prob/wiener5_lpdf_0_test.cpp b/test/unit/math/mix/prob/wiener5_lpdf_0_test.cpp new file mode 100644 index 00000000000..24901ac760b --- /dev/null +++ b/test/unit/math/mix/prob/wiener5_lpdf_0_test.cpp @@ -0,0 +1,36 @@ +#include +#include + +TEST(mathMixDouble5, wiener_lpdf) { + double y = 1.0; + double a = 2.0; + double t0 = 0.2; + double w = 0.5; + double v = 1.5; + double sv = 0.2; + stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); +} + +TEST(mathMixVar5, wiener_lpdf) { + using stan::math::var; + var y = 1.0; + var a = 2.0; + var t0 = 0.2; + var w = 0.5; + var v = 1.5; + var sv = 0.2; + stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); +} + +TEST(mathMixFVar5, wiener_lpdf) { + using stan::math::fvar; + using stan::math::var; + fvar y = 1.0; + fvar a = 2.0; + fvar t0 = 0.2; + fvar w = 0.5; + fvar v = 1.5; + fvar sv = 0.2; + double error = 1e-4; + stan::math::wiener_lpdf(y, a, t0, w, v, sv, error); +} diff --git a/test/unit/math/mix/prob/wiener5_lpdf_1_test.cpp b/test/unit/math/mix/prob/wiener5_lpdf_1_test.cpp new file mode 100644 index 00000000000..c516f0fdda7 --- /dev/null +++ b/test/unit/math/mix/prob/wiener5_lpdf_1_test.cpp @@ -0,0 +1,48 @@ +#include +#include +#include + +TEST_F(Wiener5MixArgs, wiener_lpdf_y_a_t0) { + auto f_y_a_t0 = [](const auto& w, const auto& v, const auto& sv) { + return [&w, &v, &sv](const auto& y, const auto& a, const auto& t0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_y_a_t0(w, v, sv), y, a, t0); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_y_a_w) { + auto f_y_a_w = [](const auto& t0, const auto& v, const auto& sv) { + return [&t0, &v, &sv](const auto& y, const auto& a, const auto& w) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_y_a_w(t0, v, sv), y, a, w); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_y_a_v) { + auto f_y_a_v = [](const auto& t0, const auto& w, const auto& sv) { + return [&t0, &w, &sv](const auto& y, const auto& a, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_y_a_v(t0, w, sv), y, a, v); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_y_a_sv) { + auto f_y_a_sv = [](const auto& t0, const auto& w, const auto& v) { + return [&t0, &w, &v](const auto& y, const auto& a, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_y_a_sv(t0, w, v), y, a, sv); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_y_t0_w) { + auto f_y_t0_w = [](const auto& a, const auto& v, const auto& sv) { + return [&a, &v, &sv](const auto& y, const auto& t0, const auto& w) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_y_t0_w(a, v, sv), y, t0, w); +} diff --git a/test/unit/math/mix/prob/wiener5_lpdf_2_test.cpp b/test/unit/math/mix/prob/wiener5_lpdf_2_test.cpp new file mode 100644 index 00000000000..b2d6acf2a3c --- /dev/null +++ b/test/unit/math/mix/prob/wiener5_lpdf_2_test.cpp @@ -0,0 +1,48 @@ +#include +#include +#include + +TEST_F(Wiener5MixArgs, wiener_lpdf_y_t0_v) { + auto f_y_t0_v = [](const auto& a, const auto& w, const auto& sv) { + return [&a, &w, &sv](const auto& y, const auto& t0, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_y_t0_v(a, w, sv), y, t0, v); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_y_t0_sv) { + auto f_y_t0_sv = [](const auto& a, const auto& w, const auto& v) { + return [&a, &w, &v](const auto& y, const auto& t0, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_y_t0_sv(a, w, v), y, t0, sv); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_y_w_v) { + auto f_y_w_v = [](const auto& a, const auto& t0, const auto& sv) { + return [&a, &t0, &sv](const auto& y, const auto& w, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_y_w_v(a, t0, sv), y, w, v); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_y_w_sv) { + auto f_y_w_sv = [](const auto& a, const auto& t0, const auto& v) { + return [&a, &t0, &v](const auto& y, const auto& w, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_y_w_sv(a, t0, v), y, w, sv); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_y_v_sv) { + auto f_y_v_sv = [](const auto& a, const auto& t0, const auto& w) { + return [&a, &t0, &w](const auto& y, const auto& v, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_y_v_sv(a, t0, w), y, v, sv); +} diff --git a/test/unit/math/mix/prob/wiener5_lpdf_3_test.cpp b/test/unit/math/mix/prob/wiener5_lpdf_3_test.cpp new file mode 100644 index 00000000000..d16285a95b3 --- /dev/null +++ b/test/unit/math/mix/prob/wiener5_lpdf_3_test.cpp @@ -0,0 +1,48 @@ +#include +#include +#include + +TEST_F(Wiener5MixArgs, wiener_lpdf_a_t0_w) { + auto f_a_t0_w = [](const auto& y, const auto& v, const auto& sv) { + return [&y, &v, &sv](const auto& a, const auto& t0, const auto& w) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_a_t0_w(y, v, sv), a, t0, w); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_a_t0_v) { + auto f_a_t0_v = [](const auto& y, const auto& w, const auto& sv) { + return [&y, &w, &sv](const auto& a, const auto& t0, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_a_t0_v(y, w, sv), a, t0, v); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_a_t0_sv) { + auto f_a_t0_sv = [](const auto& y, const auto& w, const auto& v) { + return [&y, &w, &v](const auto& a, const auto& t0, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_a_t0_sv(y, w, v), a, t0, sv); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_a_w_v) { + auto f_a_w_v = [](const auto& y, const auto& t0, const auto& sv) { + return [&y, &t0, &sv](const auto& a, const auto& w, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_a_w_v(y, t0, sv), a, w, v); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_a_w_sv) { + auto f_a_w_sv = [](const auto& y, const auto& t0, const auto& v) { + return [&y, &t0, &v](const auto& a, const auto& w, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_a_w_sv(y, t0, v), a, w, sv); +} diff --git a/test/unit/math/mix/prob/wiener5_lpdf_4_test.cpp b/test/unit/math/mix/prob/wiener5_lpdf_4_test.cpp new file mode 100644 index 00000000000..cb0423b9dbb --- /dev/null +++ b/test/unit/math/mix/prob/wiener5_lpdf_4_test.cpp @@ -0,0 +1,48 @@ +#include +#include +#include + +TEST_F(Wiener5MixArgs, wiener_lpdf_a_v_sv) { + auto f_a_v_sv = [](const auto& y, const auto& t0, const auto& w) { + return [&y, &t0, &w](const auto& a, const auto& v, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_a_v_sv(y, t0, w), a, v, sv); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_t0_w_v) { + auto f_t0_w_v = [](const auto& y, const auto& a, const auto& sv) { + return [&y, &a, &sv](const auto& t0, const auto& w, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_t0_w_v(y, a, sv), t0, w, v); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_t0_w_sv) { + auto f_t0_w_sv = [](const auto& y, const auto& a, const auto& v) { + return [&y, &a, &v](const auto& t0, const auto& w, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_t0_w_sv(y, a, v), t0, w, sv); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_t0_v_sv) { + auto f_t0_v_sv = [](const auto& y, const auto& a, const auto& w) { + return [&y, &a, &w](const auto& t0, const auto& v, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_t0_v_sv(y, a, w), t0, v, sv); +} + +TEST_F(Wiener5MixArgs, wiener_lpdf_w_v_sv) { + auto f_w_v_sv = [](const auto& y, const auto& a, const auto& t0) { + return [&y, &a, &t0](const auto& w, const auto& v, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, 1e-4); + }; + }; + stan::test::expect_ad(f_w_v_sv(y, a, t0), w, v, sv); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_0_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_0_test.cpp new file mode 100644 index 00000000000..374b2ed0083 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_0_test.cpp @@ -0,0 +1,53 @@ +#include +#include +#include + +TEST(mathMixDouble7, wiener_lpdf) { + double y = 1.0; + double a = 2.0; + double t0 = 0.2; + double w = 0.5; + double v = 1.5; + double sv = 0.2; + double sw = 0.2; + double st0 = 0.2; + stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0, 1e-4); +} + +TEST(mathMixVar7, wiener_lpdf) { + using stan::math::var; + var y = 1.0; + var a = 2.0; + var t0 = 0.2; + var w = 0.5; + var v = 1.5; + var sv = 0.2; + var sw = 0; + var st0 = 0.2; + stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); +} + +TEST(mathMixFVar7, wiener_lpdf) { + using stan::math::fvar; + using stan::math::var; + fvar y = 1.0; + fvar a = 2.0; + fvar t0 = 0.2; + fvar w = 0.5; + fvar v = 1.5; + fvar sv = 0.2; + fvar sw = 0; + fvar st0 = 0.2; + stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_a_t0) { + auto f_y_a_t0 = [](const auto& w, const auto& v, const auto& sv, + const auto& sw, const auto& st0) { + return + [&w, &v, &sv, &sw, &st0](const auto& y, const auto& a, const auto& t0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_a_t0(w, v, sv, sw, st0), y, a, t0); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_10_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_10_test.cpp new file mode 100644 index 00000000000..fb3c9537ec9 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_10_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_w_v_sv) { + auto f_w_v_sv = [](const auto& y, const auto& a, const auto& t0, + const auto& sw, const auto& st0) { + return + [&y, &a, &t0, &sw, &st0](const auto& w, const auto& v, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_w_v_sv(y, a, t0, sw, st0), w, v, sv); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_w_v_sw) { + auto f_w_v_sw = [](const auto& y, const auto& a, const auto& t0, + const auto& sv, const auto& st0) { + return + [&y, &a, &t0, &sv, &st0](const auto& w, const auto& v, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_w_v_sw(y, a, t0, sv, st0), w, v, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_w_v_st0) { + auto f_w_v_st0 = [](const auto& y, const auto& a, const auto& t0, + const auto& sv, const auto& sw) { + return + [&y, &a, &t0, &sv, &sw](const auto& w, const auto& v, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_w_v_st0(y, a, t0, sv, sw), w, v, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_w_sv_sw) { + auto f_w_sv_sw = [](const auto& y, const auto& a, const auto& t0, + const auto& v, const auto& st0) { + return + [&y, &a, &t0, &v, &st0](const auto& w, const auto& sv, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_w_sv_sw(y, a, t0, v, st0), w, sv, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_w_sv_st0) { + auto f_w_sv_st0 = [](const auto& y, const auto& a, const auto& t0, + const auto& v, const auto& sw) { + return + [&y, &a, &t0, &v, &sw](const auto& w, const auto& sv, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_w_sv_st0(y, a, t0, v, sw), w, sv, st0); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_11_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_11_test.cpp new file mode 100644 index 00000000000..d5a11e954d6 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_11_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_w_sw_st0) { + auto f_w_sw_st0 = [](const auto& y, const auto& a, const auto& t0, + const auto& v, const auto& sv) { + return + [&y, &a, &t0, &v, &sv](const auto& w, const auto& sw, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_w_sw_st0(y, a, t0, v, sv), w, sw, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_v_sv_sw) { + auto f_v_sv_sw = [](const auto& y, const auto& a, const auto& t0, + const auto& w, const auto& st0) { + return + [&y, &a, &t0, &w, &st0](const auto& v, const auto& sv, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_v_sv_sw(y, a, t0, w, st0), v, sv, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_v_sv_st0) { + auto f_v_sv_st0 = [](const auto& y, const auto& a, const auto& t0, + const auto& w, const auto& sw) { + return + [&y, &a, &t0, &w, &sw](const auto& v, const auto& sv, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_v_sv_st0(y, a, t0, w, sw), v, sv, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_v_sw_st0) { + auto f_v_sw_st0 = [](const auto& y, const auto& a, const auto& t0, + const auto& w, const auto& sv) { + return + [&y, &a, &t0, &w, &sv](const auto& v, const auto& sw, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_v_sw_st0(y, a, t0, w, sv), v, sw, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_sv_sw_st0) { + auto f_sv_sw_st0 = [](const auto& y, const auto& a, const auto& t0, + const auto& w, const auto& v) { + return + [&y, &a, &t0, &w, &v](const auto& sv, const auto& sw, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_sv_sw_st0(y, a, t0, w, v), sv, sw, st0); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_1_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_1_test.cpp new file mode 100644 index 00000000000..3bc1985fc14 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_1_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_a_w) { + auto f_y_a_w = [](const auto& t0, const auto& v, const auto& sv, + const auto& sw, const auto& st0) { + return + [&t0, &v, &sv, &sw, &st0](const auto& y, const auto& a, const auto& w) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_a_w(t0, v, sv, sw, st0), y, a, w); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_a_v) { + auto f_y_a_v = [](const auto& t0, const auto& w, const auto& sv, + const auto& sw, const auto& st0) { + return + [&t0, &w, &sv, &sw, &st0](const auto& y, const auto& a, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_a_v(t0, w, sv, sw, st0), y, a, v); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_a_sv) { + auto f_y_a_sv = [](const auto& t0, const auto& w, const auto& v, + const auto& sw, const auto& st0) { + return + [&t0, &w, &v, &sw, &st0](const auto& y, const auto& a, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_a_sv(t0, w, v, sw, st0), y, a, sv); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_a_sw) { + auto f_y_a_sw = [](const auto& t0, const auto& w, const auto& v, + const auto& sv, const auto& st0) { + return + [&t0, &w, &v, &sv, &st0](const auto& y, const auto& a, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_a_sw(t0, w, v, sv, st0), y, a, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_a_st0) { + auto f_y_a_st0 = [](const auto& t0, const auto& w, const auto& v, + const auto& sv, const auto& sw) { + return + [&t0, &w, &v, &sv, &sw](const auto& y, const auto& a, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_a_st0(t0, w, v, sv, sw), y, a, st0); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_2_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_2_test.cpp new file mode 100644 index 00000000000..531940c6922 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_2_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_t0_w) { + auto f_y_t0_w = [](const auto& a, const auto& v, const auto& sv, + const auto& sw, const auto& st0) { + return + [&a, &v, &sv, &sw, &st0](const auto& y, const auto& t0, const auto& w) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_t0_w(a, v, sv, sw, st0), y, t0, w); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_t0_v) { + auto f_y_t0_v = [](const auto& a, const auto& w, const auto& sv, + const auto& sw, const auto& st0) { + return + [&a, &w, &sv, &sw, &st0](const auto& y, const auto& t0, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_t0_v(a, w, sv, sw, st0), y, t0, v); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_t0_sv) { + auto f_y_t0_sv = [](const auto& a, const auto& w, const auto& v, + const auto& sw, const auto& st0) { + return + [&a, &w, &v, &sw, &st0](const auto& y, const auto& t0, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_t0_sv(a, w, v, sw, st0), y, t0, sv); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_t0_sw) { + auto f_y_t0_sw = [](const auto& a, const auto& w, const auto& v, + const auto& sv, const auto& st0) { + return + [&a, &w, &v, &sv, &st0](const auto& y, const auto& t0, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_t0_sw(a, w, v, sv, st0), y, t0, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_t0_st0) { + auto f_y_t0_st0 = [](const auto& a, const auto& w, const auto& v, + const auto& sv, const auto& sw) { + return + [&a, &w, &v, &sv, &sw](const auto& y, const auto& t0, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_t0_st0(a, w, v, sv, sw), y, t0, st0); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_3_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_3_test.cpp new file mode 100644 index 00000000000..6342b2a945c --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_3_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_w_v) { + auto f_y_w_v = [](const auto& a, const auto& t0, const auto& sv, + const auto& sw, const auto& st0) { + return + [&a, &t0, &sv, &sw, &st0](const auto& y, const auto& w, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_w_v(a, t0, sv, sw, st0), y, w, v); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_w_sv) { + auto f_y_w_sv = [](const auto& a, const auto& t0, const auto& v, + const auto& sw, const auto& st0) { + return + [&a, &t0, &v, &sw, &st0](const auto& y, const auto& w, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_w_sv(a, t0, v, sw, st0), y, w, sv); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_w_sw) { + auto f_y_w_sw = [](const auto& a, const auto& t0, const auto& v, + const auto& sv, const auto& st0) { + return + [&a, &t0, &v, &sv, &st0](const auto& y, const auto& w, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_w_sw(a, t0, v, sv, st0), y, w, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_w_st0) { + auto f_y_w_st0 = [](const auto& a, const auto& t0, const auto& v, + const auto& sv, const auto& sw) { + return + [&a, &t0, &v, &sv, &sw](const auto& y, const auto& w, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_w_st0(a, t0, v, sv, sw), y, w, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_v_sv) { + auto f_y_v_sv = [](const auto& a, const auto& t0, const auto& w, + const auto& sw, const auto& st0) { + return + [&a, &t0, &w, &sw, &st0](const auto& y, const auto& v, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_v_sv(a, t0, w, sw, st0), y, v, sv); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_4_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_4_test.cpp new file mode 100644 index 00000000000..c7462a99ca7 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_4_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_v_sw) { + auto f_y_v_sw = [](const auto& a, const auto& t0, const auto& w, + const auto& sv, const auto& st0) { + return + [&a, &t0, &w, &sv, &st0](const auto& y, const auto& v, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_v_sw(a, t0, w, sv, st0), y, v, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_v_st0) { + auto f_y_v_st0 = [](const auto& a, const auto& t0, const auto& w, + const auto& sv, const auto& sw) { + return + [&a, &t0, &w, &sv, &sw](const auto& y, const auto& v, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_v_st0(a, t0, w, sv, sw), y, v, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_sv_sw) { + auto f_y_sv_sw = [](const auto& a, const auto& t0, const auto& w, + const auto& v, const auto& st0) { + return + [&a, &t0, &w, &v, &st0](const auto& y, const auto& sv, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_sv_sw(a, t0, w, v, st0), y, sv, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_sv_st0) { + auto f_y_sv_st0 = [](const auto& a, const auto& t0, const auto& w, + const auto& v, const auto& sw) { + return + [&a, &t0, &w, &v, &sw](const auto& y, const auto& sv, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_sv_st0(a, t0, w, v, sw), y, sv, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_y_sw_st0) { + auto f_y_sw_st0 = [](const auto& a, const auto& t0, const auto& w, + const auto& v, const auto& sv) { + return + [&a, &t0, &w, &v, &sv](const auto& y, const auto& sw, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_y_sw_st0(a, t0, w, v, sv), y, sw, st0); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_5_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_5_test.cpp new file mode 100644 index 00000000000..6e26a4d9b5d --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_5_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_t0_w) { + auto f_a_t0_w = [](const auto& y, const auto& v, const auto& sv, + const auto& sw, const auto& st0) { + return + [&y, &v, &sv, &sw, &st0](const auto& a, const auto& t0, const auto& w) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_t0_w(y, v, sv, sw, st0), a, t0, w); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_t0_v) { + auto f_a_t0_v = [](const auto& y, const auto& w, const auto& sv, + const auto& sw, const auto& st0) { + return + [&y, &w, &sv, &sw, &st0](const auto& a, const auto& t0, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_t0_v(y, w, sv, sw, st0), a, t0, v); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_t0_sv) { + auto f_a_t0_sv = [](const auto& y, const auto& w, const auto& v, + const auto& sw, const auto& st0) { + return + [&y, &w, &v, &sw, &st0](const auto& a, const auto& t0, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_t0_sv(y, w, v, sw, st0), a, t0, sv); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_t0_sw) { + auto f_a_t0_sw = [](const auto& y, const auto& w, const auto& v, + const auto& sv, const auto& st0) { + return + [&y, &w, &v, &sv, &st0](const auto& a, const auto& t0, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_t0_sw(y, w, v, sv, st0), a, t0, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_t0_st0) { + auto f_a_t0_st0 = [](const auto& y, const auto& w, const auto& v, + const auto& sv, const auto& sw) { + return + [&y, &w, &v, &sv, &sw](const auto& a, const auto& t0, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_t0_st0(y, w, v, sv, sw), a, t0, st0); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_6_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_6_test.cpp new file mode 100644 index 00000000000..60e5cdc1c61 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_6_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_w_v) { + auto f_a_w_v = [](const auto& y, const auto& t0, const auto& sv, + const auto& sw, const auto& st0) { + return + [&y, &t0, &sv, &sw, &st0](const auto& a, const auto& w, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_w_v(y, t0, sv, sw, st0), a, w, v); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_w_sv) { + auto f_a_w_sv = [](const auto& y, const auto& t0, const auto& v, + const auto& sw, const auto& st0) { + return + [&y, &t0, &v, &sw, &st0](const auto& a, const auto& w, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_w_sv(y, t0, v, sw, st0), a, w, sv); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_w_sw) { + auto f_a_w_sw = [](const auto& y, const auto& t0, const auto& v, + const auto& sv, const auto& st0) { + return + [&y, &t0, &v, &sv, &st0](const auto& a, const auto& w, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_w_sw(y, t0, v, sv, st0), a, w, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_w_st0) { + auto f_a_w_st0 = [](const auto& y, const auto& t0, const auto& v, + const auto& sv, const auto& sw) { + return + [&y, &t0, &v, &sv, &sw](const auto& a, const auto& w, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_w_st0(y, t0, v, sv, sw), a, w, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_v_sv) { + auto f_a_v_sv = [](const auto& y, const auto& t0, const auto& w, + const auto& sw, const auto& st0) { + return + [&y, &t0, &w, &sw, &st0](const auto& a, const auto& v, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_v_sv(y, t0, w, sw, st0), a, v, sv); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_7_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_7_test.cpp new file mode 100644 index 00000000000..d74690692f0 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_7_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_v_sw) { + auto f_a_v_sw = [](const auto& y, const auto& t0, const auto& w, + const auto& sv, const auto& st0) { + return + [&y, &t0, &w, &sv, &st0](const auto& a, const auto& v, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_v_sw(y, t0, w, sv, st0), a, v, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_v_st0) { + auto f_a_v_st0 = [](const auto& y, const auto& t0, const auto& w, + const auto& sv, const auto& sw) { + return + [&y, &t0, &w, &sv, &sw](const auto& a, const auto& v, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_v_st0(y, t0, w, sv, sw), a, v, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_sv_sw) { + auto f_a_sv_sw = [](const auto& y, const auto& t0, const auto& w, + const auto& v, const auto& st0) { + return + [&y, &t0, &w, &v, &st0](const auto& a, const auto& sv, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_sv_sw(y, t0, w, v, st0), a, sv, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_sv_st0) { + auto f_a_sv_st0 = [](const auto& y, const auto& t0, const auto& w, + const auto& v, const auto& sw) { + return + [&y, &t0, &w, &v, &sw](const auto& a, const auto& sv, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_sv_st0(y, t0, w, v, sw), a, sv, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_a_sw_st0) { + auto f_a_sw_st0 = [](const auto& y, const auto& t0, const auto& w, + const auto& v, const auto& sv) { + return + [&y, &t0, &w, &v, &sv](const auto& a, const auto& sw, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_a_sw_st0(y, t0, w, v, sv), a, sw, st0); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_8_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_8_test.cpp new file mode 100644 index 00000000000..ace6b36c254 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_8_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_t0_w_v) { + auto f_t0_w_v = [](const auto& y, const auto& a, const auto& sv, + const auto& sw, const auto& st0) { + return + [&y, &a, &sv, &sw, &st0](const auto& t0, const auto& w, const auto& v) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_t0_w_v(y, a, sv, sw, st0), t0, w, v); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_t0_w_sv) { + auto f_t0_w_sv = [](const auto& y, const auto& a, const auto& v, + const auto& sw, const auto& st0) { + return + [&y, &a, &v, &sw, &st0](const auto& t0, const auto& w, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_t0_w_sv(y, a, v, sw, st0), t0, w, sv); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_t0_w_sw) { + auto f_t0_w_sw = [](const auto& y, const auto& a, const auto& v, + const auto& sv, const auto& st0) { + return + [&y, &a, &v, &sv, &st0](const auto& t0, const auto& w, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_t0_w_sw(y, a, v, sv, st0), t0, w, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_t0_w_st0) { + auto f_t0_w_st0 = [](const auto& y, const auto& a, const auto& v, + const auto& sv, const auto& sw) { + return + [&y, &a, &v, &sv, &sw](const auto& t0, const auto& w, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_t0_w_st0(y, a, v, sv, sw), t0, w, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_t0_v_sv) { + auto f_t0_v_sv = [](const auto& y, const auto& a, const auto& w, + const auto& sw, const auto& st0) { + return + [&y, &a, &w, &sw, &st0](const auto& t0, const auto& v, const auto& sv) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_t0_v_sv(y, a, w, sw, st0), t0, v, sv); +} diff --git a/test/unit/math/mix/prob/wiener_full_lpdf_9_test.cpp b/test/unit/math/mix/prob/wiener_full_lpdf_9_test.cpp new file mode 100644 index 00000000000..91faeb42eba --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lpdf_9_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lpdf_t0_v_sw) { + auto f_t0_v_sw = [](const auto& y, const auto& a, const auto& w, + const auto& sv, const auto& st0) { + return + [&y, &a, &w, &sv, &st0](const auto& t0, const auto& v, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_t0_v_sw(y, a, w, sv, st0), t0, v, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_t0_v_st0) { + auto f_t0_v_st0 = [](const auto& y, const auto& a, const auto& w, + const auto& sv, const auto& sw) { + return + [&y, &a, &w, &sv, &sw](const auto& t0, const auto& v, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_t0_v_st0(y, a, w, sv, sw), t0, v, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_t0_sv_sw) { + auto f_t0_sv_sw = [](const auto& y, const auto& a, const auto& w, + const auto& v, const auto& st0) { + return + [&y, &a, &w, &v, &st0](const auto& t0, const auto& sv, const auto& sw) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_t0_sv_sw(y, a, w, v, st0), t0, sv, sw); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_t0_sv_st0) { + auto f_t0_sv_st0 = [](const auto& y, const auto& a, const auto& w, + const auto& v, const auto& sw) { + return + [&y, &a, &w, &v, &sw](const auto& t0, const auto& sv, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_t0_sv_st0(y, a, w, v, sw), t0, sv, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lpdf_t0_sw_st0) { + auto f_t0_sw_st0 = [](const auto& y, const auto& a, const auto& w, + const auto& v, const auto& sv) { + return + [&y, &a, &w, &v, &sv](const auto& t0, const auto& sw, const auto& st0) { + return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + }; + }; + stan::test::expect_ad(f_t0_sw_st0(y, a, w, v, sv), t0, sw, st0); +} diff --git a/test/unit/math/prim/functor/hcubature_test.cpp b/test/unit/math/prim/functor/hcubature_test.cpp index d8f09925581..ae19259e4be 100644 --- a/test/unit/math/prim/functor/hcubature_test.cpp +++ b/test/unit/math/prim/functor/hcubature_test.cpp @@ -6,91 +6,57 @@ namespace hcubature_test { -struct my_params { - long double x; - double y; -}; - -template -stan::return_type_t f1(const T_x& x, const T_p& p) { - using T_x_ref = stan::ref_type_t; - T_x_ref x_ref = x; - stan::scalar_seq_view x_vec(x_ref); - return cos(x_vec[0]); +template +inline auto f1(const T_x& x) { + return stan::math::sum(stan::math::cos(x)); } -template -stan::return_type_t f2(const T_x& x, const T_p& p) { - using T_x_ref = stan::ref_type_t; - T_x_ref x_ref = x; - stan::scalar_seq_view x_vec(x_ref); - return cos(x_vec[0]) * cos(x_vec[1]); +template +inline auto f2(const T_x& x) { + return stan::math::prod(stan::math::cos(x)); } -template -stan::return_type_t f3(const T_x& x, const T_p& p) { - using T_x_ref = stan::ref_type_t; - T_x_ref x_ref = x; - stan::scalar_seq_view x_vec(x_ref); - my_params* pars = static_cast(p); - long double radius = (pars->x); - double result; - if (std::pow(x_vec[0], 2) + std::pow(x_vec[1], 2) < std::pow(radius, 2)) { - result = 1; - } else { - result = 0; - } - return result; +template +inline auto f3(const T_x& x, double radius) { + using stan::math::square; + return stan::return_type_t(stan::math::sum(square(x)) < square(radius)); } -template -stan::return_type_t f4(const T_x& x, const T_p& p) { - using T_x_ref = stan::ref_type_t; - using namespace stan::math; - T_x_ref x_ref = x; - stan::scalar_seq_view x_vec(x_ref); - my_params* pars = static_cast(p); - double sigma = (pars->y); - double numerator = std::pow(x_vec[0] - 0.5, 2) + std::pow(x_vec[1] - 0.5, 2); - return std::pow(TWO_OVER_SQRT_PI / (2.0 * sigma), 2) - * exp(-numerator / std::pow(sigma, 2)); +template +inline auto f4(const T_x& x, double sigma) { + using stan::math::as_array_or_scalar; + using stan::math::square; + using stan::math::sum; + using stan::math::TWO_OVER_SQRT_PI; + + auto numerator = sum(square(as_array_or_scalar(x) - 0.5)); + return square(TWO_OVER_SQRT_PI / (2.0 * sigma)) + * exp(-numerator / square(sigma)); } -template -stan::return_type_t f5(const T_x& x, const T_p& p) { - using T_x_ref = stan::ref_type_t; - using namespace stan::math; - T_x_ref x_ref = x; - stan::scalar_seq_view x_vec(x_ref); - double val = std::pow((1 - x_vec[0]) / x_vec[0], 2) - + std::pow((1 - x_vec[1]) / x_vec[1], 2) - + std::pow((1 - x_vec[2]) / x_vec[2], 2); - double scale = TWO_OVER_SQRT_PI / std::pow(x_vec[0], 2) * TWO_OVER_SQRT_PI - / std::pow(x_vec[1], 2) * TWO_OVER_SQRT_PI - / std::pow(x_vec[2], 2); +template +inline auto f5(const T_x& x) { + using stan::math::prod; + using stan::math::square; + using stan::math::TWO_OVER_SQRT_PI; + + const auto& x_arr = stan::math::as_array_or_scalar(x); + auto val = stan::math::sum(square((1 - x_arr) / x_arr)); + auto scale = prod(TWO_OVER_SQRT_PI / square(x_arr)); return exp(-val) * scale; } -template -stan::return_type_t f6(const T_x& x, const T_p& p) { - using T_x_ref = stan::ref_type_t; - T_x_ref x_ref = x; - stan::scalar_seq_view x_vec(x_ref); - return 2 * x_vec[0] * 2 * x_vec[1] * 2 * x_vec[2]; +template +inline auto f6(const T_x& x) { + const auto& x_arr = stan::math::as_array_or_scalar(x); + return stan::math::prod(2 * x_arr); } -template -stan::return_type_t f7(const T_x& x, const T_p& p) { - using T_x_ref = stan::ref_type_t; - T_x_ref x_ref = x; - stan::scalar_seq_view x_vec(x_ref); - my_params* pars = static_cast(p); - double a = (pars->y); - double result = (a / (a + 1) * std::pow((a + 1) / (a + x_vec[0]), 2)) - * (a / (a + 1) * std::pow((a + 1) / (a + x_vec[1]), 2)) - * (a / (a + 1) * std::pow((a + 1) / (a + x_vec[2]), 2)) - * (a / (a + 1) * std::pow((a + 1) / (a + x_vec[3]), 2)); - return result; +template +inline auto f7(const T_x& x, double a) { + const auto& x_arr = stan::math::as_array_or_scalar(x); + return stan::math::prod(a / (a + 1) + * stan::math::square((a + 1) / (a + x_arr))); } } // namespace hcubature_test @@ -132,11 +98,11 @@ stan::return_type_t f7(const T_x& x, const T_p& p) { * @param val correct value of integral */ -template -void test_integration(const F& f, hcubature_test::my_params* pars, int dim, - std::vector a, std::vector b, int maxEval, - double reqAbsError, std::vector reqRelError, - double val) { +template +void test_integration(const F& f, const ArgsTupleT& pars, int dim, const T_a& a, + const T_b& b, int maxEval, double reqAbsError, + const T_relerr& reqRelError, double val) { using stan::math::hcubature; for (auto tolerance : reqRelError) { @@ -151,48 +117,49 @@ 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}; - hcubature_test::my_params pars = {}; - test_integration(hcubature_test::f1, void*>, &pars, dim, - a, b, 6000, 0.0, reqRelError, 0.841471); + const Eigen::VectorXd a{{0.0}}; + const Eigen::VectorXd b{{1.0}}; + const Eigen::VectorXd reqRelError{{1e-4, 1e-6, 1e-7}}; + test_integration([](auto&& x) { return hcubature_test::f1(x); }, + 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, void*>, &pars, dim, - a, b, 6000, 0.0, reqRelError, 0.7080734); - - reqRelError = {1e-4}; - pars = {0.50124145262344534123412, 0.0}; - test_integration(hcubature_test::f3, void*>, &pars, dim, - a, b, 10000, 0.0, reqRelError, 0.1972807); + const Eigen::VectorXd a_2{{0.0, 0.0}}; + const Eigen::VectorXd b_2{{1.0, 1.0}}; + test_integration([](auto&& x) { return hcubature_test::f2(x); }, + std::make_tuple(), dim, a_2, b_2, 6000, 0.0, reqRelError, + 0.7080734); + + const Eigen::VectorXd reqRelError_2{{1e-4}}; + test_integration( + [](auto&& x, auto&& radius) { return hcubature_test::f3(x, radius); }, + 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}; - pars = {0.0, 0.1}; - test_integration(hcubature_test::f4, void*>, &pars, dim, - a, b, 6000, 0.0, reqRelError, 1); + test_integration( + [](auto&& x, auto&& sigma) { return hcubature_test::f4(x, sigma); }, + 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, void*>, &pars, dim, - a, b, 6000, 0.0, reqRelError, 1.00001); - - reqRelError = {1e-4, 1e-6, 1e-8}; - test_integration(hcubature_test::f6, void*>, &pars, dim, - a, b, 6000, 0.0, reqRelError, 1); + 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([](auto&& x) { return hcubature_test::f5(x); }, + std::make_tuple(), dim, a_3, b_3, 6000, 0.0, reqRelError_3, + 1.00001); + + const Eigen::VectorXd reqRelError_4{{1e-4, 1e-6, 1e-8}}; + test_integration([](auto&& x) { return hcubature_test::f6(x); }, + 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}; - pars = {0.0, (1 + sqrt(10.0)) / 9.0}; - test_integration(hcubature_test::f7, void*>, &pars, 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([](auto&& x, auto&& a) { return hcubature_test::f7(x, a); }, + std::make_tuple((1 + sqrt(10.0)) / 9.0), dim, a_4, b_4, + 20000, 0.0, reqRelError_3, 0.999998); } diff --git a/test/unit/math/prim/prob/wiener_full_lpdf_test.cpp b/test/unit/math/prim/prob/wiener_full_lpdf_test.cpp new file mode 100644 index 00000000000..8dedb448428 --- /dev/null +++ b/test/unit/math/prim/prob/wiener_full_lpdf_test.cpp @@ -0,0 +1,464 @@ +#include +#include +#include +#include + +#include +#include + +TEST(mathPrimScalProbWienerFullScal, valid) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, + st0 = 0.1; + EXPECT_NO_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, st0)); + rt = 5; + a = 1; + v = 1; + w = 0.5; + t0 = 0.0; + sv = 0.0; + sw = 0.0; + st0 = 0.0; + EXPECT_NO_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, st0)); +} + +// rt +TEST(mathPrimScalProbWienerFullScal, invalid_rt) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(0, a, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(-1, a, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(INFTY, a, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(-INFTY, a, t0, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lpdf(NAN, a, t0, w, v, sv, sw, st0), std::domain_error); +} + +// a +TEST(mathPrimScalProbWienerFullScal, invalid_a) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, 0, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, -1, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, INFTY, t0, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, -INFTY, t0, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, NAN, t0, w, v, sv, sw, st0), std::domain_error); +} + +// v +TEST(mathPrimScalProbWienerFullScal, invalid_v) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, + st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, INFTY, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, -INFTY, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, NAN, sv, sw, st0), std::domain_error); +} + +// w +TEST(mathPrimScalProbWienerFullScal, invalid_w) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, a, t0, -0.1, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, 0, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, 1, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, 1.1, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, INFTY, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, -INFTY, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, NAN, v, sv, sw, st0), std::domain_error); +} + +// t0 +TEST(mathPrimScalProbWienerFullScal, invalid_t0) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, a, 2, w, v, sv, sw, st0), + std::domain_error); // rt must be greater than t0 + EXPECT_THROW(wiener_lpdf(rt, a, -1, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, INFTY, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, -INFTY, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, NAN, w, v, sv, sw, st0), std::domain_error); +} + +// sv +TEST(mathPrimScalProbWienerFullScal, invalid_sv) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, -1, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, INFTY, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, -INFTY, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, NAN, sw, st0), std::domain_error); +} + +// sw +TEST(mathPrimScalProbWienerFullScal, invalid_sw) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, -1, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, 0.8, v, sv, 0.5, st0), + std::domain_error); // sw must be smaller than 2*(1-w) + EXPECT_THROW(wiener_lpdf(rt, a, t0, 0.3, v, sv, 0.7, st0), + std::domain_error); // sw must be smaller than 2*w + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, INFTY, st0), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, -INFTY, st0), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, NAN, st0), std::domain_error); +} + +// st0 +TEST(mathPrimScalProbWienerFullScal, invalid_st0) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2; + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, -1), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, INFTY), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, -INFTY), std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, NAN), std::domain_error); +} + +TEST(mathPrimScalProbWienerFullPrecScal, valid) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, + st0 = 0.1; + EXPECT_NO_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, st0, 1e-4)); + rt = 5; + a = 1; + v = 1; + w = 0.5; + t0 = 0.0; + sv = 0.0; + sw = 0.0; + st0 = 0.0; + EXPECT_NO_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, st0, 1e-4)); +} + +// rt +TEST(mathPrimScalProbWienerFullPrecScal, invalid_rt) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(0, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(-1, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(INFTY, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(-INFTY, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(NAN, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// a +TEST(mathPrimScalProbWienerFullPrecScal, invalid_a) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, 0, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, -1, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, INFTY, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, -INFTY, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, NAN, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// v +TEST(mathPrimScalProbWienerFullPrecScal, invalid_v) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, + st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, INFTY, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, -INFTY, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, NAN, sv, sw, st0, 1e-4), + std::domain_error); +} + +// w +TEST(mathPrimScalProbWienerFullPrecScal, invalid_w) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, a, t0, -0.1, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, 0, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, 1, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, 1.1, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, INFTY, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, -INFTY, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, NAN, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// t0 +TEST(mathPrimScalProbWienerFullPrecScal, invalid_t0) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, a, 2, w, v, sv, sw, st0, 1e-4), + std::domain_error); // rt must be greater than t0 + EXPECT_THROW(wiener_lpdf(rt, a, -1, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, INFTY, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, -INFTY, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, NAN, w, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// sv +TEST(mathPrimScalProbWienerFullPrecScal, invalid_sv) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, -1, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, INFTY, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, -INFTY, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, NAN, sw, st0, 1e-4), + std::domain_error); +} + +// sw +TEST(mathPrimScalProbWienerFullPrecScal, invalid_sw) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, -1, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, 0.8, v, sv, 0.5, st0, 1e-4), + std::domain_error); // sw must be smaller than 2*(1-w) + EXPECT_THROW(wiener_lpdf(rt, a, t0, 0.3, v, sv, 0.7, st0, 1e-4), + std::domain_error); // sw must be smaller than 2*w + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, INFTY, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, -INFTY, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, NAN, st0, 1e-4), + std::domain_error); +} + +// st0 +TEST(mathPrimScalProbWienerFullPrecScal, invalid_st0) { + using stan::math::INFTY; + using stan::math::wiener_lpdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2; + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, -1, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, INFTY, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, -INFTY, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lpdf(rt, a, t0, w, v, sv, sw, NAN, 1e-4), + std::domain_error); +} + +TEST(mathPrimCorrectValues, wiener_lpdf) { + // Test concrete values. True values are computed in R using the R-package + // WienR and the function WienerPDF. + std::vector y_vec = {2, 3, 4, 5, 6, 7, 8, 8.85, 8.9, 9, 1}; + std::vector a_vec + = {2.0, 2.0, 10.0, 4.0, 10.0, 1.0, 3.0, 1.7, 2.4, 11.0, 1.5}; + std::vector v_vec + = {2.0, 2.0, 4.0, 3.0, -3.0, 1.0, -1.0, -7.3, -4.9, 4.5, 3}; + std::vector w_vec = {.1, 0.5, .8, 0.7, .1, .9, .7, .92, .9, .12, 0.5}; + std::vector t0_vec + = {1e-9, 0.01, .01, .01, .01, .01, .01, .01, .01, .01, 0.1}; + std::vector sv_vec = {0, 0.2, 0, 0, .2, .2, 0, .7, 0, .7, 0.5}; + std::vector sw_vec = {0, 0, .1, 0, .1, 0, .1, .01, 0, .1, 0.2}; + std::vector st0_vec + = {0, 0, 0, 0.007, 0, .007, .007, .009, .009, .009, 0}; + + std::vector true_dens + = {-4.28564747866615, -7.52379235146909, -26.1551056209248, + -22.1939134892089, -50.0587553794834, -37.2817263586318, + -10.5428662079438, -61.5915905674246, -117.238967959795, + -12.5788594249676, -3.1448097740735}; + std::vector true_grad_y + = {-3.22509339523307, -2.91155058614589, -8.21331631900955, + -4.82948967379739, -1.50069056428102, -5.25831601347426, + -1.04831896413742, -2.67457492096193, -12.8617364931501, + -1.12047317491985, -5.68799957241344}; + std::vector true_grad_a = { + 3.25018678924105, 3.59980430191399, 0.876602303160642, 1.2215517888504, + -3.02928674030948, 67.0322498959921, 1.95334514374631, 16.4642201959135, + 5.02038145619773, 0.688439187670968, 2.63200041459657}; + std::vector true_grad_t0 + = {3.22509339523307, 2.91155058614589, 8.21331631900955, 4.82948967379739, + 1.50069056428102, 5.25831601347426, 1.04831896413742, 2.67457492096193, + 12.8617364931501, 1.12047317491985, 5.68799957241344}; + std::vector true_grad_w + = {5.67120184517318, -3.64396221090076, -38.7775057146792, + -14.1837930137393, -34.5869239580708, -10.4535345681946, + 0.679597983582904, -9.93144540834201, 2.09117200953597, + -6.0858540417876, -3.74870310978083}; + std::vector true_grad_v = { + -2.199999998, -4.44801714898178, -13.6940602985224, -13.7593709622169, + 21.5540563802381, -5.38233555673517, 8.88475440789056, 12.1280680728793, + 43.7785246930371, -5.68143495684294, -1.57639220567218}; + std::vector true_grad_sv = {0, + 3.42285198319565, + 0, + 0, + 91.9551438876654, + 4.70180879974639, + 0, + 101.80250964211, + 0, + 21.4332628706595, + 0.877556017134384}; + std::vector true_grad_sw = {0, + 0, + 10.1052188867058, + 0, + 8.72398, + 0, + -0.122807217815892, + -0.0506322723373748, + 0, + -0.0704990526706635, + 0.0827817310725268}; + std::vector true_grad_st0 = {0, + 0, + 0, + 2.42836139121338, + 0, + 2.64529825657625, + 0.524800556172613, + 1.34278261179603, + 6.55490874737353, + 0.561295838843035, + 0}; + + using stan::math::var; + double err_tol_dens = 1e-6; + double err_tol = 1e-4; + for (int i = 0; i < y_vec.size(); i++) { + var y = y_vec[i]; + var a = a_vec[i]; + var t0 = t0_vec[i]; + var w = w_vec[i]; + var v = v_vec[i]; + var sv = sv_vec[i]; + var sw = sw_vec[i]; + var st0 = st0_vec[i]; + var dens = stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0); + dens.grad(); + EXPECT_NEAR(dens.val(), true_dens[i], err_tol_dens); + EXPECT_NEAR(y.adj(), true_grad_y[i], err_tol); + EXPECT_NEAR(a.adj(), true_grad_a[i], err_tol); + EXPECT_NEAR(t0.adj(), true_grad_t0[i], err_tol); + EXPECT_NEAR(w.adj(), true_grad_w[i], err_tol); + EXPECT_NEAR(v.adj(), true_grad_v[i], err_tol); + EXPECT_NEAR(sv.adj(), true_grad_sv[i], err_tol); + EXPECT_NEAR(sw.adj(), true_grad_sw[i], err_tol); + EXPECT_NEAR(st0.adj(), true_grad_st0[i], err_tol); + } +} + +TEST(mathPrimCorrectValuesFourParamModel, wiener_lpdf) { + // Test concrete values. True values are computed in R using the R-package + // WienR and the function WienerPDF. + std::vector y_vec = {2.0}; + std::vector a_vec = {2.0}; + std::vector v_vec = {2.0}; + std::vector w_vec = {.1}; + std::vector t0_vec = {1e-9}; + + std::vector true_dens = {-4.28564747866615}; + std::vector true_grad_y = {-3.22509339523307}; + std::vector true_grad_a = {3.25018678924105}; + std::vector true_grad_t0 = {3.22509339523307}; + std::vector true_grad_w = {5.67120184517318}; + std::vector true_grad_v = {-2.199999998}; + + using stan::math::var; + double err_tol_dens = 1e-6; + double err_tol = 1e-4; + for (int i = 0; i < y_vec.size(); i++) { + var y = y_vec[i]; + var a = a_vec[i]; + var t0 = t0_vec[i]; + var w = w_vec[i]; + var v = v_vec[i]; + var dens = stan::math::wiener_lpdf(y, a, t0, w, v); + dens.grad(); + EXPECT_NEAR(dens.val(), true_dens[i], err_tol_dens); + EXPECT_NEAR(y.adj(), true_grad_y[i], err_tol); + EXPECT_NEAR(a.adj(), true_grad_a[i], err_tol); + EXPECT_NEAR(t0.adj(), true_grad_t0[i], err_tol); + EXPECT_NEAR(w.adj(), true_grad_w[i], err_tol); + EXPECT_NEAR(v.adj(), true_grad_v[i], err_tol); + } +} + +TEST(mathPrimCorrectValuesFiveParameterModel, wiener_lpdf) { + // Test concrete values. True values are computed in R using the R-package + // WienR and the function WienerPDF. + std::vector y_vec = {3.0}; + std::vector a_vec = {2.0}; + std::vector v_vec = {2.0}; + std::vector w_vec = {0.5}; + std::vector t0_vec = {0.01}; + std::vector sv_vec = {0.2}; + + std::vector true_dens = {-7.52379235146909}; + std::vector true_grad_y = {-2.91155058614589}; + std::vector true_grad_a = {3.59980430191399}; + std::vector true_grad_t0 = {2.91155058614589}; + std::vector true_grad_w = {-3.64396221090076}; + std::vector true_grad_v = {-4.44801714898178}; + std::vector true_grad_sv = {3.42285198319565}; + + using stan::math::var; + double err_tol_dens = 1e-6; + double err_tol = 1e-4; + for (int i = 0; i < y_vec.size(); i++) { + var y = y_vec[i]; + var a = a_vec[i]; + var t0 = t0_vec[i]; + var w = w_vec[i]; + var v = v_vec[i]; + var sv = sv_vec[i]; + var dens = stan::math::wiener_lpdf(y, a, t0, w, v, sv); + dens.grad(); + EXPECT_NEAR(dens.val(), true_dens[i], err_tol_dens); + EXPECT_NEAR(y.adj(), true_grad_y[i], err_tol); + EXPECT_NEAR(a.adj(), true_grad_a[i], err_tol); + EXPECT_NEAR(t0.adj(), true_grad_t0[i], err_tol); + EXPECT_NEAR(w.adj(), true_grad_w[i], err_tol); + EXPECT_NEAR(v.adj(), true_grad_v[i], err_tol); + EXPECT_NEAR(sv.adj(), true_grad_sv[i], err_tol); + } +} diff --git a/test/unit/math/rev/prob/wiener_full_test.cpp b/test/unit/math/rev/prob/wiener_full_test.cpp new file mode 100644 index 00000000000..bc36395489b --- /dev/null +++ b/test/unit/math/rev/prob/wiener_full_test.cpp @@ -0,0 +1,278 @@ +#include +#include +#include +#include + +#include +#include + +// CHECK THAT ALL VALID SCALAR TYPES ARE ACCEPTED +template +void check_scalar_types(F& f, double value, double res, double deriv) { + // - f: Function with a single parameter exposed, all others + // have to be scalars + // - value: value to be used for the parameter + // - res: expected result of calling `f` with `value` + // - deriv: expected result of partial of f with respect to + // the parameter in `value` + using stan::math::var; + double err_tol = 2e-4; + + // type double + EXPECT_NEAR(f(value), res, err_tol); + + // type var with derivative + var value_var = value; + var result_var = f(value_var); + result_var.grad(); + EXPECT_NEAR(value_of(result_var), res, err_tol); + EXPECT_NEAR(value_var.adj(), deriv, err_tol); +} + +TEST(ProbWienerFull, wiener_full_all_scalar) { + // tests all parameter types individually, with other + // parameters set to double + using stan::math::wiener_lpdf; + + std::vector rt{1, 1, 1, 1, 1, 1, 1, 1, 1}; + std::vector a{1, 1, 1, 1, 1, 1, 1, 1, 1}; + std::vector v{1, 1, 1, 1, 1, 1, 1, 1, 1}; + std::vector w{0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5}; + std::vector t0{0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0}; + std::vector sv{0.1, 0.1, 0.1, 0, 0.1, 0, 0, 0, 0}; + std::vector sw{0.1, 0.1, 0, 0.1, 0, 0.1, 0, 0, 0}; + std::vector st0{0.1, 0, 0.1, 0.1, 0, 0, 0.1, 0, 0}; + + std::vector result{-2.426204, -2.710348, -2.422505, + -2.422795, -2.70665, -2.706812, + -2.419095, -2.703112, -3.790072}; + std::vector drt{-5.437337, -5.436798, -5.437332, -5.434799, -5.436791, + -5.434801, -5.434802, -5.434802, -5.434802}; + std::vector da{5.857318, 6.395031, 5.856484, 5.858549, 6.394195, + 6.396512, 5.857722, 6.395684, 8.369604}; + std::vector dv{-0.2428443, -0.2967977, -0.2436664, + -0.2446629, -0.297619, -0.2991696, + -0.2454931, -0.3, -0.5}; + std::vector dw{-0.9891369, -0.9887305, -0.9973428, + -0.9915453, -0.9969335, -0.991674, + -0.9997794, -0.9999097, -0.9999953}; + std::vector dt0{5.437337, 5.436798, 5.437332, 5.434799, 5.436791, + 5.434801, 5.434802, 5.434802, 5.434802}; + std::vector dsv{-0.06793703, -0.07047449, -0.06797882, + 0.0, -0.07050737, 0.0, + 0.0, 0.0, 0.0}; + std::vector dsw{-0.07406705, -0.07407386, 0.0, -0.07410901, 0.0, + -0.07410686, 0.0, 0.0, 0.0}; + std::vector dst0{2.963915, 0.0, 2.963912, 2.962338, 0.0, + 0.0, 2.96234, 0.0, 0.0}; + + for (int i = 0; i < result.size(); i++) { + // rt + auto f_rt = [&](auto value) { + return wiener_lpdf(value, a[i], t0[i], w[i], v[i], sv[i], sw[i], st0[i]); + }; + check_scalar_types(f_rt, rt[i], result[i], drt[i]); + // a + auto f_a = [&](auto value) { + return wiener_lpdf(rt[i], value, t0[i], w[i], v[i], sv[i], sw[i], st0[i]); + }; + check_scalar_types(f_a, a[i], result[i], da[i]); + // v + auto f_v = [&](auto value) { + return wiener_lpdf(rt[i], a[i], t0[i], w[i], value, sv[i], sw[i], st0[i]); + }; + check_scalar_types(f_v, v[i], result[i], dv[i]); + // w + auto f_w = [&](auto value) { + return wiener_lpdf(rt[i], a[i], t0[i], value, v[i], sv[i], sw[i], st0[i]); + }; + check_scalar_types(f_w, w[i], result[i], dw[i]); + // t0 + auto f_t0 = [&](auto value) { + return wiener_lpdf(rt[i], a[i], value, w[i], v[i], sv[i], sw[i], st0[i]); + }; + check_scalar_types(f_t0, t0[i], result[i], dt0[i]); + // sv + auto f_sv = [&](auto value) { + return wiener_lpdf(rt[i], a[i], t0[i], w[i], v[i], value, sw[i], st0[i]); + }; + check_scalar_types(f_sv, sv[i], result[i], dsv[i]); + // sw + auto f_sw = [&](auto value) { + return wiener_lpdf(rt[i], a[i], t0[i], w[i], v[i], sv[i], value, st0[i]); + }; + check_scalar_types(f_sw, sw[i], result[i], dsw[i]); + // st0 + auto f_st0 = [&](auto value) { + return wiener_lpdf(rt[i], a[i], t0[i], w[i], v[i], sv[i], sw[i], value); + }; + check_scalar_types(f_st0, st0[i], result[i], dst0[i]); + } +} + +TEST(ProbWienerFullPrec, wiener_full_prec_all_scalar) { + // tests all parameter types individually, with other parameters + // set to double + using stan::math::wiener_lpdf; + double err_tol = 2e-6; + + std::vector rt{1, 1, 1, 1, 1, 1, 1, 1, 1}; + std::vector a{1, 1, 1, 1, 1, 1, 1, 1, 1}; + std::vector v{1, 1, 1, 1, 1, 1, 1, 1, 1}; + std::vector w{0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5}; + std::vector t0{0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0}; + std::vector sv{0.1, 0.1, 0.1, 0, 0.1, 0, 0, 0, 0}; + std::vector sw{0.1, 0.1, 0, 0.1, 0, 0.1, 0, 0, 0}; + std::vector st0{0.1, 0, 0.1, 0.1, 0, 0, 0.1, 0, 0}; + + std::vector result{-2.426204, -2.710348, -2.422505, + -2.422795, -2.70665, -2.706812, + -2.419095, -2.703112, -3.790072}; + std::vector drt{-5.437337, -5.436798, -5.437332, -5.434799, -5.436791, + -5.434801, -5.434802, -5.434802, -5.434802}; + std::vector da{5.857318, 6.395031, 5.856484, 5.858549, 6.394195, + 6.396512, 5.857722, 6.395684, 8.369604}; + std::vector dv{-0.2428443, -0.2967977, -0.2436664, + -0.2446629, -0.297619, -0.2991696, + -0.2454931, -0.3, -0.5}; + std::vector dw{-0.9891369, -0.9887305, -0.9973428, + -0.9915453, -0.9969335, -0.991674, + -0.9997794, -0.9999097, -0.9999953}; + std::vector dt0{5.437337, 5.436798, 5.437332, 5.434799, 5.436791, + 5.434801, 5.434802, 5.434802, 5.434802}; + std::vector dsv{-0.06793703, -0.07047449, -0.06797882, + 0.0, -0.07050737, 0.0, + 0.0, 0.0, 0.0}; + std::vector dsw{-0.07406705, -0.07407386, 0.0, -0.07410901, 0.0, + -0.07410686, 0.0, 0.0, 0.0}; + std::vector dst0{2.963915, 0.0, 2.963912, 2.962338, 0.0, + 0.0, 2.96234, 0.0, 0.0}; + + for (int i = 0; i < result.size(); i++) { + // rt + auto f_rt = [&](auto value) { + return wiener_lpdf(value, a[i], t0[i], w[i], v[i], sv[i], sw[i], st0[i], + 1e-6); + }; + check_scalar_types(f_rt, rt[i], result[i], drt[i]); + // a + auto f_a = [&](auto value) { + return wiener_lpdf(rt[i], value, t0[i], w[i], v[i], sv[i], sw[i], st0[i], + 1e-6); + }; + check_scalar_types(f_a, a[i], result[i], da[i]); + // v + auto f_v = [&](auto value) { + return wiener_lpdf(rt[i], a[i], t0[i], w[i], value, sv[i], sw[i], st0[i], + 1e-6); + }; + check_scalar_types(f_v, v[i], result[i], dv[i]); + // w + auto f_w = [&](auto value) { + return wiener_lpdf(rt[i], a[i], t0[i], value, v[i], sv[i], sw[i], st0[i], + 1e-6); + }; + check_scalar_types(f_w, w[i], result[i], dw[i]); + // t0 + auto f_t0 = [&](auto value) { + return wiener_lpdf(rt[i], a[i], value, w[i], v[i], sv[i], sw[i], st0[i], + 1e-6); + }; + check_scalar_types(f_t0, t0[i], result[i], dt0[i]); + // sv + auto f_sv = [&](auto value) { + return wiener_lpdf(rt[i], a[i], t0[i], w[i], v[i], value, sw[i], st0[i], + 1e-6); + }; + check_scalar_types(f_sv, sv[i], result[i], dsv[i]); + // sw + auto f_sw = [&](auto value) { + return wiener_lpdf(rt[i], a[i], t0[i], w[i], v[i], sv[i], value, st0[i], + 1e-6); + }; + check_scalar_types(f_sw, sw[i], result[i], dsw[i]); + // st0 + auto f_st0 = [&](auto value) { + return wiener_lpdf(rt[i], a[i], t0[i], w[i], v[i], sv[i], sw[i], value, + 1e-6); + }; + check_scalar_types(f_st0, st0[i], result[i], dst0[i]); + } +} + +// CHECK THAT ALL VALID Vector TYPES ARE ACCEPTED +template +void check_vector_types(F& f, std::vector value, double res) { + // - f: Function where all inputs are vectors + // - value: value to be used for the parameter + // - res: expected result of calling `f` with `value` + // - deriv: expected result of partial of f with respect to + // the parameter in `value` + using stan::math::var; + double err_tol = 2e-4; + + // type double + EXPECT_NEAR(f(value), res, err_tol); + + // type var with derivative + var result_var = f(value); + result_var.grad(); + EXPECT_NEAR(value_of(result_var), res, err_tol); +} + +TEST(ProbWienerFull, wiener_full_all_vector) { + // tests all parameter types individually, with other + // parameters set to std::vector + using stan::math::wiener_lpdf; + + std::vector rt{1, 1, 1, 1, 1, 1, 1, 1, 1}; + std::vector a{1, 1, 1, 1, 1, 1, 1, 1, 1}; + std::vector v{1, 1, 1, 1, 1, 1, 1, 1, 1}; + std::vector w{0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5}; + std::vector t0{0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0}; + std::vector sv{0.1, 0.1, 0.1, 0, 0.1, 0, 0, 0, 0}; + std::vector sw{0.1, 0.1, 0, 0.1, 0, 0.1, 0, 0, 0}; + std::vector st0{0.1, 0, 0.1, 0.1, 0, 0, 0.1, 0, 0}; + + double result{-24.307593}; + + // rt + auto f_rt = [&](auto value) { + return wiener_lpdf(value, a, t0, w, v, sv, sw, st0); + }; + check_vector_types(f_rt, rt, result); + // a + auto f_a = [&](auto value) { + return wiener_lpdf(rt, value, t0, w, v, sv, sw, st0); + }; + check_vector_types(f_a, a, result); + // v + auto f_v = [&](auto value) { + return wiener_lpdf(rt, a, t0, w, value, sv, sw, st0); + }; + check_vector_types(f_v, v, result); + // w + auto f_w = [&](auto value) { + return wiener_lpdf(rt, a, t0, value, v, sv, sw, st0); + }; + check_vector_types(f_w, w, result); + // t0 + auto f_t0 = [&](auto value) { + return wiener_lpdf(rt, a, value, w, v, sv, sw, st0); + }; + check_vector_types(f_t0, t0, result); + // sv + auto f_sv = [&](auto value) { + return wiener_lpdf(rt, a, t0, w, v, value, sw, st0); + }; + check_vector_types(f_sv, sv, result); + // sw + auto f_sw = [&](auto value) { + return wiener_lpdf(rt, a, t0, w, v, sv, value, st0); + }; + check_vector_types(f_sw, sw, result); + // st0 + auto f_st0 + = [&](auto value) { return wiener_lpdf(rt, a, t0, w, v, sv, sw, value); }; + check_vector_types(f_st0, st0, result); +}