####################### file: ars_cons_norm.c.r ############################## system ("R CMD SHLIB ars_cons_norm.c") dyn.load ("ars_cons_norm.so") ## a wrapper function to draw truncated normal distribution ## n is sample size, lb and up are lower and upper bounds ## dars = 1 makes the program to show ARS information ## inip is the first point sample_tnorm_ars <- function (n, lb = -Inf, ub = Inf, inip = NULL, dars = 0) { if ( is.finite (lb) & is.finite (ub)) inip <- (lb + ub) / 2 if ( is.finite (lb) & !is.finite (ub)) inip <- lb + 1 if (!is.finite (lb) & is.finite (ub)) inip <- ub - 1 if (!is.finite (lb) & !is.finite (ub)) inip <- 0 .C("R_sample_tnorm", n:n, rep (0,n), lb, ub, inip, dars, NAOK = TRUE )[[2]] } ## a direct rejection sampling for truncated normal sample_tnorm_drs <- function (n, lb = -Inf, ub = Inf) { x <- rep (0, n) for (i in 1:n) { rej <- TRUE while (rej) { x[i] <- rnorm (1) if (x[i] >= lb & x[i] <= ub) rej <- FALSE } } x } ## test the sampling procedure postscript ("ars_cons_norm-plots.ps", paper = "special", width = 10, height = 8, horiz = F) par (mfrow = c(2,2), mar = c(4,4,3,1)) ## generate random numbers with two methods n <- 5000 system.time (x1 <- sample_tnorm_ars (n, -1, Inf)) system.time (x2 <- sample_tnorm_drs (n, -1, Inf)) hist (x1, main = "ARS Samples of N(0,1)I(x>-1)") qqplot (x1,x2, xlab = "Quantiles of ARS Samples", ylab = "Quantiles of DRS Samples", main = "QQ-plot of Samples of N(0,1)I(x>-1)") ## sample from the tail of normal n <- 5000 x3 <- sample_tnorm_ars (n, 100, Inf) hist (x3, main = "Histogram of ARS Samples of N(0,1)I(x>100)") ## sample from normal constraint in a very small interval n <- 5000 x4 <- sample_tnorm_ars (n, -11, -10) hist (x4, main = "Histogram of ARS Samples of N(0,1)I(-11