Automatic Differentiation
 
Loading...
Searching...
No Matches
wiener5_lpdf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_WIENER5_LPDF_HPP
2#define STAN_MATH_PRIM_PROB_WIENER5_LPDF_HPP
3
4#include <stan/math/prim/fun.hpp>
5#include <stan/math/prim/functor/apply.hpp>
6
7namespace stan {
8namespace math {
9namespace internal {
10
11enum GradientCalc { OFF = 0, ON = 1 };
12
29template <typename T_y, typename T_a, typename T_v, typename T_w, typename T_sv>
30inline auto wiener5_compute_error_term(T_y&& y, T_a&& a, T_v&& v_value,
31 T_w&& w_value, T_sv&& sv) noexcept {
32 const auto w = 1.0 - w_value;
33 const auto v = -v_value;
34 const auto sv_sqr = square(sv);
35 const auto one_plus_svsqr_y = 1 + sv_sqr * y;
36 const auto two_avw = 2.0 * a * v * w;
37 const auto two_log_a = 2.0 * log(a);
38 return stan::math::eval((sv_sqr * square(a * w) - two_avw - square(v) * y)
39 / 2.0 / one_plus_svsqr_y
40 - two_log_a - 0.5 * log(one_plus_svsqr_y));
41}
42
59template <bool Density, GradientCalc GradW, typename T_y, typename T_a,
60 typename T_w_value, typename T_err>
61inline auto wiener5_n_terms_small_t(T_y&& y, T_a&& a, T_w_value&& w_value,
62 T_err&& error) noexcept {
63 const auto two_error = 2.0 * error;
64 const auto y_asq = y / square(a);
65 const auto two_log_a = 2 * log(a);
66 const auto log_y_asq = log(y) - two_log_a;
67 const auto w = 1.0 - w_value;
68
69 const auto n_1_factor = Density ? 2 : 3;
70 const auto n_1 = (sqrt(n_1_factor * y_asq) + w) / 2.0;
71 auto u_eps = (Density || GradW)
72 ? fmin(-1.0, LOG_TWO + LOG_PI + 2.0 * log_y_asq + two_error)
73 : fmin(-3.0, (log(8.0) - log(27.0) + LOG_PI + 4.0 * log_y_asq
74 + two_error));
75 const auto arg_mult = (Density || GradW) ? 1 : 3;
76 const auto arg = -arg_mult * y_asq * (u_eps - sqrt(-2.0 * u_eps - 2.0));
77
78 const auto n_2
79 = (arg > 0) ? GradW ? 0.5 * (sqrt(arg) + w) : 0.5 * (sqrt(arg) - w) : n_1;
80
81 return ceil(fmax(n_1, n_2));
82}
83
98template <typename T_y, typename T_a, typename T_w_value, typename T_err>
99inline auto wiener5_density_large_reaction_time_terms(T_y&& y, T_a&& a,
100 T_w_value&& w_value,
101 T_err&& error) noexcept {
102 const auto y_asq = y / square(a);
103 const auto log_y_asq = log(y) - 2 * log(a);
104 static constexpr double PI_SQUARED = pi() * pi();
105 auto n_1 = 1.0 / (pi() * sqrt(y_asq));
106 const auto two_log_piy = -2.0 * (LOG_PI + log_y_asq + error);
107 auto n_2
108 = (two_log_piy >= 0) ? sqrt(two_log_piy / (PI_SQUARED * y_asq)) : 0.0;
109 return ceil(fmax(n_1, n_2));
110}
111
127template <GradientCalc GradW, typename T_y, typename T_a, typename T_w_value,
128 typename T_err>
129inline auto wiener5_gradient_large_reaction_time_terms(T_y&& y, T_a&& a,
130 T_w_value&& w_value,
131 T_err&& error) noexcept {
132 const auto y_asq = y / square(a);
133 const auto log_y_asq = log(y) - 2 * log(a);
134 static constexpr double PI_SQUARED = pi() * pi();
135 const auto n_1_factor = GradW ? 2 : 3;
136 auto n_1 = sqrt(n_1_factor / y_asq) / pi();
137 const auto two_error = 2.0 * error;
138 const auto u_eps_arg
139 = GradW ? log(4.0) - log(9.0) + 2.0 * LOG_PI + 3.0 * log_y_asq + two_error
140 : log(3.0) - log(5.0) + LOG_PI + 2.0 * log_y_asq + error;
141 const auto u_eps = fmin(-1, u_eps_arg);
142 const auto arg_mult = GradW ? 1 : (2.0 / PI_SQUARED / y_asq);
143 const auto arg = -arg_mult * (u_eps - sqrt(-2.0 * u_eps - 2.0));
144 auto n_2 = GradW ? ((arg > 0) ? sqrt(arg / y_asq) / pi() : n_1)
145 : ((arg > 0) ? sqrt(arg) : n_1);
146 return ceil(fmax(n_1, n_2));
147}
148
167template <bool Density, GradientCalc GradW, typename T_y, typename T_a,
168 typename T_w, typename T_nsmall, typename T_nlarge>
169inline auto wiener5_log_sum_exp(T_y&& y, T_a&& a, T_w&& w_value,
170 T_nsmall&& n_terms_small_t,
171 T_nlarge&& n_terms_large_t) noexcept {
172 const auto y_asq = y / square(a);
173 const auto w = 1.0 - w_value;
174 const bool small_n_terms_small_t
175 = Density ? (2 * n_terms_small_t <= n_terms_large_t)
176 : (2 * n_terms_small_t < n_terms_large_t);
177 const auto scaling = small_n_terms_small_t ? inv(2.0 * y_asq) : y_asq / 2.0;
178 using ret_t = return_type_t<T_y, T_a, T_w, T_nsmall, T_nlarge>;
179 ret_t fplus = NEGATIVE_INFTY;
180 ret_t fminus = NEGATIVE_INFTY;
181 int current_sign;
182 if (small_n_terms_small_t) {
183 constexpr double mult = Density ? 1.0 : 3.0;
184 if (GradW) {
185 for (auto k = n_terms_small_t; k >= 1; k--) {
186 const auto w_plus_2_k = w + 2.0 * k;
187 const auto w_minus_2_k = w - 2.0 * k;
188 const auto square_w_plus_2_k_minus_offset = square(w_plus_2_k) - y_asq;
189 if (square_w_plus_2_k_minus_offset > 0) {
190 const auto summand_plus = log(square_w_plus_2_k_minus_offset)
191 - square(w_plus_2_k) * scaling;
192 fplus = log_sum_exp(fplus, summand_plus);
193 } else if (square_w_plus_2_k_minus_offset < 0) {
194 const auto summand_plus = log(-square_w_plus_2_k_minus_offset)
195 - square(w_plus_2_k) * scaling;
196 fminus = log_sum_exp(fminus, summand_plus);
197 }
198 const auto square_w_minus_2_k_minus_offset
199 = square(w_minus_2_k) - y_asq;
200 if (square_w_minus_2_k_minus_offset > 0) {
201 const auto summand_minus = log(square_w_minus_2_k_minus_offset)
202 - square(w_minus_2_k) * scaling;
203 fplus = log_sum_exp(fplus, summand_minus);
204 } else if (square_w_minus_2_k_minus_offset < 0) {
205 const auto summand_minus = log(-square_w_minus_2_k_minus_offset)
206 - square(w_minus_2_k) * scaling;
207 fminus = log_sum_exp(fminus, summand_minus);
208 }
209 }
210 const auto square_w_minus_offset = square(w) - y_asq;
211 if (square_w_minus_offset > 0) {
212 const auto new_val = log(square_w_minus_offset) - square(w) * scaling;
213 fplus = log_sum_exp(fplus, new_val);
214 } else if (square_w_minus_offset < 0) {
215 const auto new_val = log(-square_w_minus_offset) - square(w) * scaling;
216 fminus = log_sum_exp(fminus, new_val);
217 }
218 } else {
219 for (auto k = n_terms_small_t; k >= 1; k--) {
220 const auto w_plus_2_k = w + 2.0 * k;
221 const auto w_minus_2_k = w - 2.0 * k;
222 const auto summand_plus
223 = mult * log(w_plus_2_k) - square(w_plus_2_k) * scaling;
224 fplus = log_sum_exp(fplus, summand_plus);
225 const auto summand_minus
226 = mult * log(-w_minus_2_k) - square(w_minus_2_k) * scaling;
227 if (fminus <= NEGATIVE_INFTY) {
228 fminus = summand_minus;
229 } else if (summand_minus <= NEGATIVE_INFTY) {
230 continue;
231 } else if (fminus > summand_minus) {
232 fminus = fminus + log1p_exp(summand_minus - fminus);
233 } else {
234 fminus = summand_minus + log1p_exp(fminus - summand_minus);
235 }
236 }
237 const auto new_val = mult * log(w) - square(w) * scaling;
238 fplus = log_sum_exp(fplus, new_val);
239 }
240 } else { // for large t
241 constexpr double mult = (Density ? 1 : (GradW ? 2 : 3));
242 for (auto k = n_terms_large_t; k >= 1; k--) {
243 const auto pi_k = k * pi();
244 const auto check = (GradW) ? cos(pi_k * w) : sin(pi_k * w);
245 if (check > 0) {
246 fplus = log_sum_exp(
247 fplus, mult * log(k) - square(pi_k) * scaling + log(check));
248 } else if ((GradW && check < 0) || !GradW) {
249 fminus = log_sum_exp(
250 fminus, mult * log(k) - square(pi_k) * scaling + log(-check));
251 }
252 }
253 }
254 current_sign = (fplus < fminus) ? -1 : 1;
255 if (fplus == NEGATIVE_INFTY) {
256 return std::make_pair(fminus, current_sign);
257 } else if (fminus == NEGATIVE_INFTY) {
258 return std::make_pair(fplus, current_sign);
259 } else if (fplus > fminus) {
260 return std::make_pair(log_diff_exp(fplus, fminus), current_sign);
261 } else if (fplus < fminus) {
262 return std::make_pair(log_diff_exp(fminus, fplus), current_sign);
263 } else {
264 return std::make_pair(ret_t(NEGATIVE_INFTY), current_sign);
265 }
266}
267
288template <bool NaturalScale = false, typename T_y, typename T_a, typename T_w,
289 typename T_v, typename T_sv, typename T_err>
290inline auto wiener5_density(const T_y& y, const T_a& a, const T_v& v_value,
291 const T_w& w_value, const T_sv& sv,
292 T_err&& err = log(1e-12)) noexcept {
293 const auto error_term
294 = wiener5_compute_error_term(y, a, v_value, w_value, sv);
295 const auto error = (err - error_term);
296 const auto n_terms_small_t
297 = wiener5_n_terms_small_t<GradientCalc::ON, GradientCalc::OFF>(
298 y, a, w_value, error);
299 const auto n_terms_large_t
300 = wiener5_density_large_reaction_time_terms(y, a, w_value, error);
301
302 auto res = wiener5_log_sum_exp<GradientCalc::ON, GradientCalc::OFF>(
303 y, a, w_value, n_terms_small_t, n_terms_large_t)
304 .first;
305 if (2 * n_terms_small_t <= n_terms_large_t) {
306 auto log_density = error_term - 0.5 * LOG_TWO - LOG_SQRT_PI
307 - 1.5 * (log(y) - 2 * log(a)) + res;
308 return NaturalScale ? exp(log_density) : log_density;
309 } else {
310 auto log_density = error_term + res + LOG_PI;
311 return NaturalScale ? exp(log_density) : log_density;
312 }
313}
314
335template <bool WrtLog = false, typename T_y, typename T_a, typename T_w,
336 typename T_v, typename T_sv, typename T_err>
337inline auto wiener5_grad_t(const T_y& y, const T_a& a, const T_v& v_value,
338 const T_w& w_value, const T_sv& sv,
339 T_err&& err = log(1e-12)) noexcept {
340 const auto two_log_a = 2 * log(a);
341 const auto log_y_asq = log(y) - two_log_a;
342 const auto error_term
343 = wiener5_compute_error_term(y, a, v_value, w_value, sv);
344 const auto w = 1.0 - w_value;
345 const auto v = -v_value;
346 const auto sv_sqr = square(sv);
347 const auto one_plus_svsqr_y = 1 + sv_sqr * y;
348 const auto density_part_one
349 = -0.5
350 * (square(sv_sqr) * (y + square(a * w))
351 + sv_sqr * (1 - (2.0 * a * v * w)) + square(v))
352 / square(one_plus_svsqr_y);
353 const auto error = (err - error_term) + two_log_a;
354 const auto n_terms_small_t
355 = wiener5_n_terms_small_t<GradientCalc::OFF, GradientCalc::OFF>(
356 y, a, w_value, error);
357 const auto n_terms_large_t
358 = wiener5_gradient_large_reaction_time_terms<GradientCalc::OFF>(
359 y, a, w_value, error);
360 auto wiener_res = wiener5_log_sum_exp<GradientCalc::OFF, GradientCalc::OFF>(
361 y, a, w_value, n_terms_small_t, n_terms_large_t);
362 auto&& result = wiener_res.first;
363 auto&& newsign = wiener_res.second;
364 const auto error_log_density
365 = log(fmax(fabs(density_part_one - 1.5 / y), fabs(density_part_one)));
366 const auto log_density = wiener5_density<GradientCalc::OFF>(
367 y, a, v_value, w_value, sv, err - error_log_density);
368 if (2 * n_terms_small_t < n_terms_large_t) {
369 auto ans = density_part_one - 1.5 / y
370 + newsign
371 * exp(error_term - two_log_a - 1.5 * LOG_TWO - LOG_SQRT_PI
372 - 3.5 * log_y_asq + result - log_density);
373 return WrtLog ? ans * exp(log_density) : ans;
374 } else {
375 auto ans = density_part_one
376 - newsign
377 * exp(error_term - two_log_a + 3.0 * LOG_PI - LOG_TWO
378 + result - log_density);
379 return WrtLog ? ans * exp(log_density) : ans;
380 }
381}
382
403template <bool WrtLog = false, typename T_y, typename T_a, typename T_w,
404 typename T_v, typename T_sv, typename T_err>
405inline auto wiener5_grad_a(const T_y& y, const T_a& a, const T_v& v_value,
406 const T_w& w_value, const T_sv& sv,
407 T_err&& err = log(1e-12)) noexcept {
408 const auto two_log_a = 2 * log(a);
409 const auto error_term
410 = wiener5_compute_error_term(y, a, v_value, w_value, sv);
411 const auto w = 1.0 - w_value;
412 const auto v = -v_value;
413 const auto sv_sqr = square(sv);
414 const auto one_plus_svsqr_y = 1 + sv_sqr * y;
415 const auto density_part_one
416 = (-v * w + sv_sqr * square(w) * a) / one_plus_svsqr_y;
417 const auto error = err - error_term + 3 * log(a) - log(y) - LOG_TWO;
418
419 const auto n_terms_small_t
420 = wiener5_n_terms_small_t<GradientCalc::OFF, GradientCalc::OFF>(
421 y, a, w_value, error);
422 const auto n_terms_large_t
423 = wiener5_gradient_large_reaction_time_terms<GradientCalc::OFF>(
424 y, a, w_value, error);
425 auto wiener_res = wiener5_log_sum_exp<GradientCalc::OFF, GradientCalc::OFF>(
426 y, a, w_value, n_terms_small_t, n_terms_large_t);
427 auto&& result = wiener_res.first;
428 auto&& newsign = wiener_res.second;
429 const auto error_log_density = log(
430 fmax(fabs(density_part_one + 1.0 / a), fabs(density_part_one - 2.0 / a)));
431 const auto log_density = wiener5_density<GradientCalc::OFF>(
432 y, a, v_value, w_value, sv, err - error_log_density);
433 if (2 * n_terms_small_t < n_terms_large_t) {
434 auto ans
435 = density_part_one + 1.0 / a
436 - newsign
437 * exp(-0.5 * LOG_TWO - LOG_SQRT_PI - 2.5 * log(y)
438 + 2.0 * two_log_a + error_term + result - log_density);
439 return WrtLog ? ans * exp(log_density) : ans;
440 } else {
441 auto ans = density_part_one - 2.0 / a
442 + newsign
443 * exp(log(y) + error_term - 3 * (log(a) - LOG_PI) + result
444 - log_density);
445 return WrtLog ? ans * exp(log_density) : ans;
446 }
447}
448
469template <bool WrtLog = false, typename T_y, typename T_a, typename T_w,
470 typename T_v, typename T_sv, typename T_err>
471inline auto wiener5_grad_v(const T_y& y, const T_a& a, const T_v& v_value,
472 const T_w& w_value, const T_sv& sv,
473 T_err&& err = log(1e-12)) noexcept {
474 auto ans = (a * (1 - w_value) - v_value * y) / (1.0 + square(sv) * y);
475 if (WrtLog) {
476 return ans * wiener5_density<true>(y, a, v_value, w_value, sv, err);
477 } else {
478 return ans;
479 }
480}
481
502template <bool WrtLog = false, typename T_y, typename T_a, typename T_w,
503 typename T_v, typename T_sv, typename T_err>
504inline auto wiener5_grad_w(const T_y& y, const T_a& a, const T_v& v_value,
505 const T_w& w_value, const T_sv& sv,
506 T_err&& err = log(1e-12)) noexcept {
507 const auto two_log_a = 2 * log(a);
508 const auto log_y_asq = log(y) - two_log_a;
509 const auto error_term
510 = wiener5_compute_error_term(y, a, v_value, w_value, sv);
511 const auto w = 1.0 - w_value;
512 const auto v = -v_value;
513 const auto sv_sqr = square(sv);
514 const auto one_plus_svsqr_y = 1 + sv_sqr * y;
515 const auto density_part_one
516 = (-v * a + sv_sqr * square(a) * w) / one_plus_svsqr_y;
517 const auto error = (err - error_term);
518
519 const auto n_terms_small_t
520 = wiener5_n_terms_small_t<GradientCalc::OFF, GradientCalc::ON>(
521 y, a, w_value, error);
522 const auto n_terms_large_t
523 = wiener5_gradient_large_reaction_time_terms<GradientCalc::ON>(
524 y, a, w_value, error);
525 auto wiener_res = wiener5_log_sum_exp<GradientCalc::OFF, GradientCalc::ON>(
526 y, a, w_value, n_terms_small_t, n_terms_large_t);
527 auto&& result = wiener_res.first;
528 auto&& newsign = wiener_res.second;
529 const auto log_density = wiener5_density<GradientCalc::OFF>(
530 y, a, v_value, w_value, sv, err - log(fabs(density_part_one)));
531 if (2 * n_terms_small_t < n_terms_large_t) {
532 auto ans = -(density_part_one
533 - newsign
534 * exp(result - (log_density - error_term)
535 - 2.5 * log_y_asq - 0.5 * LOG_TWO - 0.5 * LOG_PI));
536 return WrtLog ? ans * exp(log_density) : ans;
537 } else {
538 auto ans
539 = -(density_part_one
540 + newsign * exp(result - (log_density - error_term) + 2 * LOG_PI));
541 return WrtLog ? ans * exp(log_density) : ans;
542 }
543}
544
565template <bool WrtLog = false, typename T_y, typename T_a, typename T_w,
566 typename T_v, typename T_sv, typename T_err>
567inline auto wiener5_grad_sv(const T_y& y, const T_a& a, const T_v& v_value,
568 const T_w& w_value, const T_sv& sv,
569 T_err&& err = log(1e-12)) noexcept {
570 const auto one_plus_svsqr_y = 1 + square(sv) * y;
571 const auto w = 1.0 - w_value;
572 const auto v = -v_value;
573 const auto t1 = -y / one_plus_svsqr_y;
574 const auto t2 = (square(a * w) + 2 * a * v * w * y + square(v * y))
575 / square(one_plus_svsqr_y);
576 const auto ans = sv * (t1 + t2);
577 return WrtLog ? ans * wiener5_density<true>(y, a, v_value, w_value, sv, err)
578 : ans;
579}
580
591template <size_t NestedIndex, typename Scalar1, typename Scalar2>
592inline void assign_err(Scalar1 arg, Scalar2 err) {
593 arg = err;
594}
595
607template <size_t NestedIndex, typename Scalar, typename... TArgs>
608inline void assign_err(std::tuple<TArgs...>& args_tuple, Scalar err) {
609 std::get<NestedIndex>(args_tuple) = err;
610}
611
629template <size_t ErrIndex, size_t NestedIndex = 0,
630 GradientCalc GradW7 = GradientCalc::OFF, bool LogResult = true,
631 typename F, typename T_err, typename... ArgsTupleT>
632inline auto estimate_with_err_check(F&& functor, T_err&& err,
633 ArgsTupleT&&... args_tuple) {
634 auto result = functor(args_tuple...);
635 auto log_fabs_result = LogResult ? log(fabs(result)) : fabs(result);
636 if (log_fabs_result < err) {
637 log_fabs_result = is_inf(log_fabs_result) ? 0 : log_fabs_result;
638 auto err_args_tuple = std::make_tuple(args_tuple...);
639 const auto new_error
640 = GradW7 ? err + log_fabs_result + LOG_TWO : err + log_fabs_result;
641 assign_err<NestedIndex>(std::get<ErrIndex>(err_args_tuple), new_error);
642 result
643 = math::apply([](auto&& func, auto&&... args) { return func(args...); },
644 err_args_tuple, functor);
645 }
646 return result;
647}
648} // namespace internal
649
672template <bool propto = false, typename T_y, typename T_a, typename T_t0,
673 typename T_w, typename T_v, typename T_sv>
674inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0,
675 const T_w& w, const T_v& v, const T_sv& sv,
676 const double& precision_derivatives = 1e-4) {
677 using T_partials_return = partials_return_t<T_y, T_a, T_t0, T_w, T_v, T_sv>;
678 using ret_t = return_type_t<T_y, T_a, T_t0, T_w, T_v, T_sv>;
679 if (!include_summand<propto, T_y, T_a, T_t0, T_w, T_v, T_sv>::value) {
680 return ret_t(0.0);
681 }
682 using T_y_ref = ref_type_t<T_y>;
683 using T_a_ref = ref_type_t<T_a>;
684 using T_t0_ref = ref_type_t<T_t0>;
685 using T_w_ref = ref_type_t<T_w>;
686 using T_v_ref = ref_type_t<T_v>;
687 using T_sv_ref = ref_type_t<T_sv>;
688
689 static constexpr const char* function_name = "wiener5_lpdf";
690
691 check_consistent_sizes(function_name, "Random variable", y,
692 "Boundary separation", a, "Drift rate", v,
693 "A-priori bias", w, "Nondecision time", t0,
694 "Inter-trial variability in drift rate", sv);
695
696 T_y_ref y_ref = y;
697 T_a_ref a_ref = a;
698 T_t0_ref t0_ref = t0;
699 T_w_ref w_ref = w;
700 T_v_ref v_ref = v;
701 T_sv_ref sv_ref = sv;
702
703 decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
704 decltype(auto) a_val = to_ref(as_value_column_array_or_scalar(a_ref));
705 decltype(auto) v_val = to_ref(as_value_column_array_or_scalar(v_ref));
706 decltype(auto) w_val = to_ref(as_value_column_array_or_scalar(w_ref));
707 decltype(auto) t0_val = to_ref(as_value_column_array_or_scalar(t0_ref));
708 decltype(auto) sv_val = to_ref(as_value_column_array_or_scalar(sv_ref));
709 check_positive_finite(function_name, "Random variable", y_val);
710 check_positive_finite(function_name, "Boundary separation", a_val);
711 check_finite(function_name, "Drift rate", v_val);
712 check_less(function_name, "A-priori bias", w_val, 1);
713 check_greater(function_name, "A-priori bias", w_val, 0);
714 check_nonnegative(function_name, "Nondecision time", t0_val);
715 check_finite(function_name, "Nondecision time", t0_val);
716 check_nonnegative(function_name, "Inter-trial variability in drift rate",
717 sv_val);
718 check_finite(function_name, "Inter-trial variability in drift rate", sv_val);
719
720 if (size_zero(y, a, t0, w, v, sv)) {
721 return ret_t(0.0);
722 }
723 const size_t N = max_size(y, a, t0, w, v, sv);
724 if (!N) {
725 return ret_t(0.0);
726 }
727
728 scalar_seq_view<T_y_ref> y_vec(y_ref);
729 scalar_seq_view<T_a_ref> a_vec(a_ref);
730 scalar_seq_view<T_t0_ref> t0_vec(t0_ref);
731 scalar_seq_view<T_w_ref> w_vec(w_ref);
732 scalar_seq_view<T_v_ref> v_vec(v_ref);
733 scalar_seq_view<T_sv_ref> sv_vec(sv_ref);
734 const size_t N_y_t0 = max_size(y, t0);
735
736 for (size_t i = 0; i < N_y_t0; ++i) {
737 if (y_vec[i] <= t0_vec[i]) {
738 std::stringstream msg;
739 msg << ", but must be greater than nondecision time = " << t0_vec[i];
740 std::string msg_str(msg.str());
741 throw_domain_error(function_name, "Random variable", y_vec[i], " = ",
742 msg_str.c_str());
743 }
744 }
745
746 const auto log_error_density = log(1e-6);
747 const auto log_error_derivative = log(precision_derivatives);
748 const double log_error_absolute_val = log(1e-12);
749 const T_partials_return log_error_absolute = log_error_absolute_val;
750 T_partials_return log_density = 0.0;
751 auto ops_partials
752 = make_partials_propagator(y_ref, a_ref, t0_ref, w_ref, v_ref, sv_ref);
753
754 static constexpr double LOG_FOUR = LOG_TWO + LOG_TWO;
755
756 // calculate density and partials
757 for (size_t i = 0; i < N; i++) {
758 // Calculate 4-parameter model without inter-trial variabilities (if
759 // sv_vec[i] == 0) or 5-parameter model with inter-trial variability in
760 // drift rate (if sv_vec[i] != 0)
761
762 const auto y_value = y_vec.val(i);
763 const auto a_value = a_vec.val(i);
764 const auto t0_value = t0_vec.val(i);
765 const auto w_value = w_vec.val(i);
766 const auto v_value = v_vec.val(i);
767 const auto sv_value = sv_vec.val(i);
768 using internal::GradientCalc;
769 auto l_density = internal::estimate_with_err_check<5, 0, GradientCalc::OFF,
770 GradientCalc::OFF>(
771 [](auto&&... args) {
772 return internal::wiener5_density<GradientCalc::OFF>(args...);
773 },
774 log_error_density - LOG_TWO, y_value - t0_value, a_value, v_value,
775 w_value, sv_value, log_error_absolute);
776
777 log_density += l_density;
778
779 const auto new_est_err = l_density + log_error_derivative - LOG_FOUR;
780
781 // computation of derivative for t and precision check in order to give
782 // the value as deriv_y to edge1 and as -deriv_y to edge5
783 const auto deriv_y
784 = internal::estimate_with_err_check<5, 0, GradientCalc::OFF,
785 GradientCalc::ON>(
786 [](auto&&... args) {
787 return internal::wiener5_grad_t<GradientCalc::OFF>(args...);
788 },
789 new_est_err, y_value - t0_value, a_value, v_value, w_value,
790 sv_value, log_error_absolute);
791
792 // computation of derivatives and precision checks
793 if (!is_constant_all<T_y>::value) {
794 partials<0>(ops_partials)[i] = deriv_y;
795 }
796 if (!is_constant_all<T_a>::value) {
797 partials<1>(ops_partials)[i]
798 = internal::estimate_with_err_check<5, 0, GradientCalc::OFF,
799 GradientCalc::ON>(
800 [](auto&&... args) {
801 return internal::wiener5_grad_a<GradientCalc::OFF>(args...);
802 },
803 new_est_err, y_value - t0_value, a_value, v_value, w_value,
804 sv_value, log_error_absolute);
805 }
806 if (!is_constant_all<T_t0>::value) {
807 partials<2>(ops_partials)[i] = -deriv_y;
808 }
809 if (!is_constant_all<T_w>::value) {
810 partials<3>(ops_partials)[i]
811 = internal::estimate_with_err_check<5, 0, GradientCalc::OFF,
812 GradientCalc::ON>(
813 [](auto&&... args) {
814 return internal::wiener5_grad_w<GradientCalc::OFF>(args...);
815 },
816 new_est_err, y_value - t0_value, a_value, v_value, w_value,
817 sv_value, log_error_absolute);
818 }
819 if (!is_constant_all<T_v>::value) {
820 partials<4>(ops_partials)[i]
821 = internal::wiener5_grad_v<GradientCalc::OFF>(
822 y_value - t0_value, a_value, v_value, w_value, sv_value,
823 log_error_absolute_val);
824 }
825 if (!is_constant_all<T_sv>::value) {
826 partials<5>(ops_partials)[i]
827 = internal::wiener5_grad_sv<GradientCalc::OFF>(
828 y_value - t0_value, a_value, v_value, w_value, sv_value,
829 log_error_absolute_val);
830 }
831 } // end for loop
832 return ops_partials.build(log_density);
833} // end wiener_lpdf
834
835// ToDo: delete old wiener_lpdf implementation to use this one
836// template <bool propto = false, typename T_y, typename T_a, typename T_t0,
837// typename T_w, typename T_v>
838// inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0,
839// const T_w& w, const T_v& v,
840// const double& precision_derivatives = 1e-4) {
841// return wiener_lpdf(y, a, t0, w, v, 0, precision_derivatives);
842//} // end wiener_lpdf
843
844} // namespace math
845} // namespace stan
846#endif
stan::scalar_seq_view
scalar_seq_view provides a uniform sequence-like wrapper around either a scalar or a sequence of scal...
Definition scalar_seq_view.hpp:18
stan::return_type_t
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
Definition return_type.hpp:218
stan::math::internal::wiener5_grad_v
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
Calculate the derivative of the wiener5 density w.r.t.
Definition wiener5_lpdf.hpp:471
stan::math::internal::wiener5_density_large_reaction_time_terms
auto wiener5_density_large_reaction_time_terms(T_y &&y, T_a &&a, T_w_value &&w_value, T_err &&error) noexcept
Calculate the 'n_terms_large_t' term for a wiener5 density.
Definition wiener5_lpdf.hpp:99
stan::math::internal::wiener5_grad_a
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
Calculate the derivative of the wiener5 density w.r.t.
Definition wiener5_lpdf.hpp:405
stan::math::internal::wiener5_grad_sv
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
Calculate the derivative of the wiener5 density w.r.t.
Definition wiener5_lpdf.hpp:567
stan::math::internal::assign_err
void assign_err(Scalar1 arg, Scalar2 err)
Utility function for replacing a value with a specified error value.
Definition wiener5_lpdf.hpp:592
stan::math::internal::wiener5_n_terms_small_t
auto wiener5_n_terms_small_t(T_y &&y, T_a &&a, T_w_value &&w_value, T_err &&error) noexcept
Calculate the 'n_terms_small_t' term for a wiener5 density or gradient.
Definition wiener5_lpdf.hpp:61
stan::math::internal::wiener5_gradient_large_reaction_time_terms
auto wiener5_gradient_large_reaction_time_terms(T_y &&y, T_a &&a, T_w_value &&w_value, T_err &&error) noexcept
Calculate the 'n_terms_large_t' term for a wiener5 gradient.
Definition wiener5_lpdf.hpp:129
stan::math::internal::wiener5_grad_t
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
Calculate the derivative of the wiener5 density w.r.t.
Definition wiener5_lpdf.hpp:337
stan::math::internal::wiener5_density
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
Calculate the wiener5 density.
Definition wiener5_lpdf.hpp:290
stan::math::internal::wiener5_grad_w
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
Calculate the derivative of the wiener5 density w.r.t.
Definition wiener5_lpdf.hpp:504
stan::math::internal::estimate_with_err_check
auto estimate_with_err_check(F &&functor, T_err &&err, ArgsTupleT &&... args_tuple)
Utility function for estimating a function with a given set of arguments, checking the result against...
Definition wiener5_lpdf.hpp:632
stan::math::internal::wiener5_compute_error_term
auto wiener5_compute_error_term(T_y &&y, T_a &&a, T_v &&v_value, T_w &&w_value, T_sv &&sv) noexcept
Calculate the 'error_term' term for a wiener5 density or gradient.
Definition wiener5_lpdf.hpp:30
stan::math::internal::wiener5_log_sum_exp
auto wiener5_log_sum_exp(T_y &&y, T_a &&a, T_w &&w_value, T_nsmall &&n_terms_small_t, T_nlarge &&n_terms_large_t) noexcept
Calculate the 'result' term and its sign for a wiener5 density or gradient.
Definition wiener5_lpdf.hpp:169
stan::math::sin
fvar< T > sin(const fvar< T > &x)
Definition sin.hpp:16
stan::math::check_nonnegative
void check_nonnegative(const char *function, const char *name, const T_y &y)
Check if y is non-negative.
Definition check_nonnegative.hpp:24
stan::math::size_zero
bool size_zero(const T &x)
Returns 1 if input is of length 0, returns 0 otherwise.
Definition size_zero.hpp:19
stan::math::fmin
fvar< T > fmin(const fvar< T > &x1, const fvar< T > &x2)
Definition fmin.hpp:14
stan::math::e
static constexpr double e()
Return the base of the natural logarithm.
Definition constants.hpp:20
stan::math::arg
fvar< T > arg(const std::complex< fvar< T > > &z)
Return the phase angle of the complex argument.
Definition arg.hpp:19
stan::math::eval
T eval(T &&arg)
Inputs which have a plain_type equal to the own time are forwarded unmodified (for Eigen expressions ...
Definition eval.hpp:20
stan::math::log
fvar< T > log(const fvar< T > &x)
Definition log.hpp:18
stan::math::NEGATIVE_INFTY
static constexpr double NEGATIVE_INFTY
Negative infinity.
Definition constants.hpp:51
stan::math::throw_domain_error
void throw_domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw a domain error with a consistently formatted message.
Definition throw_domain_error.hpp:26
stan::math::LOG_TWO
static constexpr double LOG_TWO
The natural logarithm of 2, .
Definition constants.hpp:80
stan::math::log1p_exp
fvar< T > log1p_exp(const fvar< T > &x)
Definition log1p_exp.hpp:14
stan::math::as_value_column_array_or_scalar
auto as_value_column_array_or_scalar(T &&a)
Extract the value from an object and for eigen vectors and std::vectors convert to an eigen column ar...
Definition as_value_column_array_or_scalar.hpp:22
stan::math::check_consistent_sizes
void check_consistent_sizes(const char *)
Trivial no input case, this function is a no-op.
Definition check_consistent_sizes.hpp:15
stan::math::sqrt
fvar< T > sqrt(const fvar< T > &x)
Definition sqrt.hpp:18
stan::math::LOG_SQRT_PI
static constexpr double LOG_SQRT_PI
The natural logarithm of the square root of , .
Definition constants.hpp:110
stan::math::to_ref
ref_type_t< T && > to_ref(T &&a)
This evaluates expensive Eigen expressions.
Definition to_ref.hpp:17
stan::math::LOG_PI
static constexpr double LOG_PI
The natural logarithm of , .
Definition constants.hpp:86
stan::math::check_finite
void check_finite(const char *function, const char *name, const T_y &y)
Return true if all values in y are finite.
Definition check_finite.hpp:28
stan::math::fmax
fvar< T > fmax(const fvar< T > &x1, const fvar< T > &x2)
Return the greater of the two specified arguments.
Definition fmax.hpp:23
stan::math::log_diff_exp
fvar< T > log_diff_exp(const fvar< T > &x1, const fvar< T > &x2)
Definition log_diff_exp.hpp:14
stan::math::cos
fvar< T > cos(const fvar< T > &x)
Definition cos.hpp:16
stan::math::pi
static constexpr double pi()
Return the value of pi.
Definition constants.hpp:36
stan::math::wiener_lpdf
auto wiener_lpdf(const T_y &y, const T_a &a, const T_t0 &t0, const T_w &w, const T_v &v, const T_sv &sv, const double &precision_derivatives=1e-4)
Log-density function for the 5-parameter Wiener density.
Definition wiener5_lpdf.hpp:674
stan::math::is_inf
int is_inf(const fvar< T > &x)
Returns 1 if the input's value is infinite and 0 otherwise.
Definition is_inf.hpp:21
stan::math::check_less
void check_less(const char *function, const char *name, const T_y &y, const T_high &high, Idxs... idxs)
Throw an exception if y is not strictly less than high.
Definition check_less.hpp:35
stan::math::ceil
fvar< T > ceil(const fvar< T > &x)
Definition ceil.hpp:13
stan::math::max_size
int64_t max_size(const T1 &x1, const Ts &... xs)
Calculate the size of the largest input.
Definition max_size.hpp:20
stan::math::check_greater
void check_greater(const char *function, const char *name, const T_y &y, const T_low &low, Idxs... idxs)
Throw an exception if y is not strictly greater than low.
Definition check_greater.hpp:36
stan::math::inv
fvar< T > inv(const fvar< T > &x)
Definition inv.hpp:13
stan::math::make_partials_propagator
auto make_partials_propagator(Ops &&... ops)
Construct an partials_propagator.
Definition partials_propagator.hpp:119
stan::math::apply
constexpr decltype(auto) apply(F &&f, Tuple &&t, PreArgs &&... pre_args)
Definition apply.hpp:52
stan::math::check_positive_finite
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
Definition check_positive_finite.hpp:22
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition fabs.hpp:16
stan::math::square
fvar< T > square(const fvar< T > &x)
Definition square.hpp:12
stan::math::log_sum_exp
fvar< T > log_sum_exp(const fvar< T > &x1, const fvar< T > &x2)
Definition log_sum_exp.hpp:18
stan::math::exp
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:15
stan::ref_type_t
typename ref_type_if< true, T >::type ref_type_t
Definition ref_type.hpp:55
stan::partials_return_t
typename partials_return_type< Args... >::type partials_return_t
Definition partials_return_type.hpp:44
stan
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
Definition unit_vector_constrain.hpp:15
stan::math::conjunction
Extends std::true_type when instantiated with zero or more template parameters, all of which extend t...
Definition conjunction.hpp:14
stan::math::include_summand
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...
Definition include_summand.hpp:39