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