# Modified version of rgigJD to generate random observations # from a generalized inverse Gaussian distribution in the # special case where lambda = 1. rgigJD1 <- function(n,Theta){ if(length(Theta)==2) Theta <- c(1,Theta) lambda <- 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(chi==0) stop("chi = 0, use rgamma") if(psi==0) stop("When lambda >= 0, psi must be > 0") alpha <- sqrt(psi/chi) beta <- sqrt(psi*chi) m <- abs(beta)/beta g <- function(y){ 0.5*beta*y^3 - y^2*(0.5*beta*m+lambda+1) + y*(-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)*exp(-0.25*beta*(yP+1/yP-m-1/m)) b <- (yM-m)*exp(-0.25*beta*(yM+1/yM-m-1/m)) c <- -0.25*beta*(m+1/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.25*beta*(Y+1/Y)+c){ need.value <- FALSE } } } output[i] <- Y } output/alpha }