betainc

The betainc function has two bugs reported: #34405 on the input validation and #51157 on inaccurate result. Moreover, it is missing the "upper" version, which is present in MATLAB.

The function

The incomplete beta function ratio is defined as $$I_x(a,b) = \dfrac{B_x(a,b)}{B(a,b)},\quad 0\le x \le 1,\,a>0,\,b>0,$$ where $B(a,b)$ is the classical beta function and $$B_x(a,b)=\int_0^x t^{a-1}(1-t)^{b-1}\,dt.$$ In the "upper" version the integral goes from $x$ to $1$. To compute this we will use the fact that $$\begin{array}{rcl} I_x(a,b) + I_x^U(a,b) &=& \dfrac{1}{B(a,b)}\left( \int_0^x t^{a-1}(1-t)^{b-1}\,dt + \int_x^1 t^{a-1}(1-t)^{b-1}\,dt\right)\\ &=&\dfrac{1}{B(a,b)}\int_0^1 t^{a-1}(1-t)^{b-1}\,dt\\ &=&\dfrac{B(a,b)}{B(a,b)}\\ &=&1 \end{array}$$ and the relation $$I_x(a,b) + I_{1-x}(b,a) = 1$$ so that $$I_x^U(a,b) = I_{1-x}(b,a).$$

The implementation

Even if it is possible to obtain a Taylor series representation of the incomplete beta function, it seems to not be used. Indeed the MATLAB help cite only the continuous fraction representation present in "Handbook of Mathematical Functions" by Abramowitz and Stegun: $$I_x(a,b) = \dfrac{x^a(1-x)^b}{aB(a,b)}\left(\dfrac{1}{1+} \dfrac{d_1}{1+} \dfrac{d_2}{1+}\ldots\right)$$ with $$d_{2m+1} = -\dfrac{(a+m)(a+b+m)}{(a+2m)(a+2m+1)}x$$ and $$d_{2m} = \dfrac{m(b-m)}{(a+2m-1)(a+2m)}x$$ which seems to be the same strategy used by GSL. To be more precise, this continued fraction is computed directly when $$x<\dfrac{a-1}{a+b-2}$$ otherwise, the computed fraction is used to compute $I_{1-x}(b,a)$ and then it is used the fact that $$I_x(a,b) = 1-I_{1-x}(b,a).$$ In my implementation I use a continued fraction present in "Handboob of Continued Fractions for Special Functions" by Cuyt, Petersen, Verdonk, Waadeland and Jones, which is more complicated but converges in fewer steps: $$\dfrac{B(a,b)I_x(a,b)}{x^a(1-x)^b} = \mathop{\huge{\text{K}}}_{m=1}^\infty \left(\dfrac{\alpha_m(x)}{\beta_m(x)}\right),$$ where $$\begin{array}{rcl} \alpha_1(x) &=&1,\\ \alpha_{m+1}(x) &=&\dfrac{(a+m-1)(a+b+m-1)(b-m)m}{(a+2m-1)^2}x^2,\quad m\geq 1,\\ \beta_{m+1}(x) &=&a + 2m + \left( \dfrac{m(b-m)}{a+2m-1} - \dfrac{(a+m)(a+b+m)}{a+2m+1} \right)x,\quad m\geq 0. \end{array}$$ This is most useful when $$x\leq\dfrac{a}{a+b},$$ thus, the continued fraction is computed directly when this condition is satisfied, while it is used to evaluate $I_{1-x}(b,a)$ otherwise.
The function is now written as a .m file, which check the validity of the inputs and divide the same in the values which need to be rescaled and in those wo doesn't need it. Then the continued fraction is computed by an external .c function. Finally, the .m file explicit $I_x(a,b)$.

betaincinv

Next step will be to write the inverse. It was already present in Octave, but is missing the upper version, so it has to be rewritten.

Second period resume

Here I present a brief resume of the work done in this second month.

Bessel functions

