Using Latin Hypercube Sampling with a condition that the sum of two variables should be less than one

Question

I am building an experimental design with 4 variables defined on (0,1). In notation, $x_i \in [0,1]$ with $ i=1,..., 4$). Two of these variables must satisfy the condition that $x_1 + x_2 \leq 1$. How can I perform Latin Hypercube Sampling with this condition?

I thought about rejection sampling when $x_1+x_2 > 1$, but realize that rejection sampling does not work with Latin hypercube sampling.

Answer 1

Strategy:

1. Draw $X_1, ..., X_5$ from a uniform LHS
2. Transform $X_1, X_2, X_3$ such that $X_1+X_2+X_3=1$ using the strategy I explained previously for R. The basic idea is to transform the marginal draws using the quantiles of gamma functions, then normalize those gamma quantiles. The result is a distribution like a Dirichlet distribution (although not exactly).
3. Drop $X_3$ since it is not necessary. If $X_1+X_2+X_3=1$ and $X_i > 0$ then $X_1 + X_2 < 1$.
4. Transform $X_4$ and $X_5$ to the desired distribution

require(lhs)
qdirichlet <- function(X, alpha)
{
qdirichlet is not an exact quantile function since the quantile of a
multivariate distribtion is not unique
qdirichlet is also not the quantiles of the marginal distributions since
those quantiles do not sum to one
qdirichlet is the quantile of the underlying gamma functions, normalized
This has been tested to show that qdirichlet approximates the dirichlet
distribution well and creates the correct marginal means and variances
when using a latin hypercube sample
lena <- length(alpha)
  stopifnot(is.matrix(X))
  sims <- dim(X)[1]
  stopifnot(dim(X)[2] == lena)
  if(any(is.na(alpha)) || any(is.na(X)))
    stop("NA values not allowed in qdirichlet")
Y <- matrix(0, nrow=sims, ncol=lena)
  ind <- which(alpha != 0)
  for(i in ind)
  {
    Y[,i] <- qgamma(X[,i], alpha[i], 1)
  }
  Y <- Y / rowSums(Y)
  return(Y)
}
set.seed(19753)
X <- randomLHS(500, 5)
Y <- X
transform X1, X2, X3 such that X1 + X2 + X3 =1
change the alpha parameter to change the mean of X1 and X2
Y[,1:3] <- qdirichlet(X[,1:3], rep(2,3))
transform parameter 4 and 5
Y[,4] <- qnorm(X[,4], 2, 1)
Y[,5] <- qunif(X[,5], 1, 3)
drop the unncessary X3
Y <- Y[,-3]
check that X1 + X2 < 1
stopifnot(all(Y[,1] + Y[,2] < 1.0))
plots
par(mfrow = c(2,2))
for (i in c(1,2,4,5))
  hist(X[,i], breaks = 20, main = i, xlab = "")
par(mfrow = c(2,2))
for (i in 1:4)
  hist(Y[,i], breaks = 20, main = i, xlab = "")

Answer 2

In order to implement the strategy described by @RCarnell in python, this is a translation of the function qdirichlet. The usage is similar to the one presented in the original answer

def dirichlet_ppf(X, alpha):
    # dirichlet_ppf is not an exact quantile function since the quantile of a
    #  multivariate distribtion is not unique
    # dirichlet_ppf is also not the quantiles of the marginal distributions since
    #  those quantiles do not sum to one
    # dirichlet_ppf is the quantile of the underlying gamma functions, normalized
    # This has been tested to show that dirichlet_ppf approximates the dirichlet
    #  distribution well and creates the correct marginal means and variances
    #  when using a latin hypercube sample
    #
    # Python translation of qdirichlet function by  R. Carnell
    # original: https://stats.stackexchange.com/a/476433/244679
    import numpy as np
    from scipy.stats import gamma
X = np.asarray(X)
alpha = np.asarray(alpha)

assert alpha.ndim == 1, "parameter alpha must be a vector"
assert X.ndim == 2, "parameter X must be an array with samples as rows and variables as columns"
assert X.shape[1] == alpha.shape[0], "number of variables in each row of X and length of alpha must be equal"
assert not (np.any(np.isnan(X)) or np.any(np.isnan(alpha))), "NAN values are not allowed in dirichlet_ppf"

Y = np.zeros(shape=X.shape)
for idx, a in enumerate(alpha):
    if a != 0. :
        Y[:, idx] = gamma.ppf(X[:, idx], a)

return Y / Y.sum(axis=1)[:, np.newaxis]