31 #ifndef STAN_MATH_PRIM_MAT_PROB_WIENER_LPDF_HPP 32 #define STAN_MATH_PRIM_MAT_PROB_WIENER_LPDF_HPP 68 template <
bool propto,
typename T_y,
typename T_alpha,
typename T_tau,
69 typename T_beta,
typename T_delta>
71 const T_y& y,
const T_alpha& alpha,
const T_tau& tau,
const T_beta&
beta,
72 const T_delta& delta) {
73 static const char*
function =
"wiener_lpdf";
79 static const double WIENER_ERR = 0.000001;
80 static const double PI_TIMES_WIENER_ERR =
pi() * WIENER_ERR;
81 static const double LOG_PI_LOG_WIENER_ERR =
LOG_PI +
log(WIENER_ERR);
82 static const double TWO_TIMES_SQRT_2_TIMES_SQRT_PI_TIMES_WIENER_ERR
84 static const double LOG_TWO_OVER_TWO_PLUS_LOG_SQRT_PI
86 static const double SQUARE_PI_OVER_TWO =
square(
pi()) * 0.5;
87 static const double TWO_TIMES_LOG_SQRT_PI = 2.0 *
LOG_SQRT_PI;
89 if (
size_zero(y, alpha, beta, tau, delta))
94 T_return_type lp(0.0);
110 alpha,
"A-priori bias", beta,
"Nondecision time", tau,
111 "Drift rate", delta);
124 for (
size_t i = 0; i < N_y_tau; ++i) {
125 if (y_vec[i] <= tau_vec[i]) {
126 std::stringstream msg;
127 msg <<
", but must be greater than nondecision time = " << tau_vec[i];
128 std::string msg_str(msg.str());
129 domain_error(
function,
"Random variable", y_vec[i],
" = ",
137 for (
size_t i = 0; i < N; i++) {
140 T_return_type x = (y_vec[i] - tau_vec[i]) / alpha2;
141 T_return_type kl, ks, tmp = 0;
143 T_return_type sqrt_x =
sqrt(x);
144 T_return_type log_x =
log(x);
145 T_return_type one_over_pi_times_sqrt_x = 1.0 /
pi() * sqrt_x;
149 if (PI_TIMES_WIENER_ERR * x < 1) {
151 kl =
sqrt(-2.0 *
SQRT_PI * (LOG_PI_LOG_WIENER_ERR + log_x)) / sqrt_x;
153 kl = (kl > one_over_pi_times_sqrt_x) ? kl : one_over_pi_times_sqrt_x;
155 kl = one_over_pi_times_sqrt_x;
159 T_return_type tmp_expr0
160 = TWO_TIMES_SQRT_2_TIMES_SQRT_PI_TIMES_WIENER_ERR * sqrt_x;
163 ks = 2.0 + sqrt_x *
sqrt(-2 *
log(tmp_expr0));
165 T_return_type sqrt_x_plus_one = sqrt_x + 1.0;
166 ks = (ks > sqrt_x_plus_one) ? ks : sqrt_x_plus_one;
172 T_return_type tmp_expr1 = (K - 1.0) / 2.0;
173 T_return_type tmp_expr2 =
ceil(tmp_expr1);
174 for (k = -
floor(tmp_expr1); k <= tmp_expr2; k++)
175 tmp += (one_minus_beta + 2.0 * k)
176 *
exp(-(
square(one_minus_beta + 2.0 * k)) * 0.5 / x);
177 tmp =
log(tmp) - LOG_TWO_OVER_TWO_PLUS_LOG_SQRT_PI - 1.5 * log_x;
180 for (k = 1; k <= K; ++k)
181 tmp += k *
exp(-(
square(k)) * (SQUARE_PI_OVER_TWO * x))
182 *
sin(k *
pi() * one_minus_beta);
183 tmp =
log(tmp) + TWO_TIMES_LOG_SQRT_PI;
187 lp += delta_vec[i] * alpha_vec[i] * one_minus_beta
188 -
square(delta_vec[i]) * x * alpha2 / 2.0 -
log(alpha2) + tmp;
193 template <
typename T_y,
typename T_alpha,
typename T_tau,
typename T_beta,
197 const T_beta&
beta,
const T_delta& delta) {
198 return wiener_lpdf<false>(y, alpha, tau,
beta, delta);
const double LOG_2
The natural logarithm of 2, .
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
fvar< T > sqrt(const fvar< T > &x)
void check_bounded(const char *function, const char *name, const T_y &y, const T_low &low, const T_high &high)
Check if the value is between the low and high values, inclusively.
fvar< T > log(const fvar< T > &x)
scalar_seq_view provides a uniform sequence-like wrapper around either a scalar or a sequence of scal...
bool size_zero(T &x)
Returns 1 if input is of length 0, returns 0 otherwise.
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...
fvar< T > square(const fvar< T > &x)
const double SQRT_2_TIMES_SQRT_PI
fvar< T > beta(const fvar< T > &x1, const fvar< T > &x2)
Return fvar with the beta function applied to the specified arguments and its gradient.
fvar< T > sin(const fvar< T > &x)
boost::math::tools::promote_args< double, typename scalar_type< T >::type, typename return_type< Types_pack... >::type >::type type
fvar< T > exp(const fvar< T > &x)
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
return_type< T_y, T_alpha, T_tau, T_beta, T_delta >::type wiener_lpdf(const T_y &y, const T_alpha &alpha, const T_tau &tau, const T_beta &beta, const T_delta &delta)
The log of the first passage time density function for a (Wiener) drift diffusion model for the given...
size_t max_size(const T1 &x1, const T2 &x2)
void 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.
int max(const std::vector< int > &x)
Returns the maximum coefficient in the specified column vector.
fvar< T > floor(const fvar< T > &x)
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
double pi()
Return the value of pi.
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
void check_consistent_sizes(const char *function, const char *name1, const T1 &x1, const char *name2, const T2 &x2)
Check if the dimension of x1 is consistent with x2.
fvar< T > ceil(const fvar< T > &x)