# [R] integration problem with gamma function

Leo Gürtler leog at anicca-vijja.de
Fri Sep 1 17:02:54 CEST 2006

Dear R-list members,

I have a problem with translating a mathematica script into R. The whole
script is at the end of the email (with initial values for easy
reproduction) and can be pasted directly into R. The problematic part
(which is included below of course) is

<--- Original Mathematica --->
(* p_svbar *)
UiA  = Ni (Dsi - 2Di A + A^2)/2;
UiiA = Nii (Dsii - 2Dii A + A^2)/2;
psvbar = NIntegrate[1/(UiA^(Ni/2)) 1/(UiiA^(Nii/2))
Gamma[Ni/2,UiA/(sH^2),UiA/(sL^2)]
Gamma[Nii/2,UiiA/(sH^2),UiiA/(sL^2)],{A,L,H},
MinRecursion->3];
PSVbar = psvbar/(4 Log[sH/sL]);
Print["p(sÂv|D_1D_2I)  = const. ",N[PSVbar,6]];
</--->

<--- translation to R --->
integpsvbar <- function(A)
{
# Mathematica: gamma[a,z0,z1] = gamma[a,z1] - gamma[a,z0]
UiA <- Ni*(Dsi-2*Di*A+A^2)/2
UiiA <- Nii*(Dsii-2*Dii*A+A^2)/2
1/UiA^(Ni/2) *
1/UiiA^(Nii/2) *
( pgamma(UiA/(sL^2), Ni/2) - pgamma(UiA/(sH^2), Ni/2) ) *
( pgamma(UiiA/(sL^2), Nii/2) - pgamma(UiiA/(sH^2), Nii/2) )
}

psvbar <- integrate(integpsvbar, lower=L, upper=H)
psvbar\$value
PSVbar <- psvbar\$value / (4*log(sH/sL))
PSVbar <- PSVbar * sqrt(pi) * sqrt(pi)      #note: two times pgamma in
integration function above... correct it with two times multiplied by
'sqrt(pi)' due to differences in defining incomplete gamma function
between R and mathematica
print(paste("p(sÂv|D_1D_2I)  = const. ",PSVbar,sep=""))
</--->

According to Mathematica

psvbar ~ 1.72026 * 10^20
PSVbar ~ 5.24827 * 10^19

The R part produces

psvbar ~ 1.824023 * 10^13
PSVbar ~ 17482463305549.6

The script produces proper results untill the problematic part, so it
seems quite reasonable that no other error before is responsible for the
difference of the results.
I assume that I did something wrong with the Gamma function.

The "generalized incomplete gamma function" in Mathematica (Wolfram,
p.363f.) is defined as

Gamma[a, z0, z1] := Gamma[a, z1] - Gamma[a, z0]

with

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

08-August 2006) for which Martin Mächler (08-08-2006 13:53) pointed out,
that the result has to be multiplied by "sqrt(pi)" because of (->
manpage pgamma) differences in formulating the incomplete gamma function
between Mathematica and R. However, I tried this in various ways and -
at least for me - it does not help. Thus, I assume it is another issue.

Thank you very much,

best wishes
leo gürtler

now the R script to reproduce (can be pasted directly into R):

####################
# R-Portierung aus Mathematica (Urban Studer, 90er)
# Ursprung: G.L. Bretthorst "On the difference of means"
# zuerst: 12-06-05
# zuletzt: 21-06-06

# --------------------------------------------------
# Success rates and integration bounds in the case
# of the (conservative) Bayes-Laplace prior
# --------------------------------------------------

SucRatesIntBounds <- function(Ni, Si, Nii, Sii, smin)
{
# BEGIN BLOCK
# necessary variables:
# {Nmax, Nmin}

### defintion of constants
Di <- (Si + 1) / (Ni + 2)
si <- sqrt(Di * (1 - Di) / (Ni + 3))
Dii <- (Sii + 1) / (Nii + 2)
sii <- sqrt(Dii * (1 - Dii) / (Nii + 3))

Nmax <- max(Ni, Nii)
Nmin <- min(Ni, Nii)
sL <- floor(1000 * sqrt((Nmax+1) / (Nmax+3)) / (Nmax + 2)) / 1000
sH <- ceiling(1000 / (2 * sqrt(Nmin + 3))) / 1000
L <- floor(100 * (smin + 1) / (Nmax + 2)) / 100
H <- 1 - L

return(c(Di, si, Dii, sii, sL, sH, L, H))
}
#END BLOCK
# --------------------------------------------------

Si <- 11
Ni <- 15
Sii <- 10
Nii <- 16
smin <- 0
res1 <- SucRatesIntBounds(Ni, Si, Nii, Sii, smin)
res1
#return(c(Di, si, Dii, sii, sL, sH, L, H))

names(res1) <- c("Di","si","Dii","sii","sL","sH","L","H")
res1

Di <- res1[1]
si <- res1[2]
Dii <- res1[3]
sii <- res1[4]
sL <- res1[5]
sH <- res1[6]
L <- res1[7]
H <- res1[8]

# --------------------------------------------------
# Begin: ON THE DIFFERENCE IN MEANS
# --------------------------------------------------

#DiffinMeans <- function(Ni, Di, si, Nii, Dii, sii, L, H, sL, sH)
## BEGIN BLOCK
#{

# necessary variables in the function
#{ NN, DD, Dsi, Dsii, DsD, ss,
#  dd, lownum, upnum, low, up, psv, PSV,
#  zz, lowinum, upinum, lowiinum, upiinum, psbarv, PSbarV,
#  psvbar, PSVbar, psbarvbar, PSbarVbar, cc,
#  sv, sbarv, svbar, sbarvbar,
#  samemeans, diffmeans, samevars, diffvars, diffsets },

