1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
"""
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).
"""
__all__ = ["lhs_init", "cov_init"]
from numpy import asarray, diag, empty, eye
from . import util
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
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()
|
