[R] integrate() problem {was "mathematica -> r ..."}
Martin Maechler
maechler at stat.math.ethz.ch
Tue Aug 8 09:55:50 CEST 2006
>>>>> "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
More information about the R-help
mailing list