# Function to generate random observations from a # generalized inverse Gaussian distribution using the # transformed density rejection method in special case # lambda=1. # Follows Hoermann, ACM ToMS (1995),21,182--193. rgigTDR1 <- function(n,Theta){ if(length(Theta)==2) Theta <- c(1,Theta) lambda <- Theta[1] chi <- Theta[2] psi <- Theta[3] if(chi<=0) stop("chi must be positive") if(psi<=0) stop("psi must be positive") alpha <- sqrt(psi/chi) gigBeta <- sqrt(psi*chi) lower <- 0 upper <- gigCalcRange(c(1,gigBeta,gigBeta),tol=10^(-8)) quasiDens <- function(x,gigBeta){ exp(-gigBeta*((1/x)+x)/2) } h <- function(x,gigBeta){ -gigBeta*((1/x)+x)/2 } dh <- function(x,gigBeta){ gigBeta*((1/x^2)-1)/2 } m <- 1 hm <- h(m,gigBeta) xl <- uniroot(function(x) exp(-gigBeta*((1/x)+x)/2)-exp(hm)/exp(1), lower=lower,upper=1)\$root xr <- uniroot(function(x) exp(-gigBeta*((1/x)+x)/2)-exp(hm)/exp(1), lower=1,upper=upper)\$root hxl <- h(xl,gigBeta) hxr <- h(xr,gigBeta) dhxl <- dh(xl,gigBeta) dhxr <- dh(xr,gigBeta) bl <- xl+(hm-hxl)/dhxl br <- xr+(hm-hxr)/dhxr Fhm <- exp(hm) fm <- Fhm vl <- (Fhm-exp(-dhxl*xl+hxl))/dhxl vc <- fm*(br-bl) vr <- -Fhm/dhxr output <- numeric(n) for(i in 1:n){ need.value <- TRUE while(need.value==TRUE){ U <- (vl+vc+vr)*runif(1) if (U<=vl){ X <- (log(-U*dhxl+Fhm)-hxl)/dhxl+xl lx <- exp(dhxl*(X-xl)+hxl) }else{ if(U<=(vl+vc)){ X <- ((U-vl)/vc)*(br-bl)+bl lx <- fm }else{ # Beware error in Hoermann paper here X <- (log(((U-(vl+vc)))*dhxr+Fhm)-hxr)/dhxr+xr lx <- exp(dhxr*(X-xr)+hxr) } } V <- lx*runif(1) if(X <= 0){ fX <-0 }else{ fX <- quasiDens(X,gigBeta) } if(V<=fX){ need.value <- FALSE } } output[i] <- X } return(output/alpha) }