Yup, this is a problem. When I compile and run your program gsl_invBeta.c (which should #include <gsl/gsl_cdf.h>), I get:
Observed: gsl: betainv.c:181: <http://git.savannah.gnu.org/cgit/gsl.git/tree/cdf/betainv.c#n181> ERROR: inverse failed to converge Default GSL error handler invoked. Abort trap: 6 Expected result: 1.0 Result with more precision: 0.99999999999999999999999999980178072518730245884823... For the supplied values of alpha and beta, the CDF rises very steeply near 1 and the inverse CDF is equal to 1.0 within double precision when P>0.5 (that may not be a tight bound, but is a sufficient condition). Generic workaround: #include <stdio.h> #include <gsl/gsl_sf_gamma.h> // Computes the inverse function of gsl_sf_beta_inc by binary search. double qbeta(double p, double a, double b) { double x = 0.5; for (double delta = 0.25; delta > 0; delta /= 2) { const double beta_P = gsl_sf_beta_inc(a, b, x); const double new_x = p > beta_P ? x + delta : x - delta; if (beta_P == p || x == new_x) break; x = new_x; } return x; } int main() { double a = 5.3729222; double b = 0.018122864; double p = 0.67276126; printf("%.18g\n", qbeta(p, a, b)); printf("%.7g (expected: 0.25)\n", qbeta(0.25, 1, 1)); printf("%.7g (expected: 0.9922861)\n", qbeta(0.25, 5, 0.1)); printf("%.7g (expected: 2.954261e-200)\n", qbeta(1e-50, 0.25, 0.5)); return 0; } Hope this helps, -- mj On Fri, Jan 10, 2020 at 1:30 AM Vasu Jaganath <vadie.develo...@gmail.com> wrote: > Martin and others, > > Following up, I have a particular example for which both Q and P variants > of inverse functions fail. > > where as scipy beta.ppf (equivalent to gsl_cdf_beta_Pinv) does the right > thing and converges to 1 or nearest double. > > I am attaching a very simple example with exact same data, the gsl > function fails where as scipy function does the right thing. > > Any thoughts? any workarounds? Maybe there is a way I can specify > convergence criteria? > > Thanks, > Vasu > > > On Wed, Jan 8, 2020 at 12:42 PM Vasu Jaganath <vadie.develo...@gmail.com> > wrote: > >> Thanks Martin, >> >> I will test it out. >> >> On Tue, Jan 7, 2020 at 11:16 PM Martin Jansche <mjans...@gmail.com> >> wrote: >> >>> There are many more floating point values between 0.0 and 0.001 than >>> there are between 0.999 and 1.0. The difference between 1.0 and the next >>> smaller double value is only around 1e-16, but the next larger double value >>> after 0.0 is about 1e-303. So beta_P(0.9, 1, 17) will be necessarily >>> equivalent to 1.0 due to lack of precision, whereas beta_Q(0.9, 1, 17) will >>> be 1e-17. (Haven't tried this in GSL. You may want to try and report back.) >>> >>> On Wed, Jan 8, 2020 at 1:35 AM Vasu Jaganath <vadie.develo...@gmail.com> >>> wrote: >>> >>>> Hi all, >>>> >>>> This is probably a very silly question, >>>> >>>> I don't understand why there are two separate P and Q variants for >>>> CDFs? >>>> particularly for beta distribution? >>>> >>>> >>>> https://www.gnu.org/software/gsl/doc/html/randist.html#the-beta-distribution >>>> >>>> Thanks, >>>> Vasu >>>> >>>
