Package csb :: Package statistics :: Package pdf

[frames] | no frames]

# Source Code for Package csb.statistics.pdf

      1  """
      2  Probability density functions.
      3
      4  This module defines L{AbstractDensity}: a common interface for all PDFs.
      5  Each L{AbstractDensity} describes a specific type of probability distribution,
      6  for example L{Normal} is an implementation of the Gaussian distribution:
      7
      8      >>> pdf = Normal(mu=10, sigma=1.1)
      9      >>> pdf.mu, pdf['sigma']
     10      10.0, 1.1
     11
     12  Every PDF provides an implementation of the L{AbstractDensity.evaluate}
     13  method, which evaluates the PDF for a list of input data points:
     14
     15      >>> pdf.evaluate([10, 9, 11, 12])
     16      array([ 0.3626748 ,  0.2399147 ,  0.2399147 ,  0.06945048])
     17
     18  PDF instances also behave like functions:
     19
     20      >>> pdf(data)    # the same as pdf.evaluate(data)
     21
     22  Some L{AbstractDensity} implementations may support drawing random numbers from
     23  the distribution (or raise an exception otherwise):
     24
     25      >>> pdf.random(2)
     26      array([ 9.86257083,  9.73760515])
     27
     28  Each implementation of L{AbstractDensity} may support infinite number of estimators,
     29  used to estimate and re-initialize the PDF parameters from a set of observed data
     30  points:
     31
     32      >>> pdf.estimate([5, 5, 10, 10])
     33      >>> pdf.mu, pdf.sigma
     34      (7.5, 2.5)
     35      >>> pdf.estimator
     36      <csb.statistics.pdf.GaussianMLEstimator>
     37
     38  Estimators implement the L{AbstractEstimator} interface. They are treated as
     39  pluggable tools, which can be exchanged through the L{AbstractDensity.estimator}
     40  property (you could create, initialize and plug your own estimator as well).
     41  This is a classic Strategy pattern.
     42  """
     43
     44  import numpy.random
     45  import scipy.special
     46  import csb.core
     47
     48  from abc import ABCMeta, abstractmethod
     49  from csb.core import OrderedDict
     50
     51  from csb.numeric import log, exp, psi, inv_psi, EULER_MASCHERONI
     52  from scipy.special import gammaln
     53  from numpy import array, fabs, power, sqrt, pi, mean, median, clip
     54
     55
     56 -class IncompatibleEstimatorError(TypeError):
     57      pass
     58
     59 -class ParameterNotFoundError(AttributeError):
     60      pass
     61
     62 -class ParameterValueError(ValueError):
     63
     64 -    def __init__(self, param, value):
     65
     66          self.param = param
     67          self.value = value
     68
     69          super(ParameterValueError, self).__init__(param, value)
     70
     71 -    def __str__(self):
     72          return '{0} = {1}'.format(self.param, self.value)
     73
     74 -class EstimationFailureError(ParameterValueError):
     75      pass
     76
     77 -class AbstractEstimator(object):
     78      """
     79      Density parameter estimation strategy.
     80      """
     81
     82      __metaclass__ = ABCMeta
     83
     84      @abstractmethod
     85 -    def estimate(self, context, data):
     86          """
     87          Estimate the parameters of the distribution from same {data}.
     88
     89          @param context: context distribution
     90          @type context: L{AbstractDensity}
     91          @param data: sample values
     92          @type data: array
     93
     94          @return: a new distribution, initialized with the estimated parameters
     95          @rtype: L{AbstractDensity}
     96
     97          @raise EstimationFailureError: if estimation is not possible
     98          """
     99          pass
    100
    101 -class NullEstimator(AbstractEstimator):
    102      """
    103      Does not estimate anything.
    104      """
    105 -    def estimate(self, context, data):
    106          raise NotImplementedError()
    107
    108 -class LaplaceMLEstimator(AbstractEstimator):
    109
    110 -    def estimate(self, context, data):
    111
    112          x = array(data)
    113
    114          mu = median(x)
    115          b = mean(fabs(x - mu))
    116
    117          return Laplace(mu, b)
    118
    119 -class GaussianMLEstimator(AbstractEstimator):
    120
    121 -    def estimate(self, context, data):
    122
    123          x = array(data)
    124
    125          mu = mean(x)
    126          sigma = sqrt(mean((x - mu) ** 2))
    127
    128          return Normal(mu, sigma)
    129
    130 -class InverseGaussianMLEstimator(AbstractEstimator):
    131
    132 -    def estimate(self, context, data):
    133
    134          x = array(data)
    135
    136          mu = mean(x)
    137          il = mean((1.0 / x) - (1.0 / mu))
    138
    139          if il == 0:
    140              raise EstimationFailureError('lambda', float('inf'))
    141
    142          return InverseGaussian(mu, 1.0 / il)
    143
    144 -class GammaMLEstimator(AbstractEstimator):
    145
    146 -    def __init__(self):
    147          super(GammaMLEstimator, self).__init__()
    148          self.n_iter = 1000
    149
    150
    151 -    def estimate(self, context, data):
    152
    153          mu = mean(data)
    154          logmean = mean(log(data))
    155
    156          a = 0.5 / (log(mu) - logmean)
    157
    158          for dummy in range(self.n_iter):
    159
    160              a = inv_psi(logmean - log(mu) + log(a))
    161
    162          return Gamma(a, a / mu)
    163
    164 -class GenNormalBruteForceEstimator(AbstractEstimator):
    165
    166 -    def __init__(self, minbeta=0.5, maxbeta=8.0, step=0.1):
    167
    168          self._minbeta = minbeta
    169          self._maxbeta = maxbeta
    170          self._step = step
    171
    172          super(GenNormalBruteForceEstimator, self).__init__()
    173
    174 -    def estimate(self, context, data):
    175
    176          pdf = GeneralizedNormal(1, 1, 1)
    177          data = array(data)
    178          logl = []
    179
    180          for beta in numpy.arange(self._minbeta, self._maxbeta, self._step):
    181
    182              self.update(pdf, data, beta)
    183
    184              l = pdf.log_prob(data).sum()
    185              logl.append([beta, l])
    186
    187          logl = numpy.array(logl)
    188
    189          # optimal parameters:
    190          beta = logl[ numpy.argmax(logl[:, 1]) ][0]
    191          self.update(pdf, data, beta)
    192
    193          return pdf
    194
    195 -    def estimate_with_fixed_beta(self, data, beta):
    196
    197          mu = median(data)
    198          v = mean((data - mu) ** 2)
    199          alpha = sqrt(v * exp(gammaln(1. / beta) - gammaln(3. / beta)))
    200
    201          return mu, alpha
    202
    203 -    def update(self, pdf, data, beta):
    204
    205          mu, alpha = self.estimate_with_fixed_beta(data, beta)
    206
    207          pdf.mu = mu
    208          pdf.alpha = alpha
    209          pdf.beta = beta
    210
    211          return pdf
    212
    213 -class MultivariateGaussianMLEstimator(AbstractEstimator):
    214
    215 -    def __init__(self):
    216          super(MultivariateGaussianMLEstimator, self).__init__()
    217
    218 -    def estimate(self, context, data):
    219          return MultivariateGaussian(numpy.mean(data, 0), numpy.cov(data.T))
    220
    221 -class DirichletEstimator(AbstractEstimator):
    222
    223 -    def __init__(self):
    224          super(DirichletEstimator, self).__init__()
    225          self.n_iter = 1000
    226          self.tol = 1e-5
    227
    228 -    def estimate(self, context, data):
    229
    230          log_p = numpy.mean(log(data), 0)
    231
    232          e = numpy.mean(data, 0)
    233          v = numpy.mean(data ** 2, 0)
    234          q = (e[0] - v[0]) / (v[0] - e[0] ** 2)
    235
    236          a = e * q
    237          y = a * 0
    238          k = 0
    239          while(sum(abs(y - a)) > self.tol and k < self.n_iter):
    240              y = psi(sum(a)) + log_p
    241              a = numpy.array(list(map(inv_psi, y)))
    242              k += 1
    243
    244          return Dirichlet(a)
    245
    246 -class GumbelMinMomentsEstimator(AbstractEstimator):
    247
    248 -    def estimate(self, context, data):
    249
    250          x = array(data)
    251
    252          beta = sqrt(6 * numpy.var(x)) / pi
    253          mu = mean(x) + EULER_MASCHERONI * beta
    254
    255          return GumbelMinimum(mu, beta)
    256
    257 -class GumbelMaxMomentsEstimator(AbstractEstimator):
    258
    259 -    def estimate(self, context, data):
    260
    261          x = array(data)
    262
    263          beta = sqrt(6 * numpy.var(x)) / pi
    264          mu = mean(x) - EULER_MASCHERONI * beta
    265
    266          return GumbelMaximum(mu, beta)
    267
    268
    269 -class AbstractDensity(object):
    270      """
    271      Defines the interface and common operations for all probability density
    272      functions. This is a generic class which can operate on parameters of
    273      any type (e.g. simple floats or custom parameter objects).
    274
    275      Subclasses must complete the implementation by implementing the
    276      L{AbstractDensity.log_prob} method. Subclasses could also consider--but
    277      are not obliged to--override the L{AbstractDensity.random} method. If
    278      any of the density parameters need validation, subclasses are expected to
    279      override the L{AbstractDensity._validate} method and raise
    280      L{ParameterValueError} on validation failure. Note that implementing
    281      parameter validation in property setters has almost no effect and is
    282      discouraged.
    283      """
    284
    285      __metaclass__ = ABCMeta
    286
    287
    288 -    def __init__(self):
    289
    290          self._params = OrderedDict()
    291          self._estimator = None
    292
    293          self.estimator = NullEstimator()
    294
    295 -    def __getitem__(self, param):
    296
    297          if param in self._params:
    298              return self._params[param]
    299          else:
    300              raise ParameterNotFoundError(param)
    301
    302 -    def __setitem__(self, param, value):
    303
    304          if param in self._params:
    305              self._validate(param, value)
    306              self._set(param, value)
    307          else:
    308              raise ParameterNotFoundError(param)
    309
    310 -    def _set(self, param, value):
    311          """
    312          Update the C{value} of C{param}.
    313          """
    314          self._params[param] = value
    315
    316      @property
    317 -    def estimator(self):
    318          return self._estimator
    319      @estimator.setter
    320 -    def estimator(self, strategy):
    321          if not isinstance(strategy, AbstractEstimator):
    322              raise TypeError(strategy)
    323          self._estimator = strategy
    324
    325 -    def __call__(self, x):
    326          return self.evaluate(x)
    327
    328 -    def __str__(self):
    329          name = self.__class__.__name__
    330          params = ', '.join([ '{0}={1}'.format(p, v) for p, v in self._params.items() ])
    331
    332          return '{0}({1})'.format(name, params)
    333
    334 -    def _register(self, name):
    335          """
    336          Register a new parameter name.
    337          """
    338          if name not in self._params:
    339              self._params[name] = None
    340
    341 -    def _validate(self, param, value):
    342          """
    343          Parameter value validation hook.
    344          @raise ParameterValueError: on failed validation (value not accepted)
    345          """
    346          pass
    347
    348 -    def get_params(self):
    349          return [self._params[name] for name in  self.parameters]
    350
    351 -    def set_params(self, *values, **named_params):
    352
    353          for p, v in zip(self.parameters, values):
    354              self[p] = v
    355
    356          for p in named_params:
    357              self[p] = named_params[p]
    358
    359      @property
    360 -    def parameters(self):
    361          """
    362          Get a list of all distribution parameter names.
    363          """
    364          return tuple(self._params)
    365
    366      @abstractmethod
    367 -    def log_prob(self, x):
    368          """
    369          Evaluate the logarithm of the probability of observing values C{x}.
    370
    371          @param x: values
    372          @type x: array
    373          @rtype: array
    374          """
    375          pass
    376
    377 -    def evaluate(self, x):
    378          """
    379          Evaluate the probability of observing values C{x}.
    380
    381          @param x: values
    382          @type x: array
    383          @rtype: array
    384          """
    385          x = numpy.array(x)
    386          return exp(self.log_prob(x))
    387
    388 -    def random(self, size=None):
    389          """
    390          Generate random samples from the probability distribution.
    391
    392          @param size: number of values to sample
    393          @type size: int
    394          """
    395          raise NotImplementedError()
    396
    397 -    def estimate(self, data):
    398          """
    399          Estimate and load the parameters of the distribution from sample C{data}
    400          using the current L{AbstractEstimator} strategy.
    401
    402          @param data: sample values
    403          @type data: array
    404
    405          @raise NotImplementedError: when no estimator is available for this
    406                                      distribution
    407          @raise IncompatibleEstimatorError: when the current estimator is not
    408                                             compatible with this pdf
    409          """
    410
    411          try:
    412              pdf = self.estimator.estimate(self, data)
    413
    414              for param in pdf.parameters:
    415                  self[param] = pdf[param]
    416
    417          except ParameterNotFoundError as e:
    418              raise IncompatibleEstimatorError(self.estimator)
    419
    420          except ParameterValueError as e:
    421              raise EstimationFailureError(e.param, e.value)
    422
    423
    424 -class BaseDensity(AbstractDensity):
    425      """
    426      Base abstract class for all PDFs, which operate on simple float
    427      or array-of-float parameters. Parameters of any other type will trigger
    428      TypeError-s.
    429      """
    430
    431 -    def _set(self, param, value):
    432
    433          try:
    434              if csb.core.iterable(value):
    435                  value = array(value)
    436              else:
    437                  value = float(value)
    438          except ValueError:
    439              raise TypeError(value)
    440
    441          super(BaseDensity, self)._set(param, value)
    442
    443 -class Laplace(BaseDensity):
    444
    445 -    def __init__(self, mu=0, b=1):
    446
    447          super(Laplace, self).__init__()
    448
    449          self._register('mu')
    450          self._register('b')
    451
    452          self.set_params(b=b, mu=mu)
    453          self.estimator = LaplaceMLEstimator()
    454
    455 -    def _validate(self, param, value):
    456
    457          if param == 'b' and value <= 0:
    458              raise ParameterValueError(param, value)
    459
    460      @property
    461 -    def b(self):
    462          return self['b']
    463      @b.setter
    464 -    def b(self, value):
    465          self['b'] = value
    466
    467      @property
    468 -    def mu(self):
    469          return self['mu']
    470      @mu.setter
    471 -    def mu(self, value):
    472          self['mu'] = value
    473
    474 -    def log_prob(self, x):
    475
    476          b = self.b
    477          mu = self.mu
    478
    479          return log(1 / (2. * b)) - fabs(x - mu) / b
    480
    481 -    def random(self, size=None):
    482
    483          loc = self.mu
    484          scale = self.b
    485
    486          return numpy.random.laplace(loc, scale, size)
    487
    488 -class Normal(BaseDensity):
    489
    490 -    def __init__(self, mu=0, sigma=1):
    491
    492          super(Normal, self).__init__()
    493
    494          self._register('mu')
    495          self._register('sigma')
    496
    497          self.set_params(mu=mu, sigma=sigma)
    498          self.estimator = GaussianMLEstimator()
    499
    500      @property
    501 -    def mu(self):
    502          return self['mu']
    503      @mu.setter
    504 -    def mu(self, value):
    505          self['mu'] = value
    506
    507      @property
    508 -    def sigma(self):
    509          return self['sigma']
    510      @sigma.setter
    511 -    def sigma(self, value):
    512          self['sigma'] = value
    513
    514 -    def log_prob(self, x):
    515
    516          mu = self.mu
    517          sigma = self.sigma
    518
    519          return log(1.0 / sqrt(2 * pi * sigma ** 2)) - (x - mu) ** 2 / (2 * sigma ** 2)
    520
    521 -    def random(self, size=None):
    522
    523          mu = self.mu
    524          sigma = self.sigma
    525
    526          return numpy.random.normal(mu, sigma, size)
    527
    528 -class InverseGaussian(BaseDensity):
    529
    530 -    def __init__(self, mu=1, shape=1):
    531
    532          super(InverseGaussian, self).__init__()
    533
    534          self._register('mu')
    535          self._register('shape')
    536
    537          self.set_params(mu=mu, shape=shape)
    538          self.estimator = InverseGaussianMLEstimator()
    539
    540 -    def _validate(self, param, value):
    541
    542          if value <= 0:
    543              raise ParameterValueError(param, value)
    544
    545      @property
    546 -    def mu(self):
    547          return self['mu']
    548      @mu.setter
    549 -    def mu(self, value):
    550          self['mu'] = value
    551
    552      @property
    553 -    def shape(self):
    554          return self['shape']
    555      @shape.setter
    556 -    def shape(self, value):
    557          self['shape'] = value
    558
    559 -    def log_prob(self, x):
    560
    561          mu = self.mu
    562          scale = self.shape
    563          x = numpy.array(x)
    564
    565          if numpy.min(x) <= 0:
    566              raise ValueError('InverseGaussian is defined for x > 0')
    567
    568          y = -0.5 * scale * (x - mu) ** 2 / (mu ** 2 * x)
    569          z = 0.5 * (log(scale) - log(2 * pi * x ** 3))
    570          return  z + y
    571
    572
    573 -    def random(self, size=None):
    574
    575          mu = self.mu
    576          shape = self.shape
    577
    578          mu_2l = mu / shape / 2.
    579          Y = numpy.random.standard_normal(size)
    580          Y = mu * Y ** 2
    581          X = mu + mu_2l * (Y - sqrt(4 * shape * Y + Y ** 2))
    582          U = numpy.random.random(size)
    583
    584          m = numpy.less_equal(U, mu / (mu + X))
    585
    586          return m * X + (1 - m) * mu ** 2 / X
    587
    588 -class GeneralizedNormal(BaseDensity):
    589
    590 -    def __init__(self, mu=0, alpha=1, beta=1):
    591
    592          super(GeneralizedNormal, self).__init__()
    593
    594          self._register('mu')
    595          self._register('alpha')
    596          self._register('beta')
    597
    598          self.set_params(mu=mu, alpha=alpha, beta=beta)
    599          self.estimator = GenNormalBruteForceEstimator()
    600
    601 -    def _validate(self, param, value):
    602
    603          if param in ('alpha, beta') and value <= 0:
    604              raise ParameterValueError(param, value)
    605
    606      @property
    607 -    def mu(self):
    608          return self['mu']
    609      @mu.setter
    610 -    def mu(self, value):
    611          self['mu'] = value
    612
    613      @property
    614 -    def alpha(self):
    615          return self['alpha']
    616      @alpha.setter
    617 -    def alpha(self, value):
    618          self['alpha'] = value
    619
    620      @property
    621 -    def beta(self):
    622          return self['beta']
    623      @beta.setter
    624 -    def beta(self, value):
    625          self['beta'] = value
    626
    627 -    def log_prob(self, x):
    628
    629          mu = self.mu
    630          alpha = self.alpha
    631          beta = self.beta
    632
    633          return log(beta / (2.0 * alpha)) - gammaln(1. / beta) - power(fabs(x - mu) / alpha, beta)
    634
    635 -class GeneralizedInverseGaussian(BaseDensity):
    636
    637 -    def __init__(self, a=1, b=1, p=1):
    638          super(GeneralizedInverseGaussian, self).__init__()
    639
    640          self._register('a')
    641          self._register('b')
    642          self._register('p')
    643          self.set_params(a=a, b=b, p=p)
    644
    645          self.estimator = NullEstimator()
    646
    647 -    def _validate(self, param, value):
    648
    649          if value <= 0:
    650              raise ParameterValueError(param, value)
    651
    652      @property
    653 -    def a(self):
    654          return self['a']
    655      @a.setter
    656 -    def a(self, value):
    657          self['a'] = value
    658
    659      @property
    660 -    def b(self):
    661          return self['b']
    662      @b.setter
    663 -    def b(self, value):
    664          self['b'] = value
    665
    666      @property
    667 -    def p(self):
    668          return self['p']
    669      @p.setter
    670 -    def p(self, value):
    671          self['p'] = value
    672
    673 -    def log_prob(self, x):
    674
    675          a = self['a']
    676          b = self['b']
    677          p = self['p']
    678
    679          lz = 0.5 * p * (log(a) - log(b)) - log(2 * scipy.special.kv(p, sqrt(a * b)))
    680
    681          return lz + (p - 1) * log(x) - 0.5 * (a * x + b / x)
    682
    683 -    def random(self, size=None):
    684
    685          from csb.statistics.rand import inv_gaussian
    686
    687          rvs = []
    688          burnin = 10
    689          a = self['a']
    690          b = self['b']
    691          p = self['p']
    692
    693          s = a * 0. + 1.
    694
    695          if p < 0:
    696              a, b = b, a
    697
    698          if size == None:
    699              size = 1
    700          for i in range(int(size)):
    701              for j in range(burnin):
    702
    703                  l = b + 2 * s
    704                  m = sqrt(l / a)
    705
    706                  x = inv_gaussian(m, l, shape=m.shape)
    707                  s = numpy.random.gamma(abs(p) + 0.5, x)
    708
    709              if p >= 0:
    710                  rvs.append(x)
    711              else:
    712                  rvs.append(1 / x)
    713
    714          return numpy.array(rvs)
    715
    716 -class Gamma(BaseDensity):
    717
    718 -    def __init__(self, alpha=1, beta=1):
    719          super(Gamma, self).__init__()
    720
    721          self._register('alpha')
    722          self._register('beta')
    723
    724          self.set_params(alpha=alpha, beta=beta)
    725          self.estimator = GammaMLEstimator()
    726
    727 -    def _validate(self, param, value):
    728
    729          if value <= 0:
    730              raise ParameterValueError(param, value)
    731
    732      @property
    733 -    def alpha(self):
    734          return self['alpha']
    735      @alpha.setter
    736 -    def alpha(self, value):
    737          self['alpha'] = value
    738
    739      @property
    740 -    def beta(self):
    741          return self['beta']
    742      @beta.setter
    743 -    def beta(self, value):
    744          self['beta'] = value
    745
    746 -    def log_prob(self, x):
    747
    748          a, b = self['alpha'], self['beta']
    749
    750          return a * log(b) - gammaln(clip(a, 1e-308, 1e308)) + \
    751                 (a - 1) * log(clip(x, 1e-308, 1e308)) - b * x
    752
    753 -    def random(self, size=None):
    754          return numpy.random.gamma(self['alpha'], 1 / self['beta'], size)
    755
    756 -class InverseGamma(BaseDensity):
    757
    758 -    def __init__(self, alpha=1, beta=1):
    759          super(InverseGamma, self).__init__()
    760
    761          self._register('alpha')
    762          self._register('beta')
    763
    764          self.set_params(alpha=alpha, beta=beta)
    765          self.estimator = NullEstimator()
    766
    767 -    def _validate(self, param, value):
    768
    769          if value <= 0:
    770              raise ParameterValueError(param, value)
    771
    772      @property
    773 -    def alpha(self):
    774          return self['alpha']
    775      @alpha.setter
    776 -    def alpha(self, value):
    777          self['alpha'] = value
    778
    779      @property
    780 -    def beta(self):
    781          return self['beta']
    782      @beta.setter
    783 -    def beta(self, value):
    784          self['beta'] = value
    785
    786 -    def log_prob(self, x):
    787          a, b = self['alpha'], self['beta']
    788          return a * log(b) - gammaln(a) - (a + 1) * log(x) - b / x
    789
    790 -    def random(self, size=None):
    791          return 1. / numpy.random.gamma(self['alpha'], 1 / self['beta'], size)
    792
    793 -class MultivariateGaussian(Normal):
    794
    795 -    def __init__(self, mu=numpy.zeros(2), sigma=numpy.eye(2)):
    796
    797          super(MultivariateGaussian, self).__init__(mu, sigma)
    798          self.estimator = MultivariateGaussianMLEstimator()
    799
    800 -    def random(self, size=None):
    801          return numpy.random.multivariate_normal(self.mu, self.sigma, size)
    802
    803 -    def log_prob(self, x):
    804
    805          from numpy.linalg import det
    806
    807          mu = self.mu
    808          S = self.sigma
    809          D = len(mu)
    810          q = self.__q(x)
    811          return -0.5 * (D * log(2 * pi) + log(abs(det(S)))) - 0.5 * q ** 2
    812
    813 -    def __q(self, x):
    814          from numpy import sum, dot, reshape
    815          from numpy.linalg import inv
    816
    817          mu = self.mu
    818          S = self.sigma
    819
    820          return sqrt(clip(sum(reshape((x - mu) * dot(x - mu, inv(S).T.squeeze()), (-1, len(mu))), -1), 0., 1e308))
    821
    822 -    def conditional(self, x, dims):
    823          """
    824          Return the distribution along the dimensions
    825          dims conditioned on x
    826
    827          @param x: conditional values
    828          @param dims: new dimensions
    829          """
    830          from numpy import take, dot
    831          from numpy.linalg import inv
    832
    833          dims2 = [i for i in range(self['mu'].shape[0]) if not i in dims]
    834
    835          mu1 = take(self['mu'], dims)
    836          mu2 = take(self['mu'], dims2)
    837
    838          # x1 = take(x, dims)
    839          x2 = take(x, dims2)
    840
    841          A = take(take(self['Sigma'], dims, 0), dims, 1)
    842          B = take(take(self['Sigma'], dims2, 0), dims2, 1)
    843          C = take(take(self['Sigma'], dims, 0), dims2, 1)
    844
    845          mu = mu1 + dot(C, dot(inv(B), x2 - mu2))
    846          Sigma = A - dot(C, dot(inv(B), C.T))
    847
    848          return MultivariateGaussian((mu, Sigma))
    849
    850 -class Dirichlet(BaseDensity):
    851
    852 -    def __init__(self, alpha):
    853          super(Dirichlet, self).__init__()
    854
    855          self._register('alpha')
    856
    857          self.set_params(alpha=alpha)
    858          self.estimator = DirichletEstimator()
    859
    860      @property
    861 -    def alpha(self):
    862          return self['alpha']
    863
    864      @alpha.setter
    865 -    def alpha(self, value):
    866          self['alpha'] = numpy.ravel(value)
    867
    868 -    def log_prob(self, x):
    869          #TODO check wether x is in the probability simplex
    870          alpha = self.alpha
    871          return gammaln(sum(alpha)) - sum(gammaln(alpha)) \
    872                + numpy.dot((alpha - 1).T, log(x).T)
    873
    874 -    def random(self, size=None):
    875          return numpy.random.mtrand.dirichlet(self.alpha, size)
    876
    877
    878 -class GumbelMinimum(BaseDensity):
    879
    880 -    def __init__(self, mu=0, beta=1):
    881          super(GumbelMinimum, self).__init__()
    882
    883          self._register('mu')
    884          self._register('beta')
    885
    886          self.set_params(mu=mu, beta=beta)
    887          self.estimator = GumbelMinMomentsEstimator()
    888
    889 -    def _validate(self, param, value):
    890
    891          if param == 'beta' and value <= 0:
    892              raise ParameterValueError(param, value)
    893
    894      @property
    895 -    def mu(self):
    896          return self['mu']
    897      @mu.setter
    898 -    def mu(self, value):
    899          self['mu'] = value
    900
    901      @property
    902 -    def beta(self):
    903          return self['beta']
    904      @beta.setter
    905 -    def beta(self, value):
    906          self['beta'] = value
    907
    908 -    def log_prob(self, x):
    909
    910          mu = self.mu
    911          beta = self.beta
    912
    913          z = (x - mu) / beta
    914          return log(1. / beta) + z - exp(z)
    915
    916 -    def random(self, size=None):
    917
    918          mu = self.mu
    919          beta = self.beta
    920
    921          return -numpy.random.gumbel(-mu, beta, size)
    922
    923 -class GumbelMaximum(GumbelMinimum):
    924
    925 -    def __init__(self, mu=0, beta=1):
    926
    927          super(GumbelMaximum, self).__init__(mu=mu, beta=beta)
    928          self.estimator = GumbelMaxMomentsEstimator()
    929
    930 -    def log_prob(self, x):
    931
    932          mu = self.mu
    933          beta = self.beta
    934
    935          z = (x - mu) / beta
    936          return log(1. / beta) - z - exp(-z)
    937
    938 -    def random(self, size=None):
    939
    940          mu = self.mu
    941          beta = self.beta
    942
    943          return numpy.random.gumbel(mu, beta, size)
    944