我正在尝试使用 NumPy 1.8.0 和 Python 2.7.3 将 MATLAB/Octave 程序移植到 Python。我用过 this reference作为帮助将 MATLAB 函数转换为 NumPy 方法并取得巨大成功,直到我到达我想要计算两个矩阵之间相关性的地步。
第一个矩阵是 40000x25 float ,第二个矩阵是 40000x1 整数。在 Octave 中,我使用语句 corr(a,b)
并获得一个 25x1 的浮点矩阵。在 NumPy 中尝试相应的方法 (numpy.correlate(a,b)
) 会产生错误:
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/Library/Python/2.7/site-packages/numpy-1.8.0.dev_1a9aa5a_20130415-py2.7-macosx-10.8-intel.egg/numpy/core/numeric.py", line 751, in correlate
return multiarray.correlate2(a,v,mode)
ValueError: object too deep for desired array
如果我更改代码以计算 a
的每一列的相关性,我可以让它工作,如下所示:
for i in range(25):
c2[i] = numpy.correlate(a[:,i], b)
但是,c2
数组中的值与 Octave 的输出不同。 Octave 返回一个由小于 1 的 float 组成的 25x1 矩阵。我从 NumPy 获得的值是 -270 到 900 之间的 float 。
我试图了解这两种算法在幕后做了什么,但失败得很惨。有人可以指出我的逻辑错误吗?
最佳答案
似乎存在一个 numpy.corrcoef
可以根据需要计算相关系数。但是,它的界面与Octave/Matlab corr
不同。
首先,默认情况下,该函数将行视为变量,将列视为观察值。要模仿 Octave/Matlab 的行为,您可以传递一个标志来反转它。
此外,根据 this answer ,numpy.cov
函数(我假设 corrcoef
在内部使用)返回一个 2x2 矩阵,每个矩阵都包含一个特定的协方差:
cov(a,a) cov(a,b)
cov(a,b) cov(b,b)
正如他指出的那样,[0][1]
元素是您想要的 cov(a,b)
。因此,也许这样的事情会奏效:
for i in range(25):
c2[i] = numpy.corrcoef(a[:,i], b, rowvar=0)[0][1]
作为引用,以下是您尝试过的两个功能的一些摘录。似乎他们执行的是完全不同的事情。
Octave :
— Function File: corr (x, y)
Compute matrix of correlation coefficients.
If each row of x and y is an observation and each column is a variable, then the (i, j)-th entry of corr (x, y) is the correlation between the i-th variable in x and the j-th variable in y.
corr (x,y) = cov (x,y) / (std (x) * std (y))
If called with one argument, compute corr (x, x), the correlation between the columns of x.
和 Numpy:
numpy.correlate(a, v, mode='valid', old_behavior=False)[source]
Cross-correlation of two 1-dimensional sequences.
This function computes the correlation as generally defined in signal processing texts:
z[k] = sum_n a[n] * conj(v[n+k])
with a and v sequences being zero-padded where necessary and conj being the conjugate.
关于python - MATLAB/Octave corr 和 Python numpy.correlate 有什么区别?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/16698811/