🚀 go-pugleaf

In article <fGNQ4.51301$WF.2761743@bgtnsc04-news.ops.worldnet.att.net>, 
leowood@worldnet.att.net says...
>
>Hi all --
>
>I need to compute the gamma probability density function  g(x;a,b)=C(a,b)
>x**(a-1) exp(-x/b) on the order of 10**10 calls (on many different
>occasions) and thus would like to use the FASTEST algorithm available.  It
>must also be double precision with 14 or 15 decimal places of accuracy.
>Advice/coded FORTRAN solutions greatly appreciated.  Other language
>solutions also welcome, if the speed and accuracy is achieved.
>
>Thanks, Leonard


The usual algorithm Statistics people use is the routine in Applied 
Statistics (available at Statlib http://lib.stat.cmu.edu). No idea how fast 
it's for your purpose.

-- 
T.S. Lim
lim@recursive-partitioning.com
www.Recursive-Partitioning.com
______________________________________________________________________
Get paid to write a review! http://recursive-partitioning.epinions.com

Leonard Howell wrote:
> 
> Hi all --
> 
> I need to compute the gamma probability density function  g(x;a,b)=C(a,b)
> x**(a-1) exp(-x/b) on the order of 10**10 calls (on many different
> occasions) and thus would like to use the FASTEST algorithm available.  It
> must also be double precision with 14 or 15 decimal places of accuracy.
> Advice/coded FORTRAN solutions greatly appreciated.  Other language
> solutions also welcome, if the speed and accuracy is achieved.
> 
> Thanks, Leonard
Some older refs, which may be useful as starting points for citation
searches:

DiDonato AR and Morris, AH, Computation of the incomplete gamma function
ratios and their inverses, ACM Trans MathSoftware 12 1986 pp 377 - 393

Kuki, H, Complex Gamma Function with Error Control, Comm ACM 15 (1972) pp
262-267

Spira R, Calculation of the Gamma Function by Stirling's formula, Math Comp
25 (1971) pp 317 - 322

I got these from the venerable NSWC Library of Mtahematical Subroutines
manual.

-- 
!--------------------------------------------------------------------------!
! Dr Andrej Panjkov, Analyst/Programmer                                    !
! Numerical Modelling Project, Coal Mine Engineering Group                 !
! CSIRO Division of Exploration and Mining                                 !
! Queensland Centre for Advanced Technologies        Tel:  +61 7 3212 4693 !
! PO Box 883, Kenmore (Brisbane),                    Fax:  +61 7 3212 4455 !
! Queensland, Australia, 4069                     Mobile:     0417 551 435 !
!--------------------------------------------------------------------------!
! Andrej.Panjkov@dem.csiro.au                          A.Panjkov@uq.net.au !
!--------------------------------------------------------------------------!
! This email is formatted for 76 character displays, using fixed fonts.    !
!--------------------------------------------------------------------------!

Leonard Howell asks about fast algorithms for computing the
gamma density function.

One fast computational strategy is to build a very large interpolation
table (i.e., very small increment between support points) and use
low-order interpolating polynomials. To build the table you can use
your choice of slow algorithms.

This might not help if you need values of the gamma function for
different values of the shape parameter. However, if the number of
values of the shape parameter is small enough (let's say not more
than a thousand) you might still be better off investing CPU in
building an interpolation table for each shape parameter value.

In general, given the vast number of function calls (specified as
10^10) you are justified in spending resources up front to shave
a little bit off each function call.

I don't have a good intuition about interpolating over both the
argument and the shape parameter; it depends, I suppose, on how
much the shape varies for different values of the shape parameter.

Finally, if you are doing some kind of simulation or Monte Carlo
integration, you might consider parallelizing your code, at least
to the extent of launching identical runs on separate CPU's.
Some kinds of simulation and MC integration are very amenable to
parallelization, and you can get a speed-up factor nearly equal to
the number of CPU's available; much better than anything you can
squeeze out of clever numerical analysis.

HTH,
Robert Dodier

PS. Sorry if this post is a duplicate -- the copy I posted earlier
today doesn't seem to have appeared.


Sent via Deja.com http://www.deja.com/
Before you buy.

Leonard Howell wrote:
> 
> Hi all --
> 
> I need to compute the gamma probability density function  g(x;a,b)=C(a,b)
> x**(a-1) exp(-x/b) on the order of 10**10 calls (on many different
> occasions) and thus would like to use the FASTEST algorithm available.  It
> must also be double precision with 14 or 15 decimal places of accuracy.
> Advice/coded FORTRAN solutions greatly appreciated.  Other language
> solutions also welcome, if the speed and accuracy is achieved.
> 
> Thanks, Leonard


I'd start off by trying netlib, GAMS and possibly statlib and
seeing what I could turn up there. 

netlib  - www.netlib.org
GAMS    - gams.nist.gov
statlib - lib.stat.cmu.edu

There are a number of Journals I'd look at, including 
"Communications of the ACM" (some of the older algorithms 
from which are on netlib) - this is a good place to look for 
code for something like your problem. Another place would be the 
Journal of the Royal Statistical Society - Series C (Applied
Statistics), which in the past was mostly algorithms, though 
more recently it had moved closer to its title.

If you have access to the Current Index to Statistics, I'd
give that a look.

Glen

It may also be that if you are evaluating over a grid of points,
you can save a lot of time. For example, if you are investigating
the function at a given x over a range of values of a and b, then 
note that there is a component in a only (1/Gamma(a) * x**(a-1)) and 
a component in b only (exp(-x/b)), and only a very simple component
in both (b**a). It would be silly to recompute all three things
for the whole grid - you could compute a vector of the term in a
(T(a), say), a vector for the term in b (V(b), say) and then
you only need to compute b**a each time: f = T(a)*V(b)*b**a.
Even better, you can save a lot of time computing that b**a term 
by the fact that 'a' is over a grid, so b**a(n+1) = K * b**a(n); 
errors will slowly build up, but you could do it on widely spaced
values of a (say, every 16th or 8th value). You might even check
the centre estimate in the middle of each range both from above
and below as a check on loss of accuracy, though more formal
estimation and control of accuracy would also be possible.

[If you were instead doing it over a grid of x, then the constant 
term is ... constant, and you get that the term in exp(-x/b) is a
constant multiple of the same term for the previous x (in which
case similar comments apply as for the b**a term).  If you have
a grid for x, a and b then you're really flying.]

Glen

It is safer to compute the natural logarithm of the Gamma density first,
because the Gamma function T(a) can easily overflow for moderate large value
of a. Writing g(x;a,b) = (x/b)**(a-1) exp(-x/b) / b / T(a), the logarithm of
this is ln(g(x;a,b)) = (a-1) * ln(x/b) - x/b - ln(b) - lnT(a). This formula
can be further optimized depending on how frequent x, a, and b will take new
values. Anyway, there is a site in Japan
http://momonga.t.u-tokyo.ac.jp/~ooura/
which provides a Fortran 90 codes to compute the logarithm of T(a). This is
very fast (because it uses "classical" Chebyshev approximation) and the
routine offers up to 14 to 16 significant digits.

Leonard Howell wrote:

> Hi all --
>
> I need to compute the gamma probability density function  g(x;a,b)=C(a,b)
> x**(a-1) exp(-x/b) on the order of 10**10 calls (on many different
> occasions) and thus would like to use the FASTEST algorithm available.  It
> must also be double precision with 14 or 15 decimal places of accuracy.
> Advice/coded FORTRAN solutions greatly appreciated.  Other language
> solutions also welcome, if the speed and accuracy is achieved.
>
> Thanks, Leonard