wiener_lpdf full model returns -inf for valid-support inputs in variability branches

[martonaronvarga]

Summary

A large combinatorial support sweep of the full Wiener wiener_lpdf(y, a, t0, w, v, sv, sw, st0, eps) implementation found many valid-support inputs for which the primal double-valued log density returns -inf.

This appears to be a primal numerical issue. I have discovered the issue while discussing #3322 with @Franzi2114 and @SteveBronder.

Grid Result

The sweep covered the full configured grid:

• valid support cases processed: 42,611,184 / 42,611,184
• shards completed: 512 / 512
• non-finite primal log density cases: 10,011,627
• exceptions: 0
• non-finite autodiff logp after excluding non-finite primal cases: 0
• gradient-only non-finite cases: 0
• finite-difference checks: 0 in this run

The failures were evenly distributed over supplied precision:

• eps = 1e-4: 3,337,209
• eps = 1e-6: 3,337,209
• eps = 1e-8: 3,337,209

So this does not look like the supplied precision is the problem.

Failure Locations

Failures occur in numerically difficult but valid parts of the support. They are
frequent for, among others:

• y - t0 = 1e-6, 1e-4, 10, 50
• a = 0.1, 8
• v = +/-30, +/-10

A small example from the sweep:

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

returns -inf from the primal double-valued lpdf.

Suspected Code Path

The issue appears to be in the full Wiener variability path in:

• stan/math/prim/prob/wiener_full_lpdf.hpp

The full model computes a natural-scale density through integration and then takes the log:

density = internal::wiener7_integrate(...);
log_density += log(density);

The integrand is based on:

internal::wiener5_density<GradientCalc::ON>(...)

where GradientCalc::ON == true, so this uses the natural-scale branch:

return NaturalScale ? exp(log_density) : log_density;

in stan/math/prim/prob/wiener5_lpdf.hpp.

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.

The 5-parameter path is less exposed because it accumulates the log density directly.

Expected Behavior

For valid-support inputs with positive mathematical density, the log density should
be finite. If a region is numerically too difficult, the implementation should avoid
silent natural-scale underflow to -inf or use a rescaled/log-domain integration
strategy.

Artifacts

I can attach the generated report, test code, and analysis from the grid run. The raw
failure CSVs are large, but can be sampled if useful.

Plan

I plan to submit a fix, possibly with this strategy:

1. Add scaled fixed-quadrature fallback for wiener_full when adaptive returns zero/non-finite density.
2. Make it opt-in or policy-controlled? - this is a design question that needs to be answered, since the region is quite extreme and might never come up in empirical data, but users should be notified about this behavior/policy if left unchanged.
3. Benchmark ordinary cases to confirm no measurable slowdown when fallback is not triggered.
4. Add a regression test for the failing finite-support case.