Source code for bumps.dream.initpop

Population initialization routines.

To start the analysis an initial population is required.  This will be
an array of size M x N, where M is the number of dimensions in the fitting
problem and N is the number of Markov chains.

Two functions are provided:

1. lhs_init(N, bounds) returns a latin hypercube sampling, which tests every
parameter at each of N levels.

2. cov_init(N, x, cov) returns a Gaussian sample along the ellipse
defined by the covariance matrix, cov.  Covariance defaults to
diag(dx) if dx is provided as a parameter, or to I if it is not.

Additional options are random box: rand(M, N) or random scatter: randn(M, N).

from __future__ import division, print_function

__all__ = ['lhs_init', 'cov_init']

from numpy import eye, diag, asarray, empty
from . import util

[docs]def lhs_init(N, bounds): """ Latin Hypercube Sampling Returns an array whose columns each have *N* samples from equally spaced bins between *bounds=(xmin, xmax)* for the column. DREAM bounds objects, with bounds.low and bounds.high can be used as well. Note: Indefinite ranges are not supported. """ try: xmin, xmax = bounds.low, bounds.high except AttributeError: xmin, xmax = bounds # Define the size of xmin nvar = len(xmin) # Initialize array ran with random numbers ran = util.rng.rand(N, nvar) # Initialize array s with zeros s = empty((N, nvar)) # Now fill s for j in range(nvar): # Random permutation idx = util.rng.permutation(N)+1 p = (idx-ran[:, j])/N s[:, j] = xmin[j] + p*(xmax[j]-xmin[j]) return s
[docs]def cov_init(N, x, cov=None, dx=None): """ Initialize *N* sets of random variables from a gaussian model. The center is at *x* with an uncertainty ellipse specified by the 1-sigma independent uncertainty values *dx* or the full covariance matrix uncertainty *cov*. For example, create an initial population for 20 sequences for a model with local minimum x with covariance matrix C:: pop = cov_init(cov=C, x=x, N=20) """ #return mean + dot(util.rng.randn(N, len(mean)), chol(cov)) if cov is None and dx is None: cov = eye(len(x)) elif cov is None: cov = diag(asarray(dx)**2) return util.rng.multivariate_normal(mean=x, cov=cov, size=N)
def demo(): from numpy import arange print("Three ways of calling cov_init:") print("with cov", cov_init(N=4, x=[5, 6], cov=diag([0.1, 0.001]))) print("with dx", cov_init(N=4, x=[5, 6], dx=[0.1, 0.001])) print("with nothing", cov_init(N=4, x=[5, 6])) print(""" The following array should have four columns. Column 1 should have the numbers from 10 to 19, column 2 from 20 to 29, etc. The columns are in random order with a random fractional part. """) pop = lhs_init(N=10, bounds=(arange(1, 5), arange(2, 6)))*10 print(pop) if __name__ == "__main__": demo()