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

Source Code for Package csb.statistics.pdf