# [R] integration problem with gamma function

Prof Brian Ripley ripley at stats.ox.ac.uk
Fri Sep 1 21:07:10 CEST 2006

```It is the same issue.  Martin's reply was for a specific value of 'a': the
factor is gamma(a).

Please do study the help page for pgamma (as the posting guide did ask you
to), as it has all the details.

Note that the Mathematica version can easily lead to representation
difficulties for even modest 'a'.

On Fri, 1 Sep 2006, Leo Gürtler wrote:

> 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=""))
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help