Latin hypercube sampling with python

Question

I would like to sample a distribution defined by a function in multiple dimensions (2,3,4):

f(x, y, ...) = ...

The distributions might be ugly, non standard (like a 3D spline on data, sum of gaussians ect.). To this end I would like to uniformly sample the 2..4 dimensional space, and than with an additional random number accept or reject the given point of the space into my sample.

1. Is there a ready to use python lib for this purpose?

2. Is there python lib for generating the points in this 2..4 dimensional space with latin hypercube sampling, or with other uniform sampling method? Bruteforce sampling with independent random numbers usually results in more and less dense regimes of the space.

3. if 1) and 2) doesn't exist, is there anybody who is kind enough to share his implementation for the same or similar problem.

I'll use it in a python code, but links to other solutions are also acknowledged.

Solution 1

I guess this is a late answer, but this is also for future visitors. I have just put up an implementation of multi-dimensional uniform Latin Hypercube sampling on git. It is minimal, but very easy to use. You can generate uniform random variables sampled in n dimensions using Latin Hypercube Sampling, if your variables are independent. Below is an example plot comparing Monte Carlo and Latin Hypercube Sampling with Multi-dimensional Uniformity (LHS-MDU) in two dimensions with zero correlation.

import lhsmdu
import matplotlib.pyplot as plt
import numpy

l = lhsmdu.sample(2,10) # Latin Hypercube Sampling of two variables, and 10 samples each.
k = lhsmdu.createRandomStandardUniformMatrix(2,10) # Monte Carlo Sampling

fig = plt.figure()
ax = fig.gca()
ax.set_xticks(numpy.arange(0,1,0.1))
ax.set_yticks(numpy.arange(0,1,0.1))
plt.scatter(k[0], k[1], color="b", label="LHS-MDU")
plt.scatter(l[0], l[1], color="r", label="MC")
plt.grid()
plt.show()

Solution 2

Now the pyDOE library provides a tool to generate Latin-hypercube-based samples.

https://pythonhosted.org/pyDOE/randomized.html

to generate samples over n dimensions:

lhs(n, [samples, criterion, iterations])

where n is the number of dimensions, samples as the total number of the sample space.

Solution 3

This 2-D example samples uniformly on two dimensions, chooses each point with constant probability (thus keeping a binomially distributed number of points), selects randomly and without replacement those points from the sample space, and generates a pair of vectors which you can then pass through to your function f:

import numpy as np
import random
resolution = 10
keepprob = 0.5
min1, max1 = 0., 1.
min2, max2 = 3., 11. 
keepnumber = np.random.binomial(resolution * resolution, keepprob,1)
array1,array2  = np.meshgrid(np.linspace(min1,max1,resolution),np.linspace(min2,max2,resolution))
randominixes  = random.sample(list(range(resolution * resolution)), int(keepnumber))
randominixes.sort()
vec1Sampled,vec2Sampled  = array1.flatten()[randominixes],array2.flatten()[randominixes]

Solution 4

Here is an update of Sahil M's answer for Python 3 (update from Python 2 to Python 3 and some minor code changes to match code and figure):

import lhsmdu
import matplotlib.pyplot as plt
import numpy

l = lhsmdu.sample(2,10) # Latin Hypercube Sampling of two variables, and 10 samples each.
k = lhsmdu.createRandomStandardUniformMatrix(2,10) # Monte Carlo Sampling

fig = plt.figure()
ax = fig.gca()
ax.set_xticks(numpy.arange(0,1,0.1))
ax.set_yticks(numpy.arange(0,1,0.1))
plt.scatter([k[0]], [k[1]], color="r", label="MC")
plt.scatter([l[0]], [l[1]], color="b", label="LHS-MDU")
plt.legend()
plt.grid()
plt.show()

I once encountered a Python memory error running this script. Any suggestions why this could happen or how to changes the script so it doesn't happen anymore in the future?

Solution 5

Latin Hypercube sampling is now part of SciPy since version 1.7. See the doc.

from scipy.stats.qmc import LatinHypercube

engine = LatinHypercube(d=2)
sample = engine.random(n=100)

View more solutions