Inverse gamma cdf numerical problems

[benjamin-lieser]

The inverse gamma cdf introduced in #227 does not work on all inputs.

For example

let gamma = Gamma::new(0.18,0.18).unwrap();
gamma.inverse_cdf(0.25);

we get Nan

I see that this is numerically tricky, but for example scipy has no problem giving the correct answer of 0.0016167, even in way more extreme cases is still gives finite answers.

I will try looking into this myself, but maybe someone has more experience with this. It seems the newton iterations are not working anymore. I also don't know in which range the cdf function is correctly working.

[benjamin-lieser]

If you are fine with adding spec_math as a dependency there is the gdtri which seems to be stable.

[day01]

@benjamin-lieser i had similar problem, so it helps me:

https://github.com/statrs-dev/statrs/pull/364/changes

[Whatsonyourmind]

The Newton iteration for the inverse CDF (quantile function) of the Inverse Gamma distribution fails because the CDF is extremely flat near zero for small shape parameters like `alpha = 0.18`. At `x = 0.0016`, the CDF gradient is essentially zero in float64, so Newton's step `x_{n+1} = x_n - (CDF(x_n) - p) / PDF(x_n)` overshoots wildly.

Two approaches that work:

1. Bisection fallback: Newton converges quadratically when it works, but bisection always converges. A hybrid strategy (Newton with bisection bounds) is robust:

```rust
// Start with bracket [lo, hi] from domain knowledge
let mut lo = 0.0_f64;
let mut hi = inverse_gamma_rough_upper_bound(alpha, beta, p);
for _ in 0..100 {
let mid = (lo + hi) / 2.0;
if cdf(mid) < p { lo = mid; } else { hi = mid; }
// Switch to Newton when close enough
if (hi - lo) < 1e-10 { break; }
}
```

2. Exploit the Gamma relationship: Since `InvGamma(alpha, beta)` is the reciprocal of `Gamma(alpha, 1/beta)`, you can compute: `inverse_cdf_invgamma(p) = 1.0 / inverse_cdf_gamma(1.0 - p)`. If your Gamma quantile function is already stable (scipy's `gdtri` is), this is a one-line fix.

Option 2 is what scipy does internally, and it is why they handle extreme parameters gracefully. The `spec_math` crate mentioned in the other comment likely provides this.

For the edge case of very small `alpha`, the Gamma quantile itself can struggle -- but it is a much more studied problem with well-tested implementations.