Cyt*_*rak 8 python numpy pandas rolling-computation pearson-correlation
让我先说一下，为了重现这个问题，我需要一个大数据，这是问题的一部分，我无法预测什么时候会出现这种特殊性。无论如何，数据太大（~13k 行，2 列）无法粘贴到问题中，我在帖子末尾添加了一个 pastebin 链接。
在过去的几天里，我遇到了一个奇怪的问题pandas.core.window.rolling.Rolling.corr。我有一个数据集，我试图在其中计算滚动相关性。这就是问题：
在计算
window_size=100两列 (aandb)之间的滚动 ( ) 相关性时：一些索引（一个这样的索引是12981）给出接近的0值（顺序1e-10），但理想情况下它应该返回nanorinf，（因为一列中的所有值都是常数）。但是，如果我只是计算与该索引有关的独立相关性（即包括所述索引的最后 100 行数据），或者对较少数量的行（例如 300 或 1000 而不是 13k）执行滚动计算，我得到正确的结果（即nan或inf。）
>>> df = pd.read_csv('sample_corr_data.csv') # link at the end, ## columns = ['a', 'b']
>>> df.a.tail(100).value_counts()
0.000000 86
-0.000029 3
0.000029 3
-0.000029 2
0.000029 2
-0.000029 2
0.000029 2
Name: a, dtype: int64
>>> df.b.tail(100).value_counts() # all 100 values are same
6.0 100
Name: b, dtype: int64
>>> df.a.tail(100).corr(df.b.tail(100))
nan # expected, because column 'b' has same value throughout
# Made sure of this using,
# 1. np.corrcoef, because pandas uses this internally to calculate pearson moments
>>> np.corrcoef(df.a.tail(100), df.b.tail(100))[0, 1]
nan
# 2. using custom function
>>> def pearson(a, b):
n = a.size
num = n*np.nansum(a*b) - np.nansum(a)*np.nansum(b)
den = (n*np.nansum((a**2)) - np.nansum(a)**2)*(n*np.nansum(b**2) - np.nansum(b)**2)
return num/np.sqrt(den) if den * np.isfinite(den*num) else np.nan
>>> pearson(df.a.tail(100), df.b.tail(100))
nan
>>> df.a.rolling(100).corr(df.b).tail(3)
12979 7.761921e-07
12980 5.460717e-07
12981 2.755881e-10 # This should have been NaN/inf !!
## Furthermore!!
>>> debug = df.tail(300)
>>> debug.a.rolling(100).corr(debug.b).tail(3)
12979 7.761921e-07
12980 5.460717e-07
12981 -inf # Got -inf, fine
dtype: float64
>>> debug = df.tail(3000)
>>> debug.a.rolling(100).corr(debug.b).tail(3)
12979 7.761921e-07
12980 5.460717e-07
12981 inf # Got +inf, still acceptable
dtype: float64
这继续直到9369行：
>>> debug = df.tail(9369)
>>> debug.a.rolling(100).corr(debug.b).tail(3)
12979 7.761921e-07
12980 5.460717e-07
12981 inf
dtype: float64
# then
>>> debug = df.tail(9370)
>>> debug.a.rolling(100).corr(debug.b).tail(3)
12979 7.761921e-07
12980 5.460717e-07
12981 4.719615e-10 # SPOOKY ACTION IN DISTANCE!!!
dtype: float64
>>> debug = df.tail(10000)
>>> debug.a.rolling(100).corr(debug.b).tail(3)
12979 7.761921e-07
12980 5.460717e-07
12981 1.198994e-10 # SPOOKY ACTION IN DISTANCE!!!
dtype: float64
>>> df.a.rolling(100).apply(lambda x: x.corr(df.b.reindex(x.index))).tail(3) # PREDICTABLY, VERY SLOW!
12979 7.761921e-07
12980 5.460717e-07
12981 NaN
Name: a, dtype: float64
# again this checks out using other methods,
>>> df.a.rolling(100).apply(lambda x: np.corrcoef(x, df.b.reindex(x.index))[0, 1]).tail(3)
12979 7.761921e-07
12980 5.460717e-07
12981 NaN
Name: a, dtype: float64
>>> df.a.rolling(100).apply(lambda x: pearson(x, df.b.reindex(x.index))).tail(3)
12979 7.761921e-07
12980 5.460717e-07
12981 NaN
Name: a, dtype: float64
据我了解，结果series.rolling(n).corr(other_series)应与以下内容匹配：
>>> def rolling_corr(series, other_series, n=100):
return pd.Series(
[np.nan]*(n-1) + [series[i-n: i].corr(other_series[i-n:i])
for i in range (n, series.size+1)]
)
>>> rolling_corr(df.a, df.b).tail(3)
12979 7.761921e-07
12980 5.460717e-07
12981 NaN
首先，我认为这是一个floating-point arithmetic问题（因为最初，在某些情况下，我可以通过将列 'a' 舍入到小数点后 5 位或强制转换为 来解决此问题float32，但在这种情况下，无论使用的样本数量如何，它都会存在. 因此，根据数据的大小，必须存在一些问题rolling或至少rolling会引起floating-point问题。我检查了 的源代码rolling.corr，但找不到任何可以解释这种不一致的内容。现在我很担心，有多少过去的代码受到这个问题的困扰。
这背后的原因是什么？以及如何解决这个问题？如果发生这种情况，因为可能是熊猫喜欢时速超过精度（如建议在这里），这是否意味着我永远不能可靠地使用pandas.rolling在大样本的操作？我怎么知道这种不一致会出现的大小？
测试在
注意：不同的操作系统在上述索引处返回不同的值，但都是有限的和接近的0。
如果您计算将皮尔逊公式中的总和替换为滚动总和会怎么样
def rolling_pearson(a, b, n):
a_sum = a.rolling(n).sum()
b_sum = b.rolling(n).sum()
ab_sum = (a*b).rolling(n).sum()
aa_sum = (a**2).rolling(n).sum()
bb_sum = (b**2).rolling(n).sum();
num = n * ab_sum - a_sum * b_sum;
den = (n*aa_sum - a_sum**2) * (n * bb_sum - b_sum**2)
return num / den**(0.5)
rolling_pearson(df.a, df.b, 100)
...
12977 1.109077e-06
12978 9.555249e-07
12979 7.761921e-07
12980 5.460717e-07
12981 inf
Length: 12982, dtype: float64
为了回答这个问题，我需要检查实施情况。因为确实 的最后 100 个样本的方差b为零，并且滚动相关性计算为a.cov(b) / (a.var() * b.var())**0.5。
经过一番搜索，我在这里找到了滚动方差实现，他们使用的方法是Welford 的在线算法。该算法很好，因为您可以仅使用一次乘法来添加一个样本（与累积和的方法相同），并且可以使用单个整数除法进行计算。这里用python重写一下。
...
12977 1.109077e-06
12978 9.555249e-07
12979 7.761921e-07
12980 5.460717e-07
12981 inf
Length: 12982, dtype: float64
在 pandas 实现中，他们提到了Kahan 求和，这对于在加法中获得更好的精度很重要，但结果并没有因此而改善（我没有检查它是否正确实现）。
应用 Welford 算法n=100
s = (0,0,0)
for i in range(len(df.b)):
if i >= n:
s = welford_remove(s, df.b[i-n])
s = welford_add(s, df.b[i])
finalize(s)
它给
(6.000000000000152, 4.7853099260919405e-12, 4.8336463899918594e-12)
和df.b.rolling(100).var()给予
0 NaN
1 NaN
2 NaN
3 NaN
4 NaN
...
12977 6.206061e-01
12978 4.703030e-01
12979 3.167677e-01
12980 1.600000e-01
12981 6.487273e-12
Name: b, Length: 12982, dtype: float64
误差6.4e-12略高于4.83e-12直接应用 Welford 方法给出的误差。
另一方面，(df.b**2).rolling(n).sum()-df.b.rolling(n).sum()**2/n最后一个条目的值为 0.0。
0 NaN
1 NaN
2 NaN
3 NaN
4 NaN
...
12977 61.44
12978 46.56
12979 31.36
12980 15.84
12981 0.00
Name: b, Length: 12982, dtype: float64
我希望这个解释令人满意：）