`wiener_lpdf`: faster term counts in the 4-parameter density and less repeated work in the 5-parameter gradients

未关闭
  • 哪里不对？
#3,481 0 条评论 1 个 reaction 已指派 0 人 在 GitHub 查看

维护者通常 1 天内回复

还没有人认领这个 Issue。

评估

难度
5/5
预计耗时
一周以上
新手友好度
35/100
Issue 类型
重构
描述清晰度
基本清楚
活跃度
活跃
技术栈
cpp
领域
backend, performance

调研方向

Start with stan/math/prim/prob/wiener_lpdf.hpp for the term-count changes and stan/math/prim/prob/wiener5_lpdf.hpp for the gradient and series changes; item 7 is an unprototyped larger change. The report compares against #3442, which it says should be merged before items 3–5. Benchmark with the supplied setup and validate accuracy and gradients; the example says patches 1 and 2 have not yet been run through unit tests.

由索引模型根据 Issue 内容生成。

描述

numerics performance

(assisted by Claude)

Description

I looked at the run time of the wiener_lpdf overloads and found some changes that would make them faster without losing accuracy. With per-observation a, t0, w and v, value plus gradient on develop takes:

  • about 900 ns per observation for the 4-parameter wiener_lpdf(y | a, t0, w, v);
  • about 3,600 ns per observation for the 5-parameter wiener_lpdf(y | a, t0, w, v, sv) with sv = 0, for the same values and gradients. That's about 4× slower.

This is separate from #3442, which fixes accuracy in the same files. #3442 already fixes two things we'd notice, so they aren't repeated here: the LogResult == false check in estimate_with_err_check never fires, and the gradient series use an absolute tolerance at small densities. #3442 reports the 5-parameter lpdf at 1.09–1.13× develop's run time. The 5-parameter changes below touch the same functions, so they're probably best done on top of it.

The setup, patches and benchmark code are under Example below.

4-parameter density (wiener_lpdf.hpp)

1. The term counts and the loop counter are on the autodiff stack. kl, ks, K, k and the values used to compute them are T_return_type. With var arguments, each of them, and each k++, creates a vari. They only feed ceil and floor, so their derivative is zero.

https://github.com/stan-dev/math/blob/579253050c2302549e9ca2fb6c893d5a3c96dec1/stan/math/prim/prob/wiener_lpdf.hpp#L153-L157

Computing them in double from value_of_rec(x) gives the same values and the same gradients (to 4e-15), and makes value plus gradient 1.35–1.41× faster (patch 1 below).

2. The large-time bound doesn't match Navarro & Fuss (2009). With $x = (y - \tau)/\alpha^2$ and $\epsilon = 10^{-6}$, the number of large-time terms should be

\kappa_l = \max\left(\sqrt{\frac{-2\log(\pi x \epsilon)}{\pi^2 x}},\ \frac{1}{\pi\sqrt{x}}\right)

This is the form wiener5_density_large_reaction_time_terms uses:

https://github.com/stan-dev/math/blob/579253050c2302549e9ca2fb6c893d5a3c96dec1/stan/math/prim/prob/wiener5_lpdf.hpp#L105-L110

The 4-parameter code has $\sqrt{\pi}$ where the formula has $1/\pi^2$, and $\sqrt{x}/\pi$ where it has $1/(\pi\sqrt{x})$:

https://github.com/stan-dev/math/blob/579253050c2302549e9ca2fb6c893d5a3c96dec1/stan/math/prim/prob/wiener_lpdf.hpp#L157-L168

The small-time bound matches the paper. The square-root part of the large-time bound is about 4.2× too large, so the function uses more large-time terms than it needs. It also picks the small-time series in places where the large-time series would need fewer terms:

x develop Navarro & Fuss
0.05 small time, 4 terms small time, 4 terms
0.2 small time, 5 large time, 4
0.5 small time, 6 large time, 3
1.0 large time, 7 large time, 2
2.0 large time, 5 large time, 2

On top of patch 1, the corrected bound (patch 2) makes value plus gradient 1.6–1.7× faster than develop. It is also more accurate, because the large-time series converges quickly where it is now chosen. I compared both against the 5-parameter internals with sv = 0 and a 1e-20 series tolerance, on 20,000 random points (gradient errors relative to max(|ref|, 1)):

