# Function to generate random observations from a # generalized inverse Gaussian distribution. The # algorithm is based on that given by Dagpunar (1989) rgigJD <- function(n,Theta){ lambda <- Theta[1] chi <- Theta[2] psi <- Theta[3] if(chi<0) stop("chi can not be negative") if(psi<0) stop("psi can not be negative") if((lambda>=0)&(psi==0)) stop("When lambda >= 0, psi must be > 0") if((lambda<=0)&(chi==0)) stop("When lambda <= 0, chi must be > 0") if(chi==0) stop("chi = 0, use rgamma") if(psi==0) stop("algorithm only valid for psi > 0") alpha <- sqrt(psi/chi) beta <- sqrt(psi*chi) m <- (lambda-1+sqrt((lambda-1)^2+beta^2))/beta g <- function(y){ 0.5*beta*y^3 - y^2*(0.5*beta*m+lambda+1) + y*((lambda-1)*m-0.5*beta) + 0.5*beta*m } upper <- m while(g(upper)<=0) upper <- 2*upper yM <- uniroot(g,interval=c(0,m))\$root yP <- uniroot(g,interval=c(m,upper))\$root a <- (yP-m)*(yP/m)^(0.5*(lambda-1))* exp(-0.25*beta*(yP+1/yP-m-1/m)) b <- (yM-m)*(yM/m)^(0.5*(lambda-1))* exp(-0.25*beta*(yM+1/yM-m-1/m)) c <- -0.25*beta*(m+1/m) + 0.5*(lambda-1)*log(m) output <- numeric(n) for(i in 1:n){ need.value <- TRUE while(need.value==TRUE){ R1 <- runif(1) R2 <- runif(1) Y <- m + a*R2/R1 + b*(1-R2)/R1 if(Y>0){ if(-log(R1)>=-0.5*(lambda-1)*log(Y)+0.25*beta*(Y+1/Y)+c){ need.value <- FALSE } } } output[i] <- Y } output/alpha }