dggamma <- function(beta, sigma, lambda, timet)
{

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

# consider timet below exp(-500) at exp(-500) to avoid division by zero
timet[timet<exp(-500)] <- exp(-500)

b <- lambda/sigma
dens <- (lambda>0)*((b*(timet^(b-1))*dgamma(timet^b, shape=s, rate=r)))
dens <- dens+(lambda<0)*((-b*(timet^(b-1))*dgamma(timet^b, shape=s, rate=r)))
dens
}