Merge pull request #924 from stan-dev/wiener-cdf-docs

Documentation for the wiener cdf functions

Author: WardBrian

src/bibtex/all.bib

@@ -274,7 +274,7 @@ @article{turing:1936
   volume={58}, number={345-363}, pages={5}, year={1936}
 }
 
-                  
+
 @article{kucukelbir:2015,
   title={Automatic variational inference in Stan}, author={Kucukelbir, Alp and
   Ranganath, Rajesh and Gelman, Andrew and Blei, David}, journal={Advances in
@@ -1273,7 +1273,7 @@ @article{Timonen+etal:2023:ODE-PSIS
   title={An importance sampling approach for reliable and efficient inference in
   {Bayesian} ordinary differential equation models}, author={Timonen, Juho and
   Siccha, Nikolas and Bales, Ben and L{\"a}hdesm{\"a}ki, Harri and Vehtari,
-  Aki}, journal={Stat}, year={2023}, volume = 12, number = 1, pages = {e614} 
+  Aki}, journal={Stat}, year={2023}, volume = 12, number = 1, pages = {e614}
 }
 
 @article{Vehtari+etal:2024:PSIS,
@@ -1332,3 +1332,23 @@ @article{sklar1959
   journal={The Original Document and an English Translation (March 3, 2023)},
   year={1959}
 }
+
+@article{Henrich2024,
+  author={Henrich, Franziska
+  and Hartmann, Raphael
+  and Pratz, Valentin
+  and Voss, Andreas
+  and Klauer, Karl Christoph},
+  title={The Seven-parameter Diffusion Model: an Implementation in Stan for Bayesian Analyses},
+  journal={Behavior Research Methods},
+  year={2024},
+  month={Jun},
+  day={01},
+  volume={56},
+  number={4},
+  pages={3102-3116},
+  issn={1554-3528},
+  doi={10.3758/s13428-023-02179-1},
+  url={https://doi.org/10.3758/s13428-023-02179-1}
+}
+

src/functions-reference/functions_index.qmd

@@ -157,6 +157,9 @@ For a description of argument and return types, see section
 
 ## Wiener First Passage Time Distribution
 
+For an extended explanation of how to use the `wiener_lpdf` and
+`wiener_l[c]cdf_unnorm` functions, see @Henrich2024.
+
 ### Probability density function
 
 If $\alpha \in \mathbb{R}^+$, $\tau \in \mathbb{R}^+$, $\beta \in (0, 1)$,
@@ -204,8 +207,8 @@ In the case where $s_{\delta}$, $s_{\beta}$, and $s_{\tau}$ are all $0$, this si
 \frac{\alpha}{(y-\tau)^{3/2}} \exp \! \left(- \delta \alpha \beta -
 \frac{\delta^2(y-\tau)}{2}\right) \sum_{k = - \infty}^{\infty} (2k +
 \beta) \phi \! \left(\frac{(2k + \beta)\alpha }{\sqrt{y - \tau}}\right),
-\end{equation*} where $\phi(x)$ denotes the standard normal density function;  
-
+\end{equation*}
+where $\phi(x)$ denotes the standard normal density function,
 and one representation that converges fast for large reaction-time values ("large time expansion"):
 \begin{equation*}
 \text{Wiener}(y|\alpha, \tau, \beta, \delta) =
@@ -215,6 +218,56 @@ and one representation that converges fast for large reaction-time values ("larg
 \end{equation*}
 see [@Feller1968], [@NavarroFuss2009].
 
+### Cumulative distribution function
+
+For the cumulative distribution function (cdf) there also exist two expressions
+depending on the reaction time.
+
+Let $\alpha$, $\tau$, $\beta$, $\delta$, $s_{\delta}$, $s_{\beta}$, $s_{\tau}$
+and $y$ be as above.
+
+The formula for the large-time cdf of decision times (excluding the additive
+reaction time components summarized in $\tau$ for the time being) at the upper
+boundary is stated as follows:
+
+\begin{equation}
+F(y|\alpha, \beta, \delta) = P(\alpha, \beta, \delta) -
+    \exp\left(\delta\alpha(1-\beta)-\frac{\delta^2 y}{2}\right)F_l(y|\alpha,\beta,\delta),
+\end{equation}
+where $P(\alpha,\beta,\delta)$ is the probability to hit the upper boundary, defined as
+
+\begin{equation}
+P(\alpha, \beta, \delta) =
+\begin{cases}
+    \frac{1-\exp(2\delta \alpha \beta)}{\exp(-2\delta \alpha(1-\beta)) - \exp(2\delta \alpha \beta)}, & \text{for } \delta\neq 0 \\
+    \beta, & \text{for } \delta=0,
+\end{cases}
+\end{equation}
+
+and
+
+\begin{equation}
+F_l(y|\alpha, \beta, \delta) =
+    \frac{2\pi}{\alpha^2}\sum_{k=1}^{\infty}{\frac{k\sin{k\pi(1-\beta)}}{\delta^2+(k\pi)^2/\alpha^2}}\exp(-\frac{k^2\pi^2y}{2\alpha^2}).
+\end{equation}
+
+
+The formula for the small-time cdf at the upper boundary is stated as follows:
+
+\begin{equation}
+F(y|\alpha,\beta,\delta) = \exp\left(\delta \alpha(1-\beta) -\frac{\delta^2y}{2}\right)F_s(y|\alpha, \beta,\delta),
+\end{equation} where
+
+\begin{equation}
+F_s(y|\alpha,\beta,\delta) = \sum_{k=0}^{\infty}(-1)^k\phi\left(\frac{\alpha(k+\beta^{*}_k)}
+    {\sqrt{y}} \right) \times \left( R \left(\frac{\alpha(k+\beta^{*}_k)+\delta y}{\sqrt{y}} \right)
+        + R \left(\frac{\alpha(k+\beta^*_k)-\delta y}{\sqrt{y}} \right)\right),
+\end{equation}
+
+where $\beta^*_k=(1-\beta)$ for $k$ even, $\beta^*_k=\beta$ for $k$ odd, and $R$ is Mill's ratio.
+
+The cdf for the lower boundary is $F(y|\alpha,1-\beta,-\delta)$
+
 ### Distribution statement
 
 `y ~ ` **`wiener`**`(alpha, tau, beta, delta)`
@@ -240,7 +293,7 @@ Increment target log probability density with `wiener_lupdf(y | alpha, tau, beta
 
 `real` **`wiener_lpdf`**`(reals y | reals alpha, reals tau, reals beta, reals delta)`<br>\newline
 The log of the Wiener first passage time density of `y` given boundary
-separation `alpha`, non-decision time `tau`, a-priori bias `beta`, and drift
+separation `alpha`, non-decision time `tau`, starting point `beta`, and drift
 rate `delta`.
 {{< since 2.18 >}}
 
@@ -249,7 +302,7 @@ rate `delta`.
 
 `real` **`wiener_lpdf`**`(real y | real alpha, real tau, real beta, real delta, real var_delta)`<br>\newline
 The log of the Wiener first passage time density of `y` given boundary
-separation `alpha`, non-decision time `tau`, a-priori bias `beta`, drift
+separation `alpha`, non-decision time `tau`, starting point `beta`, drift
 rate `delta`, and inter-trial drift rate variability `var_delta`.
 
 Setting `var_delta` to `0` recovers the 4-parameter signature above.
@@ -260,7 +313,7 @@ Setting `var_delta` to `0` recovers the 4-parameter signature above.
 
 `real` **`wiener_lpdf`**`(real y | real alpha, real tau, real beta, real delta, real var_delta, real var_beta, real var_tau)`<br>\newline
 The log of the Wiener first passage time density of `y` given boundary
-separation `alpha`, non-decision time `tau`, a-priori bias `beta`, drift
+separation `alpha`, non-decision time `tau`, starting point `beta`, drift
 rate `delta`, inter-trial drift rate variability `var_delta`, inter-trial
 variability of the starting point (bias) `var_beta`, and inter-trial variability
 of the non-decision time `var_tau`.
@@ -273,7 +326,7 @@ Setting `var_delta`, `var_beta`, and `var_tau` to `0` recovers the 4-parameter s
 
 `real` **`wiener_lupdf`**`(reals y | reals alpha, reals tau, reals beta, reals delta)`<br>\newline
 The log of the Wiener first passage time density of `y` given boundary
-separation `alpha`, non-decision time `tau`, a-priori bias `beta`, and drift
+separation `alpha`, non-decision time `tau`, starting point `beta`, and drift
 rate `delta`, dropping constant additive terms
 {{< since 2.25 >}}
 
@@ -283,7 +336,7 @@ rate `delta`, dropping constant additive terms
 
 `real` **`wiener_lupdf`**`(real y | real alpha, real tau, real beta, real delta, real var_delta)`<br>\newline
 The log of the Wiener first passage time density of `y` given boundary
-separation `alpha`, non-decision time `tau`, a-priori bias `beta`, drift
+separation `alpha`, non-decision time `tau`, starting point `beta`, drift
 rate `delta`, and inter-trial drift rate variability `var_delta`,
 dropping constant additive terms.
 
@@ -295,14 +348,72 @@ Setting `var_delta` to `0` recovers the 4-parameter signature above.
 
 `real` **`wiener_lupdf`**`(real y | real alpha, real tau, real beta, real delta, real var_delta, real var_beta, real var_tau)`<br>\newline
 The log of the Wiener first passage time density of `y` given boundary
-separation `alpha`, non-decision time `tau`, a-priori bias `beta`, drift
+separation `alpha`, non-decision time `tau`, starting point `beta`, drift
 rate `delta`, inter-trial drift rate variability `var_delta`, inter-trial
 variability of the starting point (bias) `var_beta`, and inter-trial variability
 of the non-decision time `var_tau`, dropping constant additive terms.
 
 Setting `var_delta`, `var_beta`, and `var_tau` to `0` recovers the 4-parameter signature above.
 {{< since 2.35 >}}
 
+---
+
+**Note:** The `lcdf` and `lccdf` functions for the `wiener` distribution are
+conditional and unnormalized, meaning that the cdf does not asymptote at 1,
+but rather at the probability to hit the upper [boundary](#boundaries).
+
+Similarly, the ccdf is defined as the probability to hit the upper boundary less
+the value of the cdf, as opposed to the more typical $1 - \textrm{cdf}$.
+
+
+<!-- real; wiener_lcdf_unnorm; (real y, real alpha, real tau, real beta, real delta); -->
+\index{{\tt \bfseries wiener\_lcdf\_unnorm }!{\tt (real y, real alpha, real tau, real beta, real delta): real}|hyperpage}
+
+`real` **`wiener_lcdf_unnorm`**`(real y, real alpha, real tau, real beta, real delta)`<br>\newline
+
+The log of the cumulative distribution function (cdf) of the Wiener distribution
+of `y` given boundary separation `alpha`, non-decision time `tau`, starting point
+`beta`, and drift rate `delta`.
+
+{{< since 2.38 >}}
+
+<!-- real; wiener_lccdf_unnorm; (real y, real alpha, real tau, real beta, real delta); -->
+\index{{\tt \bfseries wiener\_lccdf\_unnorm }!{\tt (real y, real alpha, real tau, real beta, real delta): real}|hyperpage}
+
+`real` **`wiener_lccdf_unnorm`**`(real y, real alpha, real tau, real beta, real delta)`<br>\newline
+
+The log of the complementary cumulative distribution function (ccdf) of the
+Wiener distribution of `y` given boundary separation `alpha`, non-decision time
+`tau`, starting point `beta`, and drift rate `delta`.
+
+{{< since 2.38 >}}
+
+<!-- real; wiener_lcdf_unnorm; (real y, real alpha, real tau, real beta, real delta, real var_delta, real var_beta, real var_tau); -->
+\index{{\tt \bfseries wiener\_lcdf\_unnorm }!{\tt (real y, real alpha, real tau, real beta, real delta, real var\_delta, real var\_beta, real var\_tau): real}|hyperpage}
+
+`real` **`wiener_lcdf_unnorm`**`(real y, real alpha, real tau, real beta, real delta, real var_delta, real var_beta, real var_tau)`<br>\newline
+
+The log of the cumulative distribution function (cdf) of the Wiener distribution
+of `y` given boundary separation `alpha`, non-decision time `tau`, starting point
+`beta`, drift rate `delta`, inter-trial drift rate variability `var_delta`,
+inter-trial variability of the starting point (bias) `var_beta`, and inter-trial
+variability of the non-decision time `var_tau`.
+
+{{< since 2.38 >}}
+
+
+<!-- real; wiener_lccdf_unnorm; (real y, real alpha, real tau, real beta, real delta, real var_delta, real var_beta, real var_tau); -->
+\index{{\tt \bfseries wiener\_lccdf\_unnorm }!{\tt (real y, real alpha, real tau, real beta, real delta, real var\_delta, real var\_beta, real var\_tau): real}|hyperpage}
+
+`real` **`wiener_lccdf_unnorm`**`(real y, real alpha, real tau, real beta, real delta, real var_delta, real var_beta, real var_tau)`<br>\newline
+
+The log of the complementary cumulative distribution function (ccdf) of the
+Wiener distribution of `y` given boundary separation `alpha`, non-decision time
+`tau`, starting point `beta`, drift rate `delta`, inter-trial drift rate
+variability `var_delta`, inter-trial variability of the starting point (bias)
+`var_beta`, and inter-trial variability of the non-decision time `var_tau`.
+
+{{< since 2.38 >}}
 
 ### Boundaries
 
@@ -312,3 +423,19 @@ the upper boundary only. To get the result for the lower boundary, use
 \text{Wiener}(y | \alpha, \tau, 1 - \beta, - \delta)
 \end{equation*}
 For more details, see the appendix of @Vandekerckhove-Wabersich:2014.
+
+### Vectorization
+
+The 5- and 7-argument forms of the `wiener` distribution functions (listed above
+as recieving only `real`) are implemented in such a way where they can be
+fully [vectorized](conventions_for_probability_functions.qmd#prng-vectorization),
+but currently only versions that accept all `real` and all `vector` arguments
+are exposed by Stan. If there are additional signatures that would prove useful,
+please request them by [opening an issue](https://github.com/stan-dev/stanc3/issues/new).
+
+### Tolerance tuning
+
+The 5- and 7-argument forms of the `wiener` distribution functions can also
+accept an additional `data real` argument controlling the required precision of
+the gradient calculation of the function. If omitted, this defaults to `1e-4`
+for the density and `1e-8` for the cdf functions.