Basic Latin hypercube samples and designs with package lhs

knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) require(lhs) source("VignetteCommonCode.R") graph2dLHS <- function(Alhs) { stopifnot(ncol(Alhs) == 2) sims <- nrow(Alhs) par(mar = c(4,4,2,2)) plot.default(Alhs[,1], Alhs[,2], type = "n", ylim = c(0,1), xlim = c(0,1), xlab = "Parameter 1", ylab = "Parameter 2", xaxs = "i", yaxs = "i", main = "") for (i in 1:nrow(Alhs)) { rect(floor(Alhs[i,1]*sims)/sims, floor(Alhs[i,2]*sims)/sims, ceiling(Alhs[i,1]*sims)/sims, ceiling(Alhs[i,2]*sims)/sims, col = "grey") } points(Alhs[,1], Alhs[,2], pch = 19, col = "red") abline(v = (0:sims)/sims, h = (0:sims)/sims) } # transform is a function of the kind that takes a number # transform <- function(x){r} graph2dLHSTransform <- function(Alhs, transform1, transform2, min1, max1, min2, max2) { stopifnot(ncol(Alhs) == 2) stopifnot(all(Alhs[,1] <= max1 && Alhs[,1] >= min1)) stopifnot(all(Alhs[,2] <= max2 && Alhs[,2] >= min2)) sims <- nrow(Alhs) breaks <- seq(0,1,length = sims + 1)[2:(sims)] breaksTransformed1 <- sapply(breaks, transform1) breaksTransformed2 <- sapply(breaks, transform2) par(mar = c(4,4,2,2)) plot.default(Alhs[,1], Alhs[,2], type = "n", ylim = c(min2, max2), xlim = c(min1, max1), xlab = "Parameter 1", ylab = "Parameter 2", xaxs = "i", yaxs = "i", main = "") for (si in 1:sims) { temp <- Alhs[si,] for (i in 1:sims) { if ((i == 1 && min1 <= temp[1] && breaksTransformed1[i] >= temp[1]) || (i == sims && max1 >= temp[1] && breaksTransformed1[i - 1] <= temp[1]) || (breaksTransformed1[i - 1] <= temp[1] && breaksTransformed1[i] >= temp[1])) { for (j in 1:sims) { if ((j == 1 && min2 <= temp[2] && breaksTransformed2[j] >= temp[2]) || (j == sims && max2 >= temp[2] && breaksTransformed2[j - 1] <= temp[2]) || (breaksTransformed2[j - 1] <= temp[2] && breaksTransformed2[j] >= temp[2])) { if (i == 1) { xbot <- min1 xtop <- breaksTransformed1[i] } else if (i == sims) { xbot <- breaksTransformed1[i - 1] xtop <- max1 } else { xbot <- breaksTransformed1[i - 1] xtop <- breaksTransformed1[i] } if (j == 1) { ybot <- min2 ytop <- breaksTransformed2[j] } else if (j == sims) { ybot <- breaksTransformed2[j - 1] ytop <- max2 } else { ybot <- breaksTransformed2[j - 1] ytop <- breaksTransformed2[j] } rect(xbot, ybot, xtop, ytop, col = "grey") } } } } } points(Alhs[,1], Alhs[,2], pch = 19, col = "red") abline(v = breaksTransformed1, h = breaksTransformed2) } #set.seed(1111) #A <- randomLHS(5,4) #f <- function(x){qnorm(x)} #g <- function(x){qlnorm(x, meanlog=0.5, sdlog=1)} #B <- A #B[,1] <- f(A[,1]) #B[,2] <- g(A[,2]) #graph2dLHSTransform(B[,1:2], f, g, -4, 4, 0, 8) #f <- function(x){qunif(x, 3, 5)} #B <- apply(A, 2, f) #graph2dLHSTransform(B[,1:2], f)

Theory of Latin Hypercube Sampling

For the technical basis of Latin Hypercube Sampling (LHS) and Latin Hypercube Designs (LHD) please see:

  • Stein, Michael. Large Sample Properties of Simulations Using Latin Hypercube Sampling Technometrics, Vol 28, No 2, 1987.
  • McKay, MD, et.al. A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code Technometrics, Vol 21, No 2, 1979.

This package was created to bring these designs to R and to implement many of the articles that followed on optimized sampling methods.

Create a Simple LHS

Basic LHS's are created using randomLHS.

# set the seed for reproducibility set.seed(1111) # a design with 5 samples from 4 parameters A <- randomLHS(5, 4) A

In general, the LHS is uniform on the margins until transformed (r registerFigure("X")):

r addFigureCaption("X", "Two dimensions of a Uniform random LHS with 5 samples", register=FALSE)

graph2dLHS(A[,1:2])

It is common to transform the margins of the design (the columns) into other distributions (r registerFigure("Y"))

B <- matrix(nrow = nrow(A), ncol = ncol(A)) B[,1] <- qnorm(A[,1], mean = 0, sd = 1) B[,2] <- qlnorm(A[,2], meanlog = 0.5, sdlog = 1) B[,3] <- A[,3] B[,4] <- qunif(A[,4], min = 7, max = 10) B

r addFigureCaption("Y", "Two dimensions of a transformed random LHS with 5 samples", register=FALSE)

f <- function(x){qnorm(x)} g <- function(x){qlnorm(x, meanlog = 0.5, sdlog = 1)} graph2dLHSTransform(B[,1:2], f, g, -4, 4, 0, 8)

Optimizing the Design

The LHS can be optimized using a number of methods in the lhs package. Each method attempts to improve on the random design by ensuring that the selected points are as uncorrelated and space filling as possible. r registerTable("tab1") shows some results. r registerFigure("Z"), r registerFigure("W"), and r registerFigure("G") show corresponding plots.

set.seed(101) A <- randomLHS(30, 10) A1 <- optimumLHS(30, 10, maxSweeps = 4, eps = 0.01) A2 <- maximinLHS(30, 10, dup = 5) A3 <- improvedLHS(30, 10, dup = 5) A4 <- geneticLHS(30, 10, pop = 1000, gen = 8, pMut = 0.1, criterium = "S") A5 <- geneticLHS(30, 10, pop = 1000, gen = 8, pMut = 0.1, criterium = "Maximin")

r addTableCaption("tab1", "Sample results and metrics of various LHS algorithms", register=FALSE)

Method Min Distance btwn pts Mean Distance btwn pts Max Correlation btwn pts
randomLHS r min(dist(A)) r mean(dist(A)) r max(abs(cor(A)-diag(10)))
optimumLHS r min(dist(A1)) r mean(dist(A1)) r max(abs(cor(A1)-diag(10)))
maximinLHS r min(dist(A2)) r mean(dist(A2)) r max(abs(cor(A2)-diag(10)))
improvedLHS r min(dist(A3)) r mean(dist(A3)) r max(abs(cor(A3)-diag(10)))
geneticLHS (S) r min(dist(A4)) r mean(dist(A4)) r max(abs(cor(A4)-diag(10)))
geneticLHS (Maximin) r min(dist(A5)) r mean(dist(A5)) r max(abs(cor(A5)-diag(10)))

r addFigureCaption("Z", "Pairwise margins of a randomLHS", register=FALSE)

pairs(A, pch = 19, col = "blue", cex = 0.5)

r addFigureCaption("W", "Pairwise margins of a optimumLHS", register=FALSE)

pairs(A1, pch = 19, col = "blue", cex = 0.5)

r addFigureCaption("G", "Pairwise margins of a maximinLHS", register=FALSE)

pairs(A2, pch = 19, col = "blue", cex = 0.5)