2

In numpy manual, it is said:

Instead of specifying the full covariance matrix, popular approximations include:  
    Spherical covariance (cov is a multiple of the identity matrix)

Has anybody ever specified spherical covariance? I am trying to make it work to avoid building the full covariance matrix, which is too much memory-consuming.

3
  • What is your problem exactly? are you trying to build a spherical covariance matrix that reasonably approximates your full covariance matrix? In this case, more information is needed about your specific problem. May 16, 2013 at 3:59
  • It is this: there is a 512X512 array of variables, every variable is correlated to its neighbour within the radius of 256, and the correlation coefficient is some kind of inversely proportional to the distance. So I tried to treat this array as an 1-D array with size of 512x512=262144, and build the cm(whose size is 262144x262144) for it, then make a multivariate random sample.
    – hookch
    May 16, 2013 at 16:54
  • A brute force approach is indeed impractical. What the part you quote says is what Robert Kern and Jaime have implemented in their response: it consists in completely forgetting about the correlations. If you can find a way to diagonalize analytically the full covariance matrix, or some approximation of it, you can generate independent random variables and then combine them to recover correlations… May 17, 2013 at 8:16

2 Answers 2

1

If you just have a diagonal covariance matrix, it is usually easier (and more efficient) to just scale standard normal variates yourself instead of using multivariate_normal().

>>> import numpy as np
>>> stdevs = np.array([3.0, 4.0, 5.0])
>>> x = np.random.standard_normal([100, 3])
>>> x.shape
(100, 3)
>>> x *= stdevs
>>> x.std(axis=0)
array([ 3.23973255,  3.40988788,  4.4843039 ])
2
  • Maybe it is not right to call it "spherical", since the matrix is not diagonal. Here is the whole problem:there is a 512X512 array of variables, every variable is correlated to its neighbour within the radius of 256, and the correlation coefficient is some kind of inversely proportional to the distance. So I tried to treat this array as an 1-D array with size of 512x512=262144, and build the cm(whose size is 262144x262144) for it, then make a multivariate random sample.
    – hookch
    May 16, 2013 at 16:51
  • Okay yes, this is very much not spherical in the sense that the documentation was using the term. You have a spherical covariance function in the 2-D geometric space your grid is embedded in, but that does not make a spherical covariance matrix in the 262144-D sampling space. You want to sample from a Gaussian random field. You can try using the Gaussian Process functionality in PyMC, though it won't be very efficient as it doesn't really handle the special case of grids. Googling for "python gaussian random field" might give you some leads. May 17, 2013 at 11:24
1

While @RobertKern's approach is correct, you can let numpy handle all of that for you, as np.random.normal will do broadcasting on multiple means and standard deviations:

>>> np.random.normal(0, [1,2,3])
array([ 0.83227999,  3.40954682, -0.01883329])

To get more than a single random sample, you have to give it an appropriate size:

>>> x = np.random.normal(0, [1, 2, 3], size=(1000, 3))
>>> np.std(x, axis=0)
array([ 1.00034817,  2.07868385,  3.05475583])
2
  • It seems to be working indeed, but I can't find this behavior in the documentation (docs.scipy.org/doc/numpy/reference/generated/…). Is it documented somewhere? May 17, 2013 at 7:56
  • No, it's not documented, but it is reliable. The parameters will be broadcast against each other, and each broadcasted set of parameters will be sampled independently. May 17, 2013 at 12:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.