Skip to content

Fix the underflow and gradient errors of the Wiener functions #3442

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
avehtari wants to merge 11 commits into
base: develop
Choose a base branch
from
Open

Fix the underflow and gradient errors of the Wiener functions #3442

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
ce963b4
Fix the loop condition and the list of split boxes in hcubature
avehtari Oct 11, 2026
518f9d8
Use tolerances relative to the density in the 5-parameter wiener_lpdf
avehtari Oct 11, 2026
897aa57
Compute logMill with erfcx for nonnegative arguments
avehtari Oct 11, 2026
7b255e5
Fix the signs of the small-time a and w gradients in wiener4_lcdf_unnorm
avehtari Oct 11, 2026
0b9cd2b
Compute wiener4_lccdf_unnorm without cancellation where the ccdf is s…
avehtari Oct 11, 2026
2592200
Integrate d/dsw of the 7-parameter wiener_lpdf within the limits of t…
avehtari Oct 11, 2026
fac6370
Limit the density tolerance of the 7-parameter wiener_lpdf to 0.9e-6
avehtari Oct 11, 2026
c54752e
Scale the integrands of the 7-parameter wiener_lpdf by a log shift to…
avehtari Oct 11, 2026
1a6b6fd
Compute d/dsw and d/dst0 of the 7-parameter wiener_lpdf under the int…
avehtari Oct 11, 2026
ad0b746
Fix the accuracy, the missing mass and the underflow of the 7-paramet…
avehtari Oct 11, 2026
36b88a5
Add Wiener tests for the other regions of the underflow, small-densit…
avehtari Oct 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 8 additions & 14 deletions stan/math/prim/functor/hcubature.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -433,10 +433,8 @@ inline auto hcubature(const F& integrand, const ParsTuple& pars, const int dim,
= (dim == 1) ? 15 : 1 + 4 * dim + 2 * dim * (dim - 1) + std::pow(2, dim);

auto evals_per_box = numevals;
Scalar error = err;
Scalar val = result;

if ((error <= fmax(reqRelError * fabs(val), reqAbsError))
if ((err <= fmax(reqRelError * fabs(result), reqAbsError))
|| (numevals >= maxEval)) {
return result;
}
Expand All @@ -453,8 +451,8 @@ inline auto hcubature(const F& integrand, const ParsTuple& pars, const int dim,
err_vec.reserve(numevals);
err_vec.push_back(err);
while ((numevals < maxEval)
&& (error > fmax(reqRelError * fabs(val), reqAbsError))
&& fabs(val) < INFTY) {
&& (err > fmax(reqRelError * fabs(result), reqAbsError))
&& fabs(result) < INFTY) {
auto err_idx = get_largest_box_idx(err_vec);
auto&& box = ms[err_idx];
auto w = (box.b_[box.kdiv_] - box.a_[box.kdiv_]) / 2;
Expand Down Expand Up @@ -483,23 +481,19 @@ inline auto hcubature(const F& integrand, const ParsTuple& pars, const int dim,
box_t box2(box.a_, std::move(mb), result_2, kdivide_2);
result += result_1 + result_2 - box.I_;
err += err_1 + err_2 - err_vec[err_idx];
ms[err_idx].I_ = 0;
err_vec[err_idx] = 0;
ms.push_back(std::move(box1));
// The two halves replace the box. A box that stayed in the list with
// error 0 would be split again once all other errors are 0, and its
// integral would then be counted twice.
ms[err_idx] = std::move(box1);
err_vec[err_idx] = err_1;
ms.push_back(std::move(box2));
err_vec.push_back(err_1);
err_vec.push_back(err_2);
numevals += 2 * evals_per_box;
}
result = 0.0;
err = 0.0;

for (auto&& box : ms) {
result += box.I_;
}
for (auto err_i : err_vec) {
err += err_i;
}
return result;
} // hcubature

Expand Down
Loading
Loading