Latin Hypercube Sampling from Multivariate Distribution

[iwishiwasaneagle]

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?

[Balandat]

The standard way to sample from a MVN is to draw iid normal samples and then correlate them with a root of the covariance matrix: xi = mu + Lv where LL^T = Sigma. GPyTorch has support for this on the MVN but not for LHS. We have some utility for generating Sobol qMC samples in BoTorch that could be used instead here: https://github.com/pytorch/botorch/blob/main/botorch/utils/sampling.py#L199

It’s not LHS but may be useful.

If you do need to compute CDFs of MVNs we are about to add this functionality: meta-pytorch/botorch#1394

[iwishiwasaneagle]

Thanks for your reply @Balandat. I'll keep an eye on botorch's MVN CDF. Will that support the inverse CDF as well?

[Balandat]

So since this is an MVN the CDF is not injective so the inverse CDF doesn't exist (at least in the strict sense). You probably don't want to think about inverses in this setting.