diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index be55ca6a8ef..664f97c6897 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -311,7 +311,11 @@ #include #include #include +#include +#include #include +#include +#include #include #include #include diff --git a/stan/math/prim/prob/wiener4_lccdf.hpp b/stan/math/prim/prob/wiener4_lccdf.hpp new file mode 100644 index 00000000000..115fffe8d3d --- /dev/null +++ b/stan/math/prim/prob/wiener4_lccdf.hpp @@ -0,0 +1,372 @@ +#ifndef STAN_MATH_PRIM_PROB_WIENER4_LCCDF_HPP +#define STAN_MATH_PRIM_PROB_WIENER4_LCCDF_HPP + +#include + +namespace stan { +namespace math { +namespace internal { + +/** + * Log of probability of reaching the upper bound in diffusion process + * + * @tparam T_a type of boundary + * @tparam T_w type of relative starting point + * @tparam T_v type of drift rate + * + * @param a The boundary separation + * @param w The relative starting point + * @param v The drift rate + * @return log probability to reach the upper bound + */ +template +inline auto log_wiener_prob_hit_upper(const T_a& a, const T_v& v, + const T_w& w) { + using ret_t = return_type_t; + const auto neg_v = -v; + const auto one_m_w = 1.0 - w; + if (fabs(v) == 0.0) { + return ret_t(log1m(one_m_w)); + } + const auto exponent = 2.0 * v * a * w; + // This branch is for numeric stability + if (exponent < 0) { + return ret_t(log1m_exp(exponent) + - log_diff_exp(2.0 * neg_v * a * one_m_w, exponent)); + } else { + return ret_t(log1m_exp(-exponent) - log1m_exp(2.0 * neg_v * a)); + } +} + +/** + * Calculate parts of the partial derivatives for wiener_prob_grad_a and + * wiener_prob_grad_v (on log-scale) + * + * @tparam T_a type of boundary + * @tparam T_w type of relative starting point + * @tparam T_v type of drift rate + * + * @param a The boundary separation + * @param w The relative starting point + * @param v The drift rate + * @return 'ans' term + */ +template +inline auto wiener_prob_derivative_term(const T_a& a, const T_v& v, + const T_w& w) noexcept { + using ret_t = return_type_t; + const auto exponent_m1 = log1m(1.1 * 1.0e-8); + const auto neg_v = -v; + const auto one_m_w = 1 - w; + int sign_v = neg_v < 0 ? 1 : -1; + const auto two_a_neg_v = 2.0 * a * neg_v; + const auto exponent_with_1mw = sign_v * two_a_neg_v * w; + const auto exponent = sign_v * two_a_neg_v; + const auto exponent_with_w = two_a_neg_v * one_m_w; + // truncating longer calculations, for numerical stability + if (unlikely((exponent_with_1mw >= exponent_m1) + || ((exponent_with_w >= exponent_m1) && (sign_v == 1)) + || (exponent >= exponent_m1) || neg_v == 0)) { + return ret_t(-one_m_w); + } + ret_t ans; + ret_t diff_term; + const auto log_w = log(one_m_w); + if (neg_v < 0) { + ans = LOG_TWO + exponent_with_1mw - log1m_exp(exponent_with_1mw); + diff_term = log1m_exp(exponent_with_w) - log1m_exp(exponent); + } else if (neg_v > 0) { + ans = LOG_TWO - log1m_exp(exponent_with_1mw); + diff_term = log_diff_exp(exponent_with_1mw, exponent) - log1m_exp(exponent); + } + if (log_w > diff_term) { + ans = sign_v * exp(ans + log_diff_exp(log_w, diff_term)); + } else { + ans = -sign_v * exp(ans + log_diff_exp(diff_term, log_w)); + } + if (unlikely(!is_scal_finite(ans))) { + return ret_t(NEGATIVE_INFTY); + } + return ans; +} + +/** + * Calculate wiener4 ccdf (natural-scale) + * + * @param y The reaction time in seconds + * @param a The boundary separation + * @param v The relative starting point + * @param w The drift rate + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return ccdf + */ +template +inline auto wiener4_ccdf(const T_y& y, const T_a& a, const T_v& v, const T_w& w, + T_err log_err = log(1e-12)) noexcept { + const auto prob_hit_upper = exp(log_wiener_prob_hit_upper(a, v, w)); + const auto cdf + = internal::wiener4_distribution(y, a, v, w, log_err); + return prob_hit_upper - cdf; +} + +/** + * Calculate derivative of the wiener4 ccdf w.r.t. 'a' (natural-scale) + * + * @param y The reaction time in seconds + * @param a The boundary separation + * @param v The relative starting point + * @param w The drift rate + * @param cdf The CDF value + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to a + */ +template +inline auto wiener4_ccdf_grad_a(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, T_cdf&& cdf, + T_err log_err = log(1e-12)) noexcept { + using ret_t = return_type_t; + + // derivative of the wiener probability w.r.t. 'a' (on log-scale) + auto prob_grad_a = -wiener_prob_derivative_term(a, v, w) * v; + if (!is_scal_finite(prob_grad_a)) { + prob_grad_a = ret_t(NEGATIVE_INFTY); + } + const auto log_prob_hit_upper = log_wiener_prob_hit_upper(a, v, w); + const auto cdf_grad_a = wiener4_cdf_grad_a(y, a, v, w, cdf, log_err); + return prob_grad_a * exp(log_prob_hit_upper) - cdf_grad_a; +} + +/** + * Calculate derivative of the wiener4 ccdf w.r.t. 'v' (natural-scale) + * + * @param y The reaction time in seconds + * @param a The boundary separation + * @param v The relative starting point + * @param w The drift rate + * @param cdf The CDF value + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to v + */ +template +inline auto wiener4_ccdf_grad_v(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, T_cdf&& cdf, + T_err log_err = log(1e-12)) noexcept { + using ret_t = return_type_t; + const auto log_prob_hit_upper = log_wiener_prob_hit_upper(a, v, w); + // derivative of the wiener probability w.r.t. 'v' (on log-scale) + auto prob_grad_v = -wiener_prob_derivative_term(a, v, w) * a; + if (!is_scal_finite(fabs(prob_grad_v))) { + prob_grad_v = ret_t(NEGATIVE_INFTY); + } + + const auto cdf_grad_v = wiener4_cdf_grad_v(y, a, v, w, cdf, log_err); + return prob_grad_v * exp(log_prob_hit_upper) - cdf_grad_v; +} + +/** + * Calculate derivative of the wiener4 ccdf w.r.t. 'w' (natural-scale) + * + * @param y The reaction time in seconds + * @param a The boundary separation + * @param v The relative starting point + * @param w The drift rate + * @param cdf The CDF value + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to w + */ +template +inline auto wiener4_ccdf_grad_w(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, T_cdf&& cdf, + T_err log_err = log(1e-12)) noexcept { + using ret_t = return_type_t; + const auto log_prob_hit_upper = log_wiener_prob_hit_upper(a, v, w); + // derivative of the wiener probability w.r.t. 'v' (on log-scale) + const auto exponent = -sign(v) * 2.0 * v * a * w; + auto prob_grad_w + = (v != 0) ? exp(LOG_TWO + log(fabs(v)) + log(a) - log1m_exp(exponent)) + : ret_t(1 / w); + if (v > 0) { + prob_grad_w *= exp(exponent); + } + + const auto cdf_grad_w = wiener4_cdf_grad_w(y, a, v, w, cdf, log_err); + return prob_grad_w * exp(log_prob_hit_upper) - cdf_grad_w; +} + +} // namespace internal + +/** + * Log-CCDF for the 4-parameter Wiener distribution. + * See 'wiener_full_lpdf' for more comprehensive documentation. + * + * @tparam T_y type of reaction time + * @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 + * + * @param y 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 precision_derivatives Level of precision in estimation + * @return The log of the Wiener first passage time distribution with + * the specified arguments for upper boundary responses + */ +template +inline auto wiener_lccdf(const T_y& y, const T_a& a, const T_t0& t0, + const T_w& w, const T_v& v, + const double& precision_derivatives) { + using T_partials_return = partials_return_t; + using ret_t = return_type_t; + 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 internal::GradientCalc; + + 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; + + auto y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + auto a_val = to_ref(as_value_column_array_or_scalar(a_ref)); + auto v_val = to_ref(as_value_column_array_or_scalar(v_ref)); + auto w_val = to_ref(as_value_column_array_or_scalar(w_ref)); + auto t0_val = to_ref(as_value_column_array_or_scalar(t0_ref)); + + static constexpr const char* function_name = "wiener4_lccdf"; + if (size_zero(y, a, t0, w, v)) { + return ret_t(0.0); + } + + if (!include_summand::value) { + return ret_t(0.0); + } + + check_consistent_sizes(function_name, "Random variable", y, + "Boundary separation", a, "Drift rate", v, + "A-priori bias", w, "Nondecision time", t0); + 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); + + const size_t N = max_size(y, a, t0, w, v); + + 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); + 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()); + } + } + + // for precs. 1e-6, 1e-12, see Hartmann et al. (2021), Henrich et al. (2023) + const auto log_error_cdf = log(1e-6); + const auto log_error_derivative = log(precision_derivatives); + const T_partials_return log_error_absolute = log(1e-12); + T_partials_return lccdf = 0.0; + auto ops_partials + = make_partials_propagator(y_ref, a_ref, t0_ref, w_ref, v_ref); + + const double LOG_FOUR = std::log(4.0); + + // calculate distribution and partials + for (size_t i = 0; i < N; i++) { + 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 T_partials_return cdf + = internal::estimate_with_err_check<4, 0, GradientCalc::OFF, + GradientCalc::OFF>( + [](auto&&... args) { + return internal::wiener4_distribution(args...); + }, + log_error_cdf - LOG_TWO, y_value - t0_value, a_value, v_value, + w_value, log_error_absolute); + + const auto prob_hit_upper + = exp(internal::log_wiener_prob_hit_upper(a_value, v_value, w_value)); + const auto ccdf = prob_hit_upper - cdf; + const auto log_ccdf_single_value = log(ccdf); + + lccdf += log_ccdf_single_value; + + const auto new_est_err + = log_ccdf_single_value + log_error_derivative - LOG_FOUR; + + if (!is_constant_all::value || !is_constant_all::value) { + const auto deriv_y = internal::estimate_with_err_check<5, 0>( + [](auto&&... args) { + return internal::wiener5_density(args...); + }, + new_est_err, y_value - t0_value, a_value, v_value, w_value, 0.0, + log_error_absolute); + if (!is_constant_all::value) { + partials<0>(ops_partials)[i] = -deriv_y / ccdf; + } + if (!is_constant_all::value) { + partials<2>(ops_partials)[i] = deriv_y / ccdf; + } + } + if (!is_constant_all::value) { + partials<1>(ops_partials)[i] + = internal::estimate_with_err_check<5, 0>( + [](auto&&... args) { + return internal::wiener4_ccdf_grad_a(args...); + }, + new_est_err, y_value - t0_value, a_value, v_value, w_value, cdf, + log_error_absolute) + / ccdf; + } + if (!is_constant_all::value) { + partials<3>(ops_partials)[i] + = internal::estimate_with_err_check<5, 0>( + [](auto&&... args) { + return internal::wiener4_ccdf_grad_w(args...); + }, + new_est_err, y_value - t0_value, a_value, v_value, w_value, cdf, + log_error_absolute) + / ccdf; + } + if (!is_constant_all::value) { + partials<4>(ops_partials)[i] + = internal::wiener4_ccdf_grad_v(y_value - t0_value, a_value, v_value, + w_value, cdf, log_error_absolute) + / ccdf; + } + } // for loop + return ops_partials.build(lccdf); +} +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/wiener4_lcdf.hpp b/stan/math/prim/prob/wiener4_lcdf.hpp new file mode 100644 index 00000000000..b80824fa16b --- /dev/null +++ b/stan/math/prim/prob/wiener4_lcdf.hpp @@ -0,0 +1,735 @@ +#ifndef STAN_MATH_PRIM_PROB_WIENER4_LCDF_HPP +#define STAN_MATH_PRIM_PROB_WIENER4_LCDF_HPP + +#include + +namespace stan { +namespace math { +namespace internal { + +/** + * Make the expression finite + * + * @param x Expression to test + * @return Expression or limited to maximum numeric_limit + */ +template +inline auto make_finite(const T_x& x) { + if (x < std::numeric_limits::max()) { + return x; + } else { + return std::numeric_limits::max(); + } +} + +/** + * Calculate the probability term 'P' on log scale for distribution + * + * @param a The boundary separation + * @param v The drift rate + * @param w The relative starting point + * @return 'P' term + */ +template +inline auto log_probability_distribution(const T_a& a, const T_v& v, + const T_w& w) { + using ret_t = return_type_t; + if (fabs(v) == 0.0) { + return ret_t(log1m(w)); + } + auto two_va = 2.0 * v * a; + auto minus_two_va_one_minus_w = -two_va * (1.0 - w); + // This split prevents abort errors + if (minus_two_va_one_minus_w < 0) { + const auto exp_arg = exp(minus_two_va_one_minus_w); + auto two_vaw = two_va * w; + if (two_vaw > minus_two_va_one_minus_w) { + return log1m(exp_arg) - log_diff_exp(two_vaw, minus_two_va_one_minus_w); + } else if (two_vaw < minus_two_va_one_minus_w) { + return log1m(exp_arg) - log_diff_exp(minus_two_va_one_minus_w, two_vaw); + } else { + return log1m(exp_arg) - NEGATIVE_INFTY; + } + } else { + return log1m_exp(-minus_two_va_one_minus_w) - log1m_exp(two_va); + } +} + +/** + * Calculate the probability term 'P' on log scale for grad_a and grad_v + * + * @param a The boundary separation + * @param v The drift rate + * @param w The relative starting point + * @return 'P' term + */ +template +inline auto log_probability_GradAV(const T_a& a, const T_v& v, const T_w& w) { + using ret_t = return_type_t; + if (fabs(v) == 0.0) { + return ret_t(-w); + } + auto nearly_one = ret_t(1.0 - 1.1 * 1.0e-5); + ret_t log_prob; + if (v < 0) { + const auto two_av = 2.0 * a * v; + const auto two_va_one_minus_w = (two_av * (1.0 - w)); + const auto two_avw = two_av * w; + const auto exp_two_va_one_minus_w = exp(two_va_one_minus_w); + const auto exp_two_avw = exp(two_avw); + const auto exp_two_av = exp(two_av); + if (((exp_two_va_one_minus_w >= nearly_one) || (exp_two_avw >= nearly_one)) + || (exp_two_av >= nearly_one)) { + return ret_t(-w); + } + log_prob = LOG_TWO + two_va_one_minus_w - log1m(exp_two_va_one_minus_w); + auto log_quotient = log1m(exp_two_avw) - log1m(exp_two_av); + if (log(w) > log_quotient) { + return exp(log_prob) * (w - exp(log_quotient)); + } else { + return -exp(log_prob) * (exp(log_quotient) - w); + } + } else { + const auto minus_two_av = -2.0 * a * v; + const auto minus_two_va_one_minus_w = minus_two_av * (1.0 - w); + const auto exp_minus_two_va_one_minus_w = exp(minus_two_va_one_minus_w); + const auto exp_minus_two_av = exp(minus_two_av); + if ((exp_minus_two_va_one_minus_w >= nearly_one) + || (exp_minus_two_av >= nearly_one)) { + return ret_t(-w); + } + log_prob = LOG_TWO - log1m(exp_minus_two_va_one_minus_w); + ret_t log_quotient; + if (minus_two_va_one_minus_w > minus_two_av) { + log_quotient = log_diff_exp(minus_two_va_one_minus_w, minus_two_av) + - log1m(exp_minus_two_av); + } else if (minus_two_va_one_minus_w < minus_two_av) { + log_quotient = log_diff_exp(minus_two_av, minus_two_va_one_minus_w) + - log1m(exp_minus_two_av); + } else { + log_quotient = NEGATIVE_INFTY; + } + if (log(w) > log_quotient) { + return -exp(log_prob + log_diff_exp(log(w), log_quotient)); + } else { + return exp(log_prob + log_diff_exp(log_quotient, log(w))); + } + } +} + +/** + * Log of Mill's ratio for the normal distribution + * + * @param x A scalar + * @return The natural logarithm of Mill's ratio + */ +template +inline auto logMill(T_x&& x) { + return std_normal_lcdf(-x) - std_normal_lpdf(x); +} + +/** + * Calculate the wiener4 distribution + * + * @tparam NaturalScale Whether to return the distribution on natural or + * log-scale + * + * @param y The reaction time in seconds + * @param a The boundary separation + * @param v The relative starting point + * @param w The drift rate + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return distribution + */ +template +inline auto wiener4_distribution(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, T_err log_err = log(1e-12)) { + using ret_t = return_type_t; + const auto neg_v = -v; + const auto one_m_w = 1.0 - w; + + const auto one_m_w_a_neg_v = one_m_w * a * neg_v; + + const auto K1 = 0.5 * (fabs(neg_v) / a * y - one_m_w); + const auto arg = fmax( + 0.0, fmin(1.0, exp(one_m_w_a_neg_v + square(neg_v) * y / 2.0 + log_err) + / 2.0)); + const auto K2 = (arg == 0) ? INFTY + : (arg == 1) ? NEGATIVE_INFTY + : -sqrt(y) / 2.0 / a * inv_Phi(arg); + const auto K_small_value = ceil(fmax(K1, K1 + K2)); + + const auto api = a / pi(); + const auto v_square = square(neg_v); + const auto sqrtL1 = sqrt(1.0 / y) * api; + const auto sqrtL2 = sqrt(fmax( + 1.0, -2.0 / y * square(api) + * (log_err + log(pi() * y / 2.0 * (v_square + square(pi() / a))) + + one_m_w_a_neg_v + v_square * y / 2.0))); + const auto K_large_value = ceil(fmax(sqrtL1, sqrtL2)); + + auto lg = LOG_TWO + LOG_PI - 2.0 * log(a); + + if (3 * K_small_value < K_large_value) { + const auto sqrt_y = sqrt(y); + const auto neg_vy = neg_v * y; + ret_t fplus = NEGATIVE_INFTY; + ret_t fminus = NEGATIVE_INFTY; + for (auto k = K_small_value; k >= 0; --k) { + auto rj = a * (2.0 * k + one_m_w); + auto dj = std_normal_lpdf(rj / sqrt_y); + auto pos1 = dj + logMill((rj - neg_vy) / sqrt_y); + auto pos2 = dj + logMill((rj + neg_vy) / sqrt_y); + fplus = log_sum_exp(fplus, log_sum_exp(pos1, pos2)); + rj = a * (2.0 * k + 2.0 - one_m_w); + dj = std_normal_lpdf(rj / sqrt_y); + auto neg1 = dj + logMill((rj - neg_vy) / sqrt_y); + auto neg2 = dj + logMill((rj + neg_vy) / sqrt_y); + fminus = log_sum_exp(fminus, log_sum_exp(neg1, neg2)); + } + auto ans = ret_t(0.0); + ans = fplus > fminus ? log_diff_exp(fplus, fminus) + : log_diff_exp(fminus, fplus); + ret_t log_distribution = ans - one_m_w_a_neg_v - square(neg_v) * y / 2; + return NaturalScale ? exp(log_distribution) : log_distribution; + } + const auto log_a = log(a); + const auto log_v = log(fabs(neg_v)); + ret_t fplus = NEGATIVE_INFTY; + ret_t fminus = NEGATIVE_INFTY; + for (auto k = K_large_value; k > 0; --k) { + auto log_k = log(k); + auto k_pi = k * pi(); + auto sin_k_pi_w = sin(k_pi * one_m_w); + if (sin_k_pi_w > 0) { + fplus = log_sum_exp( + fplus, log_k + - log_sum_exp(2.0 * log_v, 2.0 * (log_k + LOG_PI - log_a)) + - 0.5 * square(k_pi / a) * y + log(sin_k_pi_w)); + } else if (sin_k_pi_w < 0) { + fminus = log_sum_exp( + fminus, log_k + - log_sum_exp(2.0 * log_v, 2.0 * (log_k + LOG_PI - log_a)) + - 0.5 * square(k_pi / a) * y + log(-sin_k_pi_w)); + } + } + ret_t ans = NEGATIVE_INFTY; + ans = fplus > fminus ? log_diff_exp(fplus, fminus) + : log_diff_exp(fminus, fplus); + auto summand_1 = log_probability_distribution(a, neg_v, one_m_w); + auto summand_2 = lg + (ans - one_m_w_a_neg_v - 0.5 * square(neg_v) * y); + ret_t log_distribution = NEGATIVE_INFTY; + if (summand_1 > summand_2) { + log_distribution = log_diff_exp(summand_1, summand_2); + } else if (summand_1 < summand_2) { + log_distribution = log_diff_exp(summand_2, summand_1); + } + return NaturalScale ? exp(log_distribution) : log_distribution; +} + +/** + * Calculate derivative of the wiener4 distribution w.r.t. 'a' (natural-scale) + * + * @param y The reaction time in seconds + * @param a The boundary separation + * @param v The relative starting point + * @param w The drift rate + * @param cdf The value of the distribution + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to a + */ +template +inline auto wiener4_cdf_grad_a(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, T_cdf&& cdf, + T_err log_err = log(1e-12)) { + using ret_t = return_type_t; + const auto neg_v = -v; + const auto one_m_w = 1 - w; + + const auto one_m_w_neg_v = one_m_w * neg_v; + const auto one_m_w_a_neg_v = one_m_w_neg_v * a; + + const auto log_y = log(y); + const auto log_a = log(a); + auto C1 = ret_t( + LOG_TWO - log_sum_exp(2.0 * log(fabs(neg_v)), 2.0 * (LOG_PI - log_a))); + C1 = log_sum_exp(C1, log_y); + const auto factor = one_m_w_a_neg_v + square(neg_v) * y / 2.0 + log_err; + const auto alphK = fmin(factor + LOG_PI + log_y + log_a - LOG_TWO - C1, 0.0); + const auto K = a / pi() / sqrt(y); + const auto K_large_value + = ceil(fmax(fmax(sqrt(-2.0 * alphK / y) * a / pi(), K), ret_t(1.0))); + + const auto sqrt_y = sqrt(y); + const auto wdash = fmin(one_m_w, w); + const auto ueps + = fmin(-1.0, 2.0 * (factor + log(a) - log1p(square(neg_v) * y)) + LOG_PI); + const auto K_small + = (sqrt_y * sqrt(-(ueps - sqrt(-2.0 * ueps - 2.0))) - a * wdash) / a; + const auto K_large = sqrt_y / a - wdash; + const auto K_small_value = ceil(fmax(fmax(K_small, K_large), ret_t(1.0))); + + // Depending on the Ks use formula for small reaction times or large + // reaction times (see Navarro & Fuss, 2009) + if (K_large_value > 4 * K_small_value) { + const auto neg_vy = neg_v * y; + auto ans = ret_t(0.0); + auto F_k = ret_t(0.0); + for (auto k = K_small_value; k >= 0; --k) { + auto r_k_1 = a * (2.0 * k + one_m_w); + auto x_1 = r_k_1 - neg_vy; + auto x_over_sqrt_y_1 = x_1 / sqrt_y; + auto d_k_1 = std_normal_lpdf(r_k_1 / sqrt_y); + auto temp_1 = make_finite(exp(d_k_1 + logMill(x_over_sqrt_y_1))); + auto exp_d_k_1 = exp(d_k_1); + auto ans_1 = -temp_1 * neg_vy - sqrt_y * exp_d_k_1; + + auto x_2 = r_k_1 + neg_vy; + auto x_over_sqrt_y_2 = x_2 / sqrt_y; + auto temp_2 = make_finite(exp(d_k_1 + logMill(x_over_sqrt_y_2))); + auto ans_2 = temp_2 * neg_vy - sqrt_y * exp_d_k_1; + auto r_k_2 = a * (2.0 * k + 1.0 + w); + auto d_k_2 = std_normal_lpdf(r_k_2 / sqrt_y); + + auto x_3 = r_k_2 - neg_vy; + auto x_over_sqrt_y_3 = x_3 / sqrt_y; + auto temp_3 = make_finite(exp(d_k_2 + logMill(x_over_sqrt_y_3))); + auto exp_d_k_2 = exp(d_k_2); + auto ans_3 = -temp_3 * neg_vy - sqrt_y * exp_d_k_2; + + auto x_4 = r_k_2 + neg_vy; + auto x_over_sqrt_y_4 = x_4 / sqrt_y; + auto temp_4 = make_finite(exp(d_k_2 + logMill(x_over_sqrt_y_4))); + auto ans_4 = temp_4 * neg_vy - sqrt_y * exp_d_k_2; + + ans += (ans_1 + ans_2 + ans_3 - ans_4) * (2.0 * k + one_m_w) + + ans_3 * one_m_w; + } + F_k = make_finite(exp(one_m_w_a_neg_v + 0.5 * square(neg_v) * y)); + const auto summands_small_y = ans / (y * F_k); + return -one_m_w_neg_v * cdf + summands_small_y; + } + ret_t ans = 0.0; + for (auto k = K_large_value; k > 0; --k) { + const auto kpi = k * pi(); + const auto kpia2 = square(kpi / a); + const auto denom = square(neg_v) + kpia2; + auto last = (square(kpi) / pow(a, 3) * (y + 2.0 / denom)) * k / denom + * exp(-0.5 * kpia2 * y); + ans -= last * sin(kpi * one_m_w); + } + const ret_t prob + = make_finite(exp(log_probability_distribution(a, neg_v, one_m_w))); + const auto dav = log_probability_GradAV(a, neg_v, one_m_w); + auto dav_neg_v = dav * neg_v; + auto prob_deriv = fabs(neg_v) == 0 + ? ret_t(0.0) + : is_inf(dav_neg_v) ? NEGATIVE_INFTY : dav_neg_v * prob; + ans = (-2.0 / a - one_m_w_neg_v) * (cdf - prob) + + ans * (2.0 * pi() / square(a)) + * exp(-one_m_w_a_neg_v - 0.5 * square(neg_v) * y); + return prob_deriv + ans; +} + +/** + * Calculate derivative of the wiener4 distribution w.r.t. 'v' (natural-scale) + * + * @param y The reaction time in seconds + * @param a The boundary separation + * @param v The relative starting point + * @param w The drift rate + * @param cdf The value of the distribution + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to v + */ +template +inline auto wiener4_cdf_grad_v(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, T_cdf&& cdf, + T_err log_err = log(1e-12)) { + using ret_t = return_type_t; + const auto neg_v = -v; + const auto one_m_w = 1.0 - w; + + const auto one_m_w_a = one_m_w * a; + const auto one_m_w_a_neg_v = one_m_w_a * neg_v; + + const auto log_y = log(y); + const auto factor = one_m_w_a_neg_v + square(neg_v) * y / 2.0 + log_err; + + const auto log_a = log(a); + auto K_large_value = ret_t(1.0); + if (neg_v != 0) { + const auto temp = -make_finite(exp(log_a - LOG_PI - 0.5 * log_y)); + const auto log_v = log(fabs(neg_v)); + auto alphK_large = fmin(exp(factor + 0.5 * (7 * LOG_PI + log_y) + - 2.5 * LOG_TWO - 3 * log_a - log_v), + 1.0); + K_large_value + = fmax(ceil((alphK_large == 0) + ? ret_t(INFTY) + : (alphK_large == 1) ? ret_t(NEGATIVE_INFTY) + : temp * inv_Phi(alphK_large)), + ret_t(1.0)); + } + + const auto sqrt_y = sqrt(y); + const auto wdash = fmin(one_m_w, w); + auto K_large = fabs(neg_v) / a * y - wdash; + const auto alphK_small = factor + 0.5 * (LOG_TWO - log_y + LOG_PI); + const auto K_small + = (alphK_small < 0) ? sqrt_y * sqrt(-2.0 * alphK_small) / a - wdash : 0; + const auto K_small_value = ceil(fmax(fmax(K_small, K_large), ret_t(1.0))); + if (K_large_value > 4 * K_small_value) { + const auto sqrt_y = sqrt(y); + const auto neg_vy = neg_v * y; + auto ans = ret_t(0.0); + auto F_k = ret_t(0.0); + for (auto k = K_small_value; k >= 0; --k) { + auto r_k_1 = a * (2.0 * k + one_m_w); + auto d_k_1 = std_normal_lpdf(r_k_1 / sqrt_y); + auto x_1 = r_k_1 - neg_vy; + auto x_over_sqrt_y_1 = x_1 / sqrt_y; + auto ans_1 = make_finite(exp(d_k_1 + logMill(x_over_sqrt_y_1))); + + auto x_2 = r_k_1 + neg_vy; + auto x_over_sqrt_y_2 = x_2 / sqrt_y; + auto ans_2 = make_finite(exp(d_k_1 + logMill(x_over_sqrt_y_2))); + auto r_k_2 = a * (2.0 * k + 1.0 + w); + auto d_k_2 = std_normal_lpdf(r_k_2 / sqrt_y); + + auto x_3 = r_k_2 - neg_vy; + auto x_over_sqrt_y_3 = x_3 / sqrt_y; + auto ans_3 = make_finite(exp(d_k_2 + logMill(x_over_sqrt_y_3))); + + auto x_4 = r_k_2 + neg_vy; + auto x_over_sqrt_y_4 = x_4 / sqrt_y; + auto ans_4 = make_finite(exp(d_k_2 + logMill(x_over_sqrt_y_4))); + ans += -ans_1 * x_1 + ans_2 * x_2 + ans_3 * x_3 - ans_4 * x_4; + } + F_k = make_finite(exp(one_m_w_a_neg_v + 0.5 * square(neg_v) * y)); + const auto summands_small_y = ans / F_k; + return (one_m_w_a + neg_vy) * cdf - summands_small_y; + } + ret_t ans = 0.0; + for (auto k = K_large_value; k > 0; --k) { + const auto kpi = k * pi(); + const auto kpia2 = square(kpi / a); + const auto ekpia2y = exp(-0.5 * kpia2 * y); + const auto denom = square(neg_v) + kpia2; + const auto denomk = k / denom; + auto last = denomk * ekpia2y / denom; + ans -= last * sin(kpi * one_m_w); + } + const ret_t prob + = make_finite(exp(log_probability_distribution(a, neg_v, one_m_w))); + const auto dav = log_probability_GradAV(a, neg_v, one_m_w); + auto dav_a = dav * a; + auto prob_deriv = is_inf(dav_a) ? ret_t(NEGATIVE_INFTY) : dav_a * prob; + ans = (-one_m_w_a + v * y) * (cdf - prob) + + ans * 4.0 * v * pi() / square(a) + * exp(-one_m_w_a_neg_v - 0.5 * square(neg_v) * y); + return -(prob_deriv + ans); +} + +/** + * Calculate derivative of the wiener4 distribution w.r.t. 'w' (natural-scale) + * + * @param y The reaction time in seconds + * @param a The boundary separation + * @param v The relative starting point + * @param w The drift rate + * @param cdf The value of the distribution + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to w + */ +template +inline auto wiener4_cdf_grad_w(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, T_cdf&& cdf, + T_err log_err = log(1e-12)) { + using ret_t = return_type_t; + const auto neg_v = -v; + const auto one_m_w = 1 - w; + + const auto one_m_w_a_neg_v = one_m_w * a * neg_v; + + const auto factor = one_m_w_a_neg_v + square(neg_v) * y / 2.0 + log_err; + + const auto log_y = log(y); + const auto log_a = log(a); + const auto temp = -make_finite(exp(log_a - LOG_PI - 0.5 * log_y)); + auto alphK_large + = fmin(exp(factor + 0.5 * (LOG_PI + log_y) - 1.5 * LOG_TWO - log_a), 1.0); + alphK_large = fmax(0.0, alphK_large); + const auto K_large_value + = fmax(ceil((alphK_large == 0) + ? ret_t(INFTY) + : (alphK_large == 1) ? ret_t(NEGATIVE_INFTY) + : temp * inv_Phi(alphK_large)), + ret_t(1.0)); + + const auto sqrt_y = sqrt(y); + const auto wdash = fmin(one_m_w, w); + const auto K_large = fabs(neg_v) / a * y - wdash; + const auto lv = log1p(square(neg_v) * y); + const auto alphK_small = factor - LOG_TWO - lv; + const auto arg = fmin(exp(alphK_small), 1.0); + const auto K_small + = (arg == 0) + ? INFTY + : (arg == 1) ? NEGATIVE_INFTY : -sqrt_y / a * inv_Phi(arg) - wdash; + const auto K_small_value = ceil(fmax(fmax(K_small, K_large), ret_t(1.0))); + + if (K_large_value > 4 * K_small_value) { + const auto sqrt_y = sqrt(y); + const auto neg_vy = neg_v * y; + auto ans = ret_t(0.0); + auto F_k = ret_t(0.0); + for (auto k = K_small_value; k >= 0; --k) { + auto r_k_1 = a * (2.0 * k + one_m_w); + auto d_k_1 = std_normal_lpdf(r_k_1 / sqrt_y); + auto x_1 = r_k_1 - neg_vy; + auto x_over_sqrt_y_1 = x_1 / sqrt_y; + auto temp_1 = make_finite(exp(d_k_1 + logMill(x_over_sqrt_y_1))); + auto exp_d_k_1 = exp(d_k_1); + auto ans_1 = -temp_1 * neg_vy - sqrt_y * exp_d_k_1; + + auto x_2 = r_k_1 + neg_vy; + auto x_over_sqrt_y_2 = x_2 / sqrt_y; + auto temp_2 = make_finite(exp(d_k_1 + logMill(x_over_sqrt_y_2))); + auto ans_2 = temp_2 * neg_vy - sqrt_y * exp_d_k_1; + auto r_k_2 = a * (2.0 * k + 1.0 + w); + auto d_k_2 = std_normal_lpdf(r_k_2 / sqrt_y); + + auto x_3 = r_k_2 - neg_vy; + auto x_over_sqrt_y_3 = x_3 / sqrt_y; + auto temp_3 = make_finite(exp(d_k_2 + logMill(x_over_sqrt_y_3))); + auto exp_d_k_2 = exp(d_k_2); + auto ans_3 = -temp_3 * neg_vy - sqrt_y * exp_d_k_2; + + auto x_4 = r_k_2 + neg_vy; + auto x_over_sqrt_y_4 = x_4 / sqrt_y; + auto temp_4 = make_finite(exp(d_k_2 + logMill(x_over_sqrt_y_4))); + auto ans_4 = temp_4 * neg_vy - sqrt_y * exp_d_k_2; + + ans += (ans_1 + ans_2 + ans_3 - ans_4) * a; + } + F_k = make_finite(exp(one_m_w_a_neg_v + 0.5 * square(neg_v) * y)); + const auto summands_small_y = ans / (y * F_k); + return neg_v * a * cdf - summands_small_y; + } + ret_t ans = 0.0; + for (auto k = K_large_value; k > 0; --k) { + const auto kpi = k * pi(); + const auto kpia2 = square(kpi / a); + const auto ekpia2y = exp(-0.5 * kpia2 * y); + const auto denom = square(neg_v) + kpia2; + const auto denomk = k / denom; + auto last = kpi; + last *= denomk * ekpia2y; + ans -= last * cos(kpi * one_m_w); + } + const auto evaw = exp(-one_m_w_a_neg_v - 0.5 * square(neg_v) * y); + const ret_t prob + = make_finite(exp(log_probability_distribution(a, neg_v, one_m_w))); + + // Calculate the probability term 'P' on log scale + auto dav = ret_t(-1 / w); + if (neg_v != 0) { + auto nearly_one = ret_t(1.0 - 1.0e-6); + const auto sign_v = (neg_v < 0) ? 1 : -1; + const auto sign_two_va_one_minus_w = sign_v * (2.0 * neg_v * a * w); + const auto exp_arg = exp(sign_two_va_one_minus_w); + if (exp_arg >= nearly_one) { + dav = -1.0 / w; + } else { + auto prob = LOG_TWO + log(fabs(neg_v)) + log(a) - log1m(exp_arg); + if (neg_v < 0) { + prob += sign_two_va_one_minus_w; + } + dav = -exp(prob); + } + } + + const auto pia2 = 2.0 * pi() / square(a); + auto prob_deriv = dav * prob; + ans = v * a * (cdf - prob) + ans * pia2 * evaw; + return -(prob_deriv + ans); +} +} // namespace internal + +/** + * Log-CDF function for the 4-parameter Wiener distribution. + * See 'wiener_lpdf' for more comprehensive documentation. + * + * @tparam T_y type of reaction time + * @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 + * + * @param y 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 precision_derivatives Level of precision in estimation + * @return The log of the Wiener first passage time distribution with + * the specified arguments for upper boundary responses + */ +template +inline auto wiener_lcdf(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) { + using T_partials_return = partials_return_t; + 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 internal::GradientCalc; + using ret_t = return_type_t; + + 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; + + auto y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + auto a_val = to_ref(as_value_column_array_or_scalar(a_ref)); + auto v_val = to_ref(as_value_column_array_or_scalar(v_ref)); + auto w_val = to_ref(as_value_column_array_or_scalar(w_ref)); + auto t0_val = to_ref(as_value_column_array_or_scalar(t0_ref)); + + if (!include_summand::value) { + return ret_t(0.0); + } + + static constexpr const char* function_name = "wiener4_lcdf"; + if (size_zero(y, a, t0, w, v)) { + return ret_t(0.0); + } + + check_consistent_sizes(function_name, "Random variable", y, + "Boundary separation", a, "Drift rate", v, + "A-priori bias", w, "Nondecision time", t0); + 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); + + const size_t N = max_size(y, a, t0, w, v); + + 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); + 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()); + } + } + + // for precs. 1e-6, 1e-12, see Hartmann et al. (2021), Henrich et al. (2023) + const auto log_error_cdf = log(1e-6); + const auto log_error_derivative = log(precision_derivatives); + const T_partials_return log_error_absolute = log(1e-12); + T_partials_return lcdf = 0.0; + auto ops_partials + = make_partials_propagator(y_ref, a_ref, t0_ref, w_ref, v_ref); + + const double LOG_FOUR = std::log(4.0); + + // calculate distribution and partials + for (size_t i = 0; i < N; i++) { + 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 T_partials_return log_cdf + = internal::estimate_with_err_check<4, 0, GradientCalc::OFF, + GradientCalc::OFF>( + [](auto&&... args) { + return internal::wiener4_distribution(args...); + }, + log_error_cdf - LOG_TWO, y_value - t0_value, a_value, v_value, + w_value, log_error_absolute); + + const T_partials_return cdf = exp(log_cdf); + + lcdf += log_cdf; + + const auto new_est_err = log_cdf + log_error_derivative - LOG_FOUR; + + if (!is_constant_all::value || !is_constant_all::value) { + const auto deriv_y + = internal::estimate_with_err_check<5, 0, GradientCalc::OFF, + GradientCalc::ON>( + [](auto&&... args) { + return internal::wiener5_density(args...); + }, + new_est_err, y_value - t0_value, a_value, v_value, w_value, 0, + log_error_absolute); + + if (!is_constant_all::value) { + partials<0>(ops_partials)[i] = deriv_y / cdf; + } + if (!is_constant_all::value) { + partials<2>(ops_partials)[i] = -deriv_y / cdf; + } + } + 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::wiener4_cdf_grad_a(args...); + }, + new_est_err, y_value - t0_value, a_value, v_value, w_value, cdf, + log_error_absolute) + / cdf; + } + 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::wiener4_cdf_grad_w(args...); + }, + new_est_err, y_value - t0_value, a_value, v_value, w_value, cdf, + log_error_absolute) + / cdf; + } + if (!is_constant_all::value) { + partials<4>(ops_partials)[i] + = internal::wiener4_cdf_grad_v(y_value - t0_value, a_value, v_value, + w_value, cdf, log_error_absolute) + / cdf; + } + } // for loop + return ops_partials.build(lcdf); +} +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/wiener5_lpdf.hpp b/stan/math/prim/prob/wiener5_lpdf.hpp index 689e5c1e2fd..bd4c8beecfa 100644 --- a/stan/math/prim/prob/wiener5_lpdf.hpp +++ b/stan/math/prim/prob/wiener5_lpdf.hpp @@ -11,63 +11,64 @@ namespace internal { enum GradientCalc { OFF = 0, ON = 1 }; /** - * Calculate the 'error_term' term for a wiener5 density or gradient + * Calculate the 'log_error_term' term for a wiener5 density or gradient * - * @tparam T_y type of scalar variable + * @tparam T_y type of reaction time * @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 y The reaction time in seconds * @param a The boundary separation - * @param v_value The drift rate - * @param w_value The relative starting point + * @param v The drift rate + * @param w The relative starting point * @param sv The inter-trial variability of the drift rate - * @return 'error_term' term + * @return 'log_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; +inline auto wiener5_compute_log_error_term(T_y&& y, T_a&& a, T_v&& v, T_w&& w, + T_sv&& sv) noexcept { + const auto one_m_w = 1.0 - w; + const auto neg_v = -v; 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_avw = 2.0 * a * neg_v * one_m_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)); + return stan::math::eval( + (sv_sqr * square(a * one_m_w) - two_avw - square(neg_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 GradW Whether the calculation is for gradient with respect to 'w' + * @tparam T_y type of reaction time * @tparam T_a type of boundary separation - * @tparam T_w_value type of relative starting point + * @tparam T_w type of relative starting point * @tparam T_err type of error * - * @param y A scalar variable; the reaction time in seconds + * @param y The reaction time in seconds * @param a The boundary separation - * @param w_value The relative starting point + * @param w 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 { + typename T_w, typename T_err> +inline auto wiener5_n_terms_small_t(T_y&& y, T_a&& a, T_w&& w, + 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 two_log_a = 2.0 * log(a); const auto log_y_asq = log(y) - two_log_a; - const auto w = 1.0 - w_value; + const auto one_m_w = 1.0 - w; - const auto n_1_factor = Density ? 2 : 3; - const auto n_1 = (sqrt(n_1_factor * y_asq) + w) / 2.0; + constexpr auto n_1_factor = Density ? 2.0 : 3.0; + const auto n_1 = (sqrt(n_1_factor * y_asq) + one_m_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 @@ -75,8 +76,9 @@ inline auto wiener5_n_terms_small_t(T_y&& y, T_a&& a, T_w_value&& w_value, 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; + const auto n_2 = (arg > 0) ? GradW ? 0.5 * (sqrt(arg) + one_m_w) + : 0.5 * (sqrt(arg) - one_m_w) + : n_1; return ceil(fmax(n_1, n_2)); } @@ -84,23 +86,22 @@ inline auto wiener5_n_terms_small_t(T_y&& y, T_a&& a, T_w_value&& w_value, /** * Calculate the 'n_terms_large_t' term for a wiener5 density * - * @tparam T_y type of scalar variable + * @tparam T_y type of reaction time * @tparam T_a type of boundary separation - * @tparam T_w_value type of relative starting point + * @tparam T_w type of relative starting point * @tparam T_err type of error * - * @param y A scalar variable; the reaction time in seconds + * @param y The reaction time in seconds * @param a The boundary separation - * @param w_value The relative starting point + * @param w 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 { +template +inline auto wiener5_density_large_reaction_time_terms(T_y&& y, T_a&& a, T_w&& w, + T_err error) noexcept { const auto y_asq = y / square(a); - const auto log_y_asq = log(y) - 2 * log(a); + const auto log_y_asq = log(y) - 2.0 * 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); @@ -112,79 +113,85 @@ inline auto wiener5_density_large_reaction_time_terms(T_y&& y, T_a&& a, /** * 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 GradW Whether the calculation is for gradient with respect to 'w' + * @tparam T_y type of reaction time * @tparam T_a type of boundary separation - * @tparam T_w_value type of relative starting point + * @tparam T_w type of relative starting point * @tparam T_err type of error * - * @param y A scalar variable; the reaction time in seconds + * @param y The reaction time in seconds * @param a The boundary separation - * @param w_value The relative starting point + * @param w 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 { + T_w&& w, + T_err error) noexcept { const auto y_asq = y / square(a); - const auto log_y_asq = log(y) - 2 * log(a); + const auto log_y_asq = log(y) - 2.0 * log(a); static constexpr double PI_SQUARED = pi() * pi(); - const auto n_1_factor = GradW ? 2 : 3; + static constexpr auto n_1_factor = GradW ? 2.0 : 3.0; 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)); + if constexpr (GradW) { + const auto arg = -(u_eps - sqrt(-2.0 * u_eps - 2.0)); + auto n_2 = (arg > 0) ? sqrt(arg / y_asq) / pi() : n_1; + return ceil(fmax(n_1, n_2)); + } else { + const auto arg + = -(2.0 / PI_SQUARED / y_asq) * (u_eps - sqrt(-2.0 * u_eps - 2.0)); + auto n_2 = (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 GradW Whether the calculation is for gradient with respect to 'w' + * @tparam T_y type of reaction time * @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 y The reaction time in seconds * @param a The boundary separation - * @param w_value The relative starting point + * @param w 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, +inline auto wiener5_log_sum_exp(T_y&& y, T_a&& a, T_w&& w, T_nsmall&& n_terms_small_t, T_nlarge&& n_terms_large_t) noexcept { + using ret_t = return_type_t; + const auto y_asq = y / square(a); - const auto w = 1.0 - w_value; + const auto one_m_w = 1.0 - w; 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); + = Density ? (2.0 * n_terms_small_t <= n_terms_large_t) + : (2.0 * 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) { + if constexpr (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 w_plus_2_k = one_m_w + 2.0 * k; + const auto w_minus_2_k = one_m_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) @@ -207,18 +214,20 @@ inline auto wiener5_log_sum_exp(T_y&& y, T_a&& a, T_w&& w_value, fminus = log_sum_exp(fminus, summand_minus); } } - const auto square_w_minus_offset = square(w) - y_asq; + const auto square_w_minus_offset = square(one_m_w) - y_asq; if (square_w_minus_offset > 0) { - const auto new_val = log(square_w_minus_offset) - square(w) * scaling; + const auto new_val + = log(square_w_minus_offset) - square(one_m_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; + const auto new_val + = log(-square_w_minus_offset) - square(one_m_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 w_plus_2_k = one_m_w + 2.0 * k; + const auto w_minus_2_k = one_m_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); @@ -234,14 +243,14 @@ inline auto wiener5_log_sum_exp(T_y&& y, T_a&& a, T_w&& w_value, fminus = summand_minus + log1p_exp(fminus - summand_minus); } } - const auto new_val = mult * log(w) - square(w) * scaling; + const auto new_val = mult * log(one_m_w) - square(one_m_w) * scaling; fplus = log_sum_exp(fplus, new_val); } } else { // for large t - constexpr double mult = (Density ? 1 : (GradW ? 2 : 3)); + constexpr double mult = (Density ? 1.0 : (GradW ? 2.0 : 3.0)); 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); + const auto check = (GradW) ? cos(pi_k * one_m_w) : sin(pi_k * one_m_w); if (check > 0) { fplus = log_sum_exp( fplus, mult * log(k) - square(pi_k) * scaling + log(check)); @@ -270,44 +279,44 @@ inline auto wiener5_log_sum_exp(T_y&& y, T_a&& a, T_w&& w_value, * * @tparam NaturalScale Whether to return the density on natural (true) or * log-scale (false) - * @tparam T_y type of scalar variable + * @tparam T_y type of reaction time * @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 y The reaction time in seconds * @param a The boundary separation - * @param v_value The drift rate - * @param w_value The relative starting point + * @param v The drift rate + * @param w The relative starting point * @param sv The inter-trial variability of the drift rate - * @param err The log error tolerance + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums * @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); +inline auto wiener5_density(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, const T_sv& sv, + T_err log_err = log(1e-12)) noexcept { + const auto log_error_term = wiener5_compute_log_error_term(y, a, v, w, sv); + const auto log_error = (log_err - log_error_term); const auto n_terms_small_t - = wiener5_n_terms_small_t( - y, a, w_value, error); + = wiener5_n_terms_small_t(y, a, w, + log_error); const auto n_terms_large_t - = wiener5_density_large_reaction_time_terms(y, a, w_value, error); + = wiener5_density_large_reaction_time_terms(y, a, w, log_error); auto res = wiener5_log_sum_exp( - y, a, w_value, n_terms_small_t, n_terms_large_t) + y, a, w, 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; + auto log_density = log_error_term - 0.5 * LOG_TWO - LOG_SQRT_PI + - 1.5 * (log(y) - 2.0 * log(a)) + res; return NaturalScale ? exp(log_density) : log_density; } else { - auto log_density = error_term + res + LOG_PI; + auto log_density = log_error_term + res + LOG_PI; return NaturalScale ? exp(log_density) : log_density; } } @@ -315,66 +324,67 @@ inline auto wiener5_density(const T_y& y, const T_a& a, const T_v& v_value, /** * Calculate the derivative of the wiener5 density w.r.t. 't' * - * @tparam WrtLog Whether to return the derivative w.r.t. + * @tparam WrtLog Whether to return the derivative with respect to * the natural (true) or log-scale (false) density - * @tparam T_y type of scalar variable + * @tparam T_y type of reaction time * @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 y The reaction time in seconds * @param a The boundary separation - * @param v_value The drift rate - * @param w_value The relative starting point + * @param v The drift rate + * @param w 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 + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to 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); +inline auto wiener5_grad_t(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, const T_sv& sv, + T_err log_err = log(1e-12)) noexcept { + const auto two_log_a = 2.0 * 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 log_error_term = wiener5_compute_log_error_term(y, a, v, w, sv); + const auto one_m_w = 1.0 - w; + const auto neg_v = -v; 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(sv_sqr) * (y + square(a * one_m_w)) + + sv_sqr * (1.0 - (2.0 * a * neg_v * one_m_w)) + square(neg_v)) / square(one_plus_svsqr_y); - const auto error = (err - error_term) + two_log_a; + const auto log_error = (log_err - log_error_term) + two_log_a; const auto n_terms_small_t = wiener5_n_terms_small_t( - y, a, w_value, error); + y, a, w, log_error); const auto n_terms_large_t = wiener5_gradient_large_reaction_time_terms( - y, a, w_value, error); + y, a, w, log_error); auto wiener_res = wiener5_log_sum_exp( - y, a, w_value, n_terms_small_t, n_terms_large_t); + y, a, w, 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); + y, a, v, w, sv, log_err - error_log_density); + if (2.0 * n_terms_small_t < n_terms_large_t) { + auto ans + = density_part_one - 1.5 / y + + newsign + * exp(log_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 + * exp(log_error_term - two_log_a + 3.0 * LOG_PI - LOG_TWO + result - log_density); return WrtLog ? ans * exp(log_density) : ans; } @@ -383,65 +393,65 @@ inline auto wiener5_grad_t(const T_y& y, const T_a& a, const T_v& v_value, /** * Calculate the derivative of the wiener5 density w.r.t. 'a' * - * @tparam WrtLog Whether to return the derivative w.r.t. + * @tparam WrtLog Whether to return the derivative with respect to * the natural (true) or log-scale (false) density - * @tparam T_y type of scalar variable + * @tparam T_y type of reaction time * @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 y The reaction time in seconds * @param a The boundary separation - * @param v_value The drift rate - * @param w_value The relative starting point + * @param v The drift rate + * @param w 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 + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to 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 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; +inline auto wiener5_grad_a(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, const T_sv& sv, + T_err log_err = log(1e-12)) noexcept { + const auto two_log_a = 2.0 * log(a); + const auto log_error_term = wiener5_compute_log_error_term(y, a, v, w, sv); + const auto one_m_w = 1.0 - w; const auto sv_sqr = square(sv); - const auto one_plus_svsqr_y = 1 + sv_sqr * y; + const auto one_plus_svsqr_y = 1.0 + 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; + = (v * one_m_w + sv_sqr * square(one_m_w) * a) / one_plus_svsqr_y; + const auto log_error + = log_err - log_error_term + 3.0 * log(a) - log(y) - LOG_TWO; const auto n_terms_small_t = wiener5_n_terms_small_t( - y, a, w_value, error); + y, a, w, log_error); const auto n_terms_large_t = wiener5_gradient_large_reaction_time_terms( - y, a, w_value, error); + y, a, w, log_error); auto wiener_res = wiener5_log_sum_exp( - y, a, w_value, n_terms_small_t, n_terms_large_t); + y, a, w, n_terms_small_t, n_terms_large_t); auto&& result = wiener_res.first; auto&& newsign = wiener_res.second; - const auto error_log_density = log( + const auto log_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); + y, a, v, w, sv, log_err - log_error_log_density); + if (2.0 * 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 + log_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); + * exp(log(y) + log_error_term - 3.0 * (log(a) - LOG_PI) + + result - log_density); return WrtLog ? ans * exp(log_density) : ans; } } @@ -449,31 +459,32 @@ inline auto wiener5_grad_a(const T_y& y, const T_a& a, const T_v& v_value, /** * Calculate the derivative of the wiener5 density w.r.t. 'v' * - * @tparam WrtLog Whether to return the derivative w.r.t. + * @tparam WrtLog Whether to return the derivative with respect to * the natural (true) or log-scale (false) density - * @tparam T_y type of scalar variable + * @tparam T_y type of reaction time * @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 y The reaction time in seconds * @param a The boundary separation - * @param v_value The drift rate - * @param w_value The relative starting point + * @param v The drift rate + * @param w 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 + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to 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); +inline auto wiener5_grad_v(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, const T_sv& sv, + T_err log_err = log(1e-12)) noexcept { + auto ans = (a * (1 - w) - v * y) / (1.0 + square(sv) * y); + if constexpr (WrtLog) { + return ans * wiener5_density(y, a, v, w, sv, log_err); } else { return ans; } @@ -482,62 +493,62 @@ inline auto wiener5_grad_v(const T_y& y, const T_a& a, const T_v& v_value, /** * Calculate the derivative of the wiener5 density w.r.t. 'w' * - * @tparam WrtLog Whether to return the derivative w.r.t. + * @tparam WrtLog Whether to return the derivative with respect to * the natural (true) or log-scale (false) density - * @tparam T_y type of scalar variable + * @tparam T_y type of reaction time * @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 y The reaction time in seconds * @param a The boundary separation - * @param v_value The drift rate - * @param w_value The relative starting point + * @param v The drift rate + * @param w 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 + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to 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); +inline auto wiener5_grad_w(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, const T_sv& sv, + T_err log_err = log(1e-12)) noexcept { + const auto two_log_a = 2.0 * 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 log_error_term = wiener5_compute_log_error_term(y, a, v, w, sv); + const auto one_m_w = 1.0 - w; const auto sv_sqr = square(sv); - const auto one_plus_svsqr_y = 1 + sv_sqr * y; + const auto one_plus_svsqr_y = 1.0 + 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); + = (v * a + sv_sqr * square(a) * one_m_w) / one_plus_svsqr_y; + const auto log_error = (log_err - log_error_term); const auto n_terms_small_t - = wiener5_n_terms_small_t( - y, a, w_value, error); + = wiener5_n_terms_small_t(y, a, w, + log_error); const auto n_terms_large_t - = wiener5_gradient_large_reaction_time_terms( - y, a, w_value, error); + = wiener5_gradient_large_reaction_time_terms(y, a, w, + log_error); auto wiener_res = wiener5_log_sum_exp( - y, a, w_value, n_terms_small_t, n_terms_large_t); + y, a, w, 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) { + y, a, v, w, sv, log_err - log(fabs(density_part_one))); + if (2.0 * n_terms_small_t < n_terms_large_t) { auto ans = -(density_part_one - newsign - * exp(result - (log_density - error_term) + * exp(result - (log_density - log_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)); + auto ans = -( + density_part_one + + newsign + * exp(result - (log_density - log_error_term) + 2.0 * LOG_PI)); return WrtLog ? ans * exp(log_density) : ans; } } @@ -545,37 +556,38 @@ inline auto wiener5_grad_w(const T_y& y, const T_a& a, const T_v& v_value, /** * Calculate the derivative of the wiener5 density w.r.t. 'sv' * - * @tparam WrtLog Whether to return the derivative w.r.t. + * @tparam WrtLog Whether to return the derivative with respect to * the natural (true) or log-scale (false) density - * @tparam T_y type of scalar variable + * @tparam T_y type of reaction time * @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 y The reaction time in seconds * @param a The boundary separation - * @param v_value The drift rate - * @param w_value The relative starting point + * @param v The drift rate + * @param w 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 + * @param log_err The log error tolerance in the computation of the number + * of terms for the infinite sums + * @return Gradient with respect to 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; +inline auto wiener5_grad_sv(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, const T_sv& sv, + T_err log_err = log(1e-12)) noexcept { + const auto one_plus_svsqr_y = 1.0 + square(sv) * y; + const auto one_m_w = 1.0 - w; + const auto neg_v = -v; const auto t1 = -y / one_plus_svsqr_y; - const auto t2 = (square(a * w) + 2 * a * v * w * y + square(v * y)) + const auto t2 = (square(a * one_m_w) + 2.0 * a * neg_v * one_m_w * y + + square(neg_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; + return WrtLog ? ans * wiener5_density(y, a, v, w, sv, log_err) : ans; } /** @@ -623,22 +635,24 @@ inline void assign_err(std::tuple& args_tuple, Scalar err) { * @tparam ArgsTupleT Type of tuple of arguments for functor * * @param functor Function to apply - * @param err Error value to check against + * @param log_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, +inline auto estimate_with_err_check(F&& functor, T_err&& log_err, ArgsTupleT&&... args_tuple) { auto result = functor(args_tuple...); auto log_fabs_result = LogResult ? log(fabs(result)) : fabs(result); - if (log_fabs_result < err) { + if (log_fabs_result < log_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); + const auto new_error = GradW7 ? log_err + log_fabs_result + LOG_TWO + : log_err + log_fabs_result; + if constexpr (NestedIndex != -1) { + assign_err(std::get(err_args_tuple), new_error); + } result = math::apply([](auto&& func, auto&&... args) { return func(args...); }, err_args_tuple, functor); @@ -649,9 +663,9 @@ inline auto estimate_with_err_check(F&& functor, T_err&& err, /** * Log-density function for the 5-parameter Wiener density. - * See 'wiener_lpdf' for more comprehensive documentation + * See 'wiener_lpdf' for more comprehensive documentation. * - * @tparam T_y type of scalar + * @tparam T_y type of reaction time * @tparam T_a type of boundary * @tparam T_t0 type of non-decision time * @tparam T_w type of relative starting point @@ -659,7 +673,7 @@ inline auto estimate_with_err_check(F&& functor, T_err&& err, * @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 y The reaction time in seconds * @param a The boundary separation * @param t0 The non-decision time * @param w The relative starting point @@ -676,23 +690,13 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, const double& precision_derivatives = 1e-4) { using T_partials_return = partials_return_t; using ret_t = return_type_t; - if constexpr (!include_summand::value) { - return ret_t(0.0); - } using T_y_ref = ref_type_t; using T_a_ref = ref_type_t; using T_t0_ref = ref_type_t; using T_w_ref = ref_type_t; using T_v_ref = ref_type_t; using T_sv_ref = ref_type_t; - - 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); + using internal::GradientCalc; T_y_ref y_ref = y; T_a_ref a_ref = a; @@ -701,12 +705,24 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, 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)); + auto y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + auto a_val = to_ref(as_value_column_array_or_scalar(a_ref)); + auto v_val = to_ref(as_value_column_array_or_scalar(v_ref)); + auto w_val = to_ref(as_value_column_array_or_scalar(w_ref)); + auto t0_val = to_ref(as_value_column_array_or_scalar(t0_ref)); + auto sv_val = to_ref(as_value_column_array_or_scalar(sv_ref)); + + if constexpr (!include_summand::value) { + return ret_t(0.0); + } + + 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); 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); @@ -722,7 +738,7 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, return ret_t(0.0); } const size_t N = max_size(y, a, t0, w, v, sv); - if (!N) { + if (N == 0) { return ret_t(0.0); } @@ -744,6 +760,7 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, } } + // for precs. 1e-6, 1e-12, see Hartmann et al. (2021), Henrich et al. (2023) 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); @@ -752,7 +769,7 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, 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; + const double LOG_FOUR = std::log(4.0); // calculate density and partials for (size_t i = 0; i < N; i++) { @@ -766,7 +783,6 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, 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) { @@ -779,6 +795,7 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, const auto new_est_err = l_density + log_error_derivative - LOG_FOUR; + // computation of derivatives and precision checks // 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 @@ -833,15 +850,6 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, 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_lccdf.hpp b/stan/math/prim/prob/wiener_full_lccdf.hpp new file mode 100644 index 00000000000..c5cf883c886 --- /dev/null +++ b/stan/math/prim/prob/wiener_full_lccdf.hpp @@ -0,0 +1,407 @@ +#ifndef STAN_MATH_PRIM_PROB_WIENER_FULL_LCCDF_HPP +#define STAN_MATH_PRIM_PROB_WIENER_FULL_LCCDF_HPP + +#include +#include + +namespace stan { +namespace math { +namespace internal { + +/** + * Calculate the derivative of the wiener7 density w.r.t. 'sw' + * + * @tparam T_y type of reaction time + * @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 tolerance + * + * @param y The reaction time in seconds + * @param a The boundary separation + * @param v The drift rate + * @param w The relative starting point + * @param sw The inter-trial variability of the relative starting point + * @param wildcard This parameter space is needed for a functor. Could be + * deleted when another solution is found + * @param log_error The log error tolerance + * @return Gradient with respect to sw + */ +template +inline auto wiener7_ccdf_grad_sw(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, const T_sw& sw, + T_err log_error) { + auto low = fmax(0.0, w - sw / 2.0); + auto high = fmin(1.0, w + sw / 2.0); + + const auto lower_value = wiener4_ccdf(y, a, v, low, log_error); + const auto upper_value = wiener4_ccdf(y, a, v, high, log_error); + return 0.5 * (lower_value + upper_value) / sw; +} + +} // namespace internal + +/** \ingroup prob_dists + * Returns the log CCDF of the Wiener distribution for a + * (Wiener) drift diffusion model with up to 7 parameters. If containers + * are supplied, returns the log sum of the probabilities. + * See 'wiener_full_lpdf' for more comprehensive documentation + * As the CDF goes to the probability to hit the upper bound + * (instead of 1, as it is usually the case) when the reaction time + * goes to infinity, the CCDF is defined as + * ccdf = probability_to_hit_the_upper_bound - cdf. + * + * @tparam T_y type of reaction time + * @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 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 log probability or log sum of probabilities 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. + */ + +template +inline auto wiener_lccdf(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-8) { + using ret_t = return_type_t; + 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 internal::GradientCalc; + using T_partials_return + = partials_return_t; + + 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; + + auto y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + auto a_val = to_ref(as_value_column_array_or_scalar(a_ref)); + auto v_val = to_ref(as_value_column_array_or_scalar(v_ref)); + auto w_val = to_ref(as_value_column_array_or_scalar(w_ref)); + auto t0_val = to_ref(as_value_column_array_or_scalar(t0_ref)); + auto sv_val = to_ref(as_value_column_array_or_scalar(sv_ref)); + auto sw_val = to_ref(as_value_column_array_or_scalar(sw_ref)); + auto st0_val = to_ref(as_value_column_array_or_scalar(st0_ref)); + + if (!include_summand::value) { + return ret_t(0.0); + } + + static constexpr const char* function_name = "wiener_lccdf"; + 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); + 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.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]) { + 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)) { + std::stringstream msg; + msg << ", but must be smaller than 2*(A-priori bias) = " + << 2.0 * 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)) { + std::stringstream msg; + msg << ", but must be smaller than 2*(1-A-priori bias) = " + << 2.0 * (1.0 - 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()); + } + } + + // for precs. 1e-6, 1e-12, see Hartmann et al. (2021), Henrich et al. (2023) + const T_partials_return log_error_cdf = log(1e-6); // precision for density + const auto error_bound = precision_derivatives; // precision for + // derivatives (controllable by user) + const auto log_error_derivatives = log(error_bound); + const T_partials_return absolute_error_hcubature = 0.0; + const T_partials_return relative_error_hcubature + = .9 * error_bound; // eps_rel(Integration) + const T_partials_return log_error_absolute = log(1e-12); + const int maximal_evaluations_hcubature = 6000; + T_partials_return lccdf = 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.0; + + // calculate density and partials + for (size_t i = 0; i < N; i++) { + if (sv_vec[i] == 0 && sw_vec[i] == 0 && st0_vec[i] == 0) { + result += wiener_lccdf(y_vec[i], a_vec[i], t0_vec[i], w_vec[i], + v_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 = (sv_value != 0) + (sw_value != 0) + (st0_value != 0); + check_positive(function_name, + "(Inter-trial variability in drift rate) + " + "(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); + for (int i = 0; i < dim; i++) { + xmin[i] = 0; + xmax[i] = 1; + } + if (sv_value != 0) { + xmin[0] = -1; + xmax[0] = 1; + } + 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_cdf + LOG_TWO + 1; + + 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); + + const T_partials_return ccdf + = internal::wiener7_integrate_cdf( + [&](auto&&... args) { return internal::wiener4_ccdf(args...); }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2); + lccdf += log(ccdf); + + hcubature_err = log_error_absolute - log_error_derivatives + log(fabs(ccdf)) + + 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 auto params_dt7 = std::make_tuple( + y_value, a_value, v_value, w_value, t0_value, sv_value, sw_value, + st0_value, log_error_absolute - LOG_TWO - 9 * LOG_TWO); + // computation of derivatives and precision checks + if (!is_constant_all::value || !is_constant_all::value) { + const T_partials_return deriv_t_7 + = -internal::wiener7_integrate_cdf< + GradientCalc::OFF, GradientCalc::OFF, GradientCalc::OFF, + GradientCalc::OFF, GradientCalc::ON>( + [&](auto&&... args) { + return internal::wiener5_density(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / ccdf; + if (!is_constant_all::value) { + partials<0>(ops_partials)[i] = deriv_t_7; + } + if (!is_constant_all::value) { + partials<2>(ops_partials)[i] = -deriv_t_7; + } + } + T_partials_return deriv; + if (!is_constant_all::value) { + partials<1>(ops_partials)[i] + = internal::wiener7_integrate_cdf( + [&](auto&&... args) { + return internal::wiener4_ccdf_grad_a(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / ccdf; + } + if (!is_constant_all::value) { + partials<3>(ops_partials)[i] + = internal::wiener7_integrate_cdf( + [&](auto&&... args) { + return internal::wiener4_ccdf_grad_w(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / ccdf; + } + if (!is_constant_all::value) { + partials<4>(ops_partials)[i] + = internal::wiener7_integrate_cdf( + [&](auto&&... args) { + return internal::wiener4_ccdf_grad_v(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / ccdf; + } + if (!is_constant_all::value) { + if (sv_value == 0) { + partials<5>(ops_partials)[i] = 0.0; + } else { + partials<5>(ops_partials)[i] + = internal::wiener7_integrate_cdf< + GradientCalc::OFF, GradientCalc::OFF, GradientCalc::ON>( + [&](auto&&... args) { + return internal::wiener4_ccdf_grad_v(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / ccdf; + } + } + if (!is_constant_all::value) { + if (sw_value == 0) { + partials<6>(ops_partials)[i] = 0.0; + } else { + if (st0_value == 0 && sv_value == 0) { + deriv = internal::estimate_with_err_check<5, 0, GradientCalc::OFF, + GradientCalc::ON>( + [](auto&&... args) { + return internal::wiener7_ccdf_grad_sw(args...); + }, + hcubature_err, y_value - t0_value, a_value, v_value, w_value, + sw_value, log_error_absolute - LOG_TWO); + deriv = deriv / ccdf - 1 / sw_value; + } else { + deriv = internal::wiener7_integrate_cdf< + GradientCalc::OFF, GradientCalc::OFF, GradientCalc::OFF, + GradientCalc::ON>( + [&](auto&&... args) { + return internal::wiener4_ccdf_grad_w(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / ccdf; + } + partials<6>(ops_partials)[i] = deriv; + } + } + if (!is_constant_all::value) { + if (st0_value == 0) { + partials<7>(ops_partials)[i] = 0.0; + } else if (y_value - (t0_value + st0_value) <= 0) { + partials<7>(ops_partials)[i] = -1 / st0_value; + } else { + const auto t0_st0 = t0_value + st0_value; + if (sw_value == 0 && sv_value == 0) { + deriv = internal::estimate_with_err_check<4, 0>( + [](auto&&... args) { return internal::wiener4_ccdf(args...); }, + log_error_derivatives + log(st0_value), y_value - t0_st0, a_value, + v_value, w_value, log_error_absolute - LOG_TWO); + deriv = deriv / st0_value / ccdf - 1 / st0_value; + } else { + const int dim_st = (sv_value != 0) + (sw_value != 0); + const auto new_error = log_error_absolute - LOG_TWO; + const auto& params_st + = std::make_tuple(y_value, a_value, v_value, w_value, t0_st0, + sv_value, sw_value, 0.0, new_error); + deriv = internal::wiener7_integrate_cdf< + GradientCalc::OFF, GradientCalc::OFF, GradientCalc::OFF, + GradientCalc::OFF, GradientCalc::OFF, GradientCalc::OFF>( + [&](auto&&... args) { return internal::wiener4_ccdf(args...); }, + hcubature_err, params_st, dim_st, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2); + deriv = deriv / st0_value / ccdf - 1 / st0_value; + } + partials<7>(ops_partials)[i] = deriv; + } + } + } + return result + ops_partials.build(lccdf); +} +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/wiener_full_lcdf.hpp b/stan/math/prim/prob/wiener_full_lcdf.hpp new file mode 100644 index 00000000000..2a659ad4c0f --- /dev/null +++ b/stan/math/prim/prob/wiener_full_lcdf.hpp @@ -0,0 +1,591 @@ +#ifndef STAN_MATH_PRIM_PROB_WIENER_FULL_LCDF_HPP +#define STAN_MATH_PRIM_PROB_WIENER_FULL_LCDF_HPP + +#include +#include + +namespace stan { +namespace math { +namespace internal { + +/** + * Calculate the derivative of the wiener7 density w.r.t. 'sw' + * + * @tparam T_y type of reaction time + * @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 tolerance + * + * @param y The reaction time in seconds + * @param a The boundary separation + * @param v The drift rate + * @param w The relative starting point + * @param sw The inter-trial variability of the relative starting point + * @param wildcard This parameter space is needed for a functor. Could be + * deleted when another solution is found + * @param log_error The log error tolerance + * @return Gradient with respect to sw + */ +template +inline auto wiener7_cdf_grad_sw(const T_y& y, const T_a& a, const T_v& v, + const T_w& w, const T_sw& sw, T_err log_error) { + auto low = fmax(0.0, w - sw / 2.0); + auto high = fmin(1.0, w + sw / 2.0); + const auto lower_value + = wiener4_distribution(y, a, v, low, log_error); + const auto upper_value + = wiener4_distribution(y, a, v, high, log_error); + return 0.5 * (lower_value + upper_value) / sw; +} + +/** + * Helper function for agnostically calling wiener4 functions + * (to be integrated over) or directly calling wiener7 functions, + * accounting for the different number of arguments. + * + * Specialisation for wiener4 functions + * + * @tparam GradSW Whether the gradient of sw is computed + * @tparam F Type of Gradient/density functor + * @tparam T_y type of reaction time + * @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_cdf(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, log_error); +} + +/** + * Helper function for agnostically calling wiener4 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 reaction time + * @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_cdf(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); +} + +/** + * Implementation function for preparing arguments and functor to be passed + * to the hcubature() function for calculating wiener7 parameters via + * integration + * + * @tparam GradSW Whether the wiener7 gradient function is passed + * @tparam Wiener7FunctorT Type of functor + * @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_cdf(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&& lerr) { + using ret_t + = return_type_t; + scalar_seq_view x_vec(x); + const auto temp = (sv == 0) ? 0 : square(x_vec[0]); + const auto factor = (sv == 0) ? 0 : x_vec[0] / (1 - temp); + const auto new_v = (sv == 0) ? v : v + sv * factor; + const auto new_w + = (sw == 0) ? w : w + sw * (x_vec[sv == 0 ? 0 : 1] - 0.5); + const auto idx + = (sv == 0 && sw == 0) ? 0 : (sv != 0 && sw != 0) ? 2 : 1; + const auto new_t0 = (st0 == 0) ? t0 : t0 + st0 * x_vec[idx]; + if (y - new_t0 <= 0) { + return ret_t(0.0); + } + const auto dist = GradT ? 0 + : wiener4_distribution( + y - new_t0, a, new_v, new_w, lerr); + const auto temp2 = (sv == 0) ? 0 + : -0.5 * square(factor) - LOG_SQRT_PI + - 0.5 * LOG_TWO + log1p(temp) + - 2.0 * log1m(temp); + const auto factor_sv = GradSV ? factor : 1; + const auto factor_sw + = GradSW ? ((sv == 0) ? (x_vec[0] - 0.5) : (x_vec[1] - 0.5)) : 1; + const auto integrand + = Distribution + ? dist + : GradT + ? conditionally_grad_sw_cdf( + wiener7_functor, y - new_t0, a, v, new_w, sv, sw, + lerr) + : factor_sv * factor_sw + * conditionally_grad_sw_cdf( + wiener7_functor, y - new_t0, a, new_v, + new_w, dist, sw, lerr); + return ret_t(integrand * exp(temp2)); + }, + integration_args...); + }; + + return estimate_with_err_check<0, 8, GradW7, GradientCalc::ON>( + functor, hcubature_err, args...); +} +} // namespace internal + +/** \ingroup prob_dists + * Returns the log CDF of the Wiener distribution for a + * (Wiener) drift diffusion model with up to 7 parameters. If containers + * are supplied, returns the log sum of the probabilities. + * See 'wiener_lpdf' for more comprehensive documentation + * If the reaction time goes to infinity, the CDF goes to the probability to + * hit the upper bound (instead of 1, as it is usually the case) + * + * @tparam T_y type of reaction time + * @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 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 log probability or log sum of probabilities 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. + + + **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*. + * https://doi.org/10.3758/s13428-023-02179-1 + * - Henrich, F., & Klauer, K. C. (in press). Modeling Truncated and Censored + * Data With the Diffusion Model in Stan. *Behavior Research Methods*. + * - Linhart, J. M. (2008). Algorithm 885: Computing the Logarithm of the + * Normal Distribution. *ACM Transactions on Mathematical Software*. + * http://doi.acm.org/10.1145/1391989.1391993 + * - 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_lcdf(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-8) { + using ret_t = return_type_t; + 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 internal::GradientCalc; + using T_partials_return + = partials_return_t; + + 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; + + auto y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + auto a_val = to_ref(as_value_column_array_or_scalar(a_ref)); + auto v_val = to_ref(as_value_column_array_or_scalar(v_ref)); + auto w_val = to_ref(as_value_column_array_or_scalar(w_ref)); + auto t0_val = to_ref(as_value_column_array_or_scalar(t0_ref)); + auto sv_val = to_ref(as_value_column_array_or_scalar(sv_ref)); + auto sw_val = to_ref(as_value_column_array_or_scalar(sw_ref)); + auto st0_val = to_ref(as_value_column_array_or_scalar(st0_ref)); + + if (!include_summand::value) { + return ret_t(0); + } + + static constexpr const char* function_name = "wiener_lcdf"; + 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); + 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]) { + 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)) { + std::stringstream msg; + msg << ", but must be smaller than 2*(A-priori bias) = " + << 2.0 * 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)) { + 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()); + } + } + + // for precs. 1e-6, 1e-12, see Hartmann et al. (2021), Henrich et al. (2023) + const T_partials_return log_error_cdf = log(1e-6); // precision for density + const auto error_bound = precision_derivatives; // precision for + // derivatives (controllable by user) + const auto lerror_bound = log(error_bound); + const T_partials_return absolute_error_hcubature = 0.0; + const T_partials_return relative_error_hcubature + = .9 * error_bound; // eps_rel(Integration) + const T_partials_return log_error_absolute = log(1e-12); + const int maximal_evaluations_hcubature = 6000; + T_partials_return lcdf = 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.0; + + // calculate density and partials + for (size_t i = 0; i < N; i++) { + if (sv_vec[i] == 0 && sw_vec[i] == 0 && st0_vec[i] == 0) { + result += wiener_lcdf(y_vec[i], a_vec[i], t0_vec[i], w_vec[i], + v_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 = (sv_value != 0) + (sw_value != 0) + (st0_value != 0); + check_positive(function_name, + "(Inter-trial variability in drift rate) + " + "(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); + for (int i = 0; i < dim; i++) { + xmin[i] = 0; + xmax[i] = 1; + } + if (sv_value != 0) { + xmin[0] = -1; + xmax[0] = 1; + } + 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_cdf + LOG_TWO + 1; + + 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); + + const T_partials_return cdf + = internal::wiener7_integrate_cdf( + [&](auto&&... args) { + return internal::wiener4_distribution(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2); + lcdf += log(cdf); + + hcubature_err + = log_error_absolute - lerror_bound + log(fabs(cdf)) + 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 auto params_dt7 = std::make_tuple( + y_value, a_value, v_value, w_value, t0_value, sv_value, sw_value, + st0_value, log_error_absolute - LOG_TWO - 9 * LOG_TWO); + if (!is_constant_all::value || !is_constant_all::value) { + T_partials_return deriv_t_7 + = internal::wiener7_integrate_cdf< + GradientCalc::OFF, GradientCalc::OFF, GradientCalc::OFF, + GradientCalc::OFF, GradientCalc::ON>( + [&](auto&&... args) { + return internal::wiener5_density(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / cdf; + if (!is_constant_all::value) { + partials<0>(ops_partials)[i] = deriv_t_7; + } + if (!is_constant_all::value) { + partials<2>(ops_partials)[i] = -deriv_t_7; + } + } + T_partials_return deriv; + if (!is_constant_all::value) { + partials<1>(ops_partials)[i] + = internal::wiener7_integrate_cdf( + [&](auto&&... args) { + return internal::wiener4_cdf_grad_a(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / cdf; + } + if (!is_constant_all::value) { + partials<3>(ops_partials)[i] + = internal::wiener7_integrate_cdf( + [&](auto&&... args) { + return internal::wiener4_cdf_grad_w(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / cdf; + } + if (!is_constant_all::value) { + partials<4>(ops_partials)[i] + = internal::wiener7_integrate_cdf( + [&](auto&&... args) { + return internal::wiener4_cdf_grad_v(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / cdf; + } + if (!is_constant_all::value) { + if (sv_value == 0) { + partials<5>(ops_partials)[i] = 0; + } else { + partials<5>(ops_partials)[i] + = internal::wiener7_integrate_cdf< + GradientCalc::OFF, GradientCalc::OFF, GradientCalc::ON>( + [&](auto&&... args) { + return internal::wiener4_cdf_grad_v(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / cdf; + } + } + if (!is_constant_all::value) { + if (sw_value == 0) { + partials<6>(ops_partials)[i] = 0; + } else { + if (st0_value == 0 && sv_value == 0) { + deriv = internal::estimate_with_err_check<5, 0, GradientCalc::OFF, + GradientCalc::ON>( + [](auto&&... args) { + return internal::wiener7_cdf_grad_sw(args...); + }, + hcubature_err, y_value - t0_value, a_value, v_value, w_value, + sw_value, log_error_absolute - LOG_TWO); + deriv = deriv / cdf - 1 / sw_value; + } else { + deriv = internal::wiener7_integrate_cdf< + GradientCalc::OFF, GradientCalc::OFF, GradientCalc::OFF, + GradientCalc::ON>( + [&](auto&&... args) { + return internal::wiener4_cdf_grad_w(args...); + }, + hcubature_err, params, dim, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2) + / cdf; + } + partials<6>(ops_partials)[i] = deriv; + } + } + if (!is_constant_all::value) { + 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 && sv_value == 0) { + deriv = internal::estimate_with_err_check<4, 0, GradientCalc::OFF, + GradientCalc::ON>( + [](auto&&... args) { + return internal::wiener4_distribution( + args...); + }, + lerror_bound + log(st0_value), y_value - t0_st0, a_value, v_value, + w_value, log_error_absolute - LOG_TWO); + deriv = deriv / st0_value / cdf - 1 / st0_value; + } else { + const int dim_st = (sv_value != 0) + (sw_value != 0); + const T_partials_return new_error = log_error_absolute - LOG_TWO; + const auto& params_st + = std::make_tuple(y_value, a_value, v_value, w_value, t0_st0, + sv_value, sw_value, 0, new_error); + deriv = internal::wiener7_integrate_cdf< + GradientCalc::OFF, GradientCalc::OFF, GradientCalc::OFF, + GradientCalc::OFF, GradientCalc::OFF, GradientCalc::OFF>( + [&](auto&&... args) { + return internal::wiener4_distribution( + args...); + }, + hcubature_err, params_st, dim_st, xmin, xmax, + maximal_evaluations_hcubature, absolute_error_hcubature, + relative_error_hcubature / 2); + deriv = deriv / st0_value / cdf - 1 / st0_value; + } + partials<7>(ops_partials)[i] = deriv; + } + } + } + return result + ops_partials.build(lcdf); +} +} // 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 index 0097163650e..6d8dc650949 100644 --- a/stan/math/prim/prob/wiener_full_lpdf.hpp +++ b/stan/math/prim/prob/wiener_full_lpdf.hpp @@ -12,7 +12,7 @@ namespace internal { /** * Calculate the derivative of the wiener7 density w.r.t. 'sw' * - * @tparam T_y type of scalar variable + * @tparam T_y type of reaction time * @tparam T_a type of boundary separation * @tparam T_v type of drift rate * @tparam T_w type of relative starting point @@ -20,14 +20,14 @@ namespace internal { * @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 y 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 + * @return Gradient with respect to sw */ template @@ -52,7 +52,7 @@ inline auto wiener7_grad_sw(const T_y& y, const T_a& a, const T_v& v, * * @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_y type of reaction time * @tparam T_a type of boundary separation * @tparam T_v type of drift rate * @tparam T_w type of relative starting point @@ -76,7 +76,7 @@ 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) { + T_err log_error) { return functor(y_diff, a, v, w, sv, log_error); } @@ -89,7 +89,7 @@ inline auto conditionally_grad_sw(F&& functor, T_y&& y_diff, T_a&& a, T_v&& v, * * @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_y type of reaction time * @tparam T_a type of boundary separation * @tparam T_v type of drift rate * @tparam T_w type of relative starting point @@ -113,7 +113,7 @@ 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) { + T_err log_error) { return functor(y_diff, a, v, w, sv, sw, log_error); } @@ -212,7 +212,7 @@ inline auto wiener7_integrate(const Wiener7FunctorT& wiener7_functor, * * See \b Details below for more details on how to use \c wiener_lpdf(). * - * @tparam T_y type of scalar + * @tparam T_y type of reaction time * @tparam T_a type of boundary separation * @tparam T_t0 type of non-decision time * @tparam T_w type of relative starting point @@ -221,7 +221,7 @@ inline auto wiener7_integrate(const Wiener7FunctorT& wiener7_functor, * @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 y The reaction time in seconds * @param a The boundary separation * @param t0 The non-decision time * @param w The relative starting point @@ -322,11 +322,6 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, 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_t; using T_a_ref = ref_type_t; using T_v_ref = ref_type_t; @@ -335,18 +330,10 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, using T_sv_ref = ref_type_t; using T_sw_ref = ref_type_t; using T_st0_ref = ref_type_t; - + using internal::GradientCalc; 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; @@ -356,14 +343,27 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, 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)); + auto y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + auto a_val = to_ref(as_value_column_array_or_scalar(a_ref)); + auto v_val = to_ref(as_value_column_array_or_scalar(v_ref)); + auto w_val = to_ref(as_value_column_array_or_scalar(w_ref)); + auto t0_val = to_ref(as_value_column_array_or_scalar(t0_ref)); + auto sv_val = to_ref(as_value_column_array_or_scalar(sv_ref)); + auto sw_val = to_ref(as_value_column_array_or_scalar(sw_ref)); + auto st0_val = to_ref(as_value_column_array_or_scalar(st0_ref)); + + if (!include_summand::value) { + return ret_t(0); + } + + 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); 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); @@ -412,7 +412,7 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, [&]() STAN_COLD_PATH { std::stringstream msg; msg << ", but must be smaller than 2*(A-priori bias) = " - << 2 * w_vec[i]; + << 2.0 * w_vec[i]; std::string msg_str(msg.str()); throw_domain_error(function_name, "Inter-trial variability in A-priori bias", @@ -423,7 +423,7 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, [&]() STAN_COLD_PATH { std::stringstream msg; msg << ", but must be smaller than 2*(1-A-priori bias) = " - << 2 * (1 - w_vec[i]); + << 2.0 * (1.0 - w_vec[i]); std::string msg_str(msg.str()); throw_domain_error(function_name, "Inter-trial variability in A-priori bias", @@ -431,6 +431,8 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, }(); } } + + // for precs. 1e-6, 1e-12, see Hartmann et al. (2021), Henrich et al. (2023) // precision for density const T_partials_return log_error_density = log(1e-6); // precision for derivatives (controllable by user) @@ -478,7 +480,6 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0, 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); diff --git a/stan/math/prim/prob/wiener_lpdf.hpp b/stan/math/prim/prob/wiener_lpdf.hpp index 5c2f113a07e..05c785691cf 100644 --- a/stan/math/prim/prob/wiener_lpdf.hpp +++ b/stan/math/prim/prob/wiener_lpdf.hpp @@ -76,7 +76,7 @@ namespace math { */ template -inline return_type_t wiener_lpdf( +return_type_t wiener_lpdf( const T_y& y, const T_alpha& alpha, const T_tau& tau, const T_beta& beta, const T_delta& delta) { using T_return_type = return_type_t; diff --git a/test/unit/math/mix/prob/util.hpp b/test/unit/math/mix/prob/util.hpp index 5b1dea8ad59..e3536ca986c 100644 --- a/test/unit/math/mix/prob/util.hpp +++ b/test/unit/math/mix/prob/util.hpp @@ -45,4 +45,21 @@ class Wiener5MixArgs : public ::testing::Test { 0.2, 0.02, 9.0, 0.01, 0.2, 4, 0.2, 1e14}}; }; +class Wiener4MixArgs : 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}}; +}; + #endif diff --git a/test/unit/math/mix/prob/wiener4_lccdf_0_test.cpp b/test/unit/math/mix/prob/wiener4_lccdf_0_test.cpp new file mode 100644 index 00000000000..d32b38c7fa6 --- /dev/null +++ b/test/unit/math/mix/prob/wiener4_lccdf_0_test.cpp @@ -0,0 +1,33 @@ +#include +#include + +TEST(mathMixDouble, wiener4_lccdf) { + double y = 1.0; + double a = 2.5; + double t0 = 0.2; + double w = 0.5; + double v = 1.5; + stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); +} + +TEST(mathMixVar, wiener4_lccdf) { + using stan::math::var; + var y = 1.0; + var a = 2.5; + var t0 = 0.2; + var w = 0.5; + var v = 1.5; + stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); +} + +TEST(mathMixFVar, wiener4_lccdf) { + using stan::math::fvar; + using stan::math::var; + fvar y = 1.0; + fvar a = 2.5; + fvar t0 = 0.2; + fvar w = 0.5; + fvar v = 1.5; + double error = 1e-4; + stan::math::wiener_lccdf(y, a, t0, w, v, error); +} diff --git a/test/unit/math/mix/prob/wiener4_lccdf_1_test.cpp b/test/unit/math/mix/prob/wiener4_lccdf_1_test.cpp new file mode 100644 index 00000000000..26425508337 --- /dev/null +++ b/test/unit/math/mix/prob/wiener4_lccdf_1_test.cpp @@ -0,0 +1,48 @@ +#include +#include +#include + +TEST_F(Wiener4MixArgs, wiener_lccdf_y_a_t0) { + auto f_y_a_t0 = [](const auto& w, const auto& v) { + return [&w, &v](const auto& y, const auto& a, const auto& t0) { + return stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_a_t0(w, v), y, a, t0); +} + +TEST_F(Wiener4MixArgs, wiener_lccdf_y_a_w) { + auto f_y_a_w = [](const auto& t0, const auto& v) { + return [&t0, &v](const auto& y, const auto& a, const auto& w) { + return stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_a_w(t0, v), y, a, w); +} + +TEST_F(Wiener4MixArgs, wiener_lccdf_y_a_v) { + auto f_y_a_v = [](const auto& t0, const auto& w) { + return [&t0, &w](const auto& y, const auto& a, const auto& v) { + return stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_a_v(t0, w), y, a, v); +} + +TEST_F(Wiener4MixArgs, wiener_lccdf_y_t0_w) { + auto f_y_t0_w = [](const auto& a, const auto& v) { + return [&a, &v](const auto& y, const auto& t0, const auto& w) { + return stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_t0_w(a, v), y, t0, w); +} + +TEST_F(Wiener4MixArgs, wiener_lccdf_y_t0_v) { + auto f_y_t0_v = [](const auto& a, const auto& w) { + return [&a, &w](const auto& y, const auto& t0, const auto& v) { + return stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_t0_v(a, w), y, t0, v); +} diff --git a/test/unit/math/mix/prob/wiener4_lccdf_2_test.cpp b/test/unit/math/mix/prob/wiener4_lccdf_2_test.cpp new file mode 100644 index 00000000000..579bf94189a --- /dev/null +++ b/test/unit/math/mix/prob/wiener4_lccdf_2_test.cpp @@ -0,0 +1,48 @@ +#include +#include +#include + +TEST_F(Wiener4MixArgs, wiener_lccdf_y_w_v) { + auto f_y_w_v = [](const auto& a, const auto& t0) { + return [&a, &t0](const auto& y, const auto& w, const auto& v) { + return stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_w_v(a, t0), y, w, v); +} + +TEST_F(Wiener4MixArgs, wiener_lccdf_a_t0_w) { + auto f_a_t0_w = [](const auto& y, const auto& v) { + return [&y, &v](const auto& a, const auto& t0, const auto& w) { + return stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_a_t0_w(y, v), a, t0, w); +} + +TEST_F(Wiener4MixArgs, wiener_lccdf_a_t0_v) { + auto f_a_t0_v = [](const auto& y, const auto& w) { + return [&y, &w](const auto& a, const auto& t0, const auto& v) { + return stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_a_t0_v(y, w), a, t0, v); +} + +TEST_F(Wiener4MixArgs, wiener_lccdf_a_w_v) { + auto f_a_w_v = [](const auto& y, const auto& t0) { + return [&y, &t0](const auto& a, const auto& w, const auto& v) { + return stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_a_w_v(y, t0), a, w, v); +} + +TEST_F(Wiener4MixArgs, wiener_lccdf_t0_w_v) { + auto f_t0_w_v = [](const auto& y, const auto& a) { + return [&y, &a](const auto& t0, const auto& w, const auto& v) { + return stan::math::wiener_lccdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_t0_w_v(y, a), t0, w, v); +} diff --git a/test/unit/math/mix/prob/wiener4_lcdf_0_test.cpp b/test/unit/math/mix/prob/wiener4_lcdf_0_test.cpp new file mode 100644 index 00000000000..23a7d5d3ae1 --- /dev/null +++ b/test/unit/math/mix/prob/wiener4_lcdf_0_test.cpp @@ -0,0 +1,33 @@ +#include +#include + +TEST(mathMixDouble, wiener4_lcdf) { + double y = 1.0; + double a = 2.5; + double t0 = 0.2; + double w = 0.5; + double v = 1.5; + stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); +} + +TEST(mathMixVar, wiener4_lcdf) { + using stan::math::var; + var y = 1.0; + var a = 2.5; + var t0 = 0.2; + var w = 0.5; + var v = 1.5; + stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); +} + +TEST(mathMixFVar, wiener4_lcdf) { + using stan::math::fvar; + using stan::math::var; + fvar y = 1.0; + fvar a = 2.5; + fvar t0 = 0.2; + fvar w = 0.5; + fvar v = 1.5; + double error = 1e-4; + stan::math::wiener_lcdf(y, a, t0, w, v, error); +} diff --git a/test/unit/math/mix/prob/wiener4_lcdf_1_test.cpp b/test/unit/math/mix/prob/wiener4_lcdf_1_test.cpp new file mode 100644 index 00000000000..2cda5b5899c --- /dev/null +++ b/test/unit/math/mix/prob/wiener4_lcdf_1_test.cpp @@ -0,0 +1,48 @@ +#include +#include +#include + +TEST_F(Wiener4MixArgs, wiener_lcdf_y_a_t0) { + auto f_y_a_t0 = [](const auto& w, const auto& v) { + return [&w, &v](const auto& y, const auto& a, const auto& t0) { + return stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_a_t0(w, v), y, a, t0); +} + +TEST_F(Wiener4MixArgs, wiener_lcdf_y_a_w) { + auto f_y_a_w = [](const auto& t0, const auto& v) { + return [&t0, &v](const auto& y, const auto& a, const auto& w) { + return stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_a_w(t0, v), y, a, w); +} + +TEST_F(Wiener4MixArgs, wiener_lcdf_y_a_v) { + auto f_y_a_v = [](const auto& t0, const auto& w) { + return [&t0, &w](const auto& y, const auto& a, const auto& v) { + return stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_a_v(t0, w), y, a, v); +} + +TEST_F(Wiener4MixArgs, wiener_lcdf_y_t0_w) { + auto f_y_t0_w = [](const auto& a, const auto& v) { + return [&a, &v](const auto& y, const auto& t0, const auto& w) { + return stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_t0_w(a, v), y, t0, w); +} + +TEST_F(Wiener4MixArgs, wiener_lcdf_y_t0_v) { + auto f_y_t0_v = [](const auto& a, const auto& w) { + return [&a, &w](const auto& y, const auto& t0, const auto& v) { + return stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_t0_v(a, w), y, t0, v); +} diff --git a/test/unit/math/mix/prob/wiener4_lcdf_2_test.cpp b/test/unit/math/mix/prob/wiener4_lcdf_2_test.cpp new file mode 100644 index 00000000000..2677da3b92b --- /dev/null +++ b/test/unit/math/mix/prob/wiener4_lcdf_2_test.cpp @@ -0,0 +1,48 @@ +#include +#include +#include + +TEST_F(Wiener4MixArgs, wiener_lcdf_y_w_v) { + auto f_y_w_v = [](const auto& a, const auto& t0) { + return [&a, &t0](const auto& y, const auto& w, const auto& v) { + return stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_y_w_v(a, t0), y, w, v); +} + +TEST_F(Wiener4MixArgs, wiener_lcdf_a_t0_w) { + auto f_a_t0_w = [](const auto& y, const auto& v) { + return [&y, &v](const auto& a, const auto& t0, const auto& w) { + return stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_a_t0_w(y, v), a, t0, w); +} + +TEST_F(Wiener4MixArgs, wiener_lcdf_a_t0_v) { + auto f_a_t0_v = [](const auto& y, const auto& w) { + return [&y, &w](const auto& a, const auto& t0, const auto& v) { + return stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_a_t0_v(y, w), a, t0, v); +} + +TEST_F(Wiener4MixArgs, wiener_lcdf_a_w_v) { + auto f_a_w_v = [](const auto& y, const auto& t0) { + return [&y, &t0](const auto& a, const auto& w, const auto& v) { + return stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_a_w_v(y, t0), a, w, v); +} + +TEST_F(Wiener4MixArgs, wiener_lcdf_t0_w_v) { + auto f_t0_w_v = [](const auto& y, const auto& a) { + return [&y, &a](const auto& t0, const auto& w, const auto& v) { + return stan::math::wiener_lcdf(y, a, t0, w, v, 1e-4); + }; + }; + stan::test::expect_ad(f_t0_w_v(y, a), t0, w, v); +} diff --git a/test/unit/math/mix/prob/wiener5_lpdf_0_test.cpp b/test/unit/math/mix/prob/wiener5_lpdf_0_test.cpp index 24901ac760b..4cc6ba01367 100644 --- a/test/unit/math/mix/prob/wiener5_lpdf_0_test.cpp +++ b/test/unit/math/mix/prob/wiener5_lpdf_0_test.cpp @@ -3,7 +3,7 @@ TEST(mathMixDouble5, wiener_lpdf) { double y = 1.0; - double a = 2.0; + double a = 2.5; double t0 = 0.2; double w = 0.5; double v = 1.5; @@ -14,7 +14,7 @@ TEST(mathMixDouble5, wiener_lpdf) { TEST(mathMixVar5, wiener_lpdf) { using stan::math::var; var y = 1.0; - var a = 2.0; + var a = 2.5; var t0 = 0.2; var w = 0.5; var v = 1.5; @@ -26,7 +26,7 @@ TEST(mathMixFVar5, wiener_lpdf) { using stan::math::fvar; using stan::math::var; fvar y = 1.0; - fvar a = 2.0; + fvar a = 2.5; fvar t0 = 0.2; fvar w = 0.5; fvar v = 1.5; diff --git a/test/unit/math/mix/prob/wiener_full_lccdf_0_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_0_test.cpp new file mode 100644 index 00000000000..569c677d043 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_0_test.cpp @@ -0,0 +1,53 @@ +#include +#include +#include + +TEST(mathMixDouble, wiener_lccdf) { + 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_lccdf(y, a, t0, w, v, sv, sw, st0, 1e-4); +} + +TEST(mathMixVar, wiener_lccdf) { + 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_lccdf(y, a, t0, w, v, sv, sw, st0); +} + +TEST(mathMixFVar, wiener_lccdf) { + 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_lccdf(y, a, t0, w, v, sv, sw, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_10_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_10_test.cpp new file mode 100644 index 00000000000..e7eca727b91 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_10_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_11_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_11_test.cpp new file mode 100644 index 00000000000..11d7d7c5c27 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_11_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_1_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_1_test.cpp new file mode 100644 index 00000000000..771e20d4c57 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_1_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_2_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_2_test.cpp new file mode 100644 index 00000000000..76aabef0e32 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_2_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_3_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_3_test.cpp new file mode 100644 index 00000000000..870d294bfc0 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_3_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_4_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_4_test.cpp new file mode 100644 index 00000000000..8277f69f45b --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_4_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_5_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_5_test.cpp new file mode 100644 index 00000000000..f24cef06a49 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_5_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_6_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_6_test.cpp new file mode 100644 index 00000000000..f98faacfa6c --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_6_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_7_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_7_test.cpp new file mode 100644 index 00000000000..f7e38df10b1 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_7_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_8_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_8_test.cpp new file mode 100644 index 00000000000..39128239a62 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_8_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_9_test.cpp b/test/unit/math/mix/prob/wiener_full_lccdf_9_test.cpp new file mode 100644 index 00000000000..2a8382f3693 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lccdf_9_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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_lccdf_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_lccdf(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/mix/prob/wiener_full_lcdf_0_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_0_test.cpp new file mode 100644 index 00000000000..0dd7b4306f2 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_0_test.cpp @@ -0,0 +1,53 @@ +#include +#include +#include + +TEST(mathMixDouble, wiener_lcdf) { + 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_lcdf(y, a, t0, w, v, sv, sw, st0, 1e-4); +} + +TEST(mathMixVar, wiener_lcdf) { + 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_lcdf(y, a, t0, w, v, sv, sw, st0); +} + +TEST(mathMixFVar, wiener_lcdf) { + 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_lcdf(y, a, t0, w, v, sv, sw, st0); +} + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_10_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_10_test.cpp new file mode 100644 index 00000000000..5aa2c43df18 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_10_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_11_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_11_test.cpp new file mode 100644 index 00000000000..765f21df6dd --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_11_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_1_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_1_test.cpp new file mode 100644 index 00000000000..27e7c7b7cf3 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_1_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_2_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_2_test.cpp new file mode 100644 index 00000000000..a628fa7da98 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_2_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_3_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_3_test.cpp new file mode 100644 index 00000000000..936ea3be317 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_3_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_4_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_4_test.cpp new file mode 100644 index 00000000000..b7e0699eb47 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_4_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_5_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_5_test.cpp new file mode 100644 index 00000000000..10adf45a2bc --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_5_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_6_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_6_test.cpp new file mode 100644 index 00000000000..c9f8833d902 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_6_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_7_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_7_test.cpp new file mode 100644 index 00000000000..fe6d6061037 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_7_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_8_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_8_test.cpp new file mode 100644 index 00000000000..54350c2dc46 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_8_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_9_test.cpp b/test/unit/math/mix/prob/wiener_full_lcdf_9_test.cpp new file mode 100644 index 00000000000..566b5ff4c83 --- /dev/null +++ b/test/unit/math/mix/prob/wiener_full_lcdf_9_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST_F(Wiener7MixArgs, wiener_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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_lcdf_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_lcdf(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/prob/wiener_full_lccdf_test.cpp b/test/unit/math/prim/prob/wiener_full_lccdf_test.cpp new file mode 100644 index 00000000000..c1fdf33505d --- /dev/null +++ b/test/unit/math/prim/prob/wiener_full_lccdf_test.cpp @@ -0,0 +1,414 @@ +#include +#include +#include +#include + +#include +#include + +TEST(mathPrimScalProbWienerFullLccdfScal, valid) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + 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_lccdf(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_lccdf(rt, a, t0, w, v, sv, sw, st0)); +} + +// rt +TEST(mathPrimScalProbWienerFullLccdfScal, invalid_rt) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(0, a, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lccdf(-1, a, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lccdf(INFTY, a, t0, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(-INFTY, a, t0, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(NAN, a, t0, w, v, sv, sw, st0), std::domain_error); +} + +// a +TEST(mathPrimScalProbWienerFullLccdfScal, invalid_a) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(rt, 0, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, -1, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, INFTY, t0, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, -INFTY, t0, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, NAN, t0, w, v, sv, sw, st0), std::domain_error); +} + +// v +TEST(mathPrimScalProbWienerFullLccdfScal, invalid_v) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + 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_lccdf(rt, a, t0, w, INFTY, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, -INFTY, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, NAN, sv, sw, st0), std::domain_error); +} + +// w +TEST(mathPrimScalProbWienerFullLccdfScal, invalid_w) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, a = 1, v = -1, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(rt, a, t0, -0.1, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, 0, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, 1, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, 1.1, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, INFTY, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, -INFTY, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, NAN, v, sv, sw, st0), std::domain_error); +} + +// t0 +TEST(mathPrimScalProbWienerFullLccdfScal, invalid_t0) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, a = 1, v = -1, w = 0.5, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(rt, a, 2, w, v, sv, sw, st0), + std::domain_error); // rt must be greater than t0 + EXPECT_THROW(wiener_lccdf(rt, a, -1, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, INFTY, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, -INFTY, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, NAN, w, v, sv, sw, st0), std::domain_error); +} + +// sv +TEST(mathPrimScalProbWienerFullLccdfScal, invalid_sv) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, -1, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, INFTY, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, -INFTY, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, NAN, sw, st0), std::domain_error); +} + +// sw +TEST(mathPrimScalProbWienerFullLccdfScal, invalid_sw) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, -1, st0), std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, 0.8, v, sv, 0.5, st0), + std::domain_error); // sw must be smaller than 2*(1-w) + EXPECT_THROW(wiener_lccdf(rt, a, t0, 0.3, v, sv, 0.7, st0), + std::domain_error); // sw must be smaller than 2*w + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, INFTY, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, -INFTY, st0), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, NAN, st0), std::domain_error); +} + +// st0 +TEST(mathPrimScalProbWienerFullLccdfScal, invalid_st0) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2; + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, sw, -1), std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, sw, INFTY), std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, sw, -INFTY), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, sw, NAN), std::domain_error); +} + +TEST(mathPrimScalProbWienerFullLccdfPrecScal, valid) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + 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_lccdf(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_lccdf(rt, a, t0, w, v, sv, sw, st0, 1e-4)); +} + +// rt +TEST(mathPrimScalProbWienerFullLccdfPrecScal, invalid_rt) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(0, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(-1, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(INFTY, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(-INFTY, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(NAN, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// a +TEST(mathPrimScalProbWienerFullLccdfPrecScal, invalid_a) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(rt, 0, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, -1, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, INFTY, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, -INFTY, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, NAN, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// v +TEST(mathPrimScalProbWienerFullLccdfPrecScal, invalid_v) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + 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_lccdf(rt, a, t0, w, INFTY, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, -INFTY, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, NAN, sv, sw, st0, 1e-4), + std::domain_error); +} + +// w +TEST(mathPrimScalProbWienerFullLccdfPrecScal, invalid_w) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, a = 1, v = -1, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(rt, a, t0, -0.1, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, 0, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, 1, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, 1.1, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, INFTY, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, -INFTY, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, NAN, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// t0 +TEST(mathPrimScalProbWienerFullLccdfPrecScal, invalid_t0) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, a = 1, v = -1, w = 0.5, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(rt, a, 2, w, v, sv, sw, st0, 1e-4), + std::domain_error); // rt must be greater than t0 + EXPECT_THROW(wiener_lccdf(rt, a, -1, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, INFTY, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, -INFTY, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, NAN, w, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// sv +TEST(mathPrimScalProbWienerFullLccdfPrecScal, invalid_sv) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, -1, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, INFTY, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, -INFTY, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, NAN, sw, st0, 1e-4), + std::domain_error); +} + +// sw +TEST(mathPrimScalProbWienerFullLccdfPrecScal, invalid_sw) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, -1, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(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_lccdf(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_lccdf(rt, a, t0, w, v, sv, INFTY, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, -INFTY, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, NAN, st0, 1e-4), + std::domain_error); +} + +// st0 +TEST(mathPrimScalProbWienerFullLccdfPrecScal, invalid_st0) { + using stan::math::INFTY; + using stan::math::wiener_lccdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2; + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, sw, -1, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, sw, INFTY, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, sw, -INFTY, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lccdf(rt, a, t0, w, v, sv, sw, NAN, 1e-4), + std::domain_error); +} + +TEST(mathPrimCorrectValues, wiener_lccdf) { + /* Test concrete values. True values are computed in R using the R-package + * WienR and the function WienerCDF() with its partial derivatives: + * ccdf = cdf(big_value) - cdf + * lccdf = log(cdf(big_value) - cdf) + * lccdf' = ccdf'/ccdf = (cdf(big_value)'-cdf')/ccdf + */ + std::vector y_vec + = {.2, .3, 2.2, .1, 2, 1.7, 3, 2.5, 3.7, 1.8, .9, 1.5, 2}; + std::vector a_vec + = {3, 1.7, 2.4, 2, 1.2, 1.5, 4, 1.5, 1.7, 6, 1, 2.5, 2}; + std::vector t0_vec + = {0, 0, 0, 0, 0.2, 0.3, 0.1, 0.7, 1, 0.6, .1, .2, .2}; + std::vector w_vec + = {.7, .92, .9, .4, 0.5, 0.3, 0.45, 0.6, 0.8, 0.35, .4, .5, .6}; + std::vector v_vec + = {-1, -7.3, -4.9, 2.5, 1, 0.5, -0.7, 0.9, -1.4, 2, .8, -1, -.9}; + std::vector sv_vec = {0, 0, 0, 0, 0.5, 0, 0, 0.7, 0.8, 0, .4, .5, .2}; + std::vector sw_vec = {0, 0, 0, 0, 0, 0.1, 0, 0.2, 0, 0.3, .2, 0, .2}; + std::vector st0_vec + = {0, 0, 0, 0, 0, 0, 0.1, 0, 0.1, 0.1, .2, .3, .25}; + + std::vector true_lccdf + = {-1.92033254902358, -13.5626831740558, -33.8122958054115, + -0.0207304574971297, -7.08470781580343, -3.44486241813253, + -4.18311596420185, -4.95366166391584, -7.23109216780916, + -0.125273337362733, -3.79815495569855, -2.96111699037856, + -4.14929265642466}; + std::vector true_grad_y + = {-1.33001519743633, -31.1083518155787, -12.9346392526053, + -0.169566706752058, -3.7986585460178, -2.31784674820993, + -0.536530793068282, -2.48164917576694, -1.89222689054599, + -0.484444809563098, -5.31460285637481, -0.975151391419348, + -1.57127046017878}; + std::vector true_grad_a + = {-0.432928212420575, -0.0286062052702194, 0.308724832214765, + 0.05132866977503, 10.5080654585197, 4.37077488076664, + -0.193420148034756, 5.36345669388884, 5.14057131554871, + 0.11692460737807, 7.07166462218294, 0.365404340761721, + 1.53902298927608}; + std::vector true_grad_t0 = { + 1.33001519743633, 31.1083518155787, 12.9346392526053, 0.169566706752058, + 3.7986585460178, 2.31784674820993, 0.536530793068282, 2.48164917576694, + 1.89222689054599, 0.484444809563098, 5.31460285637481, 0.975151391419348, + 1.57127046017878}; + std::vector true_grad_w + = {4.61303457108487, 0.611152450869066, 1.93288590604027, + 0.139815038667115, -0.911757809733564, 1.54221437873658, + 3.63420918438666, -1.85952752467816, -3.65499582974638, + -1.12442042700282, 0.262989430764412, 1.3830938669116, + 0.518345884055846}; + std::vector true_grad_v + = {1.83079071623625, 2.56351255240926, 11.4026845637584, + 0.0274975992798594, -0.960807197217193, 0.132382852324283, + 5.47350811619098, -0.652456751242452, 1.57258003231325, + -0.207055671666776, -0.0823432112415456, 2.18370562291933, + 2.62134233757004}; + std::vector true_grad_sv = {0, + 0, + 0, + 0, + -0.209370036446517, + 0, + 0, + -0.430772737121245, + 1.14730002231356, + 0, + -0.293744554675593, + 1.73837609123659, + 0.997525744657476}; + std::vector true_grad_sw = {0, + 0, + 0, + 0, + 0, + -0.106375521997502, + 0, + -0.114806591325169, + 0, + -0.201319331904563, + -0.180444908443688, + 0, + -0.17718676934793}; + std::vector true_grad_st0 = {0, + 0, + 0, + 0, + 0, + 0, + 0.270517553006863, + 0, + 0.976409725660212, + 0.235042776747307, + 3.12081986216078, + 0.512275314141889, + 0.837649445725816}; + + 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 lccdf = stan::math::wiener_lccdf(y, a, t0, w, v, sv, sw, st0); + lccdf.grad(); + EXPECT_NEAR(lccdf.val(), true_lccdf[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); + } +} diff --git a/test/unit/math/prim/prob/wiener_full_lcdf_test.cpp b/test/unit/math/prim/prob/wiener_full_lcdf_test.cpp new file mode 100644 index 00000000000..bf18ff38150 --- /dev/null +++ b/test/unit/math/prim/prob/wiener_full_lcdf_test.cpp @@ -0,0 +1,405 @@ +#include +#include +#include +#include + +#include +#include + +TEST(mathPrimScalProbWienerFullLcdfScal, valid) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + 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_lcdf(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_lcdf(rt, a, t0, w, v, sv, sw, st0)); +} + +// rt +TEST(mathPrimScalProbWienerFullLcdfScal, invalid_rt) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(0, a, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(-1, a, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(INFTY, a, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(-INFTY, a, t0, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lcdf(NAN, a, t0, w, v, sv, sw, st0), std::domain_error); +} + +// a +TEST(mathPrimScalProbWienerFullLcdfScal, invalid_a) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(rt, 0, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, -1, t0, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, INFTY, t0, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, -INFTY, t0, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, NAN, t0, w, v, sv, sw, st0), std::domain_error); +} + +// v +TEST(mathPrimScalProbWienerFullLcdfScal, invalid_v) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + 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_lcdf(rt, a, t0, w, INFTY, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, -INFTY, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, NAN, sv, sw, st0), std::domain_error); +} + +// w +TEST(mathPrimScalProbWienerFullLcdfScal, invalid_w) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, a = 1, v = -1, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(rt, a, t0, -0.1, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, 0, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, 1, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, 1.1, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, INFTY, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, -INFTY, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, NAN, v, sv, sw, st0), std::domain_error); +} + +// t0 +TEST(mathPrimScalProbWienerFullLcdfScal, invalid_t0) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, a = 1, v = -1, w = 0.5, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(rt, a, 2, w, v, sv, sw, st0), + std::domain_error); // rt must be greater than t0 + EXPECT_THROW(wiener_lcdf(rt, a, -1, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, INFTY, w, v, sv, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, -INFTY, w, v, sv, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, NAN, w, v, sv, sw, st0), std::domain_error); +} + +// sv +TEST(mathPrimScalProbWienerFullLcdfScal, invalid_sv) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, -1, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, INFTY, sw, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, -INFTY, sw, st0), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, NAN, sw, st0), std::domain_error); +} + +// sw +TEST(mathPrimScalProbWienerFullLcdfScal, invalid_sw) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, -1, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, 0.8, v, sv, 0.5, st0), + std::domain_error); // sw must be smaller than 2*(1-w) + EXPECT_THROW(wiener_lcdf(rt, a, t0, 0.3, v, sv, 0.7, st0), + std::domain_error); // sw must be smaller than 2*w + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, INFTY, st0), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, -INFTY, st0), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, NAN, st0), std::domain_error); +} + +// st0 +TEST(mathPrimScalProbWienerFullLcdfScal, invalid_st0) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2; + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, sw, -1), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, sw, INFTY), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, sw, -INFTY), std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, sw, NAN), std::domain_error); +} + +TEST(mathPrimScalProbWienerFullLcdfPrecScal, valid) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + 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_lcdf(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_lcdf(rt, a, t0, w, v, sv, sw, st0, 1e-4)); +} + +// rt +TEST(mathPrimScalProbWienerFullLcdfPrecScal, invalid_rt) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(0, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(-1, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(INFTY, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(-INFTY, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(NAN, a, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// a +TEST(mathPrimScalProbWienerFullLcdfPrecScal, invalid_a) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(rt, 0, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, -1, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, INFTY, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, -INFTY, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, NAN, t0, w, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// v +TEST(mathPrimScalProbWienerFullLcdfPrecScal, invalid_v) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + 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_lcdf(rt, a, t0, w, INFTY, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, -INFTY, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, NAN, sv, sw, st0, 1e-4), + std::domain_error); +} + +// w +TEST(mathPrimScalProbWienerFullLcdfPrecScal, invalid_w) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, a = 1, v = -1, t0 = 0.1, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(rt, a, t0, -0.1, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, 0, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, 1, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, 1.1, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, INFTY, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, -INFTY, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, NAN, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// t0 +TEST(mathPrimScalProbWienerFullLcdfPrecScal, invalid_t0) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, a = 1, v = -1, w = 0.5, sv = 0.2, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(rt, a, 2, w, v, sv, sw, st0, 1e-4), + std::domain_error); // rt must be greater than t0 + EXPECT_THROW(wiener_lcdf(rt, a, -1, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, INFTY, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, -INFTY, w, v, sv, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, NAN, w, v, sv, sw, st0, 1e-4), + std::domain_error); +} + +// sv +TEST(mathPrimScalProbWienerFullLcdfPrecScal, invalid_sv) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sw = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, -1, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, INFTY, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, -INFTY, sw, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, NAN, sw, st0, 1e-4), + std::domain_error); +} + +// sw +TEST(mathPrimScalProbWienerFullLcdfPrecScal, invalid_sw) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, st0 = 0.1; + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, -1, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(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_lcdf(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_lcdf(rt, a, t0, w, v, sv, INFTY, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, -INFTY, st0, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, NAN, st0, 1e-4), + std::domain_error); +} + +// st0 +TEST(mathPrimScalProbWienerFullLcdfPrecScal, invalid_st0) { + using stan::math::INFTY; + using stan::math::wiener_lcdf; + double rt = 1, a = 1, v = -1, w = 0.5, t0 = 0.1, sv = 0.2, sw = 0.2; + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, sw, -1, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, sw, INFTY, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, sw, -INFTY, 1e-4), + std::domain_error); + EXPECT_THROW(wiener_lcdf(rt, a, t0, w, v, sv, sw, NAN, 1e-4), + std::domain_error); +} + +TEST(mathPrimCorrectValues, wiener_lcdf) { + /* Test concrete values. True values are computed in R using the R-package + * WienR and the function WienerCDF() with its partial derivatives + */ + std::vector y_vec + = {.2, .3, 2.2, .1, 2, 1.7, 3, 2.5, 3.7, 1.8, .9, 1.5, 2}; + std::vector a_vec + = {3, 1.7, 2.4, 2, 1.2, 1.5, 4, 1.5, 1.7, 6, 1, 2.5, 2}; + std::vector t0_vec + = {0, 0, 0, 0, 0.2, 0.3, 0.1, 0.7, 1, 0.6, .1, .2, .2}; + std::vector w_vec + = {.7, .92, .9, .4, 0.5, 0.3, 0.45, 0.6, 0.8, 0.35, .4, .5, .6}; + std::vector v_vec + = {-1, -7.3, -4.9, 2.5, 1, 0.5, -0.7, 0.9, -1.4, 2, .8, -1, -.9}; + std::vector sv_vec = {0, 0, 0, 0, 0.5, 0, 0, 0.7, 0.8, 0, .4, .5, .2}; + std::vector sw_vec = {0, 0, 0, 0, 0, 0.1, 0, 0.2, 0, 0.3, .2, 0, .2}; + std::vector st0_vec + = {0, 0, 0, 0, 0, 0, 0.1, 0, 0.1, 0.1, .2, .3, .25}; + + std::vector true_lcdf + = {-4.09435375436108, -1.98560937871925, -2.35200000058001, + -6.09860875887167, -0.28490338793642, -0.835235451149805, + -3.60573984770431, -0.213132053911535, -0.887413009805869, + -2.14898785081967, -0.574404607969948, -2.67530409425456, + -1.55387294861195}; + std::vector true_grad_y + = {11.6955868009794, 0.000291754616022926, 2.81009792107004e-13, + 73.9485382296513, 0.00423167900126897, 0.170505395422671, + 0.301191971384078, 0.0216747718871863, 0.00332619998061515, + 3.66549265245825, 0.211550389634303, 0.732731869409894, + 0.117239643328064}; + std::vector true_grad_a + = {-1.90909592842206, -1.1680106846017, -0.979999994944343, + -6.17343628534049, 0.189129525101159, -0.0640548615428843, + -1.01571721473935, 0.0439926894045548, -0.405445995424302, + -0.875235676490553, 0.0169546689079121, -1.10780211022062, + -0.724594716693696}; + std::vector true_grad_t0 + = {-11.6955868009794, -0.000291754616022926, -2.81009792107004e-13, + -73.9485382296513, -0.00423167900126897, -0.170505395422671, + -0.301191971384078, -0.0216747718871863, -0.00332619998061515, + -3.66549265245825, -0.211550389634303, -0.732731869409894, + -0.117239643328064}; + std::vector true_grad_w = { + 19.0909592842206, 24.8202270495257, 23.5200000150767, 20.578120951136, + 1.03356005252668, 2.72745006368943, 7.46860413650941, 0.749294067331125, + 3.99584607466979, 8.74288037387998, 1.83454127550223, 5.60378493291068, + 4.16252427817267}; + std::vector true_grad_v + = {1.04905306487445, 0.27197850836347, 0.479999997523509, + 0.977134030099708, 0.27971915432389, 0.762665629738307, + 3.35189640059343, 0.25195637817767, 0.558439280243674, + 1.59504189896887, 0.392934480883027, 1.62283843578951, + 1.22498956343152}; + std::vector true_grad_sv = {0, + 0, + 0, + 0, + -0.075059359595863, + 0, + 0, + -0.118203610656295, + 0.137718761146619, + 0, + -0.0215064672196581, + 1.05496167316742, + 0.220567350697758}; + std::vector true_grad_sw = {0, + 0, + 0, + 0, + 0, + -0.0276727265109265, + 0, + -0.0201028131785313, + 0, + 1.43818889437335, + -0.0415094057459971, + 0, + 0.25713151312856}; + std::vector true_grad_st0 = {0, + 0, + 0, + 0, + 0, + 0, + -0.151860277428223, + 0, + -0.00171635548928614, + -1.7784225450937, + -0.12422577484349, + -0.384925306866848, + -0.0625008391233113}; + + 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 lcdf = stan::math::wiener_lcdf(y, a, t0, w, v, sv, sw, st0); + lcdf.grad(); + EXPECT_NEAR(lcdf.val(), true_lcdf[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); + } +}