29 #ifndef STAN_MATH_PRIM_MAT_PROB_WIENER_LOG_HPP
30 #define STAN_MATH_PRIM_MAT_PROB_WIENER_LOG_HPP
42 #include <boost/math/distributions.hpp>
68 template <
bool propto,
69 typename T_y,
typename T_alpha,
typename T_tau,
70 typename T_beta,
typename T_delta>
71 typename return_type<T_y, T_alpha, T_tau, T_beta, T_delta>::type
72 wiener_log(
const T_y& y,
const T_alpha& alpha,
const T_tau& tau,
73 const T_beta& beta,
const T_delta& delta) {
74 static const char*
function(
"stan::math::wiener_log(%1%)");
76 using boost::math::tools::promote_args;
83 static const double WIENER_ERR = 0.000001;
84 static const double PI_TIMES_WIENER_ERR =
pi() * WIENER_ERR;
85 static const double LOG_PI_LOG_WIENER_ERR =
88 TWO_TIMES_SQRT_2_TIMES_SQRT_PI_TIMES_WIENER_ERR =
90 static const double LOG_TWO_OVER_TWO_PLUS_LOG_SQRT_PI =
92 static const double SQUARE_PI_OVER_TWO =
square(
pi()) * 0.5;
93 static const double TWO_TIMES_LOG_SQRT_PI = 2.0 *
LOG_SQRT_PI;
103 T_beta, T_delta>::type T_return_type;
104 T_return_type lp(0.0);
120 "Random variable", y,
121 "Boundary separation", alpha,
122 "A-priori bias", beta,
123 "Nondecision time", tau,
124 "Drift rate", delta);
137 T_beta, T_delta>::value) {
141 for (
size_t i = 0; i < N; i++)
142 if (y_vec[i] < tau_vec[i]) {
147 for (
size_t i = 0; i < N; i++) {
152 T_return_type x = y_vec[i];
153 T_return_type kl, ks, tmp = 0;
159 T_return_type sqrt_x =
sqrt(x);
160 T_return_type log_x =
log(x);
161 T_return_type one_over_pi_times_sqrt_x = 1.0 /
pi() * sqrt_x;
165 if (PI_TIMES_WIENER_ERR * x < 1) {
168 (LOG_PI_LOG_WIENER_ERR + log_x)) /
171 kl = (kl > one_over_pi_times_sqrt_x) ?
172 kl : one_over_pi_times_sqrt_x;
174 kl = one_over_pi_times_sqrt_x;
178 T_return_type tmp_expr0
179 = TWO_TIMES_SQRT_2_TIMES_SQRT_PI_TIMES_WIENER_ERR * sqrt_x;
182 ks = 2.0 + sqrt_x *
sqrt(-2 *
log(tmp_expr0));
184 T_return_type sqrt_x_plus_one = sqrt_x + 1.0;
185 ks = (ks > sqrt_x_plus_one) ? ks : sqrt_x_plus_one;
192 T_return_type tmp_expr1 = (K - 1.0) / 2.0;
193 T_return_type tmp_expr2 =
ceil(tmp_expr1);
194 for (k = -
floor(tmp_expr1); k <= tmp_expr2; k++)
196 tmp += (one_minus_beta + 2.0 * k) *
197 exp(-(
square(one_minus_beta + 2.0 * k)) * 0.5 / x);
200 LOG_TWO_OVER_TWO_PLUS_LOG_SQRT_PI - 1.5 * log_x;
203 for (k = 1; k <= K; k++)
206 (SQUARE_PI_OVER_TWO * x)) *
207 sin(k *
pi() * one_minus_beta);
209 TWO_TIMES_LOG_SQRT_PI;
213 lp += delta_vec[i] * alpha_vec[i] * one_minus_beta -
214 square(delta_vec[i]) * x * alpha2 / 2.0 -
221 template <
typename T_y,
typename T_alpha,
typename T_tau,
222 typename T_beta,
typename T_delta>
225 wiener_log(
const T_y& y,
const T_alpha& alpha,
const T_tau& tau,
226 const T_beta& beta,
const T_delta& delta) {
227 return wiener_log<false>(y, alpha, tau, beta, delta);
return_type< T_y, T_alpha, T_tau, T_beta, T_delta >::type wiener_log(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...
bool isfinite(const stan::math::var &v)
Checks if the given number has finite value.
fvar< T > sqrt(const fvar< T > &x)
bool check_not_nan(const char *function, const char *name, const T_y &y)
Return true if y is not NaN.
fvar< T > log(const fvar< T > &x)
size_t length(const std::vector< T > &x)
bool check_bounded(const char *function, const char *name, const T_y &y, const T_low &low, const T_high &high)
Return true if the value is between the low and high values, inclusively.
Metaprogram to calculate the base scalar return type resulting from promoting all the scalar types of...
scalar_type_helper< is_vector< T >::value, T >::type type
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...
boost::math::tools::promote_args< typename scalar_type< T1 >::type, typename scalar_type< T2 >::type, typename scalar_type< T3 >::type, typename scalar_type< T4 >::type, typename scalar_type< T5 >::type, typename scalar_type< T6 >::type >::type type
fvar< T > square(const fvar< T > &x)
const double SQRT_2_TIMES_SQRT_PI
bool isinf(const stan::math::var &v)
Checks if the given number is infinite.
fvar< T > sin(const fvar< T > &x)
fvar< T > exp(const fvar< T > &x)
bool check_positive(const char *function, const char *name, const T_y &y)
Return true if y is positive.
size_t max_size(const T1 &x1, const T2 &x2)
int max(const std::vector< int > &x)
Returns the maximum coefficient in the specified column vector.
bool check_finite(const char *function, const char *name, const T_y &y)
Return true if y is finite.
fvar< T > floor(const fvar< T > &x)
bool check_consistent_sizes(const char *function, const char *name1, const T1 &x1, const char *name2, const T2 &x2)
Return true if the dimension of x1 is consistent with x2.
double pi()
Return the value of pi.
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
VectorView is a template expression that is constructed with a container or scalar, which it then allows to be used as an array using operator[].
fvar< T > ceil(const fvar< T > &x)
double negative_infinity()
Return negative infinity.