R/rtmvnorm.R

In Laurae2/LauraeCE: Laurae's Cross-Entropy Optimization

### Converted from MATLAB by Tim J. Benham, 2014

## function [X, rho, nar, ngibbs] = rmvnrnd(mu,sigma,N,A,b,rhoThr)
## %RMVNRND Draw from the  truncated multivariate normal distribution.
## %   X = rmvnrnd(MU,SIG,N,A,B) returns in N-by-P matrix X a
## %   random sample drawn from the P-dimensional multivariate normal
## %   distribution with mean MU and covariance SIG truncated to a
## %   region bounded by the hyperplanes defined by the inequalities Ax<=B.
## %
## %   [X,RHO,NAR,NGIBBS]  = rmvnrnd(MU,SIG,N,A,B) returns the
## %   acceptance rate RHO of the accept-reject portion of the algorithm
## %   (see below), the number NAR of returned samples generated by
## %   the accept-reject algorithm, and the number NGIBBS returned by
## %   the Gibbs sampler portion of the algorithm.
## %
## %   rmvnrnd(MU,SIG,N,A,B,RHOTHR) sets the minimum acceptable
## %   acceptance rate for the accept-reject portion of the algorithm
## %   to RHOTHR. The default is the empirically identified value
## %   2.9e-4.
## %
## %   ALGORITHM
## %   The function begins by drawing samples from the untruncated MVN
## %   distribution and rejecting those which fail to satisfy the
## %   constraints. If, after a number of iterations with escalating
## %   sample sizes, the acceptance rate is less than RHOTHR it
## %   switches to a Gibbs sampler.
## %
## % ACKNOWLEDGEMENT
## %   This makes use of TruncatedGaussian by Bruno Luong (File ID:
## %   #23832) to generate draws from the one dimensional truncated normal.
## %
## %   REFERENCES
## %   Robert, C.P, "Simulation of truncated normal variables",
## %   Statistics and Computing, pp. 121-125 (1995).

## % Copyright 2011 Tim J. Benham, School of Mathematics and Physics,
## %                University of Queensland.

rtmvnorm <- function(N, mu, sigma, A, b, ..., rhoThr = NULL, maxSample = NULL) {
  
  ## Constant parameters
  defaultRhoThr <- 1e-4 # min. acceptance rate to apply accept-reject sampling.
  defaultMaxSample <- 1e6 # largest sample to draw
  
  ##
  ## Process input arguments.
  ##
  if (is.null(rhoThr) || rhoThr < 0) {
    rhoThr <- defaultRhoThr
  }
  if (is.null(maxSample) || maxSample < 0) {
    maxSample = defaultMaxSample
  }
  mu <- t(mu)
  p <- length(mu) #dimensions
  if (p < 1) {
    stop("Problem dimension must be at least 1")
  }
  if (is.null(A) || is.na(A) || !is.matrix(A) || dim(A)[2] == 0) {
    A <- matrix(rep(0, p), nrow = 1)
    b <- c(0)
  }
  A <- t(A)
  b <- t(b)
  m <- dim(A)[2] #no. constraints
  if (length(b) != m) {
    stop("A and b not conformable")
  }
  
  ##
  ## initialize return arguments
  ##
  X <- matrix(nrow = N, ncol = p)
  nar <- 0
  ngibbs <- 0
  rho <- 1
  
  ###
  ### Approach 1 : Accept/Reject
  ###
  if (rhoThr < 1) {
    
    ## Try accept-reject approach.
    n <- 0 # no. accepted
    trials <- 0
    passes <- 0
    s <- N
    
    while (n < N && (rho > rhoThr || s < maxSample)) {
      
      # cat('n:', n, a, 'rho:', rho, 's:', s, '\n')
      R <- mvrnorm(s, mu, sigma)
      YY <- R %*% A <= matrix(rep(b, s), nrow = s, byrow = TRUE)
      YY <- matrix(as.numeric(YY), nrow = s)
      R <- R[rowSums(YY) == m, , drop = FALSE]
      nr <- dim(R)[1] #no. valid proposals
      
      if (nr > 0) {
        X[(n + 1):min(N, n + nr), ] <- R[1:min(N - n, nr), ]
        nar <- nar + min(N, n + nr) - n
      }
      
      n <- n + nr
      trials <- trials + s
      rho <- n / trials
      
      if (rho > 0) {
        s <- min(maxSample, ceiling((N - n) / rho), 10 * s)
      } else {
        s <- min(maxSample, 10 * s)
      }
      
      passes <- passes + 1
      
    }
    
  }
  
  ###
  ### Approach 2: Gibbs sampler of Robert, 1995.
  ###
  if (nar < N) {
    ## % choose starting point
    if (nar > 0) {
      x <- X[nar, ]
    } else {
      x <- mu
    }
    
    ## set up inverse Sigma
    SigmaInv <- ginv(sigma)
    n <- nar
    
    while (n < N) {
      
      ## choose p new components
      for (i in 1:p) {
        
        ## Sigmai_i is the (p-1) vector derived from the i-th
        ## column of Sigma by removing the i-th row term.
        Sigmai_i <- sigma[-i, i]
        
        ## Sigma_i_iInv is the inverse of the (p-1)x(p-1)
        ## matrix derived from Sigma = (sigma(ij) ) by
        ## eliminating its i-th row and its i-th column
        
        Sigma_i_iInv <- SigmaInv[-i,-i, drop = FALSE] - SigmaInv[-i, i, drop = FALSE] %*% t(SigmaInv[-i, i, drop = FALSE]) / SigmaInv[i, i, drop = FALSE]
        
        ## x_i is the (p-1) vector of components not being updated
        ## at this iteration. /// mu_i
        x_i <- x[-i, drop = FALSE]
        mu_i <- mu[-i, drop = FALSE]
        ## mui is E(xi|x_i)
        mui <- mu[i] + t(Sigmai_i) %*% Sigma_i_iInv %*% as.matrix(x_i - mu_i, ncol = 1)
        s2i <- sigma[i, i] - t(Sigmai_i) %*% Sigma_i_iInv %*% Sigmai_i
        
        ## Find points where the line with the (p-1) components x_i
        ## fixed intersects the bounding polytope.
        ## A_i is the (p-1) x m matrix derived from A by removing
        ## the i-th row.
        A_i <- A[-i, , drop = FALSE]
        ## Ai is the i-th row of A
        Ai <- A[i, , drop = FALSE]
        c <- (b-x_i %*% A_i) / Ai
        lb <- max(c[Ai < 0])
        if (is.null(lb) || is.na(lb) || length(lb) == 0) {
          lb <- -Inf
        }
        ub <- min(c[Ai > 0])
        if (length(ub) == 0) {
          ub <- Inf
        }
        
        ## % now draw from the 1-d normal truncated to [lb, ub]
        x[i] <- rtnorm(1, mean = mui, sd = sqrt(s2i), lower = lb, upper = ub)
        
      }
      
      n <- n + 1
      X[n, ] <- x
      ngibbs <- ngibbs + 1
      
    }
  }
  
  return(list(X = X, rho = rho, nar = nar, ngibbs = ngibbs))
  
}