#================================================================================# # Goal: Sample small-shape gamma random variables via accept-reject # # Paper: arXiv:1302.1884 # # Function: rgamss # # Authors: C. Liu, R. Martin (www.math.uic.edu/~rgmartin), and N. Syring # # Version: 2nd (04/07/2015) # #================================================================================# # INPUT # n = sample size # shape = shape parameter (same for all samples) # scale = scale parameter (same for all samples; default = 1) # do.log = logical -- return values on log scale # OUTPUT # vector of length n containing gamma samples; # on log-scale depending on do.log or magnitude of shape rgamss <- function(n, shape, scale=1, do.log=TRUE) { a <- shape if(a > 0.2) oo <- rgamma(n, a, 1 / scale, log=do.log) else { e1 <- 2.71828182845905 L <- 1 / a - 1 w <- a / e1 / (1 - a) ww <- 1 / (1 + w) eta <- function(z) if(z >= 0) exp(-z) else w * L * exp(L * z) h <- function(z) exp(-z - exp(-z / a)) rh <- function(a) { repeat { U <- runif(1) z <- if(U <= ww) -log(U / ww) else log(runif(1)) / L if(h(z) / eta(z) > runif(1)) return(z) } } Z <- numeric(n) for(i in 1:n) Z[i] <- rh(a) o <- log(scale) - Z / a if(!do.log) { oo <- exp(o) if(any(oo == 0)) { oo <- o warning("Output given on log-scale since shape is small") } } else oo <- o } return(oo) } # Example: Shows inaccuracy of 'rgamma' for small shape parameter # shape <- 0.001 # log(rgamma(10, shape)) #see the "-Inf" results...? # rgamss(10, shape)