develop patch 1 + 2
log density (absolute) 3.4e-8 3.0e-9
d/da 6.5e-7 6.6e-8
d/dt0 1.9e-7 1.1e-7
d/dw 3.5e-7 5.0e-8
d/dv 8.9e-16 8.9e-16

I also truncated the standardised density at each bound on 19,551 (x, w) grid points. The largest truncation errors were 1.0e-8 for develop and 1.4e-9 for Navarro & Fuss, both within the 1e-6 target.

5-parameter density (wiener5_lpdf.hpp)

For each observation, develop evaluates six series:

  • the density three times: once directly, and again inside wiener5_grad_t and wiener5_grad_a;
  • the t, a and w gradient series, once each.

Times per point on 100k random points: wiener5_density 548 ns, wiener5_grad_t 1,323 ns, wiener5_grad_a 1,298 ns, wiener5_grad_w 415 ns.

3. The t gradient is computed even when nothing needs it. Unlike the other partials, deriv_y isn't guarded by if constexpr (is_autodiff_v<T_y> || is_autodiff_v<T_t0>):

https://github.com/stan-dev/math/blob/579253050c2302549e9ca2fb6c893d5a3c96dec1/stan/math/prim/prob/wiener5_lpdf.hpp#L866-L873

A call with all arguments double, such as log_lik in generated quantities, takes 1,807 ns per observation. About 1,300 ns of that is this gradient. Guarding it would make those calls about 3× faster.

4. wiener5_grad_a can be derived from wiener5_grad_t. The series depends on y and a only through u = y/a². Write LET for wiener5_compute_log_error_term, which includes the −2 log a. In both branches, log p = LET(y, a, v, w, sv) + h(y/a², w), so

\partial_a \log p = \partial_a \mathrm{LET} - \frac{2y}{a}\left(\partial_y \log p - \partial_y \mathrm{LET}\right)

Here $\partial_y \mathrm{LET}$ is density_part_one in wiener5_grad_t, and $\partial_a \mathrm{LET}$ is density_part_one in wiener5_grad_a minus 2/a. Once grad_t is known, this takes 4 ns instead of 1,298 ns. I checked it on 20k points (relative to max(|g|, 1)):

  • against wiener5_grad_a with a 1e-20 tolerance: 3.2e-14;
  • against a central difference of the log density, where the density is above 1e-6: 1.3e-9;
  • against wiener5_grad_a at the default tolerance, where the density is above 1e-6: 9.0e-11.

To keep the error guarantee of the a gradient, the t series needs the a gradient's tolerance when that one is tighter. In develop, grad_a already uses exactly grad_t's log tolerance minus log(2y/a):

https://github.com/stan-dev/math/blob/579253050c2302549e9ca2fb6c893d5a3c96dec1/stan/math/prim/prob/wiener5_lpdf.hpp#L362

https://github.com/stan-dev/math/blob/579253050c2302549e9ca2fb6c893d5a3c96dec1/stan/math/prim/prob/wiener5_lpdf.hpp#L426-L427

So the t series could use the smaller of the two tolerances. The re-estimate in estimate_with_err_check would then have to check both gradients.

5. The gradient functions compute the density again. wiener5_grad_t and wiener5_grad_a each call wiener5_density, at a slightly different tolerance from the outer call:

https://github.com/stan-dev/math/blob/579253050c2302549e9ca2fb6c893d5a3c96dec1/stan/math/prim/prob/wiener5_lpdf.hpp#L375-L376

https://github.com/stan-dev/math/blob/579253050c2302549e9ca2fb6c893d5a3c96dec1/stan/math/prim/prob/wiener5_lpdf.hpp#L441-L442

Computing the density once, at the tightest of the three tolerances, and passing it in would save one or two density series per observation. From the component timings, items 4 and 5 together should roughly halve the time for value plus gradient. I haven't measured that end to end.

6. The series are summed in log space, one term at a time. For each term, wiener5_log_sum_exp calls log, plus log_sum_exp, which itself costs an exp and a log1p:

https://github.com/stan-dev/math/blob/579253050c2302549e9ca2fb6c893d5a3c96dec1/stan/math/prim/prob/wiener5_lpdf.hpp#L175-L177

An alternative is to factor out the largest exponent and sum in linear space, as wiener5_grad_w has done since #3324. This gives the same log density (7.1e-15 over 20k points) and makes wiener5_density 30–40% faster (548 → 344 ns). Cancellation is the same as now, where the positive and negative sums are formed separately and combined with log_diff_exp. The scaling also protects against overflow and underflow the way the log scale does.

