qggamma <- function(beta, sigma, lambda, p)
{

s <- lambda^(-2)
r <- s * exp(( - beta * lambda)/sigma)

# consider p below 0.001 at 0.001 to avoid numerical problems
p[p<0.001] <- 0.001
# consider p above 0.999 at 0.999 to avoid numerical problems
p[p>0.999] <- 0.999

b <- lambda/sigma
quan <- (lambda>0)*((qgamma(p,shape=s, rate=r))^(1/b))
quan <- quan+(lambda<0)*((qgamma(1-p,shape=s, rate=r))^(1/b))
quan
}