计算astropy中window的变异系数

Calculate coefficient of variation of window in astropy

我有一个数组,我想使用 astropy 计算统计数据。我有的是:

from astropy.convolution import convolve
import numpy as np

x = np.random.randint(1, 10, size=(5, 5))
y = convolve(x, np.ones((3, 3)), boundary='extend', preserve_nan=True)

print(x)
print(y)

[[9 1 8 6 5]
 [4 2 1 8 4]
 [2 8 4 6 6]
 [8 4 8 5 6]
 [7 3 1 2 3]]

[[5.33333333 4.77777778 4.55555556 5.66666667 5.33333333]
 [4.55555556 4.33333333 4.88888889 5.33333333 5.55555556]
 [4.66666667 4.55555556 5.11111111 5.33333333 5.66666667]
 [5.44444444 5.         4.55555556 4.55555556 4.77777778]
 [6.         4.66666667 3.22222222 3.44444444 3.66666667]]

y 中的每个元素都是围绕 x 中该位置绘制的 3x3 框的平均值。我想要代替均值的是变异系数(标准差除以均值)。我不确定这是否可以在 astropy 中完成,或者我是否需要使用其他东西,比如

from scipy.stats import variation

Astropy 建立在 numpy 和 scipy 之上。 Numpy 是 low-level 数组库,它实现了 scipy 和 astropy 等更高级库使用的数据存储和基本操作。了解 numpy 数组的工作原理将有助于您使用 astropy。

因为你想对滚动window做统计,类似scipy.stats.variation will not help you directly. It computes statistics on elements grouped by axis, which is not how your elements are grouped. You could simulate a rolling window through axes using something like np.lib.stride_tricks.as_strided,但那很危险,容易出错,而且效率不高。相反,我建议使用类似于您已经完成的卷积方法。

首先,记住方差可以表示为

var = sum(x^2) / N - (sum(x) / N)^2

所以让我们这样做:

kernel = np.ones((3, 3))
mean_x = convolve(x, kernel, boundary='extend', preserve_nan=True)
mean_square_x = convolve(x**2, kernel, boundary='extend', preserve_nan=True)

您现在每个 window 和 sum(x^2) / N 都有 sum(x) / N。标准差就是方差的平方根:

std_x = np.sqrt(mean_square_x - mean_x**2)

你要的结果就是比例

result = std_x / mean_x

您可以更简洁地完成整个计算,如

m = convolve(x, kernel, boundary='extend', preserve_nan=True)
result = np.sqrt(convolve(x**2, kernel, boundary='extend', preserve_nan=True) - m**2) / m