Hypergeometric::pmf returns 0.0 or NaN for large populations

[samoht]

pmf divides three binomial coefficients held in f64. Once one of them passes f64::MAX the division happens between infinities, and nothing is signalled. ln_pmf, cdf and mean on the same object are correct.

let d = Hypergeometric::new(1030, 1, 515).unwrap();
d.pmf(0); // 0.0, exact value 1/2
d.pmf(1); // 0.0, exact value 1/2

let d = Hypergeometric::new(20000, 200, 300).unwrap();
d.pmf(0); // NaN, exact value 0.047931510683835526
d.pmf(3); // NaN, exact value 0.22687643066364876

C(1030, 515) exceeds f64::MAX while C(1029, 515) does not, so the first case overflows only the denominator and gives 0.0. Both coefficients overflow in the second, giving NaN.

The expected values above, computed with arbitrary-precision integers:

python3 -c "from math import comb; from fractions import Fraction as F; h=lambda N,K,n,x: F(comb(K,x)*comb(N-K,n-x),comb(N,n)); print([float(h(*a)) for a in [(1030,1,515,0),(1030,1,515,1),(20000,200,300,0),(20000,200,300,3)]])"

[RektPunk]

Hi! I was taking a look at this issue, and I think I have a working fix for the NaN and 0.0 results.
I think it can be fixed by calculating the PMF in log-space using self.ln_pmf(x).exp(), which safely avoids the overflow:

statrs/src/distribution/hypergeometric.rs

Lines 398 to 406 in 02fdab6

	fn pmf(&self, x: u64) -> f64 {
	if x > self.draws {
	0.0
	} else {
	factorial::binomial(self.successes, x)
	* factorial::binomial(self.population - self.successes, self.draws - x)
	/ factorial::binomial(self.population, self.draws)
	}
	}

Additionally, I noticed that ln_pmf lacked bounds checking, which could lead to integer underflow panics if x > self.draws. I've added the proper bounds checks (x > draws, x > successes, x < min()) so that it gracefully returns 0.0 (or f64::NEG_INFINITY for ln_pmf) instead of panicking.

I'd love to contribute this fix. Can I go ahead and open a PR for this?