R help - Aug 2006 - mathematica -> r (gamma function + integration)

Dear R-list,

I try to transform a mathematica script to R.

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

# R part
library(fOptions)

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

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


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

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

I am not that familiar with gamma functions and integration, thus I 
assume there the source of the problem can be located.
Thanks for helping me to adjust the sript.

best wishes
leo

>>>>> "Leo" == Leo G?rtler <leog at
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.

Yes.
A few remarks

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

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

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

3) But then the problem could really be with the algorithm used in
   integrate(), and indeed if you plot your integrand

      plot(integpsv, from= sL, to = sH)

   you see that indeed your integrand looks ``almost
   constant'' in the left half --- whereas that is actually not
   true but the range of the integrand varies so dramatically
   that it ``looks like'' constant 0 upto about x= .06.

   Something which   help(integrate)   warns about.

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

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

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

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

Martin Maechler, ETH Zurich

    Leo> Thanks for helping me to adjust the sript.

    Leo> best wishes
    Leo> leo