Generate samples from multivariate correlated data which have non-parametric cumulative distribution functions

Question

I have 40 samples that contain information about 6 variables (hence a 40x6 data matrix). Each variable (column) has a cumulative distribution function (marginal distribution) based on the 40 values, which is not Gaussian, whereas correlation matrix exist that relates the variable among each other (I am using Spearsman rho ). What I would like to do, is to generate 1000 samples for each variable (a 1000x6 data matrix) that satisfies the correlation values of the original sampling datasets, and follows the empirical cumulative distribution function of the original 40 samples for each variable. My understanding is that I can use copula for that purpose. So, the steps that I have in mind are:

Calculate the mean and standard deviation for each variable, and generate gaussian marginal distributions for a given correlation (Spearsman rho) matrix. The values of the gaussian variables are not restricted in [0,1]. Next, convert the marginal distribution of each variable into a unit distribution (range of values for each variable lies within [0,1]). This is done without changing the ordering of the data, hence the ranking correlation coefficient remains the same. Lastly, I am using the cumulative distribution functions of the initial dataset (40x6 data matrix) to generate simulated values for each variable, particularly by quantile-matching (for a specific value of the unit variable, which I assume that equals the cumulative distribution value (vertical axis), I find the simulated value by finding the value at the x-axis that corresponds to this cdf value). Is the above description accurate? Sorry for the long description, I am new on this and I am not very confident to jump parts of the discussion.

Answer 1

What I would like to do, is to generate 1000 samples for each variable (a 1000x6 data matrix) that satisfies the correlation values of the original sampling datasets, and follows the empirical cumulative distribution function of the original 40 samples for each variable.

You could do some fancy modeling of the margins and the copula, sure, but a possible solution is to sample from the empirical distribution itself the way bootstrapping samples with replacement from the empirical distributions. In this regard, whatever the empirical margins and empirical copula are, those become the population margins and copula, so your $iid$ samples drawn from the data follow the empirical distributions.

Example code is below.

# Load the "iris" data
#
data("iris")

# Define the sample size you want
# Unlike bootstrap, this need not coincide with the original data size
#
N <- 100 # The "iris" data set has 150 multivariate observations
         # We will take only 100

# Sample indices WITH REPLACEMENT (with replacement is crucial)
#
idx <- sample(seq(1, N, 1), N, replace = T)

# Draw iid samples from the empirical distribution by selecting the rows
# corresponding to the indices idx
# You can think of this as being like an rempirical(N) function in R-lingo or
# np.random.empirical(N) in Python/NumPy
#
S <- iris[idx, ]

# S is your sample drawn from the empirical distribution, which 
# functions as the population as far as S is concerned

While the upside is that you do not have to do fancy distribution modeling that may suffer from errors, a downside of doing this is that you only every get values that you have in your data. There is no interpolating when you do it this way (but you do have your sample following the empirical margins and correlations).