Thread Subject: betainc speed or alternate implementation

I need an implementation of the incomplete beta function (Matlab implements the regularized version which fails for the arguments I need to evaluate) and after looking through the betainc.m file I noticed that it calls the continued fraction algorithm betacore. As far as I know, this algorithm is perfectly capable of computing the values I need without divercences due to regularization.

However, if I just copy-paste that subroutine into my code it runs ~100x slower than just running betainc(z,v,w): which also uses it as a subroutine. Is the compiler capable of producing that kind of speed up for continued fraction code? If not, is there a way to access betacore without having to evaluate matlab's betainc?

As an aside, using hypergeom to evaluate betainc is several orders of magnitude slower than just calling betainc directly. Is there some problem with the way Matlab handles hypergeometic functions in general?

"Rodrigo" <guerra.remove.this@physics.harvard.edu> wrote in message
news:gjrstu$pab$1@fred.mathworks.com...
>I need an implementation of the incomplete beta function (Matlab implements
>the regularized version which fails for the arguments I need to evaluate)
>and after looking through the betainc.m file I noticed that it calls the
>continued fraction algorithm betacore. As far as I know, this algorithm is
>perfectly capable of computing the values I need without divercences due to
>regularization.

Note that BETAINC doesn't always use that subfunction; it sometimes
approximates the incomplete Beta function using ERFC and/or GAMMAINC. [I
remember that approximation comes from Abramowitz & Stegun, but I'm not at
work right now so I can't look up in which chapter it's located.]

If you have a set of parameter values for which BETAINC fails (for some
definition of "fails") can you send that set of parameter values, along with
the values you expect it to return for those parameter values, to me or to
Technical Support so we can investigate if there's a better way to evaluate
the incomplete Beta function for those parameters?

> However, if I just copy-paste that subroutine into my code it runs ~100x
> slower than just running betainc(z,v,w): which also uses it as a
> subroutine. Is the compiler capable of producing that kind of speed up
> for continued fraction code? If not, is there a way to access betacore
> without having to evaluate matlab's betainc?

BETAINC may not be using that subfunction, as I said above.

> As an aside, using hypergeom to evaluate betainc is several orders of
> magnitude slower than just calling betainc directly. Is there some
> problem with the way Matlab handles hypergeometic functions in general?

The HYPERGEOM function is part of Symbolic Math Toolbox, and uses the
symbolic kernel to perform operations symbolically rather than numerically.
Performing the calculations numerically will be faster than performing the
calculations symbolically.

--
Steve Lord
slord@mathworks.com

"Steven Lord" <slord@mathworks.com> wrote in message <gjrv72$7k5$1@fred.mathworks.com>...
>
> "Rodrigo" <guerra.remove.this@physics.harvard.edu> wrote in message
> news:gjrstu$pab$1@fred.mathworks.com...
> >I need an implementation of the incomplete beta function (Matlab implements
> >the regularized version which fails for the arguments I need to evaluate)
> >and after looking through the betainc.m file I noticed that it calls the
> >continued fraction algorithm betacore. As far as I know, this algorithm is
> >perfectly capable of computing the values I need without divercences due to
> >regularization.
>
> Note that BETAINC doesn't always use that subfunction; it sometimes
> approximates the incomplete Beta function using ERFC and/or GAMMAINC. [I
> remember that approximation comes from Abramowitz & Stegun, but I'm not at
> work right now so I can't look up in which chapter it's located.]
>
> If you have a set of parameter values for which BETAINC fails (for some
> definition of "fails") can you send that set of parameter values, along with
> the values you expect it to return for those parameter values, to me or to
> Technical Support so we can investigate if there's a better way to evaluate
> the incomplete Beta function for those parameters?
>
> > However, if I just copy-paste that subroutine into my code it runs ~100x
> > slower than just running betainc(z,v,w): which also uses it as a
> > subroutine. Is the compiler capable of producing that kind of speed up
> > for continued fraction code? If not, is there a way to access betacore
> > without having to evaluate matlab's betainc?
>
> BETAINC may not be using that subfunction, as I said above.
>
> > As an aside, using hypergeom to evaluate betainc is several orders of
> > magnitude slower than just calling betainc directly. Is there some
> > problem with the way Matlab handles hypergeometic functions in general?
>
> The HYPERGEOM function is part of Symbolic Math Toolbox, and uses the
> symbolic kernel to perform operations symbolically rather than numerically.
> Performing the calculations numerically will be faster than performing the
> calculations symbolically.
>
> --
> Steve Lord
> slord@mathworks.com
>

Hi Steve,
Thanks for the quick response. My problem with betainc is that I have to evaulate it for 0<z<~3, -2<~w<=0 , and 0<x<1. For the non-regularized version of betainc this is no problem as long as z is not too close to 1, but the powerlaw divergence as z->1 ruins the normalization.

I am not sure which algorithm is most efficient and I will try to parse through and see if gammainc/erfc are being used for these parameter values. Obviously I was comparing betacore and betainc for z,w>0 (1, .1 to be precise) with 0<x<1, so the relevant algorithm for w<0 may well be very slow.

While I have your attention, is there any numerical implementation of the F2 Appell Hypergeometric function? My using the betainc came about because I was unable to find a suitable implementation. Finding/writing one would allow me a lot more flexibility in the model these functions are being called from.

Thanks, Rodrigo

I am a bit confused now. I corrected a few unrelated bugs in the code (they multiplied the output of betacore) and the speed is not only comparable to betainc, its actually faster. I'll have to go back to make sure I didn't do any silly, but they seem to be numerically identical (up to round off in the routine they're embedded in).

I will compare the output I get from this funcion with the black-box answers I get from Mathematica. If I get a significant discrepancy I'll write back.

Thanks again, Rodrigo