7. (Larger change) One pass for all the series. The density, t and w series use the same exponentials:

  • small time: exp(-z_k^2 / (2u));
  • large time: exp(-k^2 pi^2 u / 2) and sin/cos(k pi (1 - w)).

A single loop up to the largest of the term counts could produce all three sums. That is what autodiff effectively does through the 4-parameter series, which is probably why the 4-parameter version is faster. I haven't prototyped this.

7-parameter density (wiener_full_lpdf.hpp)

The density and each gradient are separate hcubature integrations:

https://github.com/stan-dev/math/blob/579253050c2302549e9ca2fb6c893d5a3c96dec1/stan/math/prim/prob/wiener_full_lpdf.hpp#L486-L560

Each integrand evaluates the 5-parameter density again, because the WrtLog gradient functions call wiener5_density themselves. hcubature only integrates scalar functions, so the integrands can't share evaluations. A vector-valued hcubature would allow that. I've left it out of this issue.

Suggested order
  1. Items 1 and 2. They only touch wiener_lpdf.hpp, so they don't depend on #3442.
  2. After #3442 is merged: items 3, 4 and 5 in wiener5_lpdf.hpp.
  3. Items 6 and 7 as a follow-up.

Example

Setup:

  • develop at 5792530 (v5.4.0);
  • Windows 11, GCC 14.3.0 (RTools45), -O3 -march=native -DNDEBUG;
  • tables show the median of 3 runs on one laptop, and runs varied by about ±10%.

I compiled and benchmarked patches 1 and 2 but haven't run them through the unit tests.

Value plus gradient, N = 2,000 observations, per-observation a, t0, w, v (y is data):

ns/obs speed-up vs 4-param develop
4-param, develop 909 1
4-param, patch 1 643 1.35–1.41×
4-param, patch 1 + 2 542 1.59–1.74×
5-param with sv = 0, develop 3,583 0.24–0.26×
5-param, all arguments double (value only) 1,807
Patch 1: term counts in double (no change in values) 
--- a/stan/math/prim/prob/wiener_lpdf.hpp
+++ b/stan/math/prim/prob/wiener_lpdf.hpp
@@ -43,6 +43,7 @@
 #include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/sqrt.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
 #include <algorithm>
 #include <cmath>
 #include <string>
@@ -150,17 +151,20 @@
     typename scalar_type<T_beta>::type one_minus_beta = 1.0 - beta_vec[i];
     typename scalar_type<T_alpha>::type alpha2 = square(alpha_vec[i]);
     T_return_type x = (y_vec[i] - tau_vec[i]) / alpha2;
-    T_return_type kl, ks, tmp = 0;
-    T_return_type k, K;
-    T_return_type sqrt_x = sqrt(x);
-    T_return_type log_x = log(x);
-    T_return_type one_over_pi_times_sqrt_x = 1.0 / pi() * sqrt_x;
+    T_return_type tmp = 0;
+    // The number of terms has zero derivative, so compute it on values and
+    // keep it, and the loop counter, off the autodiff stack
+    const double x_val = value_of_rec(x);
+    const double sqrt_x = std::sqrt(x_val);
+    const double log_x = std::log(x_val);
+    const double one_over_pi_times_sqrt_x = 1.0 / pi() * sqrt_x;
+    double kl, ks;
 
     // calculate number of terms needed for large t:
     // if error threshold is set low enough
