# ONE-SAMPLE CONFIDENCE INTERVALS
# usage: gamma.one.ci(nr, x, cl, parameter)
#        nr = sumber of simulation runs (recommended value 100,000)
#         x = vector containing the sample data
#        cl = conf level
# parameter = "mean", "scale" or "shape"
# Example:         
# x <- c(5.1,2.4,.4,.5,2.5,.1,6.8,1.2,.5,.6,5.3,2.3,1.8,1.2,1.3,1.1,
#        .9,3.2,1.,.9,.4,.6,8,.4,2.7,.2,2,.2,.5,.8,2,2.9,.1,4)
# > gamma.one.ci(100000,x,.90,"mean")
#   lower    upper 
#  1.443292 2.553004 -> 90% CI for the gamma mean
# > gamma.one.ci(100000,x,.90,"scale")
#   lower    upper 
#  1.231482 3.076869 -> 90% CI for the scale parameter
# > gamma.one.ci(100000,x,.95,"shape")
#   lower    upper 
#  0.649768 1.52902 -> 95% CI for the shape parameter
###############################################################################
# TWO-SAMPLE CONFIDENCE INTERVALS
# usage: gamma.two.ci(nr, x1, x2, cl, parameter)
#        nr = sumber of simulation runs (recommended value 100,000)
#        x1 = vector containing the first sample data
#        x2 = vector containing the second sample data
#        cl = conf level
# parameter = "mean" or "shape" 
# Example:         
# x1 <- c(7.7,1656.0,978.0,198.6,703.4,1697.8,334.1,118.3,255.0,115.3,242.5,
#         32.7,40.6,129.6,31.4,2745.6,489.1,430.0,302.8,119.0,4.1,92.4,17.5,
#         200.7,274.7,274.7)   
# x2 <- c(26.1,26.3,87.0,95.0,1.0,372.4,17.3, 24.4, 11.5,321.2,68.5,81.2,47.3,
#         28.6,830.1,345.5,1202.6,36.6,4.9,4.9,41.1,29.0,163.0,244.3,147.8,21.7)   
# > gamma.two.ci(100000, x1, x2, .95, "mean")
#   lower  upper 
#   93.35  582.1  -> 95% CI for the difference between two gamma means
###############################################################################
# ONE-SAMPLE TESTS
# usage: gamma.one.test(nr, x, parameter,parameter0)
#        nr = sumber of simulation runs (recommended value 100,000)
#         x = vector containing the sample data
# parameter0= null value of the parameter
# parameter = "mean", "scale" or "shape"
# Example:         
# x <- c(5.1,2.4,.4,.5,2.5,.1,6.8,1.2,.5,.6,5.3,2.3,1.8,1.2,1.3,1.1,
#        .9,3.2,1.,.9,.4,.6,8,.4,2.7,.2,2,.2,.5,.8,2,2.9,.1,4)
# > gamma.one.test(0,x,"mean",2.8) # no simulation
#  MLRTstat     left_pvalue     right_pvalue
#  -2.086       .01851097       .981489 -> The mean is less than 2.8
# > gamma.one.test(0,x,"shape",.5) # no simulation
#  SLRTstat     left_pvalue     right_pvalue
#  3.198        .99857          .00143 -> The shape parameter is greater than 2.8
###############################################################################
# TWO-SAMPLE TESTS
# usage: gamma.two.test(nr, x1, x2, parameter)
#        nr = sumber of simulation runs (recommended value 100,000)
#        x1 = vector containing the first sample data
#        x2 = vector containing the second sample data
# parameter = "mean", "shape" or "scale" 
# Example:         
# x1 <- c(7.7,1656.0,978.0,198.6,703.4,1697.8,334.1,118.3,255.0,115.3,242.5,
#         32.7,40.6,129.6,31.4,2745.6,489.1,430.0,302.8,119.0,4.1,92.4,17.5,
#         200.7,274.7,274.7)   
# x2 <- c(26.1,26.3,87.0,95.0,1.0,372.4,17.3, 24.4, 11.5,321.2,68.5,81.2,47.3,
#         28.6,830.1,345.5,1202.6,36.6,4.9,4.9,41.1,29.0,163.0,244.3,147.8,21.7)   
# > gamma.two.test(10000, x1, x2, "mean")
#   SSLRTstat  left_pvalue  right_pvalue  
#   2.523      .994         .006  -> Mean 1 > Mean 2
###############################################################################

