# Function to generate random observations from a # hyperbolic distribution using the # transformed density rejection method # Follows Hoermann, ACM ToMS (1995),21,182--193. # Uses a squeeze to help with test rhypTDRsq <- function(n,Theta){ hyperbPi <- Theta[1] zeta <- Theta[2] delta <- Theta[3] mu <- Theta[4] if(zeta<=0) stop("zeta must be positive") if(delta<=0) stop("delta must be positive") # restrict to delta = 1, mu = 0 lowerupper <- calculate.range(c(hyperbPi,zeta,1,0)) lower <- lowerupper[1] upper <- lowerupper[2] alpha <- zeta * sqrt(1 + hyperbPi^2) hyperbBeta <- zeta * hyperbPi quasiDens <- function(x,alpha,hyperbBeta){ quasiDens <- exp(-alpha*sqrt(1+x^2)+hyperbBeta*x) } h <- function(x,alpha,hyperbBeta){ h <- -alpha*sqrt(1+x^2)+hyperbBeta*x } dh <- function(x,alpha,hyperbBeta){ h <- -alpha*x/sqrt(1+x^2)+hyperbBeta } xl <- uniroot(function(x) -alpha*sqrt(1+x^2)+hyperbBeta*x+zeta+1, lower=lower,upper=hyperbPi)\$root xr <- uniroot(function(x) -alpha*sqrt(1+x^2)+hyperbBeta*x+zeta+1, lower=hyperbPi,upper=upper)\$root hxl <- h(xl,alpha,hyperbBeta) hxr <- h(xr,alpha,hyperbBeta) dhxl <- dh(xl,alpha,hyperbBeta) dhxr <- dh(xr,alpha,hyperbBeta) bl <- xl+(-zeta-hxl)/dhxl br <- xr+(-zeta-hxr)/dhxr Fhm <- exp(-zeta) fm <- Fhm vl <- Fhm/dhxl vc <- fm*(br-bl) vr <- -Fhm/dhxr # define squeeze function sl <- (-zeta-hxl)/(hyperbPi-xl) sr <- (-zeta-hxr)/(hyperbPi-xr) 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 here X <- (log(((U-(vl+vc)))*dhxr+Fhm)-hxr)/dhxr+xr lx <- exp(dhxr*(X-xr)+hxr) } } V <- lx*runif(1) fX <- quasiDens(X,alpha,hyperbBeta) if(Xxl)&(V<=exp(-zeta-(hyperbPi-X)*sl))){ need.value <- FALSE }else{ if(V<=fX){ need.value <- FALSE } } }else{ #X>=hyperbPi if((X