precision of gammap()

[paulvonhippel]

I have a question about the precision of Stata's statistical functions, in particular gammap() and gammaptail(). The results are not quite lining up with Mathematica for me. The differences are small but noticeable. For example, with k=.8078858 and L=3.35227, Stata's

gammaptail(1/k, L)

returns .0544, while the equivalent function in Mathematica returns .0495.

That's about a 10% difference. I tend to think of Mathematica as a gold standard for accuracy. Does Stata use an approximation that could be off by 10% with these input?

[Kreshna Gopal (StataCorp)]

Given k = 0.8078858 and L = 3.35227,

di %20.16g gammaptail(1/k,L)

gives 0.05441228319270654.

gammaptail(a,x) is the reverse (upper tail or survivor) cumulative gamma distribution with shape parameter a. I don't know which function in Mathematica Paul used to get 0.0495. To generate the result of gammaptail(a,x) in Mathematica, one option is to use the two-parameter formulation of Mathematica's Gamma distribution, with shape parameter a and scale parameter 1, as follows:

1.0 - CDF[GammaDistribution[a,1], x]

This indeed gives the same result as Stata's gammaptail(a,x). The same equation can be derived for gammaptail() from Stata's and Mathematica's documentation. Stata documents gammaptail() in [D] Data Management manual.

We have compared Stata's computation of gammaptail(1/k,L), with k and L as defined above, to Mathematica, and two other sources: the multiprecision Boost library and the website http://keisan.casio.com. Mathematica 9.0 (Home Edition) was used to compute the expression above; it does not give more than 7 decimal places, even if we try the option to increase precision. This is how they compare:

Stata: 0.05441228319270654
Keisan: 0.05441228319270653
Boost: 0.05441228319270656
Mathematica: 0.0544123

Stata, Boost, and Keisan match up to 16 decimal places.

Note that gammaptail() is also known as the complement to the incomplete gamma function. Mathematica also defines this function as the upper incomplete gamma function, with a difference: Stata includes the normalizing constant (which is the gamma function of a) whereas Mathematica does not. This may be a source of confusion.

--Kreshna