#PB CI FOR THE MEAN
gamma.one.ci = function(nr, x, cl, parameter){
if(parameter == "mean"){
pb.ci.mean(nr,x,cl)}
else if(parameter == "shape"){
pb.ci.shape(nr,x,cl)}
else if(parameter == "scale"){
pb.ci.scale(nr, x, cl)
}
else{
return(c("parameter is mean, scale or shape"))
}
}
###############################################################################
# PB CI FOR A GAMMA MEAN
pb.ci.mean = function(nr, x, cl){
pivot = seq(1:nr)
al = (1-cl)/2; ll = floor(nr*al); 
n = length(x)
xb = mean(x); xd = mean(log(x))
ml =  gam.mles(n, xb, xd)
ah = ml[1]; bh = ml[2]
for(i in 1:nr){
x = rgamma(n,ah)*bh
xbs = mean(x); xds = mean(log(x))
ml = gam.mles(n, xbs, xds)
ahs = ml[1]; bhs = ml[2]
pivot[i] = sqrt(n*ahs)*(xbs -xb)/xbs
}
pivot = sort(pivot); tlow = pivot[ll]; tupp = pivot[nr-ll]
el = xb-tupp*xb/sqrt(n*ah)
ul = xb-tlow*xb/sqrt(n*ah)
return(list(lower=el,upper=ul))
}
#PB CI FOR THE SHAPE PARAMETER
pb.ci.shape = function(nr, x, cl){
pivot = seq(1:nr)
al = (1-cl)/2; ll = floor(nr*al)
n = length(x)
xb = mean(x); xd = mean(log(x))
ml = gam.mles(n, xb, xd)
ah = ml[1]; bh = ml[2]
sd0 = sqrt(ah/(n*(ah*trigamma(ah)-1)))
for(i in 1:nr){
x = rgamma(n, ah)*bh
xbs = mean(x); xds = mean(log(x))
ml = gam.mles(n, xbs, xds)
ahs = ml[1]; bhs = ml[2]
sds = sqrt(ahs/(n*(ahs*trigamma(ahs)-1)))
pivot[i] = (ahs -ah)/sds
}
pivot = sort(pivot)
tlow = pivot[ll]; tupp = pivot[nr-ll]
el = ah-tupp*sd0
ul = ah-tlow*sd0
return(list(lower=el,upper=ul))
}
#PB CI FOR THE SCLAE PARAMETER
pb.ci.scale = function(nr, x, cl){
pivot = seq(1:nr)
al = (1-cl)/2; ll = floor(nr*al)
n = length(x)
xb = mean(x); xd = mean(log(x))
ml = gam.mles(n, xb, xd)
ah = ml[1]; bh = ml[2]
sd0 = sqrt(bh**2*trigamma(ah)/(n*(ah*trigamma(ah)-1.0)))
for(i in 1:nr){
x = rgamma(n, ah)*bh
xbs = mean(x); xds = mean(log(x))
ml = gam.mles(n, xbs, xds)
ahs = ml[1]; bhs = ml[2]
sds = sqrt(bhs**2*trigamma(ahs)/(n*(ahs*trigamma(ahs)-1)))
pivot[i] = (bhs -bh)/sds
}
pivot = sort(pivot)
tlow = pivot[ll]; tupp = pivot[nr-ll]
el = bh-tupp*sd0
ul = bh-tlow*sd0
return(list(lower=el,upper=ul))
}
#########################################################################
# TWO-SAMPLE CIS
gamma.two.ci = function(nr, x1, x2, cl, parameter){
if(parameter == "shape"){
pb.ci.dif.shapes(nr, x1, x2, cl)
}
else if(parameter == "mean"){
pb.ci.dif.means(nr, x1, x2, cl)}
else{
return(c("parameter is mean or shape "))
}
}
#CIS FOR THE DIFFERENCE BETWEEN TWO MEANS
pb.ci.dif.means = function(nr, x1, x2, cl){
pivot = seq(1:nr)
al = (1-cl)/2; ll = floor(nr*al)
n1 = length(x1); n2 = length(x2)
xb1 = mean(x1); xd1 = mean(log(x1))
xb2 = mean(x2); xd2 = mean(log(x2))
ml =  gam.mles(n1, xb1, xd1)
ah1 = ml[1]; bh1 = ml[2]
ml =  gam.mles(n2, xb2, xd2)
ah2 = ml[1]; bh2 = ml[2]
sd0 = sqrt(xb1**2/n1/ah1+xb2**2/n2/ah2)
t0 = (xb1-xb2)/sd0
for(i in 1:nr){
x1 = rgamma(n1,ah1)*bh1
xbs1 = mean(x1); xds1 = mean(log(x1))
ml = gam.mles(n1, xbs1, xds1)
ahs1 = ml[1]; bhs1 = ml[2]
x2 = rgamma(n2,ah2)*bh2
xbs2 = mean(x2); xds2 = mean(log(x2))
ml = gam.mles(n2, xbs2, xds2)
ahs2 = ml[1]; bhs2 = ml[2]
sds = sqrt(xbs1**2/n1/ahs1+xbs2**2/n2/ahs2)
pivot[i] = (xbs1-xbs2-(xb1-xb2))/sds
}
pivot = sort(pivot)
tlow = pivot[ll]; tupp = pivot[nr-ll]
el = xb1-xb2-tupp*sd0
ul = xb1-xb2-tlow*sd0
return(list(lower=el,upper=ul))
}
#PB CI FOR THE DIFFERENCE BETWEEN TWO SHAPE PARAMETERS
pb.ci.dif.shapes = function(nr, x1, x2, cl){
pivot = seq(1:nr)
n1 = length(x1); n2 = length(x2)
al = (1-cl)/2; ll = floor(nr*al)
ml = gam.mles.data(x1)
ah1 = ml[1]; bh1 = ml[2]
ml = gam.mles.data(x2)
ah2 = ml[1]; bh2 = ml[2]
v10 = ah1/(n1*(ah1*trigamma(ah1)-1))
v20 = ah2/(n2*(ah2*trigamma(ah2)-1))
sd0 = sqrt(v10+v20)
for(i in 1:nr){
x1 = rgamma(n1,ah1)*bh1
xbs1 = mean(x1); xds1 = mean(log(x1))
ml = gam.mles(n1, xbs1, xds1)
ahs1 = ml[1]; bhs1 = ml[2]
x2 = rgamma(n2,ah2)*bh2
xbs2 = mean(x2); xds2 = mean(log(x2))
ml = gam.mles(n2, xbs2, xds2)
ahs2 = ml[1]; bhs2 = ml[2]
v1s = ahs1/(n1*(ahs1*trigamma(ahs1)-1))
v2s = ahs2/(n2*(ahs2*trigamma(ahs2)-1))
sds = sqrt(v1s+v2s)
pivot[i] = (ahs1-ahs2-(ah1-ah2))/sds
}
pivot = sort(pivot)
tlow = pivot[ll]; tupp = pivot[nr-ll]
el = ah1-ah2-tupp*sd0
ul = ah1-ah2-tlow*sd0
return(list(lower=el,upper=ul))
}
#################################################################
#TESTS: ONE-SAMPLE CASE
gamma.one.test = function(nr, x, parameter, parameter0){
if(parameter == "shape"){
gamma.test.shape(nr, x, parameter0)}
else if(parameter == "scale"){
gamma.test.scale(nr, x, parameter0)}
else if(parameter == "mean"){
gamma.test.mean(x, parameter0)} # no simulation, nr is not used
else{
return(c("parameter is shape, scale or mean"))
}
}
#TESTS FOR THE SHAPE PARAMETER
gamma.test.shape = function(nr,x,a0){
n = length(x)
stat0 = gamma.shape.stat(x,a0)
bhc = mean(x)/a0
pval = 0
for(i in 1:nr){
x <- rgamma(n,a0)*bhc
st = gamma.shape.stat(x,a0)
if(st > stat0){pval = pval+1}
}
rpval = pval/nr; lpval = 1-rpval
return(list(SLRTstat = stat0,left_pvalue=lpval,right_pvalue=rpval))
}
gamma.shape.stat = function(x,a0){
n = length(x)
xb = mean(x); xd = mean(log(x))
ml = gam.mles(n, xb, xd)
ah = ml[1]; bh = ml[2]
sl = lgamma(a0)-lgamma(ah)+(ah-a0)*(xd-log(xb)-1)+(ah*log(ah)-a0*log(a0))
stat = sign(ah-a0)*sqrt(2*n*sl)
return(stat)
}
#TESTS FOR THE SCALE PARAMETER
gamma.test.scale = function(nr,x,b0){
n = length(x)
stat0 = gamma.scale.stat(x,b0)
pval = 0
for(i in 1:nr){
x = rgamma(n,1)
stats = gamma.scale.stat(x,1)
if(stats > stat0){pval = pval +1}
}
rpval = pval/nr; lpval = 1-rpval
return(list(SLRTstat = stat0, left_pvalue=lpval, right_pvalue=rpval))
}
gamma.scale.stat = function(x,b0){
n = length(x)
xb = mean(x); xd = mean(log(x))
ml = gam.mles(n, xb, xd)
ah = ml[1]; bh = ml[2]
sc = xd-log(b0)
ahc = digaminv(sc)
sl = lgamma(ahc)-lgamma(ah)+ah*(xd-log(bh)-1)-ahc*(xd-log(b0))+xb/b0
stat = sign(bh-b0)*sqrt(2*n*sl)
return(stat)
}
#TESTS FOR THE MEAN
gamma.test.mean <- function(x, mu0){
n <- length(x)
xb <- mean(x)
xd <- mean(log(x))
ml <- gam.mles(n,xb,xd)
ah <- ml[1]; bh <- ml[2]; muh <- xb
ml <- cons.mles.mean(n, xb, xd, mu0)
ahc <- ml[1]; bhc <- ml[2]
sl <- lgamma(ahc)-lgamma(ah)+ahc*(log(mu0)-xd-log(ahc)+xb/mu0)-
      ah*(log(xb)-xd-log(ah)+1.0)
qty <- sqrt(n*ah)*(muh/mu0-1)
qty <- qty*sqrt(trigamma(ah)-1/ah)/sqrt(trigamma(ahc)-1/ahc)
r <- sign(muh-mu0)*sqrt(2.0*n*sl)
stat <- r-log(r/qty)/r
lpval = pnorm(stat); rpval = 1-lpval;
return(list(MLRTstat = stat, left_pvalue = lpval, right_pvalue = rpval))
}

