Asymptotic Approximations to Betainc Function

In this report, we investigate the four approximations, namely A&S 26.5.17-20, to the betainc function.

Contents

Test AS 26.5.19

This approximation is best when a = b, or at least a ~ b. Any order of magnitude difference between a and b are discouraged. Also, it is better to have a < b rather than b < a. The following code compares this approximation to the continued fraction expansion (done by m2betainc, since a=1e-6<1e-7, note that betainc would use approximations when b > 1e-7).

clear
blist=logspace(3,10,40);
a=1e6;
for k=1:numel(blist);

    b=blist(k);
    x=a./(a+b)+ (-10/sqrt(a+b):1/10/sqrt(a+b):10/sqrt(a+b));
    x=union(x,b./(a+b)+ (-10/sqrt(a+b):1/10/sqrt(a+b):10/sqrt(a+b)));
    x=x(x<1 & x>0); x=sort(x);

    P = @(z) 0.5*erfc(-z/sqrt(2));
    Z = @(z) 1/sqrt(2*pi).* exp(-z.^2/2);

    a1=2/3*(b-a)./sqrt((a+b-2).*(a-1).*(b-1));
    a2=1/12*(1/(a-1)+1/(b-1)-13/(a+b-1));
    a3=-8/15*(a1.*(a2+3/(a+b-2)));

    y=sign(x-a/(a+b)).*sqrt(2 * ( (a+b-1).*log1p(1/(a+b-2))+(a-1)*log1p(((a-1)*(1-x)-b*x)/(a+b-1)./x) ...
                 +(b-1)*log1p((x*(b-1)-a*(1-x))/(a+b-1)./(1-x))));

    p1=P(y)-Z(y).*(a1+a2*(y-a1)/(1+a2) + a3*(1+y.^2/2)/(1+a2));

    p2= m2betainc(x,a,b);

    er=max(abs(p1-p2));
    imp=find(p1>1e-290 & p2>1e-290);
    p1=p1(imp);p2=p2(imp);x=x(imp);
    err(k)=max(abs((p1-p2)./p2));
end
loglog(blist,err), xlabel('b'), ylabel('Relative difference'),title(['a = ' num2str(a)])

Test AS 26.5.18

This approximation is best when a << b, and it seems like an order of magnitude difference or so is enough to get 4 digit accuracy. The following code tests it against m2betainc.

blist=logspace(4,9,12);
a=1e6;
for k=1:numel(blist);

    b=blist(k);

    x=a./(a+b)+ (-10/sqrt(a*b):1/100/sqrt(a*b):10/sqrt(a*b));
    x=union(x,b./(a+b)+ (-10/sqrt(a+b):1/10/sqrt(a+b):10/sqrt(a+b)));
    x=x(x>0 & x<1); x=sort(x);

    w=b*(x./(1-x));

    p1=gammainc(w,a) +  exp(-w+a*log(w)-gammaln(a)).*((a-1-w)/2/b+ ...
         1/(2*b)^2*(a.^3/2-5/3*a.^2+1.5*a-1/3-w.*(1.5*a.^2-11/6*a+1/3)+ ...
         w.^2*(1.5*a-1/6)-0.5*w.^3));

    p2=m2betainc(x,a,b);
    p0=betainc(x,a,b);

    imp=(p0>1e-290 & p1 >1e-290 & p2>1e-290);
    p0=p0(imp);p1=p1(imp);p2=p2(imp);x=x(imp);
    err0=max(abs((p0-p2)./p2));err1=max(abs((p1-p2)./p2));

    erlist(k)=err1;
end
loglog(blist,erlist); xlabel('b'), ylabel('Relative Difference')
title(['a = ' num2str(a)])

Test AS 26.5.17

This approximation is best when b<<a, and preferably at least three orders of magnitude difference exists between a and b. The following code demonstrates this by comparing it to m2betainc.

alist=logspace(7.5,12,16);
b=1e6;
for k=1:numel(alist);

    a=alist(k);
    x=a./(a+b)+ (-10/sqrt(a+b):1/10/sqrt(a+b):10/sqrt(a+b));
    x=union(x,b./(a+b)+ (-10/sqrt(a+b):1/10/sqrt(a+b):10/sqrt(a+b)));
    x=x(x<1 & x>0); x=sort(x);

    N=a+b/2-1/2;
    y=-N*log(x);

    p1=gammainc(y,b)...
        +exp(b.*log(y)-y-gammaln(b-1)).*(-1/24/N.^2.*(b+1+y)+ ...
        1/5760/N.^4.*((b-3).*(b-2).*(5*b+7).*(b+1+y)-(5*b-7).*(b+3+y).*y.^2));

    p2=m2betainc(1-x,b,a);

    er=max(abs(p1-p2));
    imp=find(p1>1e-290 & p2>1e-290);
    p1=p1(imp);p2=p2(imp);x=x(imp);

    rer(k)=max(abs((p1-p2)./p2));

end

loglog(alist,rer);xlabel('a'),ylabel('Relative Difference')
title(['b = ' num2str(b)])

Test AS 26.5.20

This is the approximation originally used by betainc, the function abetainc is a modified version of betainc which always uses this approximation and never does the continued fractions. We are going to compare it to m2betainc.

Although claimed to be good over all the ranges, this approximation seems to be doing best when a ~ b. When a<10, and b is large(or vice versa), it is reported to have large relative errors (see figure below, also the report on Geck289905).

clear
alist=logspace(3,10,30);
b=1e6;
for k=1:numel(alist);

    a=alist(k);

    x=[a./(a+b)+ (-10/sqrt(a*b):1/10/sqrt(a*b):10/sqrt(a*b))];
    x=union(x,b./(a+b)+ (-10/sqrt(a*b):1/10/sqrt(a*b):10/sqrt(a*b)));
    x=x(x<1 & x>0); x=sort(x);

    p1=abetainc(x,a,b);
    p2=m2betainc(x,a,b);

    er=max(abs(p1-p2));
    imp=find(p1>1e-290 & p2>1e-290);
    p1=p1(imp);p2=p2(imp);x=x(imp);

    rer(k)=max(abs((p1-p2)./p2));

end

loglog(alist,rer);xlabel('a'),ylabel('Relative Difference')
title(['b = ' num2str(b)])

Some Remainining Questions

In this report, we have investigated three asymptotic expansions of the betainc function. We have compared them with the continued fractions expansion. There are several issues:

Published with MATLAB® 7.4