wiener_lpdf.hpp
Go to the documentation of this file.
1 // Original code from which Stan's code is derived:
2 // Copyright (c) 2013, Joachim Vandekerckhove.
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted
7 // provided that the following conditions are met:
8 //
9 // * Redistributions of source code must retain the above copyright notice,
10 // * this list of conditions and the following disclaimer.
11 // * Redistributions in binary form must reproduce the above copyright notice,
12 // * this list of conditions and the following disclaimer in the
13 // * documentation and/or other materials provided with the distribution.
14 // * Neither the name of the University of California, Irvine nor the names
15 // * of its contributors may be used to endorse or promote products derived
16 // * from this software without specific prior written permission.
17 //
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
22 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
28 // THE POSSIBILITY OF SUCH DAMAGE.
29 
30 #ifndef STAN_MATH_PRIM_MAT_PROB_WIENER_LPDF_HPP
31 #define STAN_MATH_PRIM_MAT_PROB_WIENER_LPDF_HPP
32 
45 #include <boost/math/distributions.hpp>
46 #include <algorithm>
47 #include <cmath>
48 #include <string>
49 
50 namespace stan {
51 namespace math {
52 
53 /**
54  * The log of the first passage time density function for a (Wiener)
55  * drift diffusion model for the given \f$y\f$,
56  * boundary separation \f$\alpha\f$, nondecision time \f$\tau\f$,
57  * relative bias \f$\beta\f$, and drift rate \f$\delta\f$.
58  * \f$\alpha\f$ and \f$\tau\f$ must be greater than 0, and
59  * \f$\beta\f$ must be between 0 and 1. \f$y\f$ should contain
60  * reaction times in seconds (strictly positive) with
61  * upper-boundary responses.
62  *
63  * @param y A scalar variate.
64  * @param alpha The boundary separation.
65  * @param tau The nondecision time.
66  * @param beta The relative bias.
67  * @param delta The drift rate.
68  * @return The log of the Wiener first passage time density of
69  * the specified arguments.
70  */
71 template <bool propto, typename T_y, typename T_alpha, typename T_tau,
72  typename T_beta, typename T_delta>
74  const T_y& y, const T_alpha& alpha, const T_tau& tau, const T_beta& beta,
75  const T_delta& delta) {
76  static const char* function = "wiener_lpdf";
77 
78  using std::exp;
79  using std::log;
80  using std::pow;
81 
82  static const double WIENER_ERR = 0.000001;
83  static const double PI_TIMES_WIENER_ERR = pi() * WIENER_ERR;
84  static const double LOG_PI_LOG_WIENER_ERR = LOG_PI + log(WIENER_ERR);
85  static const double TWO_TIMES_SQRT_2_TIMES_SQRT_PI_TIMES_WIENER_ERR
86  = 2.0 * SQRT_2_TIMES_SQRT_PI * WIENER_ERR;
87  static const double LOG_TWO_OVER_TWO_PLUS_LOG_SQRT_PI
88  = LOG_TWO / 2 + LOG_SQRT_PI;
89  static const double SQUARE_PI_OVER_TWO = square(pi()) * 0.5;
90  static const double TWO_TIMES_LOG_SQRT_PI = 2.0 * LOG_SQRT_PI;
91 
92  if (size_zero(y, alpha, beta, tau, delta))
93  return 0.0;
94 
96  T_return_type;
97  T_return_type lp(0.0);
98 
99  check_not_nan(function, "Random variable", y);
100  check_not_nan(function, "Boundary separation", alpha);
101  check_not_nan(function, "A-priori bias", beta);
102  check_not_nan(function, "Nondecision time", tau);
103  check_not_nan(function, "Drift rate", delta);
104  check_finite(function, "Boundary separation", alpha);
105  check_finite(function, "A-priori bias", beta);
106  check_finite(function, "Nondecision time", tau);
107  check_finite(function, "Drift rate", delta);
108  check_positive(function, "Random variable", y);
109  check_positive(function, "Boundary separation", alpha);
110  check_positive(function, "Nondecision time", tau);
111  check_bounded(function, "A-priori bias", beta, 0, 1);
112  check_consistent_sizes(function, "Random variable", y, "Boundary separation",
113  alpha, "A-priori bias", beta, "Nondecision time", tau,
114  "Drift rate", delta);
115 
116  size_t N = std::max(max_size(y, alpha, beta), max_size(tau, delta));
117  if (!N)
118  return 0.0;
119 
120  scalar_seq_view<T_y> y_vec(y);
121  scalar_seq_view<T_alpha> alpha_vec(alpha);
122  scalar_seq_view<T_beta> beta_vec(beta);
123  scalar_seq_view<T_tau> tau_vec(tau);
124  scalar_seq_view<T_delta> delta_vec(delta);
125 
126  size_t N_y_tau = max_size(y, tau);
127  for (size_t i = 0; i < N_y_tau; ++i) {
128  if (y_vec[i] <= tau_vec[i]) {
129  std::stringstream msg;
130  msg << ", but must be greater than nondecision time = " << tau_vec[i];
131  std::string msg_str(msg.str());
132  domain_error(function, "Random variable", y_vec[i], " = ",
133  msg_str.c_str());
134  }
135  }
136 
138  return 0;
139 
140  for (size_t i = 0; i < N; i++) {
141  typename scalar_type<T_beta>::type one_minus_beta = 1.0 - beta_vec[i];
142  typename scalar_type<T_alpha>::type alpha2 = square(alpha_vec[i]);
143  T_return_type x = (y_vec[i] - tau_vec[i]) / alpha2;
144  T_return_type kl, ks, tmp = 0;
145  T_return_type k, K;
146  T_return_type sqrt_x = sqrt(x);
147  T_return_type log_x = log(x);
148  T_return_type one_over_pi_times_sqrt_x = 1.0 / pi() * sqrt_x;
149 
150  // calculate number of terms needed for large t:
151  // if error threshold is set low enough
152  if (PI_TIMES_WIENER_ERR * x < 1) {
153  // compute bound
154  kl = sqrt(-2.0 * SQRT_PI * (LOG_PI_LOG_WIENER_ERR + log_x)) / sqrt_x;
155  // ensure boundary conditions met
156  kl = (kl > one_over_pi_times_sqrt_x) ? kl : one_over_pi_times_sqrt_x;
157  } else {
158  kl = one_over_pi_times_sqrt_x; // set to boundary condition
159  }
160  // calculate number of terms needed for small t:
161  // if error threshold is set low enough
162  T_return_type tmp_expr0
163  = TWO_TIMES_SQRT_2_TIMES_SQRT_PI_TIMES_WIENER_ERR * sqrt_x;
164  if (tmp_expr0 < 1) {
165  // compute bound
166  ks = 2.0 + sqrt_x * sqrt(-2 * log(tmp_expr0));
167  // ensure boundary conditions are met
168  T_return_type sqrt_x_plus_one = sqrt_x + 1.0;
169  ks = (ks > sqrt_x_plus_one) ? ks : sqrt_x_plus_one;
170  } else { // if error threshold was set too high
171  ks = 2.0; // minimal kappa for that case
172  }
173  if (ks < kl) { // small t
174  K = ceil(ks); // round to smallest integer meeting error
175  T_return_type tmp_expr1 = (K - 1.0) / 2.0;
176  T_return_type tmp_expr2 = ceil(tmp_expr1);
177  for (k = -floor(tmp_expr1); k <= tmp_expr2; k++)
178  tmp += (one_minus_beta + 2.0 * k)
179  * exp(-(square(one_minus_beta + 2.0 * k)) * 0.5 / x);
180  tmp = log(tmp) - LOG_TWO_OVER_TWO_PLUS_LOG_SQRT_PI - 1.5 * log_x;
181  } else { // if large t is better...
182  K = ceil(kl); // round to smallest integer meeting error
183  for (k = 1; k <= K; ++k)
184  tmp += k * exp(-(square(k)) * (SQUARE_PI_OVER_TWO * x))
185  * sin(k * pi() * one_minus_beta);
186  tmp = log(tmp) + TWO_TIMES_LOG_SQRT_PI;
187  }
188 
189  // convert to f(t|v,a,w) and return result
190  lp += delta_vec[i] * alpha_vec[i] * one_minus_beta
191  - square(delta_vec[i]) * x * alpha2 / 2.0 - log(alpha2) + tmp;
192  }
193  return lp;
194 }
195 
196 template <typename T_y, typename T_alpha, typename T_tau, typename T_beta,
197  typename T_delta>
199 wiener_lpdf(const T_y& y, const T_alpha& alpha, const T_tau& tau,
200  const T_beta& beta, const T_delta& delta) {
201  return wiener_lpdf<false>(y, alpha, tau, beta, delta);
202 }
203 
204 } // namespace math
205 } // namespace stan
206 #endif
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
Definition: return_type.hpp:20
T max(const caf::Proxy< T > &a, T b)
void check_finite(const char *function, const char *name, const T_y &y)
fvar< T > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:12
void check_bounded(const char *function, const char *name, const T_y &y, const T_low &low, const T_high &high)
double delta
Definition: runWimpSim.h:98
const double LOG_PI
Definition: constants.hpp:145
constexpr T pow(T x)
Definition: pow.h:75
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:12
const double LOG_SQRT_PI
Definition: constants.hpp:147
Float_t tmp
Definition: plot.C:36
Double_t beta
bool size_zero(T &x)
Definition: size_zero.hpp:18
Double_t K
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:12
const double LOG_TWO
Definition: constants.hpp:151
const double SQRT_2_TIMES_SQRT_PI
Definition: constants.hpp:136
fvar< T > sin(const fvar< T > &x)
Definition: sin.hpp:10
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
void check_not_nan(const char *function, const char *name, const T_y &y)
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)
Definition: wiener_lpdf.hpp:73
size_t max_size(const T1 &x1, const T2 &x2)
Definition: max_size.hpp:9
void domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
void check_positive(const char *function, const char *name, const T_y &y)
double pi()
Definition: constants.hpp:82
void check_consistent_sizes(const char *function, const char *name1, const T1 &x1, const char *name2, const T2 &x2)
const double SQRT_PI
Definition: constants.hpp:134
fvar< T > ceil(const fvar< T > &x)
Definition: ceil.hpp:11