[martonaronvarga]

Summary

Fixes #3329.

Adds an automatic log-scaled fallback for the full Wiener density when
natural-scale integration produces a zero, subnormal, or non-finite result.
The public api is unchanged.

Tests

• Prim Wiener full tests: 23/23 passed
• Rev Wiener full tests: 5/5 passed
• Standalone header compilation passed
• cpplint passed
• Benchmarks found no measurable slowdown in ordinary support regions

Side Effects

None expected.

Release notes

Fixed underflow in the full Wiener log density and its gradients for extreme
valid parameter values.

Checklist

• Copyright holder: Marton A. Varga

• the basic tests are passing

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[martonaronvarga]

@avehtari This is the PR fixing underflows mentioned in #3329. Unfortunately, I don't have permissions to assign reviewers, so please assign yourself for the time being.

[avehtari]

This doesn't compile for forward or mixed mode.** CI job "mix fun and prob tests 5" fails on all three platforms with 88 errors from three lines:

• line 852: std::isfinite(density) with var or fvar;
• line 854: wiener7_log_params{y_value, …}, which narrows T_partials_return to double;
• line 863: xmin/xmax passed as const Eigen::VectorXd&.

[martonaronvarga]

Sorry, forgot about mixed modes. Fixed in bd47838.

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.04	0.04	1.02	1.85% faster
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.16	0.17	0.99	-0.97% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.86	0.86	1.01	0.7% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	1.0	-0.22% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.01	0.85% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	7.97	8.07	0.99	-1.32% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.74	3.74	1.0	-0.15% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	32.27	33.12	0.97	-2.62% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.38	1.01	1.34% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.95	6.95	1.0	-0.0% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	135.29	129.82	1.04	4.05% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.21	3.22	1.0	-0.24% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.6	0.61	0.99	-1.37% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.75	25.88	0.99	-0.5% slower
performance.compilation	407.17	399.05	1.02	1.99% faster
Mean result: 1.0025071235394414

Jenkins Console Log
Jenkins Build Stages
Commit hash: 44a50912ff3f2e71023ea29a6f91adf7c6b1cf87

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           52 bits physical, 57 bits virtual
CPU(s):                                  192
On-line CPU(s) list:                     0-191
Thread(s) per core:                      2
Core(s) per socket:                      48
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              25
Model:                                   17
Model name:                              AMD EPYC 9474F 48-Core Processor
Stepping:                                1
Frequency boost:                         enabled
CPU MHz:                                 1497.776
CPU max MHz:                             4114.4229
CPU min MHz:                             1500.0000
BogoMIPS:                                7189.39
Virtualization:                          AMD-V
L1d cache:                               3 MiB
L1i cache:                               3 MiB
L2 cache:                                96 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-47,96-143
NUMA node1 CPU(s):                       48-95,144-191
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Not affected
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Enhanced / Automatic IBRS; IBPB conditional; STIBP always-on; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Mitigation; Clear CPU buffers
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good amd_lbr_v2 nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 pcid sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba perfmon_v2 ibrs ibpb stibp ibrs_enhanced vmmcall fsgsbase bmi1 avx2 smep bmi2 erms invpcid cqm rdt_a avx512f avx512dq rdseed adx smap avx512ifma clflushopt clwb avx512cd sha_ni avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local user_shstk avx512_bf16 clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin cppc arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif x2avic v_spec_ctrl vnmi avx512vbmi umip pku ospke avx512_vbmi2 gfni vaes vpclmulqdq avx512_vnni avx512_bitalg avx512_vpopcntdq la57 rdpid overflow_recov succor smca fsrm flush_l1d debug_swap

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

[avehtari]

Thanks for working on this. While testing this PR, I (with help from Claude) tried an alternative: keep the existing integrator and scale the integrand by a per-observation log shift, which is a simpler change. In addition, the tests found several other defects in develop, so I propose a separate PR with all fixes. I did test also adding the additional fixes on top of this PR, but building on top the log shift approach did produce better results. Comparison with this PR (head bd47838; for double and var it gives the same values as c8b33c4, where I measured):

