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
>
> 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=""))
>
> ______________________________________________
> R-help at 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.
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595