Latin Hypercube Sampling with input correlation matrix in Python

Question

I am attempting to take into account the interaction among 5 different parameters in a Latin hypercube design in Python 3.8. However, I can't understand how, probably also because - not being my field - I never used these methods. I initially posted this on Stack Overflow, but it seems more relevant here.

I have an initial ensemble of 500 simulations, where I perturb 5 parameters. The 500 members are initially created by generating a Latin hypercube with 5 dimensions, sampled 500 times from uniform distributions.

The parameters are input model parameters, and I have limited knowledge of them. The outputs are forecasts of a plume. I have observations available for comparing the forecasts to the real event, satellite data that detected the plume. However, the observations do not provide any information about these parameters.

Among the 5 parameters, I am expecting correlation especially among X, Y and Z as all the three parameters are used to calculate the source strength in the model.

The comparison between forecasts and real observations (satellite data), however, will allow me to see which forecast is closer to the observations. Therefore the ensemble members closer to the observations will help in constraining my input model parameters ranges.

I use the parameters ranges from these retained simulations to fit new distributions and resample each parameter in a Latin hypercube for a new ensemble. The data are skewed and Gamma is a good fit. So, for now without considering interaction, I am operating like this:

```python    
# Fit the parameters with a gamma distribution:
a_X, loc_X, scale_X = stats.gamma.fit(X)
a_Y, loc_Y, scale_Y = stats.gamma.fit(Y)
a_Z, loc_Z, scale_Z = stats.gamma.fit(Z)
a_R, loc_R, scale_R = stats.gamma.fit(R)
a_T, loc_T, scale_T = stats.gamma.fit(T)

# Create my new LHS design using the LHS function from pyDOE:
parameter_sets = lhs(5, samples=500)

# Transform it to be gamma distributed:
parameter_sets[:, 0] = gamma(a=a_X, loc=loc_X, scale=scale_X).ppf(parameter_sets[:, 0])
parameter_sets[:, 1] = gamma(a=a_Y, loc=loc_Y, scale=scale_Y).ppf(parameter_sets[:, 1])
parameter_sets[:, 2] = gamma(a=a_Z, loc=loc_Z, scale=scale_Z).ppf(parameter_sets[:, 2])
parameter_sets[:, 3] = gamma(a=a_R, loc=loc_R, scale=scale_R).ppf(parameter_sets[:, 3])
parameter_sets[:, 4] = gamma(a=a_T, loc=loc_T, scale=scale_T).ppf(parameter_sets[:, 4])

If I run these new 500 simulations with this new parameter set, and compare them with the observations, the number of simulations closer to the observations increase considerably. However, I now want to improve the method, and take into account interaction among the parameters. I created a correlation matrix, and it seems that the major interactions are among X, Y and Z (I was expecting it). The correlation matrix I get, is something like this

           X         Y         Z         R         T
X       1.000000 -0.477159 -0.135367  0.074064 -0.045986
Y      -0.477159  1.000000 -0.501240 -0.127851 -0.139549
Z      -0.135367 -0.501240  1.000000 -0.078035 -0.026585
R       0.074064 -0.127851 -0.078035  1.000000 -0.004523
T      -0.045986 -0.139549 -0.026585 -0.004523  1.000000

However, I can’t figure out how to integrate the correlation matrix in the LHS.. By calculating the means and the covariance matrix I could use a multivariate normal, as I often read in these cases. But it assumes that all inputs are normally distributed that is not my case. And also in this case, I can't figure out then how to sample the multivariate from the LHS. Any help is greatly appreciated!

EDIT:

after looking into using multivariate copulas (still need to understand the concept properly, but I wanted to give it a go already), I modified the method like this (using this package specifically: Copulas):

from copulas.multivariate import GaussianMultivariate
from copulas.univariate import GammaUnivariate

# Here I specify which distribution to use, gamma in my case     
dist = GaussianMultivariate(distribution={

        "X": GammaUnivariate,
        "Y": GammaUnivariate,
        "Z": GammaUnivariate,
        "R": GammaUnivariate,
        "T": GammaUnivariate})
    copula = GaussianMultivariate()
    copula.fit(data)

    # Create a new LH
    test_lhs = lhs(5, samples=500)

    # Transform test_lhs based on the copula
    test_lhs[:, 0] = copula.univariates[0].ppf(test_lhs[:, 0])
    test_lhs[:, 1] = copula.univariates[1].ppf(test_lhs[:, 1])
    test_lhs[:, 2] = copula.univariates[2].ppf(test_lhs[:, 2])
    test_lhs[:, 3] = copula.univariates[3].ppf(test_lhs[:, 3])
    test_lhs[:, 4] = copula.univariates[4].ppf(test_lhs[:, 4])

I'll keep looking into this, but I wonder if I am on the right direction or not?

Answer 1

I don't think you are on the right path. If I understand your question correctly, you have the following situation:

1. You have empirical observations from some process with features ${X, Y, Z, R, T}$. You don't state it, but I assume there is some dependent variable as well.
2. You next draw a Latin hypercube sample with marginal distributions that match ${X, Y, Z, R, T}$. A Latin hypercube is normally drawn for the independent parameters, not for dependent parameters.
3. You then compare the observations to the Latin hypercube and keep the ones that are the closest.

That process is backward from the purpose of Latin hypercube sampling. You have two choices:

Choice A: Stick with observations that you have from your experiment and conduct the analysis. The Latin hypercube design is not adding any information or statistical power to the pre-existing observations.

Choice B: Start over. Draw the parameters from the Latin hypercube first. Run your process at the settings described by the Latin hypercube parameters. Analyze the output of the process.

A Latin hypercube sample is an experimental design sample that should be drawn before any experiments are conducted. See this paper to get a better understanding of why Latin hypercubes are used: https://www.jstor.org/stable/1268522

If I misunderstood your situation, I'm happy to edit my answer to be more helpful.

EDIT: (Second Attempt)

Observed Data: Plume detections from satellite data. No observations of $X,Y,Z,R,T$ Simulation Data: 500 Latin hypercube samples from $X, Y, Z, R, T$ that lead to plume simulations. Select the closest simulations that match the observed data. (How many are retained?) Fit gamma distributions to the parameters that match the observed data. Resample a Latin hypercube from the fitted gamma distributions. Desire: want to have a correlated sample on the second Latin hypercube.

Here is how I would attack this problem. I assume that your goal is to determine the $X, Y, Z, R, T$ that created the observed plume detection.

Observed Data: Plume detections from satellite data. No observations of $X,Y,Z,R,T$ Simulation Data: 500 Latin hypercube samples from $X, Y, Z, R, T$ that lead to plume simulations. Model: Create a GLM or other type of model using the simulation data as training data. Prediction: Using the model, find predicted the set of $X, Y, Z, R, T$ that lead to the observed data (with 95% confidence)

The reason I like this method better than continuing with another round of Latin hypercube sampling is that it allows for a measure of statistical uncertainty and it doesn't rely on a decision of how close is close enough to subselect the first lhs. Finally, once you have the predictions from the model, you can validate the predictions with another round of simulations.