Latin hypercube sampler¶

Welcome to the lhs documentation. This package implements Latin hypercube sampling in order to draw near-random samples of parameter values from multi-dimensional distributions. This package is primarily intended for scenario modelling.

License¶

The code is distributed under the terms of the BSD 3-Clause license (see LICENSE), and the documentation is distributed under the terms of the Creative Commons BY-SA 4.0 license.

A simple example¶

>>> import lhs
>>> import numpy as np
>>> import scipy.stats
>>> # Define normal distributions for alpha, beta, and gamma.
>>> alpha = {'name': 'norm', 'args': {'loc': 5, 'scale': 1}}
>>> beta = {'name': 'norm', 'args': {'loc': 10, 'scale': 2}}
>>> gamma = {'distribution': scipy.stats.norm(loc=0, scale=0.1)}
>>> # Collect the distributions in a dictionary.
>>> distributions = {'alpha': alpha, 'beta': beta, 'gamma': gamma}
>>> # Define the number of samples to draw.
>>> num_samples = 100
>>> # Create a pseudo-random number generator.
>>> seed = 12345
>>> rand = np.random.default_rng(seed=seed)
>>> # Draw samples for all three parameters.
>>> # This returns a dictionary of 1-D NumPy arrays.
>>> samples = lhs.draw(rand, num_samples, distributions)
>>> # Print the first sampled value for alpha.
>>> print(samples['alpha'][0])
4.538220677728063
>>> # Collect the samples for each parameter into a 2-D array.
>>> table = np.column_stack([values for values in samples.values()])
>>> # Print the first five samples.
>>> print(table[:5])
[[ 4.53822068e+00  1.14024721e+01  6.26798628e-02]
 [ 3.48555342e+00  1.15887999e+01 -8.20221380e-02]
 [ 6.91828587e+00  1.34933849e+01  3.75658088e-04]
 [ 6.02953073e+00  1.44184693e+01 -1.67071425e-01]
 [ 2.08898887e+00  9.06176782e+00  2.83558279e-02]]
>>> # Save the samples table to a text file.
>>> # np.savetxt('samples.out', table)

A more detailed example¶

Define a distribution for each independent parameter (in this example, R0 and eta):

def independent_params():
    return {
        # Basic reproduction number.
        'R0': {
            'name': 'constant',
            'args': {'value': 2.53},
            'shape': 1,
        },
        # Proportion of infections that require hospitalisation ("severe").
        # See doi:10.1371/journal.pone.0014505 for details.
        'eta': {
            'distribution': scipy.stats.loguniform(a=0.05, b=0.1),
            'shape': 1,
        },
    }

Note

See Supported distributions for details about defining distributions.

Define the name and size of each dependent parameter (if any; in this example, alpha_m):

def dependent_params():
    return {
        # Proportion of non-severe infections that seek clinical care.
        # See doi:10.1371/journal.pone.0014505 for details.
        'alpha_m': {'shape': 1},
    }

Define the distribution for each dependent parameter (if any), given the sampled values for every independent parameter, in a separate function:

def dependent_dists(indep_values, dep_params):
    # Scale eta from [0.001, 0.1] to [0.2, 1.0]
    eta = indep_values['eta']
    scale = (100 * eta - 0.1) * 0.8 / 9.9 + 0.2
    # Calculate the mean and bounds for alpha_m.
    mean_val = scale * 0.5
    min_val = scale * 0.25
    max_val = scale * 0.75
    # Calculate the mean and variance for alpha_m.
    varn = scale * 0.1
    span = max_val - min_val
    scaled_mean = (mean_val - min_val) / span
    scaled_var = varn / span
    # Calculate the shape parameters for alpha_m.
    a = np.square(scaled_mean / scaled_var) * (1 - scaled_mean) - scaled_mean
    b = a * (1 - scaled_mean) / scaled_mean

    dist = scipy.stats.beta(a=a, b=b, loc=min_val, scale=span)

    return {
        'alpha_m': {'distribution': dist},
    }

Draw samples for all of the model parameters:

def draw_samples(seed, num_samples):
    rng = np.random.default_rng(seed=seed)
    samples = lhs.draw(rng, num_samples, independent_params(),
                       dependent_params(), dependent_dists)
    return samples

Note

If there are no dependent parameters, you can omit the dependent parameter arguments and simply call:

samples = lhs.draw(rng, num_samples, independent_params())

See Drawing samples for further details about drawing samples.

You can then use these samples to, e.g., define initial model parameters.

def test_draw_samples():
    samples = draw_samples(seed=12345, num_samples=1000)
    # Ensure that R0 is constant.
    assert(all(samples['R0'] == 2.53))
    # Ensure that all eta values lie within the interval bounds.
    assert(all(samples['eta'] > 0.05))
    assert(all(samples['eta'] < 0.10))
    # Ensure that between 5% and 75% of infections lead to non-severe clinical
    # presentations.
    assert(all(samples['alpha_m'] > 0.05))
    assert(all(samples['alpha_m'] < 0.75))

User Documentation

Development