R/sample.R

In Tinflex: A Universal Non-Uniform Random Number Generator

#############################################################################
##
##  Sampling routine for generator.
##  (R version)
##
#############################################################################


## --------------------------------------------------------------------------

Tinflex.sample <- function(gen, n=1) {
    ## ------------------------------------------------------------------------
    ## Draw a sample of size 'n'.
    ## (C version)
    ## ------------------------------------------------------------------------
    ##   gen ... S3 object generated by function 'Tinflex.setup'
    ##   n   ... sample size
    ## ------------------------------------------------------------------------
    ## Return: random sample.
    ## ------------------------------------------------------------------------
    
    switch(class(gen),
           "Tinflex"  = .Call(C_Tinflex_RC_sample, gen, n),
           "TinflexC" = .Call(C_Tinflex_C_sample, gen$Cgen, n),
           stop("Argument 'gen' is invalid")
           )
}  ## -- Tinflex.sample() -- ##

## --------------------------------------------------------------------------

Tinflex.sample.R <- function(gen, n=1) {
    ## ------------------------------------------------------------------------
    ## Draw a sample of size 'n'.
    ## (R version)
    ## ------------------------------------------------------------------------
    ##   gen ... S3 object generated by function 'Tinflex.setup'
    ##   n   ... sample size
    ## ------------------------------------------------------------------------
    ## Return: random sample.
    ## ------------------------------------------------------------------------
    
    ## Get parameters for hat and squeeze.
    ivs <- gen$ivs
    n.ivs <- ncol(ivs)-1
    areas <- gen$ivs["A.ht", 1:n.ivs]
    lpdf <- gen$lpdf
    
    ## Allocate memory for storing random sample.
    rv <- numeric(n)
    
    ## Create sample of size 'n'.
    for (k in 1:n) {
        
        while (TRUE) {
            ## Acceptance-rejection loop.
            
            ## Draw interval.
            ## We simply could use 
            ##   i <- sample.int(n.ivs, size=1, prob=areas)
            ## For compatibility with C code, however, we use inversion method.
            
            ## Sequential search.
            ##      U <- gen$A.ht.tot * runif(1)
            ##      sum <- 0
            ##      for (j in 1:n.ivs) {
            ##        i <- j
            ##        sum <- sum + ivs["A.ht",i]
            ##        if (sum >= U) break
            ##      }
            
            ## Use guide table.
            U <- runif(1)
            j = 1+gen$gt[1+floor(U*n.ivs)];
            U <- U * gen$A.ht.tot;
            for (i in j:n.ivs) {
                if (gen$Acum[i] >= U) break;
            }
            
            ## Get parameters for hat in interval.
            a <- as.numeric(ivs["ht.a",i])
            b <- as.numeric(ivs["ht.b",i])
            y <- as.numeric(ivs["ht.y",i])
            x0 <- as.numeric(ivs["x",i])
            cT <- as.numeric(ivs["c",i])
            
            ## Need a uniform random number in interval (0, ivs["A.ht",i]).
            ##    U <- ivs["A.ht",i] * runif(1)
            ## However, we can recycle U for this task.
            U <- ivs["A.ht",i] + U - gen$Acum[i]
            
            ## Generate from "hat distribution":
            ##    X = y + ( FTinv(cT, FT(cT,t) + b*U) - a ) / b;
            
            ## For numerical reasons we have to distinguish
            ## between different values of 'cT'.
            
            if (cT == 0) {
                ## Case: T(x)=log(x)
                
                z <- U * b / exp(a+b*(x0-y))
                if (isTRUE(abs(z) > 1.e-6)) {
                    X <- y + (log(exp(a + b*(x0-y))+b*U)-a)/b
                    ## or for |z| < Inf:
                    ##  X <- x0 + U / exp(a+b*(x0-y)) * 1/z * log(1+z)
                } else {
                    ## We need approximation by Taylor polynomial to avoid
                    ## severe round-off errors.
                    X <- x0 + U / exp(a+b*(x0-y)) * (1 - z/2 + z*z/3)
                }
            }
            
            else if (cT == -0.5) {
                ## Case: T(x) = -1/sqrt(x)
                
                z <- U * b * (a+b*(x0-y))
                if (isTRUE(abs(z) > 1.e-6)) {
                    X <- y + (FTinv(cT, FT(cT, a + b*(x0-y))+b*U)-a)/b
                } else {
                    X <- x0 + U * (a+b*(x0-y))^2 * (1 + z + z*z)
                }
            }
            
            else if (cT == 1) {
                ## Case: T(x) = x
                
                z <- U * b / (a+b*(x0-y))^2
                if (isTRUE(abs(z) > 1.e-6)) {
                    X <- y + (FTinv(cT, FT(cT, a + b*(x0-y))+b*U)-a)/b
                } else {
                    X <- x0 + U / (a+b*(x0-y)) * (1 - z/2 + z*z/2)
                }
            }
            
            else {
                ## Case: T(x)=sign(c)*x^c
                ## For all other cases we only use a rough approximation in
                ## case of numerical errors.
                if (isTRUE(abs(b)>1e-10)) {
                    X <- y + (FTinv(cT, FT(cT, a + b*(x0-y))+b*U)-a)/b
                } else {
                    U <- U / ivs["A.ht",i]
                    X <- (1-U)*ivs["x",i] + U*ivs["x",i+1]
                }
            }
            
            ## Compute hat and squeeze at X.
            hx <- Tinv(cT, a + b * (X-y))
            
            if (ivs["A.sq",i] > 0)
                sx <- Tinv(cT, ivs["sq.a",i]+ivs["sq.b",i]*(X-ivs["sq.y",i]))
            else
                sx <- 0
            
            ## Accept or reject.
            V <- hx * runif(1)
            
            if (V <= sx) break
            if (V <= exp(lpdf(X))) break
        }
        
        ## Store random point.
        rv[k] <- X
    }
    
    ## Return random sample.
    return (rv)

}  ## -- Tinflex.sample.R() -- ##

## --------------------------------------------------------------------------