cons.mles.mean <- function(n, xb, xd, mu){
s <- xb/mu-xd-1.0+log(mu)  
a0 <- (3.0-s+sqrt((s-3.0)^2+24*s))/12/s
l <- 1
repeat{
ans <- (log(a0)-digamma(a0)-s)
a1 <- a0-ans/(1/a0-trigamma(a0))
if(abs(a1-a0) <= 1.0e-7 | l >= 30){break}
a0 <- a1
l <- l+1}
ahc <- a1
bhc <- xb/ahc
return(c(ahc,bhc))
}
#################################################################
# TWO-SAMPLE TESTS
#################################################################
gamma.two.test = function(nr, x1, x2, parameter){
if(parameter == "shape"){
gamma.test.diff.shapes(x1,x2)}
else if(parameter == "scale"){
gamma.test.diff.scales(nr,x1,x2)}
else if(parameter == "mean"){
gamma.test.diff.means(nr, x1, x2)}
else{
return(c("parameter is mean, scale, shape"))
}
}
# TEST FOR THE DIFFERENCE BTEWEEN TWO SHAPE PARAMETERS
gamma.test.diff.shapes = function(x1,x2){
n1 = length(x1); n2 = length(x2)
xb1 = mean(x1); xd1 = mean(log(x1))
xb2 = mean(x2); xd2 = mean(log(x2))
s1 = log(xb1)-xd1
s2 = log(xb2)-xd2
ml = gam.mles(n1,xb1,xd1)
ah1 = ml[1]
ml = gam.mles(n2,xb2,xd2)
ah2 = ml[1]
df1 = (1+1/(1+4.3*ah1)**2)*(n1-1.0)
df2 = (1+1/(1+4.3*ah2)**2)*(n2-1.0)
stat = (n1*s1/(n1-1))/(n2*s2/(n2-1))
print(c(df1,df2))
rpval = 1-pf(stat, df1, df2)
lpval = 1-rpval
return(list(SBEstat=stat,left_pvalue=lpval,right_pvalue=rpval))
}

