bug: `Beta` with `shape_a < 1`: `pdf(0.0)` returns `inf` or `0.0` depending on `shape_b`, `ln_pdf(0.0)` returns `-inf`

[FBruzzesi]

Same context and premise as in #432

The issue

For Beta(a, b) with a < 1 the density diverges at x = 0. Evaluating at that endpoint, I get:

• ln_pdf(0.0) is f64::NEG_INFINITY for every b I tried;
• pdf(0.0) is inf for b up to and including 80, and 0.0 from b = 81 onwards.

So pdf(0.0) returns a different answer depending on b, and pdf(0.0) and exp(ln_pdf(0.0)) do not agree in either case. The same thing happens at x = 1 when b < 1.

Environment: statrs 0.19.0 (default-features = false, features = ["std"]), rustc 1.97.1, aarch64-apple-darwin, --release.

Reproducer

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

fn main() {
    for b in [1.0, 50.0, 80.0, 81.0, 100.0, 500.0] {
        let d = Beta::new(0.1, b).unwrap();
        println!(
            "Beta(0.1, {b:>5}): pdf(0) = {:>5e}   ln_pdf(0) = {:>5e}   exp(ln_pdf(0)) = {:e}",
            d.pdf(0.0),
            d.ln_pdf(0.0),
            d.ln_pdf(0.0).exp()
        );
    }
}

Beta(0.1,     1): pdf(0) =   inf   ln_pdf(0) =  -inf   exp(ln_pdf(0)) = 0e0
Beta(0.1,    50): pdf(0) =   inf   ln_pdf(0) =  -inf   exp(ln_pdf(0)) = 0e0
Beta(0.1,    80): pdf(0) =   inf   ln_pdf(0) =  -inf   exp(ln_pdf(0)) = 0e0
Beta(0.1,    81): pdf(0) =   0e0   ln_pdf(0) =  -inf   exp(ln_pdf(0)) = 0e0
Beta(0.1,   100): pdf(0) =   0e0   ln_pdf(0) =  -inf   exp(ln_pdf(0)) = 0e0
Beta(0.1,   500): pdf(0) =   0e0   ln_pdf(0) =  -inf   exp(ln_pdf(0)) = 0e0

The upper endpoint behaves the same way:

Beta::new(1.0,   0.1).unwrap().pdf(1.0)  // inf
Beta::new(100.0, 0.1).unwrap().pdf(1.0)  // 0.0
// ln_pdf(1.0) is -inf in both cases

Expected

The density is x^(a-1) * (1-x)^(b-1) / B(a, b). For a < 1 the exponent a - 1 is negative, so x^(a-1) diverges as x -> 0, which puts the limit at +inf for both pdf and ln_pdf, independently of b. For comparison, that's what scipy returns:

>>> from scipy.stats import beta
>>> [beta.pdf(0.0, 0.1, b) for b in (1.0, 80.0, 100.0, 500.0)]
[inf, inf, inf, inf]
>>> [beta.logpdf(0.0, 0.1, b) for b in (1.0, 80.0, 100.0, 500.0)]
[inf, inf, inf, inf]

---

I would be interested in trying opening a PR 🙏🏼

[day01]

fixed in #438