• wiener_full_lpdf.hpp: +162/−35 (log shift, gradients for small sw and st0, two develop fixes), against +426 in this PR.
• Integration: the existing hcubature, with a fix for two bugs in it (+8/−14), against a second adaptive integrator of about 355 lines. With the two bugs, hcubature runs to the 6000-evaluation limit once the first box fails, and can count a split box twice.
• Speed of the 7-parameter lpdf with gradients (300 ordinary points, default precision, one core): 14 times faster with sw and st0, 11 times with st0 only, 0.95 times with sw only; about 1.2 times faster at the wiener_lpdf full model returns -inf for valid-support inputs in variability branches #3329 point and scan (precision 1e-8). Most of the gain in the ordinary region comes from the hcubature fix.
• wiener_lpdf full model returns -inf for valid-support inputs in variability branches #3329 point and scan (y − t0 = 3e-6 to 3e-5, sw = 1e-8, st0 = 1e-3): both give the log density to 5e-11. d/dsw: error at most 5e-10, against −1e8 at 3 of 16 scan points (develop's integration limits) and errors of 1e-5 to 4e-5 at 7 points (develop's end-point formula). Your fallback is accurate at the other 6 points.
• Extreme grid (128 points around the wiener_lpdf full model returns -inf for valid-support inputs in variability branches #3329 locations): finite at 108 points, −inf at 20 points (true log density −4.5e4 to −3.2e7), no exceptions. Errors below 1e-7 at 104 points; at the other 4 (a narrow 2-D peak) errors of 6.4e-4, with the same log-density error in this PR. This PR: finite at 112 points, but 24 more have log-density errors (0.13 to 0.47 at 8 points, 1e-5 to 4e-5 at 16) and 12 have d/dsw errors (develop's integration limits and series tolerance); exceptions at 16 points. The log-density errors come from the absolute tolerance log(1e-12) that the fallback passes to wiener5_density and wiener5_grad_*. With a tolerance relative to the factor in front of the series, params.log_error + fmin(0.0, wiener5_compute_log_error_term(reaction_time, params.a, params.v, relative_start, params.sv)) in wiener7_log_node_value, all 24 errors are below 2e-12 and the gradients there below 2e-9 (tested on this PR's head).
• The gradients for small sw and st0 use the idea of your fallback: differentiation under the integral. The end-point formulas multiply the integration error by 1/sw (1/st0); in develop the d/dsw error at sw = 1e-8 is up to 1.1 at the default precision.
• Other develop defects fixed: 5-parameter lpdf series truncated where the density is small (value errors up to 1e-2 at ordinary points); cancellation in P − F in the 4-parameter lccdf (−inf at 8 of 300 ordinary points) and in the 7-parameter lccdf (−inf below about 1e-16 of P; now finite down to about 1e-308); sign errors in the 4-parameter lcdf gradients (also in the 7-parameter lcdf and lccdf gradients, errors up to 3e-3); 7-parameter lccdf missing the probability mass of trials with t0 + τ > y (log-ccdf errors up to 1.4); d/dsw with an absolute series tolerance when st0 = 0 (relative errors 2e3).

EDIT: new PR is #3442

[martonaronvarga]

Thanks for the detailed comparison. I integrated the directly applicable fallback improvement in 33fbfab: each fallback node now makes the Wiener-5 series tolerance relative to its prefactor. I added a regression where the previous fallback returned -38.6717878 instead of the 50-digit reference -38.699512995; the corrected node matches the reference value and active scores to 1e-10. Local checks pass: mixed 4/4, prim 24/24, rev 5/5.\n\nI did not duplicate the broader hcubature, small-sw/st0, lcdf, and lccdf changes from #3442 here. Those form a coherent alternative implementation and are better reviewed in their existing PR. Given its smaller log-shift design and broader accuracy/performance results, #3442 appears to be the preferred merge candidate; this PR can remain available as an independently tested fallback comparison and then be closed as superseded if maintainers agree.

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.04	0.04	0.99	-0.66% slower
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.17	0.16	1.03	2.98% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.86	0.85	1.0	0.22% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	1.0	0.17% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	0.99	-1.43% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	8.02	8.1	0.99	-1.07% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.78	3.7	1.02	1.98% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	33.17	32.75	1.01	1.26% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.39	0.97	-3.45% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.94	6.97	1.0	-0.46% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	137.19	131.22	1.05	4.35% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.2	3.22	0.99	-0.55% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.61	0.6	1.01	0.87% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.86	25.7	1.01	0.62% faster
performance.compilation	389.11	385.77	1.01	0.86% faster
Mean result: 1.004145656666831

Jenkins Console Log
Jenkins Build Stages
Commit hash: 44a50912ff3f2e71023ea29a6f91adf7c6b1cf87

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           52 bits physical, 57 bits virtual
CPU(s):                                  192
On-line CPU(s) list:                     0-191
Thread(s) per core:                      2
Core(s) per socket:                      48
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              25
Model:                                   17
Model name:                              AMD EPYC 9474F 48-Core Processor
Stepping:                                1
Frequency boost:                         enabled
CPU MHz:                                 1493.824
CPU max MHz:                             4114.4229
CPU min MHz:                             1500.0000
BogoMIPS:                                7189.39
Virtualization:                          AMD-V
L1d cache:                               3 MiB
L1i cache:                               3 MiB
L2 cache:                                96 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-47,96-143
NUMA node1 CPU(s):                       48-95,144-191
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Not affected
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Enhanced / Automatic IBRS; IBPB conditional; STIBP always-on; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Mitigation; Clear CPU buffers
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good amd_lbr_v2 nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 pcid sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba perfmon_v2 ibrs ibpb stibp ibrs_enhanced vmmcall fsgsbase bmi1 avx2 smep bmi2 erms invpcid cqm rdt_a avx512f avx512dq rdseed adx smap avx512ifma clflushopt clwb avx512cd sha_ni avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local user_shstk avx512_bf16 clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin cppc arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif x2avic v_spec_ctrl vnmi avx512vbmi umip pku ospke avx512_vbmi2 gfni vaes vpclmulqdq avx512_vnni avx512_bitalg avx512_vpopcntdq la57 rdpid overflow_recov succor smca fsrm flush_l1d debug_swap

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

[stevieb-bot]

Action: Not Ready for Human Review

Commit checked: 33fbfabf424c

Summary

Thank you for this PR and for following up on the review comments.

The PR adds a log-scaled fallback integration to the 7-parameter wiener_lpdf, which runs when the natural-scale integral of the density is zero, subnormal, or non-finite.
The description and comments claim that this fixes the -inf results of #3329 for the value and the gradients, and that it compiles in forward and mixed mode after bd47838.
They also say that the tests and the repository checks pass and that ordinary inputs are not slower.

I read the code and the tests at 33fbfab.
In a container without network access I ran the prim and rev Wiener full tests, the 12 mix test files for wiener_full_lpdf, the standalone header check of wiener_full_lpdf.hpp, and cpplint on the three changed files, all at the PR head.
I also copied the PR's two test files into the merge base with develop (43bc568) and ran them there to see which tests fail on develop.

The output of each run is in this gist.
Its README.md gives the command and commit behind each file.

• Claims and Evidence: 8 Confirmed, 2 Partially Confirmed, 3 Not Verified.
• Description vs Code: 4 Confirmed, 1 Partially Confirmed, 1 Contradicted.

Claims and Evidence

Each section below lists facts from the PR summary and discussion that an agent was able to confirm, partially confirm, or not confirm at all, or that a run of the code in the PR contradicted.
Please review the "Partially Confirmed", "Not Verified", and "Contradicted" facts for facts missing evidence found by this agent.

Confirmed Facts

FACT: The #3329 example returns -inf on develop

Quote:

returns -inf from the primal double-valued lpdf.

(issue #3329 by @martonaronvarga)

STATUS: Confirmed

The PR's rev regression test runs this example with var arguments.
For var, the integration in wiener_lpdf runs on the double values, so the value it returns is the double-valued lpdf.
On the merge base with the PR's test copied in, the test reports log_density.val() as -inf.

Evidence:

The rev regression test at the issue's example.
Code: wiener_full_test.cpp:407-449.

What I read: On develop, wiener_full_lpdf.hpp lines 334-335 set T_partials_return to partials_return_t of the argument types, which is double for var.
Lines 461-494 take the values with .val(i) and integrate them as T_partials_return.

What I ran: python3 runTests.py -j8 test/unit/math/rev/prob/wiener_full_test.cpp on the merge base with the PR's test copied in, where the value is -inf against the expected -4998.0685384084609.
The same test passes at the PR head.
Output in the gist: wiener_full_rev-base.out, wiener_full_rev-head.out.

---

FACT: On develop, small log densities become 0 in the integrand, so the lpdf returns log(0)

Quote:

This means very small but finite log-density values can be converted to 0.0 before or during integration, after which the public lpdf returns log(0) = -inf.

(issue #3329 by @martonaronvarga)

STATUS: Confirmed

On develop, wiener_lpdf integrates wiener5_density<GradientCalc::ON>.
Because GradientCalc::ON is 1, wiener5_density returns exp(log_density).
The integral is then passed to log.
A log density low enough that exp underflows gives 0, so the integrand is 0 and the result is log(0).

Evidence:

The code on develop that the issue names.
Code: wiener5_lpdf.hpp:298-322 on develop, wiener_full_lpdf.hpp:486-494 on develop.

What I read: At the merge base 43bc568, wiener5_lpdf.hpp line 11 defines ON = 1.
Lines 317 and 320 return NaturalScale ? exp(log_density) : log_density.
In wiener_full_lpdf.hpp, lines 486-492 integrate wiener5_density<GradientCalc::ON> and line 494 adds log(density).

---

FACT: The prim and rev Wiener full tests pass

Quote:

Prim Wiener full tests: 23/23 passed [...]
Rev Wiener full tests: 5/5 passed

(PR description, Tests)

STATUS: Confirmed

At the PR head, wiener_full_lpdf_test.cpp has 24 tests and all 24 pass.
wiener_full_test.cpp in rev has 5 tests and all 5 pass.
The prim file has 23 tests at c8b33c4, and 33fbfabf42 added the 24th, logScaledNodeUsesRelativeSeriesTolerance.

Evidence:

The prim and rev test files for the full Wiener density.
Code: wiener_full_lpdf_test.cpp, wiener_full_test.cpp.

What I ran: python3 runTests.py -j8 on each of the two files at the PR head.
Output in the gist: wiener_full_prim-head.out, wiener_full_rev-head.out.

---

FACT: The standalone header compilation passes

Quote:

Standalone header compilation passed

(PR description, Tests)

STATUS: Confirmed

The repository's header check for the changed header, make stan/math/prim/prob/wiener_full_lpdf.hpp-test, compiles at the PR head.

Evidence:

The repository's header test target for the changed header.
Code: wiener_full_lpdf.hpp.

What I ran: make -j8 stan/math/prim/prob/wiener_full_lpdf.hpp-test at the PR head, which exits with status 0.
Output in the gist: header_test-head.out.

---

FACT: cpplint passes

Quote:

cpplint passed

(PR description, Tests)

STATUS: Confirmed

With the flags of the repository's make cpplint target, cpplint does not report any errors in the three files the PR changes.

Evidence:

The repository's cpplint configuration, run on the changed files.
Code: make/cpplint.

What I ran: lib/cpplint_2.0.2/cpplint.py with the filters from make/cpplint on wiener_full_lpdf.hpp and the two test files, at the PR head.
It prints only "Done processing" for each file and exits with status 0.
Output in the gist: cpplint-head.out.

---

FACT: bd47838 fixes the forward and mixed mode compilation

Quote:

Sorry, forgot about mixed modes.
Fixed in bd47838.

(comment by @martonaronvarga)

STATUS: Confirmed

avehtari reported that c8b33c4 does not compile in forward or mixed mode, at the std::isfinite(density) check, the wiener7_log_params initializer, and the xmin/xmax arguments. bd47838 wraps all three in value_of_rec.
At the PR head all 12 mix test files for wiener_full_lpdf compile and pass.
Together they run expect_ad over all 56 triples of the eight arguments.
expect_ad also runs fvar<double>, fvar<var>, fvar<fvar<double>>, and fvar<fvar<var>>.

Evidence:

The repository's mix tests for wiener_full_lpdf.
Code: bd47838399, wiener_full_lpdf.hpp:856-879, wiener_full_lpdf_1_test.cpp.

What I ran: python3 runTests.py -j8 on test/unit/math/mix/prob/wiener_full_lpdf_0_test.cpp through wiener_full_lpdf_11_test.cpp at the PR head.
Output in the gist: wiener_full_lpdf_mix-head.out.

---

FACT: The mixed, prim, and rev tests pass at 33fbfab

Quote:

Local checks pass: mixed 4/4, prim 24/24, rev 5/5.

(comment by @martonaronvarga)

STATUS: Confirmed

At the PR head the prim file passes 24 of 24 tests and the rev file 5 of 5.
The comment does not say which mix tests make up "4/4".
wiener_full_lpdf_0_test.cpp is the only one of the 12 mix files for wiener_full_lpdf with 4 tests.
It passes, as do the other 11.

Evidence:

The prim, rev, and mix test files for the full Wiener density.
Code: wiener_full_lpdf_test.cpp, wiener_full_test.cpp, wiener_full_lpdf_0_test.cpp.

What I ran: python3 runTests.py -j8 on the prim file, the rev file, and the 12 mix files at the PR head.
Output in the gist: wiener_full_prim-head.out, wiener_full_rev-head.out, wiener_full_lpdf_mix-head.out.

---

FACT: The corrected node matches the reference value and the active scores to 1e-10

Quote:

the corrected node matches the reference value and active scores to 1e-10.

(comment by @martonaronvarga)

STATUS: Confirmed

logScaledNodeUsesRelativeSeriesTolerance compares the log density and all eight scores of wiener7_log_node_value with the test's reference values at an absolute tolerance of 1e-10.
The test passes at the PR head.
The source of the reference values is covered by the fact check "The node test's reference value is a 50-digit reference", under Not Verified.

Evidence:

The node test in the prim test file.
Code: wiener_full_lpdf_test.cpp:488-512.

What I ran: python3 runTests.py -j8 test/unit/math/prim/prob/wiener_full_lpdf_test.cpp at the PR head, where logScaledNodeUsesRelativeSeriesTolerance passes.
Output in the gist: wiener_full_prim-head.out.

---

Partially Confirmed Facts

FACT: The PR fixes #3329

Quote:

Fixes #3329.

(PR description, Summary)

STATUS: Partially Confirmed

#3329 reports -inf from the double-valued 7-parameter lpdf at many valid inputs and gives one full example, y=1e-6, a=0.1, t0=0, w=0.0001, v=-10, sv=0.1, sw=1e-8, st0=0.001, eps=1e-8.
The PR tests this example in two places.
In the prim test, underflowRegression3329 runs it with three variability settings (sw only, st0 only, and both) and expects finite values near -4982.6 and -4998.07.
In the rev test, ProbWienerFull_underflowRegression3329 runs the case with both and checks the value and all eight gradients.
At the PR head both tests pass.
On the merge base with the PR's tests copied in, the rev test fails because the value is -inf and the gradients are NaN.
The prim test file does not compile on the merge base, because it also tests the new internal functions wiener7_log_node_value and wiener7_log_integrate.
The sw-only and st0-only cases were therefore only run at the head.
The issue also names other regions where its sweep found -inf, such as y - t0 of 1e-4, 10, and 50, a = 8, and v = +/-30.
The PR does not test any point from those regions.
Neither test states where its expected values come from.
The tests show that the value is finite and matches those numbers to 1e-8.

Evidence:

The prim and rev regression tests at the #3329 example.
Tests for the other regions the issue names, and the source of the expected values, are missing.
Code: wiener_full_lpdf_test.cpp:467-486, wiener_full_test.cpp:407-449.

What I ran: python3 runTests.py -j8 on the prim and rev test files at the PR head (24 of 24 and 5 of 5 pass), and on the merge base with the PR's two test files copied in.
On the merge base the rev test ProbWienerFull_underflowRegression3329 fails with log_density.val() equal to -inf and NaN gradients, and the prim test file does not compile.
Output in the gist: wiener_full_prim-head.out, wiener_full_rev-head.out, wiener_full_rev-base.out, wiener_full_prim-base.out.

Request: Could you add regression points from the other failure regions that #3329 lists, for example one with y - t0 = 10 or 50 and one with a = 8 or v = +/-30?
It would also help to say how the expected values in the two regression tests were computed, for example with a script or a reference implementation.
Moving the tests of the internal functions to their own file would let the issue's three cases run against develop as well.

---

FACT: The gradients no longer underflow for extreme valid parameter values

Quote:

Fixed underflow in the full Wiener log density and its gradients for extreme valid parameter values.

(PR description, Release notes)

STATUS: Partially Confirmed

The rev regression test computes all eight gradients at the #3329 example with both sw and st0 nonzero.
It compares them with finite differences of the double-valued wiener_lpdf, with a tolerance of 5e-5 times the larger of 1 and the finite difference (5e-3 for sw).
At the PR head the test passes.
On the merge base with the test copied in, the gradients are NaN for seven arguments and -1000 for st0, and the finite differences are NaN.
The test covers one point.
It does not cover the sw-only and st0-only cases, which use a one-dimensional integral (line 831), or any of the other regions the issue names.
The log-density part of this sentence is covered by the fact check of Fixes #3329.

Evidence:

The rev regression test at one point.
Gradient tests at other points and for one-dimensional integrals are missing.
Code: wiener_full_test.cpp:407-449.

What I ran: python3 runTests.py -j8 test/unit/math/rev/prob/wiener_full_test.cpp at the PR head (5 of 5 pass) and on the merge base with the PR's test copied in (the regression test fails with NaN gradients).
Output in the gist: wiener_full_rev-head.out, wiener_full_rev-base.out.

Request: Could you add gradient checks for the sw-only and st0-only cases of the #3329 example, and for at least one point from another region the issue lists?

---

Not Verified Facts

FACT: The sweep found 10,011,627 non-finite log densities among 42,611,184 valid inputs

Quote:

non-finite primal log density cases: 10,011,627

(issue #3329 by @martonaronvarga)

STATUS: Not Verified

The issue gives the counts of a sweep over a grid of valid inputs and the split of those counts over eps.
The issue says the test code and the report can be attached, but neither the issue nor the PR includes the sweep code or the grid.
The PR's tests cover one point of this grid.

Evidence:

Counts only.
The sweep code and the grid definition are missing.

Request: Could you attach the sweep code and the grid it ran, for example as a gist linked from the issue?
Running it at the PR head would also show how many of the 10,011,627 cases the fallback makes finite.

---

FACT: Benchmarks found no measurable slowdown in ordinary support regions

Quote:

Benchmarks found no measurable slowdown in ordinary support regions

(PR description, Tests)

STATUS: Not Verified

The PR does not include or link a benchmark, and the description does not name the inputs, the machine, or the number of repetitions.
Lines 856-858 take the fallback only when the hcubature result is below DBL_MIN or not finite.
For ordinary inputs, the added work is the two checks in that if.
A benchmark is still needed to back a statement about measured timings.
The stan-buildbot end-to-end table in the PR conversation lists only the stat_comp_benchmarks models, which do not use the Wiener density.

Evidence:

A sentence only.
The benchmark code and the inputs that count as ordinary support regions are missing.

Request: Could you add the benchmark that found no slowdown, with the inputs and the results at develop and at the PR head?
A google-benchmark file or a script in a gist linked from the description would be enough.

---

FACT: The node test's reference value is a 50-digit reference

Quote:

instead of the 50-digit reference -38.699512995

(comment by @martonaronvarga)

STATUS: Not Verified

logScaledNodeUsesRelativeSeriesTolerance checks the node's log density against -38.699512995035121 and the eight scores against hard-coded values.
Six scores are active, and three of the expected values are 0 (the sv score, because sv = 0, and the two inactive scores).
The PR does not include the code that computed these values at 50 digits, and the comment does not say which tool or formula was used.

Evidence:

Hard-coded values in the test.
The code that computed the references is missing.

Request: Could you add the script that computed the reference log density and scores (for example an mpmath script), or link it from the test?

---

Description vs Code

Each section below lists facts from the PR summary and discussion that an agent was able to confirm, partially confirm, or not confirm at all, or that a run of the code in the PR contradicted.
Please review the "Partially Confirmed", "Not Verified", and "Contradicted" facts for facts missing evidence found by this agent.

Confirmed Facts

FACT: A log-scaled fallback runs when the natural-scale density is zero, subnormal, or non-finite

Quote:

Adds an automatic log-scaled fallback for the full Wiener density when natural-scale integration produces a zero, subnormal, or non-finite result.

(PR description, Summary)

STATUS: Confirmed

After the existing hcubature call, lines 856-858 take the fallback unless the density is at least std::numeric_limits<double>::min() and finite.
This covers zero, subnormal, NaN, and infinite results.
The fallback calls wiener7_log_integrate (lines 417-525), which integrates in each box relative to that box's largest log density and returns the log of the total.
The fallback does not need a flag or argument to turn on.

---

FACT: The public API is unchanged

Quote:

The public api is unchanged.

(PR description, Summary)

STATUS: Confirmed

The signature of wiener_lpdf is the same as on develop (line 320 on develop and line 682 at the head), and the PR's diff to wiener_full_lpdf.hpp only adds lines.
The new types and functions (wiener7_log_params, wiener7_log_node_value, wiener7_log_integrate, and the others) are in the stan::math::internal namespace, at lines 174-526.

---

FACT: Each fallback node makes the series tolerance relative to its prefactor

Quote:

each fallback node now makes the Wiener-5 series tolerance relative to its prefactor.

(comment by @martonaronvarga)

STATUS: Confirmed

In wiener7_log_node_value, lines 250-255 add fmin(0.0, wiener5_compute_log_error_term(...)) to the log error, and the density and all wiener5_grad_* calls in that function use this log_error.
33fbfab made this change.

---

FACT: The fallback uses one log scale for the density and the derivative integrals

Quote:

Retry full-model integration with a shared log-density scale when the natural-scale density underflows.
Apply the same scale to derivative integrals so reverse-mode ratios remain finite and consistent.

(commit c8b33c42)

STATUS: Confirmed

Each box of the fallback scales its nodes by the box's largest log density (line 363), and the derivative sums at line 376 use the same scaled density.
The boxes are then added on the scale of the largest box (lines 441-453), and each partial is the ratio of a derivative sum to the density sum.
As the comment at line 349 says, each box has its own scale.
The density and derivative integrals in a box share that scale.

---

Partially Confirmed Facts

FACT: Primal and reverse-mode regressions cover the sw and st0 branches

Quote:

Add primal and reverse-mode regressions for issue #3329 across starting-point and nondecision-time variability branches.

(commit c8b33c42)

STATUS: Partially Confirmed

The prim test underflowRegression3329 runs the #3329 example with sw only, st0 only, and both (lines 476-477).
The rev test ProbWienerFull_underflowRegression3329 runs only the case with both sw and st0.
Reverse mode is not tested for the sw-only and st0-only branches.

Request: Could you add the sw-only and st0-only cases to the rev regression test, as in the prim test?

---

Contradicted Facts

FACT: The PR has no side effects

Quote:

None expected.

(PR description, Side Effects)

STATUS: Contradicted

Results only change for observations whose natural-scale density is below DBL_MIN or is not finite.
For those observations the fallback either returns a finite log density or, when its integration does not converge, throws std::domain_error with the message "log-scaled integration failed" (lines 879-897).
On develop the same inputs return -inf, NaN, +inf, or a finite value from the log of a subnormal density.
logScaledIntegrationReportsNonconvergence checks that wiener7_log_integrate reports converged == false, but the tests do not check that wiener_lpdf throws.
A model can therefore see an exception where it saw -inf before.
avehtari's comment reports exceptions at 16 of the 128 points of his extreme grid at bd47838.
Because the code adds this new behavior, the statement that no side effects are expected does not match the code.

Request: Could the Side Effects section mention that wiener_lpdf can now throw std::domain_error when the log-scaled fallback does not converge, and that subnormal densities now go through the fallback?

---

Not Checked

• The measurements in avehtari's comment of 2026-10-11 (the speed ratios, the wiener_lpdf full model returns -inf for valid-support inputs in variability branches #3329 scan, and the extreme grid of 128 points), which compare this PR with Fix the underflow and gradient errors of the Wiener functions #3442 and develop.
  They are a reviewer's measurements, not claims by the author of this PR, and the comment does not include the code behind them.
• The stan-buildbot end-to-end benchmark table.
  It lists only the stat_comp_benchmarks models, which do not use the Wiener density.
• The sentence in wiener_lpdf full model returns -inf for valid-support inputs in variability branches #3329 that the 5-parameter path is less exposed because it accumulates the log density directly.
  It is part of the issue's diagnosis, and the PR does not change the 5-parameter code.
• The checklist items in the description (idiomatic C++, doxygen, tested), which are the PR template's checkboxes.
  The new internal functions, such as wiener7_log_node_value and wiener7_log_integrate, do not have doxygen blocks.
  The testing of the change is covered by the fact checks above.
• The author's closing remark that Fix the underflow and gradient errors of the Wiener functions #3442 appears to be the preferred merge candidate.
  This is a suggestion to the maintainers about which PR to merge.