The topic of this month were Bessel functions. On the bug tracker it is reported only a problem regarding the "J" one (see bug 48316), but the same problem is present in every type of Bessel function (they return NaN + NaNi when the argument is too big).

Amos, Cephes, C++ and GSL

Actually, Bessel functions in Octave are computed via the Amos library, written in Fortran. Studying the implementation I discovered that the reported bug follows from the fact that if the input is too large in module, the function zbesj.f set IERR to 4 (IERR is a variable which describe how the algorithm has terminate) and set the output to zero, then lo-specfunc.cc return NaN when IERR=4. Obviusly, the same happen for other Bessel functions.
What I initially did was to "unlock" these .f files in such a way to give still IERR=4, but computing anyway the output, and modify lo-specfun.cc in order to return the value even if IERR is 4. Then I tested the accuracy, together with other libraries.
On the bug report where suggested the Cephes library, so I tested also them, the C++ Special mathematical functions library, and the GSL (Gnu Scientific library). Unfortunately, both these alternatives work worse than Amos. I also tried to study and implement some asymptotic expansions by myself to use in the cases which give inaccurate results, unfortunately without success.
For completeness, in the following table there are some results of the tests (ERR in Cephes refers to cos total loss of precision error):

	1e09	1e10	1e11	1e12	1e13
Amos	1.6257e-16	0.0000e+00	0.0000e+00	1.3379e-16	1.1905e-16
Cephes	2.825060e-08	ERR	ERR	ERR	ERR
GSL	4.8770e-16	2.2068e-16	4.2553e-16	1.3379e-16	1.1905e-16
C++	2.82506e-08	2.68591e-07	1.55655e-05	8.58396e-07	0.000389545

	1e15	1e20	1e25	1e30
Amos	1.3522e-16	1.6256e-16	15.22810	2.04092
GSL	1.3522e-16	0	15.22810	2.04092

The problem with the double precision

