[frames] | no frames]

# 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.
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
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Thu Sep 22 04:25:14 2011 http://epydoc.sourceforge.net