Generate a Monte Carlo sample with Sobol' LDQR sequence

monte-carlo-sample.r

# Generate a Monte Carlo sample using Sobol' low-discrepancy quasi-random sequences
# James Keirstead
# 3 February 2012
#
# Random sampling with R's standard methods is inefficient for Monte Carlo analysis as
# the sampled values do not cover the parameter space evenly.  This Gist allows users
# to create parameter samples using Sobol' sequences to get around this problem.

# makeMCSample
#    Makes a Monte Carlo sample using Sobol' sequences
#
# Parameters:
#     n = number of samples to draw
#  vals = list describing distributions, each consisting of:
#         - a variable name 'var', 
#         - the distribution name, 'dist', e.g. 'unif' for Uniform (see R's distribution names)
#         - the distribution parameters, 'params' (names vary by distribution)
#
# Returns:
#   a data frame with an index column "n" and each distribution in a named column
#
makeMCSample <- function(n, vals) {
  # Packages to generate quasi-random sequences
  # and rearrange the data
  require(randtoolbox)
  require(plyr)

  # Generate a Sobol' sequence
  sob <- sobol(n, length(vals))

  # Fill a matrix with the values
  # inverted from uniform values to
  # distributions of choice
  samp <- matrix(rep(0,n*(length(vals)+1)), nrow=n)
  samp[,1] <- 1:n
  for (i in 1:length(vals)) {
    l <- vals[[i]]
    dist <- l$dist
    params <- l$params
    fname <- paste("q",dist,sep="")
    samp[,i+1] <- do.call(fname,c(list(p=sob[,i]),params))
  }

  # Convert matrix to data frame and add labels
  samp <- as.data.frame(samp)
  names(samp) <- c("n",laply(vals, function(l) l$var))
  return(samp)
}


# An example 
n <- 1000  # number of simulations to run

# List described the distribution of each variable
library(triangle)   # for triangular distributed values
vals <- list(list(var="Uniform",
               dist="unif",
               params=list(min=0,max=1)),
          list(var="Normal",
               dist="norm",
               params=list(mean=0,sd=1)),
          list(var="Weibull",
               dist="weibull",
               params=list(shape=2,scale=1)),
          list(var="Triangular",
               dist="triangle",
               params=list(a=0,b=2,c=0.5)))

# Generate the sample
samp <- makeMCSample(n,vals)