-    if (PI_TIMES_WIENER_ERR * x < 1) {
+    if (PI_TIMES_WIENER_ERR * x_val < 1) {
       // compute bound
-      kl = sqrt(-2.0 * SQRT_PI * (LOG_PI_LOG_WIENER_ERR + log_x)) / sqrt_x;
+      kl = std::sqrt(-2.0 * SQRT_PI * (LOG_PI_LOG_WIENER_ERR + log_x)) / sqrt_x;
       // ensure boundary conditions met
       kl = (kl > one_over_pi_times_sqrt_x) ? kl : one_over_pi_times_sqrt_x;
     } else {
@@ -168,28 +172,28 @@
     }
     // calculate number of terms needed for small t:
     // if error threshold is set low enough
-    T_return_type tmp_expr0 = TWO_TIMES_SQRT_TWO_PI_TIMES_WIENER_ERR * sqrt_x;
+    const double tmp_expr0 = TWO_TIMES_SQRT_TWO_PI_TIMES_WIENER_ERR * sqrt_x;
     if (tmp_expr0 < 1) {
       // compute bound
-      ks = 2.0 + sqrt_x * sqrt(-2 * log(tmp_expr0));
+      ks = 2.0 + sqrt_x * std::sqrt(-2 * std::log(tmp_expr0));
       // ensure boundary conditions are met
-      T_return_type sqrt_x_plus_one = sqrt_x + 1.0;
+      const double sqrt_x_plus_one = sqrt_x + 1.0;
       ks = (ks > sqrt_x_plus_one) ? ks : sqrt_x_plus_one;
     } else {     // if error threshold was set too high
       ks = 2.0;  // minimal kappa for that case
     }
-    if (ks < kl) {   // small t
-      K = ceil(ks);  // round to smallest integer meeting error
-      T_return_type tmp_expr1 = (K - 1.0) / 2.0;
-      T_return_type tmp_expr2 = ceil(tmp_expr1);
-      for (k = -floor(tmp_expr1); k <= tmp_expr2; k++) {
-        tmp += (one_minus_beta + 2.0 * k)
-               * exp(-(square(one_minus_beta + 2.0 * k)) * 0.5 / x);
+    if (ks < kl) {                     // small t
+      const double K = std::ceil(ks);  // round to smallest integer meeting error
+      const double tmp_expr1 = (K - 1.0) / 2.0;
+      const double tmp_expr2 = std::ceil(tmp_expr1);
+      for (double k = -std::floor(tmp_expr1); k <= tmp_expr2; k++) {
+        const auto w_plus_2_k = one_minus_beta + 2.0 * k;
+        tmp += w_plus_2_k * exp(-(square(w_plus_2_k)) * 0.5 / x);
       }
-      tmp = log(tmp) - LOG_TWO_OVER_TWO_PLUS_LOG_SQRT_PI - 1.5 * log_x;
-    } else {         // if large t is better...
-      K = ceil(kl);  // round to smallest integer meeting error
-      for (k = 1; k <= K; ++k) {
+      tmp = log(tmp) - LOG_TWO_OVER_TWO_PLUS_LOG_SQRT_PI - 1.5 * log(x);
+    } else {                           // if large t is better...
+      const double K = std::ceil(kl);  // round to smallest integer meeting error
+      for (double k = 1; k <= K; ++k) {
         tmp += k * exp(-(square(k)) * (SQUARE_PI_OVER_TWO * x))
                * sin(k * pi() * one_minus_beta);
       }
Patch 2, on top of patch 1: Navarro & Fuss large-time bound 
--- a/stan/math/prim/prob/wiener_lpdf.hpp
+++ b/stan/math/prim/prob/wiener_lpdf.hpp
@@ -157,14 +157,14 @@
     const double x_val = value_of_rec(x);
     const double sqrt_x = std::sqrt(x_val);
     const double log_x = std::log(x_val);
-    const double one_over_pi_times_sqrt_x = 1.0 / pi() * sqrt_x;
+    const double one_over_pi_times_sqrt_x = 1.0 / (pi() * sqrt_x);
     double kl, ks;
 
     // calculate number of terms needed for large t:
     // if error threshold is set low enough
     if (PI_TIMES_WIENER_ERR * x_val < 1) {
       // compute bound
-      kl = std::sqrt(-2.0 * SQRT_PI * (LOG_PI_LOG_WIENER_ERR + log_x)) / sqrt_x;
+      kl = std::sqrt(-2.0 * (LOG_PI_LOG_WIENER_ERR + log_x)) / (pi() * sqrt_x);
       // ensure boundary conditions met
       kl = (kl > one_over_pi_times_sqrt_x) ? kl : one_over_pi_times_sqrt_x;
     } else {
Benchmark and checks for items 4 and 6

wiener_lpdf_b.hpp and wiener_lpdf_c.hpp are copies of wiener_lpdf.hpp with patch 1, and with patches 1 + 2. In each copy, the include guard and the function are renamed (wiener_lpdf_b, wiener_lpdf_c).

#include <stan/math/rev.hpp>
#include "wiener_lpdf_b.hpp"  // patch 1
#include "wiener_lpdf_c.hpp"  // patch 1 + 2
#include <chrono>
#include <cstdio>
#include <random>
#include <vector>
using namespace stan::math;
using internal::GradientCalc;

template <typename F>
double best_ns(F&& f, int reps, int n) {
  double best = 1e300;
  for (int r = 0; r < reps; ++r) {
    auto t0 = std::chrono::steady_clock::now();
    f();
    auto t1 = std::chrono::steady_clock::now();
    best = std::min(best, std::chrono::duration<double, std::nano>(t1 - t0).count() / n);
  }
  return best;
}

// d/da log p from d/dy log p: the series depends on y and a only through y/a^2
double grad_a_from_t(double y, double a, double v, double w, double sv, double gt) {
  const double q = 1.0 - w, sv2 = sv * sv, opsy = 1.0 + sv2 * y;
  const double dlet_da = (v * q + sv2 * q * q * a) / opsy - 2.0 / a;
  const double dlet_dy = -0.5 * (sv2 * sv2 * (y + a * a * q * q) + sv2 * (1.0 + 2.0 * a * v * q) + v * v)
                         / (opsy * opsy);
  return dlet_da - (2.0 * y / a) * (gt - dlet_dy);
}

// wiener5_density with the series summed in linear space after factoring out
// the largest exponent (as wiener5_grad_w already does)
double density_linear(double y, double a, double v, double w, double sv, double log_err) {
  const double let = internal::wiener5_compute_log_error_term(y, a, v, w, sv);
  const double n_small
      = internal::wiener5_n_terms_small_t<true, GradientCalc::OFF>(y, a, w, log_err - let);
  const double n_large = internal::wiener5_density_large_reaction_time_terms(y, a, w, log_err - let);
  const double u = y / (a * a), q = 1.0 - w;
  if (2 * n_small <= n_large) {
    const double m = -q * q / (2.0 * u);
    double s = q;
    for (int k = 1; k <= static_cast<int>(n_small); ++k) {
      const double zp = q + 2.0 * k, zm = q - 2.0 * k;
      s += zp * std::exp(-zp * zp / (2.0 * u) - m) + zm * std::exp(-zm * zm / (2.0 * u) - m);
    }
    return let - 0.5 * LOG_TWO - LOG_SQRT_PI - 1.5 * std::log(u) + std::log(std::fabs(s)) + m;
  }
  const double c = 0.5 * pi() * pi() * u;
  double s = 0.0;
  for (int k = 1; k <= static_cast<int>(n_large); ++k)
    s += k * std::sin(k * pi() * q) * std::exp(-(k * k - 1.0) * c);
  return let + std::log(std::fabs(s)) - c + LOG_PI;
}

int main() {
  // Part 1: public functions, value + gradient, reverse mode
  const int N = 2000, reps = 200;
  std::mt19937 rng(1);
  std::uniform_real_distribution<double> ua(0.8, 2.0), ut0(0.15, 0.35), uw(0.3, 0.7), uv(-1.5, 2.5),
      udt(0.15, 1.5);
  std::vector<double> y(N), a0(N), t00(N), w0(N), v0(N);
  for (int i = 0; i < N; ++i) {
    a0[i] = ua(rng); t00[i] = ut0(rng); w0[i] = uw(rng); v0[i] = uv(rng); y[i] = t00[i] + udt(rng);
  }
  std::vector<double> g_ref;
  double lp_ref = 0, t_ref = 0;
  auto run = [&](const char* name, auto&& f) {
    double lp_val = 0;
    std::vector<double> g;
    double t = best_ns([&] {
      std::vector<var> a(a0.begin(), a0.end()), t0(t00.begin(), t00.end()), w(w0.begin(), w0.end()),
          v(v0.begin(), v0.end());
      var lp = f(a, t0, w, v);
      lp.grad();
      lp_val = lp.val();
      g.clear();
      for (int i = 0; i < N; ++i) {
        g.push_back(a[i].adj()); g.push_back(t0[i].adj()); g.push_back(w[i].adj()); g.push_back(v[i].adj());
      }
      recover_memory();
    }, reps, N);
    if (g_ref.empty()) { g_ref = g; lp_ref = lp_val; t_ref = t; }
    double md = 0;
    for (size_t i = 0; i < g.size(); ++i) md = std::max(md, std::fabs(g[i] - g_ref[i]));
    std::printf("| %-46s | %7.0f | %5.2fx | %.1e | %.1e |\n", name, t, t_ref / t, std::fabs(lp_val - lp_ref), md);
  };
  run("4-param, develop", [&](auto& a, auto& t0, auto& w, auto& v) { return wiener_lpdf(y, a, t0, w, v); });
  run("4-param, patch 1", [&](auto& a, auto& t0, auto& w, auto& v) { return wiener_lpdf_b(y, a, t0, w, v); });
  run("4-param, patch 1 + 2", [&](auto& a, auto& t0, auto& w, auto& v) { return wiener_lpdf_c(y, a, t0, w, v); });
  run("5-param with sv = 0", [&](auto& a, auto& t0, auto& w, auto& v) { return wiener_lpdf(y, a, t0, w, v, 0.0); });
  volatile double sink = 0;
  double t_dbl = best_ns([&] { sink = sink + wiener_lpdf(y, a0, t00, w0, v0, 0.0); }, reps, N);
  std::printf("5-param, all arguments double (value only): %.0f ns/obs\n\n", t_dbl);

  // Part 2: 5-parameter internals on 100k random points
  const double le = std::log(1e-12);
  std::mt19937 rng2(42);
  std::uniform_real_distribution<double> uy(0.05, 2.5), ua2(0.5, 3.0), uv2(-4, 4), uw2(0.1, 0.9), usv(0.0, 1.0);
  struct Pt { double y, a, v, w, sv; };
  std::vector<Pt> pts(100000);
  for (auto& p : pts) p = {uy(rng2), ua2(rng2), uv2(rng2), uw2(rng2), (rng2() % 2) ? usv(rng2) : 0.0};
  auto per_point = [&](const char* name, auto&& f) {
    double t = best_ns([&] {
      double acc = 0;
      for (auto& p : pts) acc += f(p);
      sink = sink + acc;
    }, 5, pts.size());
    std::printf("| %-40s | %7.0f |\n", name, t);
  };
  per_point("wiener5_density", [&](const Pt& p) { return internal::wiener5_density<false>(p.y, p.a, p.v, p.w, p.sv, le); });
  per_point("density, linear-space sum", [&](const Pt& p) { return density_linear(p.y, p.a, p.v, p.w, p.sv, le); });
  per_point("wiener5_grad_t", [&](const Pt& p) { return internal::wiener5_grad_t<false>(p.y, p.a, p.v, p.w, p.sv, le); });
  per_point("wiener5_grad_a", [&](const Pt& p) { return internal::wiener5_grad_a<false>(p.y, p.a, p.v, p.w, p.sv, le); });
  per_point("grad_a from grad_t (identity)", [&](const Pt& p) { return grad_a_from_t(p.y, p.a, p.v, p.w, p.sv, 0.3); });
  per_point("wiener5_grad_w", [&](const Pt& p) { return internal::wiener5_grad_w<false>(p.y, p.a, p.v, p.w, p.sv, le); });

  // Part 3: accuracy of the identity and of the linear-space sum
  const double lt = std::log(1e-20);
  double id_tight = 0, id_fd = 0, id_def = 0, lin = 0;
  for (size_t i = 0; i < 20000; ++i) {
    const auto& p = pts[i];
    const double ga_t = internal::wiener5_grad_a<false>(p.y, p.a, p.v, p.w, p.sv, lt);
    const double gi_t = grad_a_from_t(p.y, p.a, p.v, p.w, p.sv,
                                      internal::wiener5_grad_t<false>(p.y, p.a, p.v, p.w, p.sv, lt));
    const double h = 1e-5 * p.a;
    const double fd = (internal::wiener5_density<false>(p.y, p.a + h, p.v, p.w, p.sv, lt)
                       - internal::wiener5_density<false>(p.y, p.a - h, p.v, p.w, p.sv, lt))
                      / (2 * h);
    id_tight = std::max(id_tight, std::fabs(ga_t - gi_t) / std::max(1.0, std::fabs(ga_t)));
    const double lp = internal::wiener5_density<false>(p.y, p.a, p.v, p.w, p.sv, le);
    if (lp > std::log(1e-6)) {
      id_fd = std::max(id_fd, std::fabs(gi_t - fd) / std::max(1.0, std::fabs(fd)));
      const double ga = internal::wiener5_grad_a<false>(p.y, p.a, p.v, p.w, p.sv, le);
      const double gi = grad_a_from_t(p.y, p.a, p.v, p.w, p.sv,
                                      internal::wiener5_grad_t<false>(p.y, p.a, p.v, p.w, p.sv, le));
      id_def = std::max(id_def, std::fabs(ga - gi) / std::max(1.0, std::fabs(ga)));
    }
    lin = std::max(lin, std::fabs(lp - density_linear(p.y, p.a, p.v, p.w, p.sv, le)));
  }
  std::printf("identity vs wiener5_grad_a, 1e-20 tolerance: %.1e\n", id_tight);
  std::printf("identity vs wiener5_grad_a, default tolerance, density > 1e-6: %.1e\n", id_def);
  std::printf("identity vs central difference, density > 1e-6: %.1e\n", id_fd);
  std::printf("linear-space vs log-space log density: %.1e\n", lin);
}
Accuracy check for patches 1 + 2 
// Accuracy of the 4-parameter wiener_lpdf (develop, patch 1 + 2) against the
// 5-parameter internals with sv = 0 and a 1e-20 series tolerance
#include <stan/math/rev.hpp>
#include "wiener_lpdf_c.hpp"
#include <cstdio>
#include <random>
using namespace stan::math;

struct Res { double lp, ga, gt0, gw, gv; };

template <typename F>
Res eval(F&& f, double y, double a, double t0, double w, double v) {
  var av = a, t0v = t0, wv = w, vv = v;
  var lp = f(y, av, t0v, wv, vv);
  lp.grad();
  Res r{lp.val(), av.adj(), t0v.adj(), wv.adj(), vv.adj()};
  recover_memory();
  return r;
}

int main() {
  const double lt = std::log(1e-20);
  std::mt19937 rng(7);
  std::uniform_real_distribution<double> ua(0.5, 3.0), uw(0.2, 0.8), uv(-4, 4), udt(0.05, 3.0);
  double e_dev[5] = {0}, e_new[5] = {0};
  for (int i = 0; i < 20000; ++i) {
    const double a = ua(rng), w = uw(rng), v = uv(rng), dt = udt(rng), t0 = 0.2, y = t0 + dt;
    const double ref[5] = {internal::wiener5_density<false>(dt, a, v, w, 0.0, lt),
                           internal::wiener5_grad_a<false>(dt, a, v, w, 0.0, lt),
                           -internal::wiener5_grad_t<false>(dt, a, v, w, 0.0, lt),
                           internal::wiener5_grad_w<false>(dt, a, v, w, 0.0, lt),
                           internal::wiener5_grad_v<false>(dt, a, v, w, 0.0, lt)};
    const Res d = eval([](auto&&... x) { return wiener_lpdf(x...); }, y, a, t0, w, v);
    const Res c = eval([](auto&&... x) { return wiener_lpdf_c(x...); }, y, a, t0, w, v);
    const double dv[5] = {d.lp, d.ga, d.gt0, d.gw, d.gv}, cv[5] = {c.lp, c.ga, c.gt0, c.gw, c.gv};
    for (int k = 0; k < 5; ++k) {
      const double scale = k == 0 ? 1.0 : std::max(1.0, std::fabs(ref[k]));
      e_dev[k] = std::max(e_dev[k], std::fabs(dv[k] - ref[k]) / scale);
      e_new[k] = std::max(e_new[k], std::fabs(cv[k] - ref[k]) / scale);
    }
  }
  const char* names[5] = {"log density (abs)", "d/da", "d/dt0", "d/dw", "d/dv"};
  for (int k = 0; k < 5; ++k) std::printf("| %-18s | %9.1e | %13.1e |\n", names[k], e_dev[k], e_new[k]);
}

Expected Output

  • Items 1, 3, 4, 5 and 6 should give the same values to floating-point precision.
  • Item 2 gives slightly different values, and they are more accurate (table above).
  • Expected speed-ups:
    • 4-parameter value plus gradient: about 1.4× with patch 1, about 1.65× with patches 1 + 2.
    • 5-parameter value with all arguments double: about 3× (item 3).
    • 5-parameter value plus gradient: roughly half the time (items 4 and 5, projected).
    • wiener5_density: a further 30–40% faster (item 6).
Current Version:

v5.4.0 (develop at 5792530)

主要语言
C++
星标
845
派生
225
平均合并
3 天 9 小时
30 天内合并 PR
12

环境准备

从这里开始

  1. 先读完整个 Issue，再读项目的贡献指南。
  2. 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
  3. Fork 仓库，在一个分支上完成修改。
  4. 提交 Pull Request，并在描述里引用这个 Issue 编号。

stan-dev/math 的其他 Issue

查看 stan-dev/math 的全部 Issue

相似的 Issue

更多 C++ Issue

把新 issue 发到你的邮箱

精选适合新手参与的 GitHub issue 摘要。