2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ After some research, I came to realize that the editions by @rcarnell describe my question much more clearly than I could express initially. Thanks for the support $\endgroup$ – rms Jul 9 at 17:57
2
$\begingroup$

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 = "")
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ thank you for the answer, @RCarnell ! Just a small additional question: Is it possible to get a uniform distribution by changing the line Y[,1:3] <- qdirichlet(X[,1:3], rep(2,3)) to Y[,1:3] <- qdirichlet(X[,1:3], rep(1,3))? I understand that this would be the case with regular Dirichlet distribution, as answered in this question $\endgroup$ – rms Jul 11 at 1:32
  • $\begingroup$ No. With the regular Dirichlet distribution and with this qdirichlet function, if you specify $\alpha = (1,1)$ you get two uniform distributions. If you specify $\alpha = (1,1,1)$, you get "flat" distributions, but not uniform. They are decreasing in density from left to right. I don't agree with the answer that was given to the question you cited, I'll take a look at responding. Please accept my answer if it was helpful. $\endgroup$ – R Carnell Jul 12 at 19:08
  • $\begingroup$ could you indicate reference(s) I can use to support the usage of this approach? The objective is both self-learning and writing a methodology section. $\endgroup$ – rms Jul 16 at 22:10
  • 1
    $\begingroup$ For Latin hypercubes: Large Sample Properties of Simulations Using Latin Hypercube Sampling, Stein, 1987, Technometrics. For the fact that a Dirichlet distribution can be simulated by drawing gamma distributions and dividing by the sum, web.archive.org/web/20150219021331/https://…. For the qdirichlet function, just me! $\endgroup$ – R Carnell Jul 17 at 1:39
0
$\begingroup$

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]
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ @RCarnell, I've translated your R code to Python. Please check if the credits are ok. $\endgroup$ – rms Jul 16 at 22:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.