### defintion of constants
NN <- Ni+Nii
DD <- (Ni * Di + Nii * Dii) / NN
Dsi <- (Ni-1) / Ni * si^2 + Di^2
Dsii <- (Nii-1) / Nii * sii^2 + Dii^2
DsD <- (Ni * Dsi + Nii * Dsii) / NN
ss <- sqrt(NN * (DsD - DD^2) / (NN-1))

### descriptive statistics
print(" \n------------- Data ---------------------------------------\n")
print(paste("N_1 = ",Ni ," :     Mean_1 ± s_1    = ", Di," ± ",si, sep=""))
print(paste("N_2 = ",Nii," :     Mean_2 ± s_2    = ", Dii," ± ",sii,
sep=""))
print(paste("N   = ", NN ," :  Mean_comb ± s_comb = ", DD," ± ",ss, sep=""))
print(" ")
print(paste("s_L = ",sL,", s_H = ",sH,";  L = ",L, ", H = ",H,"  (smin =
",smin,")",sep=""))
if(L < DD)
{
print(paste("L - Mean_comb < 0 : ", (L<DD), "\n", "             (->
'+'-sign between Gamma-fcts o.k.)", sep=""))
} else print(paste("L - Mean_comb < 0 : ", (L<DD), "\n", "
(-> '+'-sign between Gamma-fcts false!)", sep=""))
print(paste(" \n ", sep=""))

#(* p_sv *)
dd <- NN * (DsD - DD^2)
lownum <- NN * (L-DD)^2
upnum  <- NN * (H-DD)^2

#low = lownum / (2*s^2)
#up  = upnum / (2*s^2)
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) )
}
psv <- integrate(integpsv, lower=sL, upper=sH)
#???    (* + if (L-DD) < 0 *)
PSV <- psv\$value / sqrt(2*NN) * sqrt(pi)    # normalizing factor due to R
# Gamma(1/2)= sqrt(pi)
# psv ~ 1.39848 * 10^20
# ~ 3.44715 * 10 ^20
# ~ 5.24827 * 10 ^20
# ~ 1.52788 * 10 ^20
print("------------- Results ------------------------------------\n")
print(paste("p(sv|D_1D_2I)   = const. ",PSV, sep=""))

### p_sbarv
zz <- Ni * (Dsi - Di^2) + Nii * (Dsii - Dii^2)
lowinum <- Ni * (L - Di)^2
upinum <- Ni * (H - Di)^2
lowiinum <- Nii * (L - Dii)^2
upiinum <- Nii * (H - Dii)^2

#lowi <- lowinum / (2*s^2)
#upi  <- upinum / (2*s^2)
#lowii <- lowiinum / (2*s^2)
#upii <- upiinum / (2*s^2)
integpsbarv <- function(s)
{
1/(s^(NN-1)) * exp(-zz/(2*s^2)) *
( pgamma(upinum/(2*s^2), 1/2) + pgamma(lowinum/(2*s^2), 1/2) ) *
( pgamma(upiinum/(2*s^2), 1/2) + pgamma(lowiinum/(2*s^2), 1/2) )
}
psbarv <- integrate(integpsbarv, lower=sL, upper=sH)
#??? (* + if (L-DD) < 0 *)
PSbarV <- psbarv\$value / (2*(H-L) * sqrt(Ni*Nii)) * sqrt(pi) * sqrt(pi)
# 2mal mit sqrt(pi) mal nehmen = "* pi", weil zei pgamma-Ausdrücke
enthalten sind in der Integration
print(paste("p(Âsv|D_1D_2I)  = const. ",PSbarV, sep=""))

###############ALL ABOVE PRODUCES PROPER RESULTS COMPARED TO
###############MATHEMATICA SCRIPT

###############PROBLEMATIC PART STARTS HERE#################
### p_svbar
###
# Mathematica
#(* p_svbar *)
#UiA  = Ni (Dsi - 2Di A + A^2)/2;
#UiiA = Nii (Dsii - 2Dii A + A^2)/2;
#psvbar = NIntegrate[1/(UiA^(Ni/2)) 1/(UiiA^(Nii/2))
#   Gamma[Ni/2,UiA/(sH^2),UiA/(sL^2)]
#   Gamma[Nii/2,UiiA/(sH^2),UiiA/(sL^2)],{A,L,H},
#    MinRecursion->3];
#PSVbar = psvbar/(4 Log[sH/sL]);
#Print["p(sÂv|D_1D_2I)  = const. ",N[PSVbar,6]];

###R trial
integpsvbar <- function(A)
{
# Mathematica: gamma[a,z0,z1] = gamma[a,z1] - gamma[a,z0]
UiA <- Ni*(Dsi-2*Di*A+A^2)/2
UiiA <- Nii*(Dsii-2*Dii*A+A^2)/2
1/UiA^(Ni/2) *
1/UiiA^(Nii/2) *
( pgamma(UiA/(sL^2), Ni/2) - pgamma(UiA/(sH^2), Ni/2) ) *
( pgamma(UiiA/(sL^2), Nii/2) - pgamma(UiiA/(sH^2), Nii/2) )
}

psvbar <- integrate(integpsvbar, lower=L, upper=H)
psvbar\$value
PSVbar <- psvbar\$value / (4*log(sH/sL)) * sqrt(pi) * sqrt(pi) # two
times sqrt(pi) to correct for differences between R and mathematica
print(paste("p(sÂv|D_1D_2I)  = const. ",PSVbar,sep=""))