Source Code for Module gmisclib.multivariate_normal

 1  import Num 
 2  import types 
 3   
 4 -class multivariate_normal: 
 5          __doc__ = """multivariate_normal(mean, cov) or multivariate_normal(mean, cov, [m, n, ...]) 
 6            returns an array containing multivariate normally distributed random numbers 
 7            with specified mean and covariance. 
 8   
 9            mean must be a 1 dimensional array. cov must be a square two dimensional 
10            array with the same number of rows and columns as mean has elements. 
11   
12            The first form returns a single 1-D array containing a multivariate 
13            normal. 
14   
15            The second form returns an array of shape (m, n, ..., cov.shape[0]). 
16            In this case, output[i,j,...,:] is a 1-D array containing a multivariate 
17            normal.""" 
18   
19 -        def __init__(self, mean, cov): 
20                  cov = Num.asarray(cov) 
21                  mean = Num.array(mean, copy=True) 
22                  # Check preconditions on arguments 
23                  if len(mean.shape) != 1: 
24                          raise TypeError, "mean must be 1 dimensional." 
25                  if (len(cov.shape) != 2) or (cov.shape[0] != cov.shape[1]): 
26                          raise TypeError, "cov must be 2 dimensional and square." 
27                  if mean.shape[0] != cov.shape[0]: 
28                          raise TypeError, "mean and cov must have same length." 
29                  u, s, self.v = Num.LA.singular_value_decomposition(cov) 
30                  self.ss = Num.sqrt(s) 
31                  self.mean = Num.array(mean, copy=True)  # Copy, so no one changes it 
32                                                  # underneath you. 
33   
34   
35 -        def sample(self, shape=[]): 
36                  # Compute shape of output 
37                  if isinstance(shape, types.IntType): 
38                          final_shape = [shape] 
39                  else: 
40                          final_shape = list(shape[:]) 
41                  final_shape.append(self.mean.shape[0]) 
42                  # Create a matrix of independent standard normally distributed random 
43                  # numbers. The matrix has rows with the same length as mean and as 
44                  # many rows are necessary to form a matrix of shape final_shape. 
45                  x = Num.RA.standard_normal(Num.multiply.reduce(final_shape)) 
46                  x.shape = (Num.multiply.reduce(final_shape[0:-1]), final_shape[-1]) 
47                  # Transform matrix of standard normals into matrix where each row 
48                  # contains multivariate normals with the desired covariance. 
49                  # Compute A such that matrixmultiply(transpose(A),A) == cov. 
50                  # Then the matrix products of the rows of x and A has the desired 
51                  # covariance. Note that sqrt(s)*v where (u,s,v) is the singular value 
52                  # decomposition of cov is such an A. 
53                  x = Num.matrixmultiply(x*self.ss, self.v) 
54                  # The rows of x now have the correct covariance but mean 0. Add 
55                  # mean to each row. Then each row will have mean mean. 
56                  Num.add(self.mean, x, x) 
57                  x.shape = final_shape 
58                  return x 
59   
60 -def test(): 
61          import sys 
62          N = 30000 
63          cov = Num.array([[1, 0.1, -0.2], [0.1, 1.1, 0.2], [-0.2, 0.2, 2.0]]) 
64          x = multivariate_normal(Num.zeros((3,), Num.Float), 
65                                  cov) 
66          s = Num.zeros((3,3), Num.Float) 
67          for i in range(N): 
68                  z = x.sample() 
69                  if i%300 == 0: 
70                          sys.stdout.writelines('.') 
71                  s += Num.outerproduct(z, z) 
72          print 
73          print cov 
74          print s/N 
75   
76  if __name__ == '__main__': 
77          test() 
78