Optimizing Sampling Strategy to Enhance Uniformity Under Conditional Constraints

Question

I am facing a challenge in a project that involves sampling from a design space defined by 10 variables. I use Latin Hypercube Sampling (LHS) and/or Sobol sequences, and initially, the samples are uniformly distributed. The complication arises when I introduce a condition based on a derived variable (DV), which is the product of four of these variables.

Here's the setup:

• I sample 10 variables using LHS or Sobol sequences.
• I calculate DV from four of these variables.
• I retain a sample only if DV falls within a specific range, discarding others.

This filtering condition aims to ensure that samples meet specific criteria but appears to bias the distribution of DV, affecting sample quality and uniformity. I am considering integrating a resampling strategy that targets regions where DV is underrepresented. The idea is to resample these specific areas more intensively to improve the uniformity of DV across the sampled data.

Questions:

1. Are there established best practices or modifications to the sampling process that can maintain a uniform distribution while adhering to conditional constraints on DV?
2. Does anyone have experience with or suggestions on using resampling to enhance data uniformity in the presence of such conditional constraints?

Any references or pointers on balancing uniform sampling with conditional filtering would be greatly appreciated. This scenario may resonate with many who manage constrained sampling spaces where maintaining a representative and unbiased sample distribution is crucial.

Thank you for any insights or references you can provide!

Answer 1

This procedure below is not guaranteed to work, but it will find a solution for some problems depending on the marginal distributions (not specified in the question), the tightness of bounds on the product of interest, the number of Latin hypercube samples required, and the maximum number of iterations allowed.

There are general references for optimizing latin hypercubes for a variety of criteria, but none that I know of that are directly on point for this problem 1. I used these algorithms in developing my R package lhs

Process:

• Create a simple random sample using rejection sampling
• Create empirical inverse cumulative distribution functions from the rejection sample
• Draw a Latin hypercube
• Transform the margins of the hypercube using the empirical inverse CDFs
• Iterate through the hypercube swapping one margin of two points until a hypercube is found that meets the criteria

Example below:

set.seed(3803)  
# Assumptions
#
# v1 ~ runif(1, 4)
# v2 ~ runif(1E-6, 2)
# v3 ~ runif(2, 6)
# v4 ~ runif(1E-6, 0.1)
# vx ~ runif(0, 1)
#
# lower < v1*v2*v3*v4 < upper

N <- 100000 # size of simple random sample
upper <- 1
lower <- 0.3
Nlhs <- 50 # size of Latin hypercube
klhs <- 10

# create a simple random sample for rejection sampling
reject_samp <- data.frame(
  v1 = runif(N, 1, 4),
  v2 = runif(N, 1E-6, 2),
  v3 = runif(N, 2, 6),
  v4 = runif(N, 1E-6, 0.1)
)
p <- with(reject_samp, v1*v2*v3*v4)
ind <- which(p < upper & p > lower)
reject_samp <- reject_samp[ind,]

# visualize marginal pdfs
par(mfrow=c(2,2))  
for (i in 1:4) {
  hist(reject_samp[,i], main = paste0("rejection sampled v", i), breaks = 30)
}

# create empirical distribution functions based on the rejection sample
F1 <- ecdf(reject_samp[,1])
Finv1 <- function(p) quantile(reject_samp[,1], probs = p)
F2 <- ecdf(reject_samp[,2])
Finv2 <- function(p) quantile(reject_samp[,2], probs = p)
F3 <- ecdf(reject_samp[,3])
Finv3 <- function(p) quantile(reject_samp[,3], probs = p)
F4 <- ecdf(reject_samp[,4])
Finv4 <- function(p) quantile(reject_samp[,4], probs = p)

# create the Latin hypercube
X <- as.data.frame(lhs::randomLHS(Nlhs, klhs))
names(X) <- paste0("v", 1:klhs)

# transform with the empirical cdf  
transform <- function(W) {
  W$v1 <- Finv1(W$v1)
  W$v2 <- Finv2(W$v2)
  W$v3 <- Finv3(W$v3)
  W$v4 <- Finv4(W$v4)
  return(W)
}

X_t <- transform(X)

# check if the products are in range
p <- with(X_t, v1*v2*v3*v4)

ind <- which(p < lower | p > upper)
length(ind) # number out of range
#> [1] 27

# products are not in range, so use the Column-wise Pairwise swapping algorithm to find a suitable solution
pp <- length(which(p > lower & p < upper)) / nrow(X)

print(paste("starting at ", pp))
#> [1] "starting at  0.46"

swap <- function(i1, i2, j1, j2, dat) {
  temp <- dat[i1, j1]
  dat[i1, j1] <- dat[i2, j2]
  dat[i2, j2] <- temp
  return(dat)
}

maxiter <- 20000 # change as needed
iter <- 1
ppbest <- pp
while (pp < 1 & iter < maxiter) {
  for (j in 1:4) {
    for (i1 in 1:(Nlhs-1)) {
      for (i2 in i1:Nlhs) {
        if (pp < 1) {
          X <- swap(i1, i2, j, j, X)
          X_t <- transform(X)
          p <- with(X_t, v1*v2*v3*v4)
          pp <- length(which(p > lower & p < upper)) / nrow(X)
          if (pp <= ppbest) {
            # swap back
            X <- swap(i2, i1, j, j, X)
          } else {
            # stay here and update best
            ppbest <- pp
            break
          }
          iter <- iter + 1
        }
      }
      if (pp == 1 | iter > maxiter) break
    }
    if (pp == 1 | iter > maxiter) break
  }
}

print(paste("ending at ", pp))
#> [1] "ending at  1"

# check latin hypercube properties
all(apply(X, 2, function(x) sum(floor(x*Nlhs)+1)) == Nlhs*(Nlhs+1)/2)
#> [1] TRUE
all(X > 0 & X < 1)
#> [1] TRUE

# check all products are in range
X_t <- transform(X)
p <- with(X_t, v1*v2*v3*v4)
all(p > lower & p < upper)
#> [1] TRUE

par(mfrow=c(2,2))  
pairs(X[,1:4])

pairs(X_t[,1:4])

^{Created on 2024-04-28 with reprex v2.1.0}