<p>经过一番思考和几次测试,我想出了以下简单的代码:</p>
<pre><code>def var_covar(arr, means=None, axis=1, normalize=True, bias=None):
"""
Computes the covariance of array arr. Optionally it takes means
separately instead of computing it from the data.
PARAMETERS
:arr: array_like; A NumPy array or other data type that can be converted
into an NumPy array.
:means: Optional, array_like; A vector of means for each observation in
arr. Mostly would be population means.
:axis: Optional, int, default 1; Operates either by row or by column.
Default assumes experiments are along rows, and observations are
along columns.
:normalize: Optional, bool, default True; Normalizes the data by 1/(N-bias),
where N is number of observations.
:bias: Optional, float; Bias is subtracted from N, to account for sample
mean. If means are passed, it is assumed they are population means
and so bias is set to 0. Else to 1. Passing a value overrides the
default value of bias.
RETURNS
:cov: array_likep; A 2-D NumPy array if arr is a matrix. A float if arr
is a vector
NOTES:
Maximum number of supported dimensions is 2
It is assumed that passed means are population means, and not sample
means. This assumption affects bias. Pass a value for bias to override
the default bias.
It is recommended that length of means equals number of observations in
arr. If that is not the case, means will be broadcast, but it may throw
an error.
axis must be either 0 or 1. Higher dimensions are not supported.
"""
arr = np.asarray(arr)
tbias = bias
if arr.ndim > 2:
raise RuntimeError('Cannot handle arr larger than a 2x2 matrix')
if axis not in [0,1]:
raise ValueError('axis must be either 0 or 1')
elif axis == 0:
arr = arr.T
if means is None:
means = np.average(arr, axis=1)
bias = 1
else:
bias = 0
means = np.asarray(means)
arr = np.subtract(arr.T, means).T
cov = np.dot(arr, arr.T)
if normalize:
if tbias is not None:
bias = tbias
factor = 1/(arr.shape[0] - bias)
cov = np.multiply(cov, factor)
return cov
</code></pre>
<p>它比较快,完全矢量化。我把这个分享给其他需要它的人。在</p>
<p>如果你认为你做得比我好,请随时发表答案。我很乐意接受它作为答案。在</p>