Cross-covariance function
Greetings, I recently needed to calculate the cross-covariance of two random vectors, (e.g. I have two matricies, X and Y, the columns of which are observations of one variable, and I wish to generate a matrix pairing each value of X and Y) and so I wrote a small utility function to do so, and I'd like to try and get it integrated into numpy core, if it is deemed useful. I have never submitted a patch to numpy before, so I'm not sure as to the protocol; do I ask someone on this list to review the code? Are there conventions I should be aware of? Etc... Thank you all, -E
Hi Eliot, Le 19/01/2012 07:50, Elliot Saba a écrit :
I recently needed to calculate the cross-covariance of two random vectors, (e.g. I have two matricies, X and Y, the columns of which are observations of one variable, and I wish to generate a matrix pairing each value of X and Y) I don't see how does your function relates to numpy.cov [1]. Is it an "extended case" function or is there a difference in the underlying math ?
Best, Pierre [1] numpy.cov docstring : http://docs.scipy.org/doc/numpy/reference/generated/numpy.cov.html
Den 20.01.2012 13:39, skrev Pierre Haessig:
I don't see how does your function relates to numpy.cov [1]. Is it an "extended case" function or is there a difference in the underlying math ?
If X is rank n x p, then np.cov(X, rowvar=False) is equal to S after cX = X - X.mean(axis=0)[np.newaxis,:] S = np.dot(cX.T, cX)/(n-1.) If we also have Y rank n x p, then the upper p x p quadrant of np.cov(X, y=Y, rowvar=False) is equal to S after XY = np.hstack(X,Y) cXY = XY - XY.mean(axis=0)[np.newaxis,:] S = np.dot(cXY.T, cXY)/(n-1.) Thus we can see thatthe total cocariance is composed of four parts: S[:p,:p] = np.dot(cX.T, cX)/(n-1.) # upper left S[:p,p:] = np.dot(cXY.T, cYY)/(n-1.) # upper right S[p:,:p] = np.dot(cY.T, cX)/(n-1.) # lower left S[p:,:p] = np.dot(cYX.T, cYX)/(n-1.) # lower right Often we just want the upper-right p x p quadrant. Thus we can save 75 % of the cpu time by not computing the rest. Sturla
Le 20/01/2012 16:30, Sturla Molden a écrit :
Often we just want the upper-right p x p quadrant. Thanks for the explanation. If I understood it correctly, you're interested in the *cross*-covariance block of the matrix (and now I understand better Elliot's message). Actually, I thought that was the behavior of the np.cor function ! But you're right it's not ! [source code] The seconde 'y' argument just gets concatenated with the first one 'm'.
I would go further and ask why it so. People around may have use cases in mind, because I have not. Otherwise, I feel that the default behavior of cov when called with two arguments should be what Sturla and Elliot just described. Best, Pierre (that is something like this : def cov(X, Y=None): if Y is None: Y = X else: assert Y.shape == X.shape # or something like that # [...jumping to the end of the existing code...] if not rowvar: return (dot(X.T, Y.conj()) / fact).squeeze() else: return (dot(X, Y.T.conj()) / fact).squeeze() ) [source code] https://github.com/numpy/numpy/blob/master/numpy/lib/function_base.py
participants (3)
-
Elliot Saba -
Pierre Haessig -
Sturla Molden