Package csb :: Package statistics :: Module rand |
[frames] | no frames]

# Source Code for Module csb.statistics.rand

      1  """
      2  Random number generators
      3  """
      4
      5 -def probability_transform(shape, inv_cum, cum_min=0., cum_max=1.):
      6      """
      7      Generic sampler based on the probability transform.
      8
      9      @param shape: shape of the random sample
     10      @param inv_cum: inversion of the cumulative density function from which one seeks to sample
     11      @param cum_min: lower value of the cumulative distribution
     12      @param cum_max: upper value of the cumulative distribution
     13      @return: random variates of the PDF implied by the inverse cumulative distribution
     14      """
     15      from numpy.random import random
     16
     17      return inv_cum(cum_min + random(shape) * (cum_max - cum_min))
     18
     19 -def truncated_gamma(shape=None, alpha=1., beta=1., x_min=None, x_max=None):
     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
     51 -def truncated_normal(shape=None, mu=0., sigma=1., x_min=None, x_max=None):
     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
     83 -def sample_dirichlet(alpha, n_samples=1):
     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
    100 -def sample_sphere3d(radius=1., n_samples=1):
    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
    126 -def sample_from_histogram(p, n_samples=1):
    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
    146 -def gen_inv_gaussian(a, b, p, burnin=10):
    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
    172 -def inv_gaussian(mu=1., _lambda=1., shape=None):
    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
    188 -def random_rotation(A, n_iter=10, initial_values=None):
    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          ## sample alpha and gamma
    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          ## sample beta
    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