Re: [R] Incomplete gamma function (was "integrate()", {was "mathematica -> r ..."})

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Tue 08 Aug 2006 - 21:53:55 EST

>>>>> "MM" == Martin Maechler <maechler@stat.math.ethz.ch> >>>>> on Tue, 8 Aug 2006 09:55:50 +0200 writes:

 >>>>> "Leo" == Leo Gürtler <leog@anicca-vijja.de>  >>>>> on Tue, 08 Aug 2006 00:13:19 +0200 writes:

     Leo> Dear R-list,
     Leo> I try to transform a mathematica script to R.

     Leo> #######relevant part of the Mathematica script
     Leo> (* p_sv *)
     Leo> dd = NN (DsD - DD^2);
     Leo> lownum = NN (L-DD)^2;
     Leo> upnum  = NN (H-DD)^2;
     Leo> low = lownum/(2s^2);
     Leo> up  = upnum/(2s^2);
     Leo> psv = NIntegrate[1/(s^NN) Exp[-dd/(2s^2)]
     Leo>    (Gamma[1/2,0,up] + Gamma[1/2,0,low]),{s,sL,sH}, MinRecursion-> 3];
     Leo>    PSV = psv/Sqrt[2NN];
     Leo> Print["------------- Results ------------------------------------"];
     Leo> Print[" "];
     Leo> Print["p(sv|D_1D_2I)   = const. ",N[PSV,6]];
     Leo> ########

     Leo> # R part
     Leo> library(fOptions)

     Leo> ###raw values for reproduction
     Leo> NN <- 58
     Leo> dd <- 0.411769
     Leo> lownum <- 20.81512
     Leo> upnum <- 6.741643
     Leo> sL <- 0.029
     Leo> sH <- 0.092
     Leo> ###

     Leo> integpsv <- function(s) { 1 / (s^NN) * exp(-dd / (2 * s^2)) *
     Leo>    ( (igamma((upnum/(2*s^2)),1/2) - igamma(0,1/2) ) +
     Leo>    (igamma((lownum/(2*s^2)),1/2) - igamma(0,1/2) ) )
     Leo> }
     Leo> psv <- integrate(integpsv, lower=sL, upper=sH)
     Leo> PSV <- psv$value / sqrt(2*NN)
     Leo> print("------------- Results ------------------------------------\n")
     Leo> print(paste("p(sv|D_1D_2I)   = const. ",PSV, sep=""))


     Leo> The results of variable "PSV" are not the same.

     Leo> In mathematica -> PSV ~ 2.67223e+47
     Leo> with rounding errors due to the initial values, in R -> PSV ~ 1.5e+47

     Leo> I am not that familiar with gamma functions and integration, thus I 
     Leo> assume there the source of the problem can be located.

    MM> Yes.
    MM> A few remarks

    MM> 1) No need to use package "fOptions" and igamma(); 
    MM> standard R's  pgamma() is all you need
    MM> {igamma() has added value only for *complex* arguments!}

    MM> 2) igamma(0, 1/2) == pgamma(0, 1/2) == 0 , so you can really     MM> drop them from your integrand.

    MM> integpsv <- function(s) { 
    MM>    1 / (s^NN) * exp(-dd / (2 * s^2)) *
    MM>    ( pgamma(upnum/(2*s^2), 1/2) + pgamma(lownum/(2*s^2), 1/2) )
    MM> }

    

    [............]

    MM> However, if I experiment, using integrate() in two parts, or using many other
    MM> numerical integration approximators,
    MM> all methods give ( your 'psv', not PSV )

    MM> integrate(integpsv, lower=sL, upper=sH)

    MM> a value of 1.623779e+48 (which leads to your PSV of 1.5076e+47)

    MM> Could it be that you are not using the same definition of     MM> incomplete gamma in Mathematica and R ?

Offlist, Leo sent me Mathematica's definition of

   Gamma[a, z0, z1] := integral_z0^z1 t^(a-1) exp(-t) dt

Now if you compare this with what ?pgamma (not ?gamma !) tells you, namely that R uses (Abramowitz and Stegun's definition of the incomplete gamma function)

   pgamma(x, a) = 1/ Gamma(a) * integral_0^x t^(a-1) exp(-t) dt

which has a normalizing factor: In your case above, it is Gamma(1/2) with its well-known value of sqrt(pi).

And indeed, if you multiply the current result by sqrt(pi), you get what you want -- and did get from Mathematica:

> (1.623779e+48 / sqrt(2*NN)) * sqrt(pi)
[1] 2.672224e+47

Regards,
Martin Maechler, ETH Zurich



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Tue Aug 08 22:00:07 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 08 Aug 2006 - 22:17:40 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.