Inverse cdf inconsistency with cdf for the gamma distibution

[RonBarakBackal]

First of all, thank you for making this crate! I have been working in R and MATLAB and searched for a Rust crate.

When I run the cdf function of gamma with parameters shape = 3, rate = 0.5, to x = 20.5567 in R, MATLAB, or Rust (using statrs) the answer is identical:
0.997796914589247

But inversing it in R and Matlab results in very accurate 20.5567000...
while running this code :

use statrs::distribution::{Gamma, Continuous,ContinuousCDF};
use statrs::statistics::Distribution;
let gamma_dist = Gamma::new(3.0,0.5).unwrap();
let x = gamma_dist.cdf(20.5567);
let t = gamma_dist.inverse_cdf(x);
println!("{}",x);
println!("{}",t);

gives for the inverse:
20.556948....
Which is relatively high error

Thanks a lot and have a great day!
Ron

[YeungOnion]

Thinking we could implement up to convergence with a higher upper iteration bound, or perhaps we could change to another algorithm:

• brent bisection
• newton raphson as we know the derivative, just need a way to handle not going out of $p \in (0,1)$ and potentially causing a logical loop, perhaps bisect a few steps first?

ref: scipy implementation