Latin Hypercube Sampling of a Gaussian Process I have a small set ( $n \approx 200$ ) $(x,y)$ data points representing the find location of a person on a landscape. From this, I generate a $20\times 20$ histogram and fit the fairly basic model
class GPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super(GPModelWithDerivatives, self).__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ZeroMean()
        self.base_kernel = gpytorch.kernels.MaternKernel(ard_num_dims=2)
        self.covar_module = gpytorch.kernels.ScaleKernel(self.base_kernel)
        
    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

with good results

If I were to use the "true" histogram to sample points, it would be contained to the few cells with values which is why I've used a GP. In order to "sample" a position from the GP, I do the following:
predictions = likelihood(model(test_x))
mean_pred_np = predictions.mean.numpy()
p = (mean_pred_np+np.abs(mean_pred_np.min()))
p /= p.sum()
ind = np.random.choice(np.arange(mean_pred_np.size), p = p,size=1000)
pts = test_x[ind]

(where pts are the white dots on the image above, larger = more samples at that point)
Ideally, I would like to sample this using Latin Hypercube Sampling as the output of this is the input to another Monte-Carlo simulation.  However, the crux for me is that gpytorch.distributions.MultivariateNormal has no CDF implementation, and as far as I am aware the Multivariate Normal also has no inverse CDF. I can probably do some hacky gradient descent to fit the $\frac{1}{n^2}$ sized areas but before I go down that route I thought I'd ask here for some help.
Is LHS using GPytorch feasible or is there some other solution that I should be exploring?
 A: Basic desired analysis pipeline:

*

*2D data (200 pts)

*Fit a bivariate normal distribution to the data

*Sample from the bivariate normal with a Latin hypercube sample

*Continue with another Monte Carlo simulation

Here are some options:

*

*Instead of fitting the bivariate normal, fit two univariate normals to the margins and use those to transform the Latin hypercube.  Issue: loss of the original covariance structure

*Instead of fitting the bivariate normal distribution, just bootstrap sample the original 200 data points in the last Monte Carlo simulation.  Issue: If you believe that the other parts of the 2D space are valid, but didn't show up in the original sample, you might want to sample from more than the simple data locations.

*Instead of a 20x20 histogram of the original data, create a 10x10 histogram of the original data and then take that forward to the final Monte Carlo sample (sample a grid cell, then sample inside the grid cell).  This will "spread-out" the original data points, but will stay more true to the original data.  Issue: still might have missing areas from the 2D space in the last sample.

*Create a 2D histogram directly from the density of the bivariate normal, then sample the grid cell with the histogram probability and then sample uniformly inside the grid cell.  Issue:  Not really a Latin hypercube.

