[R] integration problem with gamma function

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

Please note: I asked a quite similar question (thread from 07-August/ 
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=""))