[R] mathematica -> r (gamma function + integration)
Leo Gürtler
leog at anicca-vijja.de
Tue Aug 8 00:13:19 CEST 2006
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
More information about the R-help
mailing list