#TESTS FOR THE DIFFERENCE BETWEEN TWO SCALE PARAMETERS
gamma.test.diff.scales <- function(nr,x1,x2){
n1 <- length(x1); n2 <- length(x2)
slrt0 <- slrt.two.scales(x1,x2)
pval <- 0
for(i in 1:nr){
y1 <- rgamma(n1,1)
y2 <- rgamma(n2,1)
stats <- slrt.two.scales(y1,y2)
if(stats > slrt0){pval <- pval + 1}
}
rpval <- pval/nr; lpval = 1-rpval;
return(list(SLRTstat = slrt0, left_pvalue = lpval, right_pvalue = rpval))
}
slrt.two.scales <- function(x1,x2){
n1 <- length(x1); n2 <- length(x2)
sn <- n1+n2
xb1 <- mean(x1); xd1 <- mean(log(x1))
xb2 <- mean(x2); xd2 <- mean(log(x2))
xbb = n1*xb1+n2*xb2
h1 = n1*xb1/xbb; h2 = 1- h1
ml <- gam.mles(n1,xb1,xd1)
ah1 = ml[1]; bh1 = ml[2]
ml <- gam.mles(n2,xb2,xd2)
ah2 = ml[1]; bh2 = ml[2]
ml <- mles.cons.shapes.nr(x1,x2)
ah1c = ml[1]; ah2c = ml[2]; bhc = ml[3]
ulf = -n1*lgamma(ah1)-n1*ah1*log(bh1)-n1*xb1/bh1+(ah1-1)*n1*xd1
ulf = ulf-n2*lgamma(ah2)-n2*ah2*log(bh2)-n2*xb2/bh2+(ah2-1)*n2*xd2
clf = -n1*lgamma(ah1c)-n1*ah1c*log(xbb)+n1*ah1c*log(n1*ah1c+n2*ah2c)
clf = clf-h1*(n1*ah1c+n2*ah2c)+(ah1c-1)*n1*xd1
clf = clf-n2*lgamma(ah2c)-n2*ah2c*log(xbb)+n2*ah2c*log(n1*ah1c+n2*ah2c)
clf = clf-h2*(n1*ah1c+n2*ah2c)+(ah2c-1)*n2*xd2
sl = ulf-clf
slrt = sign(bh1-bh2)*sqrt(2*sl)
return(slrt)
}
# this prog calculates constrained mles of the shape parameters when H0: b1 = b2
mles.cons.shapes.nr <- function(x1,x2){
f <- matrix(c(1:4),2,2)
n1 <- length(x1); n2 <- length(x2)
sn <- n1+n2
xb1 <- mean(x1); xd1 <- mean(log(x1))
xb2 <- mean(x2); xd2 <- mean(log(x2))
xbb = log((n1*xb1+n2*xb2))
s1 = xbb-xd1
s2 = xbb-xd2
a10 = (3.0-s1+sqrt((s1-3.0)^2+24.0*s1))/12.0/s1
a20 = (3.0-s2+sqrt((s2-3.0)^2+24.0*s2))/12.0/s2
l = 1
repeat{
anb = n1*a10+n2*a20
f1 = log(anb)-digamma(a10)-s1
f2 = log(anb)-digamma(a20)-s2
f[1,1] = n1/anb-trigamma(a10)
f[1,2] = n2/anb
f[2,1] = n1/anb
f[2,2] = n2/anb-trigamma(a20)
fi = solve(f)
a1n = a10-fi[1,1]*f1-fi[1,2]*f2
a2n = a20-fi[2,1]*f1-fi[2,2]*f2
err = max(abs(f1),abs(f2))
if(err <= 1.0e-4 | l > 300){break}
a10 = a1n; a20 = a2n
l = l +1
}
bh = (n1*xb1+n2*xb2)/(n1*a1n+n2*a2n)
return(c(a1n,a2n, bh))
}
#SSLRT FOR TESTING THE DIFFERENCE BETWEEN TWO GAMMA MEANS
gamma.test.diff.means = function(nr,x1,x2){
stat = seq(1:nr)
n1 = length(x1)
n2 = length(x2) 
xb1 = mean(x1); xd1 = mean(log(x1))
xb2 = mean(x2); xd2 = mean(log(x2))
ml <- gam.mles(n1,xb1,xd1)
ah1 = ml[1]; bh1 = ml[2]
ml <- gam.mles(n2,xb2,xd2)
ah2 = ml[1]; bh2 = ml[2]
stat0 = slrt.dif.means(n1, xb1, xd1, n2, xb2, xd2, ah1, bh1, ah2, bh2)
ml = cons.mles.nr(n1, n2, xb1, xd1, xb2, xd2, ah1, ah2)
uh = ml[1]; ah1c = ml[2]; ah2c = ml[3]
for(i in 1:nr){
x1 = rgamma(n1,ah1c)*uh/ah1c
x2 = rgamma(n2,ah2c)*uh/ah2c
xb1 = mean(x1); xd1 = mean(log(x1))
xb2 = mean(x2); xd2 = mean(log(x2))
ml <- gam.mles(n1,xb1,xd1)
ah1 = ml[1]; bh1 = ml[2]
ml <- gam.mles(n2,xb2,xd2)
ah2 = ml[1]; bh2 = ml[2]
stat[i] = slrt.dif.means(n1, xb1, xd1, n2, xb2, xd2, ah1, bh1, ah2, bh2)
}
xb = mean(stat); s = sd(stat)
slrt0 = (stat0-xb)/s
lpval = pnorm(slrt0); rpval = 1-lpval
return(list(SLRTstat = slrt0, left_pvalue = lpval, right_pvalue = rpval))
}
cons.mles.nr = function(n1, n2, xb1, xd1, xb2, xd2, ah1, ah2){
f = matrix(c(1:4),2,2)
a10 = ah1
a20 = ah2
uh = (n1*xb1+n2*xb2)/(n1+n2)
s1 = log(uh)-xd1+xb1/uh-1
s2 = log(uh)-xd2+xb2/uh-1
a10 = (3.0-s1+sqrt((s1-3.0)^2+24.0*s1))/12.0/s1
a20 = (3.0-s2+sqrt((s2-3.0)^2+24.0*s2))/12.0/s2
l = 1
repeat{
abar = n1*a10+n2*a20
uh = (n1*a10*xb1+n2*a20*xb2)/abar
dig10 = digamma(a10); dig20 = digamma(a20)
tri10 = trigamma(a10); tri20 = trigamma(a20)
f1 = n1*(log(a10)-dig10-log(uh)+xd1-xb1/uh+1)
f2 = n2*(log(a20)-dig20-log(uh)+xd2-xb2/uh+1)
f[1,1] = n1*(1/a10-tri10+n1*(xb1-uh)**2/abar/uh**2)
f[2,2] = n2*(1/a20-tri20+n2*(xb2-uh)**2/abar/uh**2)
f[1,2] = n2*n1*(xb2-uh)*(xb1-uh)/abar/uh**2
f[2,1] = f[1.2]
fi = solve(f)
a11 = a10-fi[1,1]*f1-fi[1,2]*f2
a21 = a20-fi[2,1]*f1-fi[2,2]*f2
err = sqrt((a11-a10)^2+(a21-a20)^2)
if(err <= 1.0e-5 | l > 300 | a10 <= 0 | a20 <= 0){break}
a10 = a11
a20 = a21
l = l + 1
}
return(c(uh,a11,a21))
}
slrt.dif.means = function(n1, xb1, xd1, n2, xb2, xd2, ah1, bh1, ah2, bh2){
ml = cons.mles.nr(n1, n2, xb1, xd1, xb2, xd2, ah1, ah2)
uh = ml[1]; ah1c = ml[2]; ah2c = ml[3]
lu1 = n1*(-lgamma(ah1)-ah1*log(bh1)-xb1/bh1+(ah1-1.0)*xd1)
lu2 = n2*(-lgamma(ah2)-ah2*log(bh2)-xb2/bh2+(ah2-1.0)*xd2)
lu = lu1 + lu2
lc1 = n1*(-lgamma(ah1c)-ah1c*log(uh/ah1c)-xb1*ah1c/uh+(ah1c-1.0)*xd1)
lc2 = n2*(-lgamma(ah2c)-ah2c*log(uh/ah2c)-xb2*ah2c/uh+(ah2c-1.0)*xd2)
lc = lc1 + lc2
slrt = sign(xb1-xb2)*sqrt(2.0*(lu-lc))
if(is.nan(slrt)){slrt = 0}
return(slrt)
}
# CALCULATION OF THE MLES GIVEN AM AND LOG(GM) OF THE SAMPLE
gam.mles = function(n, xb, xd){
s = log(xb)-xd
a0 = (3.0-s+sqrt((s-3.0)^2+24.0*s))/12.0/s
l = 1
repeat{
ans = (log(a0)-digamma(a0)-s)
a1 = a0-ans/(1.0/a0-trigamma(a0))
if(abs(ans) <= 1.0e-7 | l >= 30){break}
a0 = a1
l = l+1}
ah = a1; bh = xb/a1
return(c(ah,bh))
}

