# Modified version of rgig to generate random observations # from a generalized inverse Gaussian distribution in the # special case where lambda = 1. rgig1 <- 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") 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){ needValue <- TRUE while(needValue == 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){ needValue <- FALSE } } } output[i] <- Y } output/alpha } ## End of rgig1