1 """
2 Random number generators
3 """
4
18
20 """
21 Generate random variates from a lower-and upper-bounded gamma distribution.
22
23 @param shape: shape of the random sample
24 @param alpha: shape parameter (alpha > 0.)
25 @param beta: scale parameter (beta >= 0.)
26 @param x_min: lower bound of variate
27 @param x_max: upper bound of variate
28 @return: random variates of lower-bounded gamma distribution
29 """
30 from scipy.special import gammainc, gammaincinv
31 from numpy.random import gamma
32 from numpy import inf
33
34 if x_min is None and x_max is None:
35 return gamma(alpha, 1 / beta, shape)
36 elif x_min is None:
37 x_min = 0.
38 elif x_max is None:
39 x_max = inf
40
41 x_min = max(0., x_min)
42 x_max = min(1e300, x_max)
43
44 a = gammainc(alpha, beta * x_min)
45 b = gammainc(alpha, beta * x_max)
46
47 return probability_transform(shape,
48 lambda x, alpha=alpha: gammaincinv(alpha, x),
49 a, b) / beta
50
52 """
53 Generates random variates from a lower-and upper-bounded normal distribution
54
55 @param shape: shape of the random sample
56 @param mu: location parameter
57 @param sigma: width of the distribution (sigma >= 0.)
58 @param x_min: lower bound of variate
59 @param x_max: upper bound of variate
60 @return: random variates of lower-bounded normal distribution
61 """
62 from scipy.special import erf, erfinv
63 from numpy.random import standard_normal
64 from numpy import inf, sqrt
65
66 if x_min is None and x_max is None:
67 return standard_normal(shape) * sigma + mu
68 elif x_min is None:
69 x_min = -inf
70 elif x_max is None:
71 x_max = inf
72
73 x_min = max(-1e300, x_min)
74 x_max = min(+1e300, x_max)
75 var = sigma ** 2 + 1e-300
76 sigma = sqrt(2 * var)
77
78 a = erf((x_min - mu) / sigma)
79 b = erf((x_max - mu) / sigma)
80
81 return probability_transform(shape, erfinv, a, b) * sigma + mu
82
84 """
85 Sample points from a dirichlet distribution with parameter alpha.
86
87 @param alpha: alpha parameter of a dirichlet distribution
88 @type alpha: array
89 """
90 from numpy import array, sum, transpose, ones
91 from numpy.random import gamma
92
93 alpha = array(alpha, ndmin=1)
94 X = gamma(alpha,
95 ones(len(alpha)),
96 [n_samples, len(alpha)])
97
98 return transpose(transpose(X) / sum(X, -1))
99
101 """
102 Sample points from 3D sphere.
103
104 @param radius: radius of the sphere
105 @type radius: float
106
107 @param n_samples: number of samples to return
108 @type n_samples: int
109
110 @return: n_samples times random cartesian coordinates inside the sphere
111 @rtype: numpy array
112 """
113 from numpy.random import random
114 from numpy import arccos, transpose, cos, sin, pi, power
115
116 r = radius * power(random(n_samples), 1 / 3.)
117 theta = arccos(2. * (random(n_samples) - 0.5))
118 phi = 2 * pi * random(n_samples)
119
120 x = cos(phi) * sin(theta) * r
121 y = sin(phi) * sin(theta) * r
122 z = cos(theta) * r
123
124 return transpose([x, y, z])
125
127 """
128 returns the indice of bin according to the histogram p
129
130 @param p: histogram
131 @type p: numpy.array
132 @param n_samples: number of samples to generate
133 @type n_samples: integer
134 """
135
136 from numpy import add, less, argsort, take, arange
137 from numpy.random import random
138
139 indices = argsort(p)
140 indices = take(indices, arange(len(p) - 1, -1, -1))
141
142 c = add.accumulate(take(p, indices)) / add.reduce(p)
143
144 return indices[add.reduce(less.outer(c, random(n_samples)), 0)]
145
147 """
148 Sampler based on Gibbs sampling.
149 Assumes scalar p.
150 """
151 from numpy.random import gamma
152 from numpy import sqrt
153
154 s = a * 0. + 1.
155
156 if p < 0:
157 a, b = b, a
158
159 for i in range(burnin):
160
161 l = b + 2 * s
162 m = sqrt(l / a)
163
164 x = inv_gaussian(m, l, shape=m.shape)
165 s = gamma(abs(p) + 0.5, x)
166
167 if p >= 0:
168 return x
169 else:
170 return 1 / x
171
173 """
174 Generate random samples from inverse gaussian.
175 """
176 from numpy.random import standard_normal, random
177 from numpy import sqrt, less_equal, clip
178
179 mu_2l = mu / _lambda / 2.
180 Y = mu * standard_normal(shape) ** 2
181 X = mu + mu_2l * (Y - sqrt(4 * _lambda * Y + Y ** 2))
182 U = random(shape)
183
184 m = less_equal(U, mu / (mu + X))
185
186 return clip(m * X + (1 - m) * mu ** 2 / X, 1e-308, 1e308)
187
189 """
190 Generation of three-dimensional random rotations in
191 fitting and matching problems, Habeck 2009.
192
193 Generate random rotation R from::
194
195 exp(trace(dot(transpose(A), R)))
196
197 @param A: generating parameter
198 @type A: 3 x 3 numpy array
199
200 @param n_iter: number of gibbs sampling steps
201 @type n_iter: integer
202
203 @param initial_values: initial euler angles alpha, beta and gamma
204 @type initial_values: tuple
205
206 @rtype: 3 x 3 numpy array
207 """
208 from numpy import cos, sin, dot, pi, clip
209 from numpy.linalg import svd, det
210 from random import vonmisesvariate, randint
211 from csb.numeric import euler
212
213
214 def sample_beta(kappa, n=1):
215 from numpy import arccos
216 from csb.numeric import log, exp
217 from numpy.random import random
218
219 u = random(n)
220
221 if kappa != 0.:
222 x = clip(1 + 2 * log(u + (1 - u) * exp(-kappa)) / kappa, -1., 1.)
223 else:
224 x = 2 * u - 1
225
226 if n == 1:
227 return arccos(x)[0]
228 else:
229 return arccos(x)
230
231
232 U, L, V = svd(A)
233
234 if det(U) < 0:
235 L[2] *= -1
236 U[:, 2] *= -1
237 if det(V) < 0:
238 L[2] *= -1
239 V[2] *= -1
240
241 if initial_values is None:
242 beta = 0.
243 else:
244 alpha, beta, gamma = initial_values
245
246 for _i in range(n_iter):
247
248
249 phi = vonmisesvariate(0., clip(cos(beta / 2) ** 2 * (L[0] + L[1]), 1e-308, 1e10))
250 psi = vonmisesvariate(pi, sin(beta / 2) ** 2 * (L[0] - L[1]))
251 u = randint(0, 1)
252
253 alpha = 0.5 * (phi + psi) + pi * u
254 gamma = 0.5 * (phi - psi) + pi * u
255
256
257 kappa = cos(phi) * (L[0] + L[1]) + cos(psi) * (L[0] - L[1]) + 2 * L[2]
258 beta = sample_beta(kappa)
259
260 return dot(U, dot(euler(alpha, beta, gamma), V))
261