# CALCULATION OF THE MLES GIVEN THE SAMPLE
gam.mles.data = function(x){
n = length(x)
xb = mean(x)
xd = mean(log(x))
s = log(xb)-xd
a0 = (3.0-s+sqrt((s-3.0)^2+24.0*s))/12.0/s
l = 1
repeat{
ans = (log(a0)-digamma(a0)-s)
a1 = a0-ans/(1.0/a0-trigamma(a0))
if(abs(ans) <= 1.0e-7 | l >= 30){break}
a0 = a1
l = l+1}
ah = a1; bh = xb/a1
return(c(ah,bh))
}

#AUXIALLARY PROGRAM TO CALCULATE INVERSE DIGAMMA FUNCTION
digaminv = function(c){
fn = function(x,c){
return(digamma(x)-c)
}
x = exp(c)
y = exp(c)+1
fx = fn(x,c)
fy = fn(y,c)
c0 = (x*fy-y*fx)/(fy-fx)
l = 1
repeat{
c1 = (x*fy-y*fx)/(fy-fx)
if(fn(c1,c)*fx < 0){
y = c1
fy = fn(y,c)
if(fn(c0,c)*fn(c1,c) > 0){fx = fx/2}
}
else{
x = c1
fx = fn(x,c)
if(fn(c0,c)*fn(c1,c) > 0){fy = fy/2}
}
c0 = c1
if(l > 100 || abs(fn(c1,c)) <= 1.0e-7){break}
l = l + 1
}
return(c1)
}