# 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
}