1# -*- coding: utf-8 -*- 2"""Multivariate Normal and t distributions 3 4 5 6Created on Sat May 28 15:38:23 2011 7 8@author: Josef Perktold 9 10TODO: 11* renaming, 12 - after adding t distribution, cov does not make sense for Sigma DONE 13 - should mean also be renamed to mu, if there will be distributions 14 with mean != mu 15* not sure about corner cases 16 - behavior with (almost) singular sigma or transforms 17 - df <= 2, is everything correct if variance is not finite or defined ? 18* check to return possibly univariate distribution for marginals or conditional 19 distributions, does univariate special case work? seems ok for conditional 20* are all the extra transformation methods useful outside of testing ? 21 - looks like I have some mixup in definitions of standardize, normalize 22* new methods marginal, conditional, ... just added, typos ? 23 - largely tested for MVNormal, not yet for MVT DONE 24* conditional: reusing, vectorizing, should we reuse a projection matrix or 25 allow for a vectorized, conditional_mean similar to OLS.predict 26* add additional things similar to LikelihoodModelResults? quadratic forms, 27 F distribution, and others ??? 28* add Delta method for nonlinear functions here, current function is hidden 29 somewhere in miscmodels 30* raise ValueErrors for wrong input shapes, currently only partially checked 31 32* quantile method (ppf for equal bounds for multiple testing) is missing 33 http://svitsrv25.epfl.ch/R-doc/library/mvtnorm/html/qmvt.html seems to use 34 just a root finder for inversion of cdf 35 36* normalize has ambiguous definition, and mixing it up in different versions 37 std from sigma or std from cov ? 38 I would like to get what I need for mvt-cdf, or not 39 univariate standard t distribution has scale=1 but std>1 40 FIXED: add std_sigma, and normalize uses std_sigma 41 42* more work: bivariate distributions, 43 inherit from multivariate but overwrite some methods for better efficiency, 44 e.g. cdf and expect 45 46I kept the original MVNormal0 class as reference, can be deleted 47 48 49See Also 50-------- 51sandbox/examples/ex_mvelliptical.py 52 53Examples 54-------- 55 56Note, several parts of these examples are random and the numbers will not be 57(exactly) the same. 58 59>>> import numpy as np 60>>> import statsmodels.sandbox.distributions.mv_normal as mvd 61>>> 62>>> from numpy.testing import assert_array_almost_equal 63>>> 64>>> cov3 = np.array([[ 1. , 0.5 , 0.75], 65... [ 0.5 , 1.5 , 0.6 ], 66... [ 0.75, 0.6 , 2. ]]) 67 68>>> mu = np.array([-1, 0.0, 2.0]) 69 70multivariate normal distribution 71-------------------------------- 72 73>>> mvn3 = mvd.MVNormal(mu, cov3) 74>>> mvn3.rvs(size=3) 75array([[-0.08559948, -1.0319881 , 1.76073533], 76 [ 0.30079522, 0.55859618, 4.16538667], 77 [-1.36540091, -1.50152847, 3.87571161]]) 78 79>>> mvn3.std 80array([ 1. , 1.22474487, 1.41421356]) 81>>> a = [0.0, 1.0, 1.5] 82>>> mvn3.pdf(a) 830.013867410439318712 84>>> mvn3.cdf(a) 850.31163181123730122 86 87Monte Carlo integration 88 89>>> mvn3.expect_mc(lambda x: (x<a).all(-1), size=100000) 900.30958999999999998 91>>> mvn3.expect_mc(lambda x: (x<a).all(-1), size=1000000) 920.31197399999999997 93 94multivariate t distribution 95--------------------------- 96 97>>> mvt3 = mvd.MVT(mu, cov3, 4) 98>>> mvt3.rvs(size=4) 99array([[-0.94185437, 0.3933273 , 2.40005487], 100 [ 0.07563648, 0.06655433, 7.90752238], 101 [ 1.06596474, 0.32701158, 2.03482886], 102 [ 3.80529746, 7.0192967 , 8.41899229]]) 103 104>>> mvt3.pdf(a) 1050.010402959362646937 106>>> mvt3.cdf(a) 1070.30269483623249821 108>>> mvt3.expect_mc(lambda x: (x<a).all(-1), size=1000000) 1090.30271199999999998 110 111>>> mvt3.cov 112array([[ 2. , 1. , 1.5], 113 [ 1. , 3. , 1.2], 114 [ 1.5, 1.2, 4. ]]) 115>>> mvt3.corr 116array([[ 1. , 0.40824829, 0.53033009], 117 [ 0.40824829, 1. , 0.34641016], 118 [ 0.53033009, 0.34641016, 1. ]]) 119 120get normalized distribution 121 122>>> mvt3n = mvt3.normalized() 123>>> mvt3n.sigma 124array([[ 1. , 0.40824829, 0.53033009], 125 [ 0.40824829, 1. , 0.34641016], 126 [ 0.53033009, 0.34641016, 1. ]]) 127>>> mvt3n.cov 128array([[ 2. , 0.81649658, 1.06066017], 129 [ 0.81649658, 2. , 0.69282032], 130 [ 1.06066017, 0.69282032, 2. ]]) 131 132What's currently there? 133 134>>> [i for i in dir(mvn3) if not i[0]=='_'] 135['affine_transformed', 'cdf', 'cholsigmainv', 'conditional', 'corr', 'cov', 136'expect_mc', 'extra_args', 'logdetsigma', 'logpdf', 'marginal', 'mean', 137'normalize', 'normalized', 'normalized2', 'nvars', 'pdf', 'rvs', 'sigma', 138'sigmainv', 'standardize', 'standardized', 'std', 'std_sigma', 'whiten'] 139 140>>> [i for i in dir(mvt3) if not i[0]=='_'] 141['affine_transformed', 'cdf', 'cholsigmainv', 'corr', 'cov', 'df', 'expect_mc', 142'extra_args', 'logdetsigma', 'logpdf', 'marginal', 'mean', 'normalize', 143'normalized', 'normalized2', 'nvars', 'pdf', 'rvs', 'sigma', 'sigmainv', 144'standardize', 'standardized', 'std', 'std_sigma', 'whiten'] 145 146""" 147import numpy as np 148from scipy import special 149 150from statsmodels.sandbox.distributions.multivariate import mvstdtprob 151from .extras import mvnormcdf 152 153 154def expect_mc(dist, func=lambda x: 1, size=50000): 155 '''calculate expected value of function by Monte Carlo integration 156 157 Parameters 158 ---------- 159 dist : distribution instance 160 needs to have rvs defined as a method for drawing random numbers 161 func : callable 162 function for which expectation is calculated, this function needs to 163 be vectorized, integration is over axis=0 164 size : int 165 number of random samples to use in the Monte Carlo integration, 166 167 168 Notes 169 ----- 170 this does not batch 171 172 Returns 173 ------- 174 expected value : ndarray 175 return of function func integrated over axis=0 by MonteCarlo, this will 176 have the same shape as the return of func without axis=0 177 178 Examples 179 -------- 180 181 integrate probability that both observations are negative 182 183 >>> mvn = mve.MVNormal([0,0],2.) 184 >>> mve.expect_mc(mvn, lambda x: (x<np.array([0,0])).all(-1), size=100000) 185 0.25306000000000001 186 187 get tail probabilities of marginal distribution (should be 0.1) 188 189 >>> c = stats.norm.isf(0.05, scale=np.sqrt(2.)) 190 >>> expect_mc(mvn, lambda x: (np.abs(x)>np.array([c, c])), size=100000) 191 array([ 0.09969, 0.0986 ]) 192 193 or calling the method 194 195 >>> mvn.expect_mc(lambda x: (np.abs(x)>np.array([c, c])), size=100000) 196 array([ 0.09937, 0.10075]) 197 198 199 ''' 200 def fun(x): 201 return func(x) # * dist.pdf(x) 202 rvs = dist.rvs(size=size) 203 return fun(rvs).mean(0) 204 205def expect_mc_bounds(dist, func=lambda x: 1, size=50000, lower=None, upper=None, 206 conditional=False, overfact=1.2): 207 '''calculate expected value of function by Monte Carlo integration 208 209 Parameters 210 ---------- 211 dist : distribution instance 212 needs to have rvs defined as a method for drawing random numbers 213 func : callable 214 function for which expectation is calculated, this function needs to 215 be vectorized, integration is over axis=0 216 size : int 217 minimum number of random samples to use in the Monte Carlo integration, 218 the actual number used can be larger because of oversampling. 219 lower : None or array_like 220 lower integration bounds, if None, then it is set to -inf 221 upper : None or array_like 222 upper integration bounds, if None, then it is set to +inf 223 conditional : bool 224 If true, then the expectation is conditional on being in within 225 [lower, upper] bounds, otherwise it is unconditional 226 overfact : float 227 oversampling factor, the actual number of random variables drawn in 228 each attempt are overfact * remaining draws. Extra draws are also 229 used in the integration. 230 231 232 Notes 233 ----- 234 this does not batch 235 236 Returns 237 ------- 238 expected value : ndarray 239 return of function func integrated over axis=0 by MonteCarlo, this will 240 have the same shape as the return of func without axis=0 241 242 Examples 243 -------- 244 >>> mvn = mve.MVNormal([0,0],2.) 245 >>> mve.expect_mc_bounds(mvn, lambda x: np.ones(x.shape[0]), 246 lower=[-10,-10],upper=[0,0]) 247 0.24990416666666668 248 249 get 3 marginal moments with one integration 250 251 >>> mvn = mve.MVNormal([0,0],1.) 252 >>> mve.expect_mc_bounds(mvn, lambda x: np.dstack([x, x**2, x**3, x**4]), 253 lower=[-np.inf,-np.inf], upper=[np.inf,np.inf]) 254 array([[ 2.88629497e-03, 9.96706297e-01, -2.51005344e-03, 255 2.95240921e+00], 256 [ -5.48020088e-03, 9.96004409e-01, -2.23803072e-02, 257 2.96289203e+00]]) 258 >>> from scipy import stats 259 >>> [stats.norm.moment(i) for i in [1,2,3,4]] 260 [0.0, 1.0, 0.0, 3.0] 261 262 263 ''' 264 #call rvs once to find length of random vector 265 rvsdim = dist.rvs(size=1).shape[-1] 266 if lower is None: 267 lower = -np.inf * np.ones(rvsdim) 268 else: 269 lower = np.asarray(lower) 270 if upper is None: 271 upper = np.inf * np.ones(rvsdim) 272 else: 273 upper = np.asarray(upper) 274 275 def fun(x): 276 return func(x) # * dist.pdf(x) 277 278 rvsli = [] 279 used = 0 #remain = size #inplace changes size 280 total = 0 281 while True: 282 remain = size - used #just a temp variable 283 rvs = dist.rvs(size=int(remain * overfact)) 284 total += int(size * overfact) 285 286 rvsok = rvs[((rvs >= lower) & (rvs <= upper)).all(-1)] 287 #if rvsok.ndim == 1: #possible shape problems if only 1 random vector 288 rvsok = np.atleast_2d(rvsok) 289 used += rvsok.shape[0] 290 291 rvsli.append(rvsok) #[:remain]) use extras instead 292 print(used) 293 if used >= size: 294 break 295 rvs = np.vstack(rvsli) 296 print(rvs.shape) 297 assert used == rvs.shape[0] #saftey check 298 mean_conditional = fun(rvs).mean(0) 299 if conditional: 300 return mean_conditional 301 else: 302 return mean_conditional * (used * 1. / total) 303 304 305def bivariate_normal(x, mu, cov): 306 """ 307 Bivariate Gaussian distribution for equal shape *X*, *Y*. 308 309 See `bivariate normal 310 <http://mathworld.wolfram.com/BivariateNormalDistribution.html>`_ 311 at mathworld. 312 """ 313 X, Y = np.transpose(x) 314 mux, muy = mu 315 sigmax, sigmaxy, tmp, sigmay = np.ravel(cov) 316 sigmax, sigmay = np.sqrt(sigmax), np.sqrt(sigmay) 317 Xmu = X-mux 318 Ymu = Y-muy 319 320 rho = sigmaxy/(sigmax*sigmay) 321 z = Xmu**2/sigmax**2 + Ymu**2/sigmay**2 - 2*rho*Xmu*Ymu/(sigmax*sigmay) 322 denom = 2*np.pi*sigmax*sigmay*np.sqrt(1-rho**2) 323 return np.exp( -z/(2*(1-rho**2))) / denom 324 325 326 327class BivariateNormal(object): 328 329 330 #TODO: make integration limits more flexible 331 # or normalize before integration 332 333 def __init__(self, mean, cov): 334 self.mean = mu 335 self.cov = cov 336 self.sigmax, self.sigmaxy, tmp, self.sigmay = np.ravel(cov) 337 self.nvars = 2 338 339 def rvs(self, size=1): 340 return np.random.multivariate_normal(self.mean, self.cov, size=size) 341 342 def pdf(self, x): 343 return bivariate_normal(x, self.mean, self.cov) 344 345 def logpdf(self, x): 346 #TODO: replace this 347 return np.log(self.pdf(x)) 348 349 def cdf(self, x): 350 return self.expect(upper=x) 351 352 def expect(self, func=lambda x: 1, lower=(-10,-10), upper=(10,10)): 353 def fun(x, y): 354 x = np.column_stack((x,y)) 355 return func(x) * self.pdf(x) 356 from scipy.integrate import dblquad 357 return dblquad(fun, lower[0], upper[0], lambda y: lower[1], 358 lambda y: upper[1]) 359 360 def kl(self, other): 361 '''Kullback-Leibler divergence between this and another distribution 362 363 int f(x) (log f(x) - log g(x)) dx 364 365 where f is the pdf of self, and g is the pdf of other 366 367 uses double integration with scipy.integrate.dblquad 368 369 limits currently hardcoded 370 371 ''' 372 fun = lambda x : self.logpdf(x) - other.logpdf(x) 373 return self.expect(fun) 374 375 def kl_mc(self, other, size=500000): 376 fun = lambda x : self.logpdf(x) - other.logpdf(x) 377 rvs = self.rvs(size=size) 378 return fun(rvs).mean() 379 380class MVElliptical(object): 381 '''Base Class for multivariate elliptical distributions, normal and t 382 383 contains common initialization, and some common methods 384 subclass needs to implement at least rvs and logpdf methods 385 386 ''' 387 #getting common things between normal and t distribution 388 389 390 def __init__(self, mean, sigma, *args, **kwds): 391 '''initialize instance 392 393 Parameters 394 ---------- 395 mean : array_like 396 parameter mu (might be renamed), for symmetric distributions this 397 is the mean 398 sigma : array_like, 2d 399 dispersion matrix, covariance matrix in normal distribution, but 400 only proportional to covariance matrix in t distribution 401 args : list 402 distribution specific arguments, e.g. df for t distribution 403 kwds : dict 404 currently not used 405 406 ''' 407 408 self.extra_args = [] 409 self.mean = np.asarray(mean) 410 self.sigma = sigma = np.asarray(sigma) 411 sigma = np.squeeze(sigma) 412 self.nvars = nvars = len(mean) 413 #self.covchol = np.linalg.cholesky(sigma) 414 415 416 #in the following sigma is original, self.sigma is full matrix 417 if sigma.shape == (): 418 #iid 419 self.sigma = np.eye(nvars) * sigma 420 self.sigmainv = np.eye(nvars) / sigma 421 self.cholsigmainv = np.eye(nvars) / np.sqrt(sigma) 422 elif (sigma.ndim == 1) and (len(sigma) == nvars): 423 #independent heteroskedastic 424 self.sigma = np.diag(sigma) 425 self.sigmainv = np.diag(1. / sigma) 426 self.cholsigmainv = np.diag( 1. / np.sqrt(sigma)) 427 elif sigma.shape == (nvars, nvars): #python tuple comparison 428 #general 429 self.sigmainv = np.linalg.pinv(sigma) 430 self.cholsigmainv = np.linalg.cholesky(self.sigmainv).T 431 else: 432 raise ValueError('sigma has invalid shape') 433 434 #store logdetsigma for logpdf 435 self.logdetsigma = np.log(np.linalg.det(self.sigma)) 436 437 def rvs(self, size=1): 438 '''random variable 439 440 Parameters 441 ---------- 442 size : int or tuple 443 the number and shape of random variables to draw. 444 445 Returns 446 ------- 447 rvs : ndarray 448 the returned random variables with shape given by size and the 449 dimension of the multivariate random vector as additional last 450 dimension 451 452 453 ''' 454 raise NotImplementedError 455 456 def logpdf(self, x): 457 '''logarithm of probability density function 458 459 Parameters 460 ---------- 461 x : array_like 462 can be 1d or 2d, if 2d, then each row is taken as independent 463 multivariate random vector 464 465 Returns 466 ------- 467 logpdf : float or array 468 probability density value of each random vector 469 470 471 this should be made to work with 2d x, 472 with multivariate normal vector in each row and iid across rows 473 does not work now because of dot in whiten 474 475 ''' 476 477 478 raise NotImplementedError 479 480 def cdf(self, x, **kwds): 481 '''cumulative distribution function 482 483 Parameters 484 ---------- 485 x : array_like 486 can be 1d or 2d, if 2d, then each row is taken as independent 487 multivariate random vector 488 kwds : dict 489 contains options for the numerical calculation of the cdf 490 491 Returns 492 ------- 493 cdf : float or array 494 probability density value of each random vector 495 496 ''' 497 raise NotImplementedError 498 499 500 def affine_transformed(self, shift, scale_matrix): 501 '''affine transformation define in subclass because of distribution 502 specific restrictions''' 503 #implemented in subclass at least for now 504 raise NotImplementedError 505 506 def whiten(self, x): 507 """ 508 whiten the data by linear transformation 509 510 Parameters 511 ---------- 512 x : array_like, 1d or 2d 513 Data to be whitened, if 2d then each row contains an independent 514 sample of the multivariate random vector 515 516 Returns 517 ------- 518 np.dot(x, self.cholsigmainv.T) 519 520 Notes 521 ----- 522 This only does rescaling, it does not subtract the mean, use standardize 523 for this instead 524 525 See Also 526 -------- 527 standardize : subtract mean and rescale to standardized random variable. 528 """ 529 x = np.asarray(x) 530 return np.dot(x, self.cholsigmainv.T) 531 532 def pdf(self, x): 533 '''probability density function 534 535 Parameters 536 ---------- 537 x : array_like 538 can be 1d or 2d, if 2d, then each row is taken as independent 539 multivariate random vector 540 541 Returns 542 ------- 543 pdf : float or array 544 probability density value of each random vector 545 546 ''' 547 return np.exp(self.logpdf(x)) 548 549 def standardize(self, x): 550 '''standardize the random variable, i.e. subtract mean and whiten 551 552 Parameters 553 ---------- 554 x : array_like, 1d or 2d 555 Data to be whitened, if 2d then each row contains an independent 556 sample of the multivariate random vector 557 558 Returns 559 ------- 560 np.dot(x - self.mean, self.cholsigmainv.T) 561 562 Notes 563 ----- 564 565 566 See Also 567 -------- 568 whiten : rescale random variable, standardize without subtracting mean. 569 570 571 ''' 572 return self.whiten(x - self.mean) 573 574 def standardized(self): 575 '''return new standardized MVNormal instance 576 ''' 577 return self.affine_transformed(-self.mean, self.cholsigmainv) 578 579 580 def normalize(self, x): 581 '''normalize the random variable, i.e. subtract mean and rescale 582 583 The distribution will have zero mean and sigma equal to correlation 584 585 Parameters 586 ---------- 587 x : array_like, 1d or 2d 588 Data to be whitened, if 2d then each row contains an independent 589 sample of the multivariate random vector 590 591 Returns 592 ------- 593 (x - self.mean)/std_sigma 594 595 Notes 596 ----- 597 598 599 See Also 600 -------- 601 whiten : rescale random variable, standardize without subtracting mean. 602 603 604 ''' 605 std_ = np.atleast_2d(self.std_sigma) 606 return (x - self.mean)/std_ #/std_.T 607 608 def normalized(self, demeaned=True): 609 '''return a normalized distribution where sigma=corr 610 611 if demeaned is True, then mean will be set to zero 612 613 ''' 614 if demeaned: 615 mean_new = np.zeros_like(self.mean) 616 else: 617 mean_new = self.mean / self.std_sigma 618 sigma_new = self.corr 619 args = [getattr(self, ea) for ea in self.extra_args] 620 return self.__class__(mean_new, sigma_new, *args) 621 622 def normalized2(self, demeaned=True): 623 '''return a normalized distribution where sigma=corr 624 625 626 627 second implementation for testing affine transformation 628 ''' 629 if demeaned: 630 shift = -self.mean 631 else: 632 shift = self.mean * (1. / self.std_sigma - 1.) 633 return self.affine_transformed(shift, np.diag(1. / self.std_sigma)) 634 #the following "standardizes" cov instead 635 #return self.affine_transformed(shift, self.cholsigmainv) 636 637 638 639 @property 640 def std(self): 641 '''standard deviation, square root of diagonal elements of cov 642 ''' 643 return np.sqrt(np.diag(self.cov)) 644 645 @property 646 def std_sigma(self): 647 '''standard deviation, square root of diagonal elements of sigma 648 ''' 649 return np.sqrt(np.diag(self.sigma)) 650 651 652 @property 653 def corr(self): 654 '''correlation matrix''' 655 return self.cov / np.outer(self.std, self.std) 656 657 expect_mc = expect_mc 658 659 def marginal(self, indices): 660 '''return marginal distribution for variables given by indices 661 662 this should be correct for normal and t distribution 663 664 Parameters 665 ---------- 666 indices : array_like, int 667 list of indices of variables in the marginal distribution 668 669 Returns 670 ------- 671 mvdist : instance 672 new instance of the same multivariate distribution class that 673 contains the marginal distribution of the variables given in 674 indices 675 676 ''' 677 indices = np.asarray(indices) 678 mean_new = self.mean[indices] 679 sigma_new = self.sigma[indices[:,None], indices] 680 args = [getattr(self, ea) for ea in self.extra_args] 681 return self.__class__(mean_new, sigma_new, *args) 682 683 684#parts taken from linear_model, but heavy adjustments 685class MVNormal0(object): 686 '''Class for Multivariate Normal Distribution 687 688 original full version, kept for testing, new version inherits from 689 MVElliptical 690 691 uses Cholesky decomposition of covariance matrix for the transformation 692 of the data 693 694 ''' 695 696 697 def __init__(self, mean, cov): 698 self.mean = mean 699 self.cov = cov = np.asarray(cov) 700 cov = np.squeeze(cov) 701 self.nvars = nvars = len(mean) 702 703 704 #in the following cov is original, self.cov is full matrix 705 if cov.shape == (): 706 #iid 707 self.cov = np.eye(nvars) * cov 708 self.covinv = np.eye(nvars) / cov 709 self.cholcovinv = np.eye(nvars) / np.sqrt(cov) 710 elif (cov.ndim == 1) and (len(cov) == nvars): 711 #independent heteroskedastic 712 self.cov = np.diag(cov) 713 self.covinv = np.diag(1. / cov) 714 self.cholcovinv = np.diag( 1. / np.sqrt(cov)) 715 elif cov.shape == (nvars, nvars): #python tuple comparison 716 #general 717 self.covinv = np.linalg.pinv(cov) 718 self.cholcovinv = np.linalg.cholesky(self.covinv).T 719 else: 720 raise ValueError('cov has invalid shape') 721 722 #store logdetcov for logpdf 723 self.logdetcov = np.log(np.linalg.det(self.cov)) 724 725 def whiten(self, x): 726 """ 727 whiten the data by linear transformation 728 729 Parameters 730 ---------- 731 X : array_like, 1d or 2d 732 Data to be whitened, if 2d then each row contains an independent 733 sample of the multivariate random vector 734 735 Returns 736 ------- 737 np.dot(x, self.cholcovinv.T) 738 739 Notes 740 ----- 741 This only does rescaling, it does not subtract the mean, use standardize 742 for this instead 743 744 See Also 745 -------- 746 standardize : subtract mean and rescale to standardized random variable. 747 """ 748 x = np.asarray(x) 749 if np.any(self.cov): 750 #return np.dot(self.cholcovinv, x) 751 return np.dot(x, self.cholcovinv.T) 752 else: 753 return x 754 755 def rvs(self, size=1): 756 '''random variable 757 758 Parameters 759 ---------- 760 size : int or tuple 761 the number and shape of random variables to draw. 762 763 Returns 764 ------- 765 rvs : ndarray 766 the returned random variables with shape given by size and the 767 dimension of the multivariate random vector as additional last 768 dimension 769 770 Notes 771 ----- 772 uses numpy.random.multivariate_normal directly 773 774 ''' 775 return np.random.multivariate_normal(self.mean, self.cov, size=size) 776 777 def pdf(self, x): 778 '''probability density function 779 780 Parameters 781 ---------- 782 x : array_like 783 can be 1d or 2d, if 2d, then each row is taken as independent 784 multivariate random vector 785 786 Returns 787 ------- 788 pdf : float or array 789 probability density value of each random vector 790 791 ''' 792 793 return np.exp(self.logpdf(x)) 794 795 def logpdf(self, x): 796 '''logarithm of probability density function 797 798 Parameters 799 ---------- 800 x : array_like 801 can be 1d or 2d, if 2d, then each row is taken as independent 802 multivariate random vector 803 804 Returns 805 ------- 806 logpdf : float or array 807 probability density value of each random vector 808 809 810 this should be made to work with 2d x, 811 with multivariate normal vector in each row and iid across rows 812 does not work now because of dot in whiten 813 814 ''' 815 x = np.asarray(x) 816 x_whitened = self.whiten(x - self.mean) 817 SSR = np.sum(x_whitened**2, -1) 818 llf = -SSR 819 llf -= self.nvars * np.log(2. * np.pi) 820 llf -= self.logdetcov 821 llf *= 0.5 822 return llf 823 824 expect_mc = expect_mc 825 826 827class MVNormal(MVElliptical): 828 '''Class for Multivariate Normal Distribution 829 830 uses Cholesky decomposition of covariance matrix for the transformation 831 of the data 832 833 ''' 834 __name__ == 'Multivariate Normal Distribution' 835 836 837 def rvs(self, size=1): 838 '''random variable 839 840 Parameters 841 ---------- 842 size : int or tuple 843 the number and shape of random variables to draw. 844 845 Returns 846 ------- 847 rvs : ndarray 848 the returned random variables with shape given by size and the 849 dimension of the multivariate random vector as additional last 850 dimension 851 852 Notes 853 ----- 854 uses numpy.random.multivariate_normal directly 855 856 ''' 857 return np.random.multivariate_normal(self.mean, self.sigma, size=size) 858 859 def logpdf(self, x): 860 '''logarithm of probability density function 861 862 Parameters 863 ---------- 864 x : array_like 865 can be 1d or 2d, if 2d, then each row is taken as independent 866 multivariate random vector 867 868 Returns 869 ------- 870 logpdf : float or array 871 probability density value of each random vector 872 873 874 this should be made to work with 2d x, 875 with multivariate normal vector in each row and iid across rows 876 does not work now because of dot in whiten 877 878 ''' 879 x = np.asarray(x) 880 x_whitened = self.whiten(x - self.mean) 881 SSR = np.sum(x_whitened**2, -1) 882 llf = -SSR 883 llf -= self.nvars * np.log(2. * np.pi) 884 llf -= self.logdetsigma 885 llf *= 0.5 886 return llf 887 888 def cdf(self, x, **kwds): 889 '''cumulative distribution function 890 891 Parameters 892 ---------- 893 x : array_like 894 can be 1d or 2d, if 2d, then each row is taken as independent 895 multivariate random vector 896 kwds : dict 897 contains options for the numerical calculation of the cdf 898 899 Returns 900 ------- 901 cdf : float or array 902 probability density value of each random vector 903 904 ''' 905 #lower = -np.inf * np.ones_like(x) 906 #return mvstdnormcdf(lower, self.standardize(x), self.corr, **kwds) 907 return mvnormcdf(x, self.mean, self.cov, **kwds) 908 909 @property 910 def cov(self): 911 '''covariance matrix''' 912 return self.sigma 913 914 def affine_transformed(self, shift, scale_matrix): 915 '''return distribution of an affine transform 916 917 for full rank scale_matrix only 918 919 Parameters 920 ---------- 921 shift : array_like 922 shift of mean 923 scale_matrix : array_like 924 linear transformation matrix 925 926 Returns 927 ------- 928 mvt : instance of MVNormal 929 instance of multivariate normal distribution given by affine 930 transformation 931 932 Notes 933 ----- 934 the affine transformation is defined by 935 y = a + B x 936 937 where a is shift, 938 B is a scale matrix for the linear transformation 939 940 Notes 941 ----- 942 This should also work to select marginal distributions, but not 943 tested for this case yet. 944 945 currently only tested because it's called by standardized 946 947 ''' 948 B = scale_matrix #tmp variable 949 mean_new = np.dot(B, self.mean) + shift 950 sigma_new = np.dot(np.dot(B, self.sigma), B.T) 951 return MVNormal(mean_new, sigma_new) 952 953 def conditional(self, indices, values): 954 r'''return conditional distribution 955 956 indices are the variables to keep, the complement is the conditioning 957 set 958 values are the values of the conditioning variables 959 960 \bar{\mu} = \mu_1 + \Sigma_{12} \Sigma_{22}^{-1} \left( a - \mu_2 \right) 961 962 and covariance matrix 963 964 \overline{\Sigma} = \Sigma_{11} - \Sigma_{12} \Sigma_{22}^{-1} \Sigma_{21}.T 965 966 Parameters 967 ---------- 968 indices : array_like, int 969 list of indices of variables in the marginal distribution 970 given : array_like 971 values of the conditioning variables 972 973 Returns 974 ------- 975 mvn : instance of MVNormal 976 new instance of the MVNormal class that contains the conditional 977 distribution of the variables given in indices for given 978 values of the excluded variables. 979 980 981 ''' 982 #indices need to be nd arrays for broadcasting 983 keep = np.asarray(indices) 984 given = np.asarray([i for i in range(self.nvars) if i not in keep]) 985 sigmakk = self.sigma[keep[:, None], keep] 986 sigmagg = self.sigma[given[:, None], given] 987 sigmakg = self.sigma[keep[:, None], given] 988 sigmagk = self.sigma[given[:, None], keep] 989 990 991 sigma_new = sigmakk - np.dot(sigmakg, np.linalg.solve(sigmagg, sigmagk)) 992 mean_new = self.mean[keep] + \ 993 np.dot(sigmakg, np.linalg.solve(sigmagg, values-self.mean[given])) 994 995# #or 996# sig = np.linalg.solve(sigmagg, sigmagk).T 997# mean_new = self.mean[keep] + np.dot(sigmakg, values-self.mean[given]) 998# sigma_new = sigmakk - np.dot(sigmakg, sig) 999 return MVNormal(mean_new, sigma_new) 1000 1001 1002#redefine some shortcuts 1003np_log = np.log 1004np_pi = np.pi 1005sps_gamln = special.gammaln 1006 1007class MVT(MVElliptical): 1008 1009 __name__ == 'Multivariate Student T Distribution' 1010 1011 def __init__(self, mean, sigma, df): 1012 '''initialize instance 1013 1014 Parameters 1015 ---------- 1016 mean : array_like 1017 parameter mu (might be renamed), for symmetric distributions this 1018 is the mean 1019 sigma : array_like, 2d 1020 dispersion matrix, covariance matrix in normal distribution, but 1021 only proportional to covariance matrix in t distribution 1022 args : list 1023 distribution specific arguments, e.g. df for t distribution 1024 kwds : dict 1025 currently not used 1026 1027 ''' 1028 super(MVT, self).__init__(mean, sigma) 1029 self.extra_args = ['df'] #overwrites extra_args of super 1030 self.df = df 1031 1032 def rvs(self, size=1): 1033 '''random variables with Student T distribution 1034 1035 Parameters 1036 ---------- 1037 size : int or tuple 1038 the number and shape of random variables to draw. 1039 1040 Returns 1041 ------- 1042 rvs : ndarray 1043 the returned random variables with shape given by size and the 1044 dimension of the multivariate random vector as additional last 1045 dimension 1046 - TODO: Not sure if this works for size tuples with len>1. 1047 1048 Notes 1049 ----- 1050 generated as a chi-square mixture of multivariate normal random 1051 variables. 1052 does this require df>2 ? 1053 1054 1055 ''' 1056 from .multivariate import multivariate_t_rvs 1057 return multivariate_t_rvs(self.mean, self.sigma, df=self.df, n=size) 1058 1059 1060 def logpdf(self, x): 1061 '''logarithm of probability density function 1062 1063 Parameters 1064 ---------- 1065 x : array_like 1066 can be 1d or 2d, if 2d, then each row is taken as independent 1067 multivariate random vector 1068 1069 Returns 1070 ------- 1071 logpdf : float or array 1072 probability density value of each random vector 1073 1074 ''' 1075 1076 x = np.asarray(x) 1077 1078 df = self.df 1079 nvars = self.nvars 1080 1081 x_whitened = self.whiten(x - self.mean) #should be float 1082 1083 llf = - nvars * np_log(df * np_pi) 1084 llf -= self.logdetsigma 1085 llf -= (df + nvars) * np_log(1 + np.sum(x_whitened**2,-1) / df) 1086 llf *= 0.5 1087 llf += sps_gamln((df + nvars) / 2.) - sps_gamln(df / 2.) 1088 1089 return llf 1090 1091 def cdf(self, x, **kwds): 1092 '''cumulative distribution function 1093 1094 Parameters 1095 ---------- 1096 x : array_like 1097 can be 1d or 2d, if 2d, then each row is taken as independent 1098 multivariate random vector 1099 kwds : dict 1100 contains options for the numerical calculation of the cdf 1101 1102 Returns 1103 ------- 1104 cdf : float or array 1105 probability density value of each random vector 1106 1107 ''' 1108 lower = -np.inf * np.ones_like(x) 1109 #std_sigma = np.sqrt(np.diag(self.sigma)) 1110 upper = (x - self.mean)/self.std_sigma 1111 return mvstdtprob(lower, upper, self.corr, self.df, **kwds) 1112 #mvstdtcdf does not exist yet 1113 #return mvstdtcdf(lower, x, self.corr, df, **kwds) 1114 1115 @property 1116 def cov(self): 1117 '''covariance matrix 1118 1119 The covariance matrix for the t distribution does not exist for df<=2, 1120 and is equal to sigma * df/(df-2) for df>2 1121 1122 ''' 1123 if self.df <= 2: 1124 return np.nan * np.ones_like(self.sigma) 1125 else: 1126 return self.df / (self.df - 2.) * self.sigma 1127 1128 def affine_transformed(self, shift, scale_matrix): 1129 '''return distribution of a full rank affine transform 1130 1131 for full rank scale_matrix only 1132 1133 Parameters 1134 ---------- 1135 shift : array_like 1136 shift of mean 1137 scale_matrix : array_like 1138 linear transformation matrix 1139 1140 Returns 1141 ------- 1142 mvt : instance of MVT 1143 instance of multivariate t distribution given by affine 1144 transformation 1145 1146 1147 Notes 1148 ----- 1149 1150 This checks for eigvals<=0, so there are possible problems for cases 1151 with positive eigenvalues close to zero. 1152 1153 see: http://www.statlect.com/mcdstu1.htm 1154 1155 I'm not sure about general case, non-full rank transformation are not 1156 multivariate t distributed. 1157 1158 y = a + B x 1159 1160 where a is shift, 1161 B is full rank scale matrix with same dimension as sigma 1162 1163 ''' 1164 #full rank method could also be in elliptical and called with super 1165 #after the rank check 1166 B = scale_matrix #tmp variable as shorthand 1167 if not B.shape == (self.nvars, self.nvars): 1168 if (np.linalg.eigvals(B) <= 0).any(): 1169 raise ValueError('affine transform has to be full rank') 1170 1171 mean_new = np.dot(B, self.mean) + shift 1172 sigma_new = np.dot(np.dot(B, self.sigma), B.T) 1173 return MVT(mean_new, sigma_new, self.df) 1174 1175 1176def quad2d(func=lambda x: 1, lower=(-10,-10), upper=(10,10)): 1177 def fun(x, y): 1178 x = np.column_stack((x,y)) 1179 return func(x) 1180 from scipy.integrate import dblquad 1181 return dblquad(fun, lower[0], upper[0], lambda y: lower[1], 1182 lambda y: upper[1]) 1183 1184if __name__ == '__main__': 1185 1186 from numpy.testing import assert_almost_equal, assert_array_almost_equal 1187 1188 examples = ['mvn'] 1189 1190 mu = (0,0) 1191 covx = np.array([[1.0, 0.5], [0.5, 1.0]]) 1192 mu3 = [-1, 0., 2.] 1193 cov3 = np.array([[ 1. , 0.5 , 0.75], 1194 [ 0.5 , 1.5 , 0.6 ], 1195 [ 0.75, 0.6 , 2. ]]) 1196 1197 1198 if 'mvn' in examples: 1199 bvn = BivariateNormal(mu, covx) 1200 rvs = bvn.rvs(size=1000) 1201 print(rvs.mean(0)) 1202 print(np.cov(rvs, rowvar=0)) 1203 print(bvn.expect()) 1204 print(bvn.cdf([0,0])) 1205 bvn1 = BivariateNormal(mu, np.eye(2)) 1206 bvn2 = BivariateNormal(mu, 4*np.eye(2)) 1207 fun = lambda x : np.log(bvn1.pdf(x)) - np.log(bvn.pdf(x)) 1208 print(bvn1.expect(fun)) 1209 print(bvn1.kl(bvn2), bvn1.kl_mc(bvn2)) 1210 print(bvn2.kl(bvn1), bvn2.kl_mc(bvn1)) 1211 print(bvn1.kl(bvn), bvn1.kl_mc(bvn)) 1212 mvn = MVNormal(mu, covx) 1213 mvn.pdf([0,0]) 1214 mvn.pdf(np.zeros((2,2))) 1215 #np.dot(mvn.cholcovinv.T, mvn.cholcovinv) - mvn.covinv 1216 1217 cov3 = np.array([[ 1. , 0.5 , 0.75], 1218 [ 0.5 , 1.5 , 0.6 ], 1219 [ 0.75, 0.6 , 2. ]]) 1220 mu3 = [-1, 0., 2.] 1221 mvn3 = MVNormal(mu3, cov3) 1222 mvn3.pdf((0., 2., 3.)) 1223 mvn3.logpdf((0., 2., 3.)) 1224 #comparisons with R mvtnorm::dmvnorm 1225 #decimal=14 1226# mvn3.logpdf(cov3) - [-7.667977543898155, -6.917977543898155, -5.167977543898155] 1227# #decimal 18 1228# mvn3.pdf(cov3) - [0.000467562492721686, 0.000989829804859273, 0.005696077243833402] 1229# #cheating new mean, same cov 1230# mvn3.mean = np.array([0,0,0]) 1231# #decimal= 16 1232# mvn3.pdf(cov3) - [0.02914269740502042, 0.02269635555984291, 0.01767593948287269] 1233 1234 #as asserts 1235 r_val = [-7.667977543898155, -6.917977543898155, -5.167977543898155] 1236 assert_array_almost_equal( mvn3.logpdf(cov3), r_val, decimal = 14) 1237 #decimal 18 1238 r_val = [0.000467562492721686, 0.000989829804859273, 0.005696077243833402] 1239 assert_array_almost_equal( mvn3.pdf(cov3), r_val, decimal = 17) 1240 #cheating new mean, same cov, too dangerous, got wrong instance in tests 1241 #mvn3.mean = np.array([0,0,0]) 1242 mvn3c = MVNormal(np.array([0,0,0]), cov3) 1243 r_val = [0.02914269740502042, 0.02269635555984291, 0.01767593948287269] 1244 assert_array_almost_equal( mvn3c.pdf(cov3), r_val, decimal = 16) 1245 1246 mvn3b = MVNormal((0,0,0), 1) 1247 fun = lambda x : np.log(mvn3.pdf(x)) - np.log(mvn3b.pdf(x)) 1248 print(mvn3.expect_mc(fun)) 1249 print(mvn3.expect_mc(fun, size=200000)) 1250 1251 1252 mvt = MVT((0,0), 1, 5) 1253 assert_almost_equal(mvt.logpdf(np.array([0.,0.])), -1.837877066409345, 1254 decimal=15) 1255 assert_almost_equal(mvt.pdf(np.array([0.,0.])), 0.1591549430918953, 1256 decimal=15) 1257 1258 mvt.logpdf(np.array([1.,1.]))-(-3.01552989458359) 1259 1260 mvt1 = MVT((0,0), 1, 1) 1261 mvt1.logpdf(np.array([1.,1.]))-(-3.48579549941151) #decimal=16 1262 1263 rvs = mvt.rvs(100000) 1264 assert_almost_equal(np.cov(rvs, rowvar=0), mvt.cov, decimal=1) 1265 1266 mvt31 = MVT(mu3, cov3, 1) 1267 assert_almost_equal(mvt31.pdf(cov3), 1268 [0.0007276818698165781, 0.0009980625182293658, 0.0027661422056214652], 1269 decimal=18) 1270 1271 mvt = MVT(mu3, cov3, 3) 1272 assert_almost_equal(mvt.pdf(cov3), 1273 [0.000863777424247410, 0.001277510788307594, 0.004156314279452241], 1274 decimal=17) 1275