# Basic Latin hypercube samples and designs with package lhs" In lhs: Latin Hypercube Samples

```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){return(qnorm(x,mean=0, std=1))}
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 && breaksTransformed1[i] >= temp) ||
(i == sims && max1 >= temp && breaksTransformed1[i - 1] <= temp) ||
(breaksTransformed1[i - 1] <= temp && breaksTransformed1[i] >= temp))
{
for (j in 1:sims)
{
if ((j == 1 && min2 <= temp && breaksTransformed2[j] >= temp) ||
(j == sims && max2 >= temp && breaksTransformed2[j - 1] <= temp) ||
(breaksTransformed2[j - 1] <= temp && breaksTransformed2[j] >= temp))
{
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)
```

## Try the lhs package in your browser

Any scripts or data that you put into this service are public.

lhs documentation built on Dec. 28, 2022, 2:59 a.m.