Search code examples
gnuplot

Gnuplot: Incomplete beta function


Gnuplot's ibeta function seems to be implemented in a way that fails on regular inputs close to the boundary of the domain, but not consistently for closer approaches.

For example, I get the following results (on Gnuplot 5.2 patchlevel 8):

gnuplot> print ibeta(.1,1e-3,.93), ibeta(.1,1e-3,.94), ibeta(.1,1e-3,.95)   
0.0123146537078641 0.0124763579765399
gnuplot> print ibeta(.1,1e-3,.93), ibeta(.1,1e-3,.94), ibeta(.1,1e-3,.95) 
                                                       ^
         undefined value

or, even weirder,

gnuplot> print ibeta(.1,5e-3,.93)   
0.0591743782874705
gnuplot> print ibeta(.1,5e-3,.95)   

gnuplot> print ibeta(.1,5e-3,.95) 
               ^
         undefined value

gnuplot> print ibeta(.1,5e-3,.99)   
0.0685440281786021

Is there a way to change some sensitivity threshold or whatnot to get gnuplot to compute closer to the boundaries?


Solution

  • I'm afraid not. The current gnuplot implementation of ibeta uses a continued fraction approximation with known limitation in the valid domain. From the current documentation:

    gnuplot> help ibeta
     The `ibeta(p,q,x)` function returns the incomplete beta function of the real
     parts of its arguments. p, q > 0 and x in [0:1].  If the arguments are
     complex, the imaginary components are ignored.  The function is approximated by
     the method of continued fractions (Abramowitz and Stegun, 1964).
     The approximation is only accurate in the region x < (p-1)/(p+q-2). 
    

    The parameters you show are outside this restricted domain.

    There is now a development branch of gnuplot whose focus is on updating support for complex functions, newer algorithms with higher precision, and an expanded set of special functions. I don't expect this to land in a new release version any time soon, but I will add ibeta() to the list of candidate functions for updating if it isn't already there.