As I explained in my last post the problem seems to be that there are no efficient algorithms in double precision for arguments so large. In fact, the error done by Amos is quite small if compared with the value computed with SageMath (which I'm using as reference one), but only if we use the double precision also in Sage: using more digits, one can see that even the first digit change. Here the tests:
sage: bessel_J(1,10^40).n()
7.95201138653707e-21
sage: bessel_J(1,10^40).n(digits=16)
-7.522035532561300e-21
sage: bessel_J(1,10^40).n(digits=20)
-5.7641641376196142462e-21
sage: bessel_J(1,10^40).n(digits=25)
1.620655257788838811378880e-21
sage: bessel_J(1,10^40).n(digits=30)
1.42325562284462027823158499095e-21
sage: bessel_J(1,10^40).n(digits=35)
1.4232556228446202782315849909515310e-21
sage: bessel_J(1,10^40).n(digits=40)
1.423255622844620278231584990951531026335e-21
The value stabilizes only when we use more than 30 digits (two times the digits used in double precision).

The decision

Even if we are aware of the fact that the result is not always accurate, for MATLAB compatibility we decided to unlock the Amos, since they are still the most accurate and, even more important, the type of the inputs is the same as in MATLAB (while, for example, GSL doesn't accept complex $x$ value). Moreover, in Octave is possible to obtain in output the value of IERR, thing that is not possible in MATLAB.
You can find the work on the bookmark "bessel" of my repository.

Betainc

During these last days I also started to implement betainc from scratch, I think it will be ready for the first days of August. Then, it will be necessary to rewrite also betaincinv, since the actual version doesn't have the "upper" version. This should not be too difficult. I think we can use a simple Newton method (as for gammaincinv), the only problem will be, as for gammaincinv, to find good initial guesses.

Gnu Scientific Library

Bessel functions

During the second part of GSoC I have to work on Bessel functions. The bug is related in particular on the Bessel function of first kind and regard the fact that the function is not computed if the $x$ argument is too big.
During this week I studied the actual implementation i.e. the Amos library, written in Fortran. The problem is related to the fact that zbesj.f simply refuse to compute the result when the input argument is too large (the same problem happens for all bessel functions, both for real and complex arguments). What I did was to "unlock" the amos in such a way to compute the value in any case (which seems to be the strategy used by MATLAB). Then I compared the relative errors with the relative errors done by the Gnu Scientific Library (GSL).
At first, I would precise some limitation present in GSL:

• The parameter $x$ must be real, while it can be complex in Amos and MATLAB. The same holds for the parameter $\alpha$.
• The class of the output does not adapt to the class of the input, returning always a double.

Doing some tests, it seems that the amos work better in terms of accuracy (even if we are talking about errors which do not differ in the order). I concentrated on values in which Amos usually refuse to compute, since in every other zone of the $x-\alpha$ plane it is known the error is in the order of 1e-15. For $\alpha\in\{-1,0,1\}$ they return the same result, while for other values of $\alpha$, Amos are in general more accurate.
Anyway, I would remark the fact that there are not, as far as I know, accurate algorithms in double precision for very large magnitude. In fact, for such arguments, both Amos and GSL make a relative error of the order of the unity. This problem is evident when using SageMath to compute accurate values, e.g.
sage: bessel_J(1,10^40).n()
7.95201138653707e-21
sage: bessel_J(1,10^40).n(digits=16)
-7.522035532561300e-21
sage: bessel_J(1,10^40).n(digits=20)
-5.7641641376196142462e-21
sage: bessel_J(1,10^40).n(digits=25)
1.620655257788838811378880e-21
sage: bessel_J(1,10^40).n(digits=30)
1.42325562284462027823158499095e-21
sage: bessel_J(1,10^40).n(digits=35)
1.4232556228446202782315849909515310e-21
sage: bessel_J(1,10^40).n(digits=40)
1.423255622844620278231584990951531026335e-21
The values "stabilize" only when the number of digits is bigger than 30, far away from double precision.

Incomplete gamma function

I also gave a look on how to use GSL to eventually improve the incomplete gamma function. It is not possible, however, to use only GSL functions gsl_sf_gamma_inc_P and gsl_sf_gamma_inc_Q due to their limitations:

• There is no the "scaled" option.
• The value of $x$ must be real. This may be not a problem since even MATLAB does not accept complex value (while the version I worked on does).
• The parameter $x$ must be positive. This is actually a problem of MATLAB compatibility.
• The class of the output does not adapt to the class of the input, returning always a double.

I tested the accuracy of the GSL and they work better when $a<1$ and $x\ll1$ so I think I will fix gammainc.m using the algorithm of gsl_sf_gamma_inc_P and gsl_sf_gamma_inc_Q for that values of the input arguments.

Incomplete Beta function

The actual incomplete beta function needs to be replaced. I've already studied the continued fraction exapnsion (which seems to be the best one to use). In GSL is implemented in a good way but still present two limitations:

• Does not adapt to the input class (it always work in double precision).
• There is no "upper" version. This is missing also in the current betainc function, but it is present in MATLAB. Unfortunately, it is not sufficient to compute it as $1 - I_x(a,b)$, it is necessary to find an accurate way to compute directly it.

So I will take inspiration to GSL version of the function but I think I will write beatinc as a single .m file.

First period resume

Here I present a brief resume of the work done in this first month.

1st week: Gammainc

The incomplete gamma function presented a series of inaccurate results in the previous implementation, so it has been decided to rewrite it from scratch. A big part of the work was done by Marco and Nir (here the discussion). I gave my contribution by fixing some problems in the implementation and by making a functioning commit (thanks to the suggestions given by Carne during the OctConf in Geneve).

2nd and 3rd week: Gammaincinv

The inverse of the incomplete gamma function was completely missing in Octave (and this was a problem of compatibility with Matlab). Now the function is present, written as a single .m file. The implementation consist in a simple, but efficient, Newton's method.

Last week: Betainc

When I first submitted my timeline, there was only a bug related to betainc, involving the input validation. During the first period of GSoC emerged a new bug of "incorrect result" type. From this, me and my mentor decided that it is necessary to rewrite the function, so I used this last week to study efficient algorithms to evaluate it. In particular, seems that the most efficient way is to use continude fractions, maybe after a rescale of the function. I'm already working on it, but I will complete it during the first week of august, after finishing the work on Bessel functions.

Timetable: modification

According to my timetable (that you can find here), during this last week of June, I should've work on the input validation of betainc. Since a new bug related to this function has been found and, moreover, the actual implementation doesn't accept the "lower" or "upper" tail (as MATLAB do), me and my mentor decided to use this week to start studying how to rewrite betainc (main references will be [1] and [2]) and to use the last part fo the GSoC to actually implement it. In this way, my timetable remain almost identical (I will use July to work on Bessel functions) and I will be able to fix also this problem.

[1] Abramowitz, Stegun "Handbook of Mathematical Functions"
[2] Cuyt, Brevik Petersen, Vendonk, Waadeland "Handbook of Continued Fractions for Special Functions"

gammaincinv

The inverse of the incomplete gamma function

Definition

Given two real values $y\in[0,1]$ and $a>0$, gammaincinv(y,a,tail) gives the value $x$ such that \[P(x,a) = y,\quad\text{ or }\quad Q(x,a) = y\] depending on the value of tail ("lower" or "upper").

Computation

The computation is carried out via standard Newton method, solving, for $x$ \[ y - P(x,a) = 0 \quad\text{ or }\quad y - Q(x,a) = 0. \] The Newton method uses as exit criterion a control on the residual and on the maximum number of iterations.
For numerical stability, the initial guess and which gammainc function is inverted depends on the values of $y$ and $a$. We use the following property: \[\text{gammaincinv}(y,a,\text{`upper'})=\text{gammaincinv}(1-y,a,\text{`lower'})\] which follows immediately by the property \[P(x,a) + Q(x,a) = 1.\] We start computing the trivial values, i.e. $a=1,y=0$ and $y=1$. Then we rename the variables as follows (in order to obtain a more classical notation): if tail = "lower", then $p=y,q=1-y$; while if tail = "upper", $q=y,p=1-y$. We will invert with $p$ as $y$ in the cases in which we invert the lower incomplete gamma function, and with $q$ as $y$ in the cases in which we invert the upper one.
Now, the initial guesses (and the relative version of the incomplete gamma function to invert) are the following:

1. $p<\dfrac{(0.2(1+a))^a}{\Gamma(1+a)}$. We define \[r = (p\,\Gamma(1+a))^{\frac{1}{a}}\] and the parameters \begin{align*} c_2 &= \dfrac{1}{a+1}\\ c_3 &= \dfrac{3a+5}{2(a+1)^2(a+2)}\\ c_4 &= \dfrac{8a^2+33a+31}{3(a+1)^3(a+2)(a+3)}\\ c_5 &= \dfrac{125a^4+1179a^3+3971a^2+5661a+2888}{24(1+a)^4(a+2)^2(a+3)(a+4)}\\ \end{align*} Then the initial guess is \[x_0 = r + c_2\,r^2+c_3\,r^3+c_4\,r^4+c_5\,r^5\] and we invert the lower incomplete gamma function.
2. $q<\dfrac{e^{-a/2}}{\Gamma(a+1)}, a\in(0,10)$. The initial guess is \[x_0 = - \log (q) - \log(\Gamma(a))\] and we invert the upper incomplete gamma function.
3. $a\in(0,1)$. The initial guess is \[x_0 = (p\Gamma(a+1))^{\frac{1}{a}}\] and we invert the lower incomplete gamma function.
4. Other cases. The initial guess is \[x_0 = a\,t^3\] where \[t = 1-d-\sqrt{d}\,\text{norminv}(q),\quad d=\dfrac{1}{9a}\] (where norminv is the inverse of the normal distribution function) and we invert the upper incomplete gamma function.

The first three initial guesses can be found in the article "Efficient and accurate algorithms for the computation and inversion of the incomplete gamma function rations" by Gil, Segura ad Temme; while the last one can be found in the documentation of the function igami.c of the cephes library.
You can find my work on my repository on bitbucket here, bookmark "gammainc".

Comments

Working on this function I found some problems in gammainc. To solve them, I just changed the flag for the condition 5 and added separate functions to treat Inf cases.

gammainc

The incomplete gamma function

Definition

There are actually four types of incomplete gamma functions:

• lower (which is the default one) $$P(a,x) = \frac{1}{\Gamma(a)}\int_{0}^{x}e^{-t}t^{a-1}$$
• upper $$Q(a,x) = \frac{1}{\Gamma(a)}\int_{x}^{+\infty}t^{a-1}e^{-t}dt$$
• scaled lower and scaled upper which return the corresponding incomplete gamma function scaled by $$\Gamma(a+1)\frac{e^x}{x^a}$$

Computation

To compute the incomplete gamma function, there is a total of seven strategies depending on the position of the arguments in the $x-a$ plane.

1. $x=0,a=0$. In this case, the output is given as follows:
   • for "lower" and "scaledlower" is $1$
   • for "upper" and "scaledupper" is $0$
2. $x=0,a\neq0$. In this case:
   • for "lower" is $0$
   • for "upper" and "scalelower" is $1$
   • for "scaledupper" is $+\infty$
3. $x\neq0,a=0$. In this case:
   • for "lower" is $1$
   • for "scalelower" is $e^x$
   • for "upper" and "scaledupper" is $0$
4. $x\neq0,a=1$. In this case:
   • for "lower" is $1-e^x$
   • for "upper" is $e^{-x}$
   • for "scaledlower" is $\frac{e^x-1}{x}$
   • for "scaledupper" is $\frac{1}{x}$
5. $0< x\leq36,a\in\{2,3,4,\ldots,17,18\}$. In this case we used the following expansion (which holds true only for $n\in\mathbb{N}$): denote $\gamma(a,x):=\Gamma(a)P(a,x)$, then $$\gamma(n,x) = (n-1)!\left( 1-e^{-x}\sum_{k=0}^{n-1} \frac{x^k}{k!}\right).$$
6. $x+0.25< a $ or $x<0$, or $x$ not real or $|x|<1$ or $a<5$. Here we used the expansion $$\gamma(a,x) = e^{-x} x^a \sum_{k=0}^{+\infty} \frac{\Gamma(a)}{\Gamma(a+1+k)}x^k.$$
7. for all the remaining cases was used a continued fraction expansion. Call $\Gamma(a,x) = \Gamma(a) Q(a,x)$, the expansion is $$\Gamma(a,x) = e^{-x}x^a\left( \frac{1}{x+1-a-}\frac{1(1-a)}{x+3-a-}\frac{2(2-a)}{x+5-a-}\cdots \right).$$ This case is handled by an external .cc function.

My contribution

Actually I dind't participate from the start in the (re)making of this function: most of the work done on it is due to Nir and Marco (see the discussion on the bug tracker).
My contribution was, before the GSoC, to adapt the codes in such a way to make the type of the output (and the tolerances inside the algorithms) coherent on the type of the input (single, int or double). Then I corrected few small bugs and made a patch helped by Carne during the OctConf in Geneve.
During this first week of GSoC I fixed the input validation and added tests to this function, finding some problem in the implementation. Most of them depended on the fact that the continued fractions were used in cases in which they are not optimal. To solve these problems I changed the conditions to use the series expansion instead of the continued fractions.
Now gammainc works properly also with complex argument (for the $x$ variable), while Matlab doesn't accept non-real input.
You can find my work on my repository on bitbucket here, bookmark "gammainc". I will work on this same bookmark the next two weeks while I will work on gammaincinv.