bug: `LogNormal::pdf` in the left tail: 0.0 where the density is finite, then values up to 2x off

[FBruzzesi]

In the left tail, LogNormal::pdf(x) and ln_pdf(x).exp() on the same object disagree over two adjacent bands of x:

1. a band where pdf returns 0.0 and ln_pdf(x).exp() returns an ordinary f64,
2. immediately above it, a band where pdf returns a finite value that differs from the density by up to a factor of two.

Both bands sit well above the point where the density underflows, both widen with scale, and both move with location. scipy.stats.lognorm and mpmath agree with ln_pdf(x).exp() throughout.

Environment: statrs main on 128c9ae6d52c2be4574b6d8669136e344b87f65c, rustc 1.97.1, aarch64-apple-darwin.

Reproducer

use statrs::distribution::{Continuous, LogNormal};

fn main() {
    // 1. pdf returns 0.0 where the density is an ordinary double
    let d = LogNormal::new(0.0, 15.0).unwrap();
    println!("{:e}\t{:e}", d.pdf(1e-252), d.ln_pdf(1e-252).exp());

    // 2. just above that band, pdf returns a finite value about 1.9x the density
    let d = LogNormal::new(0.0, 10.0).unwrap();
    println!("{:e}\t{:e}", d.pdf(2.24e-168), d.ln_pdf(2.24e-168).exp());

    // 3. the case this came out of
    let d = LogNormal::new(0.0, 5.0).unwrap();
    println!("{:e}\t{:e}", d.pdf(1.38e-87), d.ln_pdf(1.38e-87).exp());
}

0e0                     3.047968563730997e-75
8.79927122410328e-158   4.605967253151374e-158
0e0                     2.0733810257862521e-262

Expected

case	pdf	scipy.stats.lognorm(s=sigma).pdf	mpmath, mp.dps = 60
LogNormal(0, 15), x = 1e-252	0e0	3.047969e-75	3.047969e-75
LogNormal(0, 10), x = 2.24e-168	8.799271e-158	4.605967e-158	4.605967e-158
LogNormal(0, 5), x = 1.38e-87	0e0	2.073381e-262	2.073381e-262

---

As those above seems like "lucky" findings, I asked Claude to run a full sweep. What follows is generated by it.

Claude Opus 5.0 analysis

How far the two bands reach

LogNormal(0, sigma), scanning x at 200 points per power of 10. "Width" is the span of band 1 in powers of 10: 1e-96 to 1e-84 is 12 of them.

sigma	band 1: pdf returns 0.0	width	largest density it returns 0.0 for	band 2: worst relative error
1	6.4e-18 to 1.7e-17	0.4	3.886588e-308	0.91 at x = 1.72e-17
2	5.3e-36 to 2.9e-34	1.7	1.412756e-291	0.91 at x = 2.95e-34
3	5.2e-55 to 5.0e-51	4.0	5.904035e-275	0.91 at x = 5.07e-51
5	5.8e-96 to 1.5e-84	11.4	1.273613e-241	0.91 at x = 1.50e-84
7	4.7e-141 to 4.4e-118	23.0	3.161836e-208	0.91 at x = 4.42e-118
10	5.4e-217 to 2.2e-168	48.6	4.449485e-158	0.91 at x = 2.24e-168
15	5e-324 to 3.3e-252	71.8 or more	1.975989e-74	0.97 at x = 3.31e-252
19	5e-324 to 2.9e-319	4.8 or more	1.813587e-07	0.96 at x = 2.88e-319

Reading the sigma = 10 row: pdf returns 0.0 over 48.6 powers of 10, the largest density it calls zero is 4.4e-158, and the first x above that band, 2.24e-168, is the second line of the reproducer, where pdf returns 8.80e-158 against a true 4.61e-158.

Four notes on that scan:

• Band 2 is wide too. Measured out to where pdf comes back within 1e-9 relative, it spans 0.2 powers of 10 at sigma = 1, 1.1 at sigma = 5, 2.3 at sigma = 10 and 3.4 at sigma = 15, starting at the top of band 1 and easing back to ordinary 1e-14 agreement above it.
• The sigma = 15 and sigma = 19 bands run off the bottom of the f64 range, so those two widths are floors rather than measurements.
• Below band 1 the density really does underflow (LogNormal(0, 5).pdf(1e-100) is about 2.4e-362), so 0.0 is the right answer there.
• Stepping sigma by 0.01, the largest one with a band is 19.28, and there is none at 19.5 or above. Every sigma I sampled from 0.05 upward has one, but below sigma around 0.6 the affected densities are all under 1e-313, so that end looks immaterial in practice.

location moves both bands

(location, scale)	band 1: pdf returns 0.0	width	largest density it returns 0.0 for
(0, 5)	5.8e-96 to 1.5e-84	11.4	1.273613e-241
(-100, 5)	8.2e-145 to 5.5e-128	16.8	3.394508e-198
(100, 5)	9.4e-47 to 4.0e-41	5.6	4.778570e-285

Band 2 moves with it, with worst relative errors of 0.93 and 0.90 on the two shifted rows. With location = 100 the first x that returns 0.0 is 4.0e-41, much closer to ordinary magnitudes than any location = 0 case.

Where the same scan finds nothing

LogNormal(0, 5).pdf(1e84) and LogNormal(0, 10).pdf(1e200) both return 0.0, and mpmath puts both densities below the f64 range, so the right tail is fine. Normal::pdf for std_dev in {1e-3, 1, 100}, over the same grid across both tails, agrees with ln_pdf().exp() to 4.9e-14 or better everywhere.

Two more corners of the same tail, with subnormal x

Less likely to matter than the two bands, but they came out of the same scan:

call	returns	mpmath, mp.dps = 60
LogNormal::new(0.0, 0.1).unwrap().pdf(5e-324)	NaN	density is 0 (log density -2.77e7)
LogNormal::new(0.0, 0.1).unwrap().ln_pdf(1e-323)	inf	-2.765723e7
LogNormal::new(0.0, 20.0).unwrap().pdf(5e-324)	4.736269e20	5.668494e20

The NaN and the inf show up for scale <= 0.5, and scipy.stats.lognorm returns inf for those same two calls, so mpmath is the only reference I have there. The third row is 16% off with no band 1 anywhere at that scale (same at scale = 50: 9.856067e272 against 1.179601e273), and ln_pdf().exp() matches mpmath to 2e-13 on both.

[YeungOnion]

Likely the exp(d^2) underflow.

We can probably use one of

• ln_pdf(x).exp() in the regimes where d<-300
• replace whole interval as backed by ln_pdf
• rewrite using exponent and mantissa and simplify so you can get the difference between the exponent of x and of d^2 (scaled by ln2) before exponentiating