>>>>> "Leo" == Leo G?rtler <[EMAIL PROTECTED]> >>>>> 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
______________________________________________ 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.