# Function for the density with delta=1 and mu=0
g <- function(x,Theta){
  hyperbPi <- Theta[1]
  zeta <- Theta[2]
  exp(-zeta*(sqrt(1+hyperbPi^2)*sqrt(1+x^2)-hyperbPi*x))
}

# Function to generate random observations from a
# hyperbolic distribution using a simple "slice"
# sampler.
rhypSS <- function(n,Theta,burn=0,nth=1){
  hyperbPi <- Theta[1]
  zeta <- Theta[2]
  delta <- Theta[3]
  mu <- Theta[4]

  # The total number of observations actually generated
  N <- burn + nth*(n-1) + 1

  z <- numeric(N)
  y <- numeric(N)

  if(burn==0){
    y[1] <- runif(1,0,g(rghyp(1,
       as.numeric(hyperbChangePars(1,2,c(hyperbPi,zeta,1,0)))),
       Theta))
  } else{
    y[1] <- runif(1,0,g(hyperbPi,Theta))
  }

  for(i in 1:N){
    a <- 1
    b <- 2*hyperbPi/zeta*log(y[i])
    c <- 1 + hyperbPi^2 - (1/zeta)^2*(log(y[i]))^2
    z[i] <- runif(1,(-b-sqrt(b^2-4*a*c))/(2*a),
                  (-b+sqrt(b^2-4*a*c))/(2*a))
    if(i<N) y[i+1] <- runif(1,0,g(z[i],Theta))
  }

  x <- z[(burn+1):N][(0:(n-1))*nth+1]*delta + mu
  x
}