Pandas扩展/滚动窗口相关性计算及p值

Question

假设我有一个DataFrame，我想在其中计算两列之间的滚动或扩展Pearson相关性。

import numpy as np
import pandas as pd
import scipy.stats as st


df = pd.DataFrame({'x': np.random.rand(10000), 'y': np.random.rand(10000)})

利用内置的 pandas 功能，可以快速计算此内容。

expanding_corr = df['x'].expanding(50).corr(df['y'])
rolling_corr = df['x'].rolling(50).corr(df['y'])

然而，如果我希望获得与这些相关性相关的p值，我最好定义一个自定义滚动函数并将apply传递给groupby对象。

def custom_roll(df, w, **kwargs):

    v = df.values
    d0, d1 = v.shape
    s0, s1 = v.strides
    a = np.lib.stride_tricks.as_strided(v, (d0 - (w - 1), w, d1), (s0, s0, s1))
    rolled_df = pd.concat({
        row: pd.DataFrame(values, columns=df.columns)
        for row, values in zip(df.index[(w-1):], a)
    })
    return rolled_df.groupby(level=0, **kwargs)

c_df = custom_roll(df, 50).apply(lambda df: st.pearsonr(df['x'], df['y']))

c_df现在包含适当的相关性以及重要的p值。

然而，与内置的pandas方法相比，这种方法非常慢，这意味着它不适用于实际上我在优化过程中计算这些相关性数千次。此外，我不确定如何扩展custom_roll函数以适用于扩展窗口。

有人可以指导我如何利用numpy以向量化速度获取扩展窗口上的p值吗？

- mch56

Answer 1

我无法想到一种聪明的方式来直接使用rolling 在pandas中完成这个任务，但请注意，您可以根据相关系数计算p值。

皮尔逊相关系数遵循Student's t-distribution，您可以通过将其插入由不完全beta函数定义的累积分布函数来获取p值，scipy.special.betainc。听起来很复杂，但只需要几行代码即可完成。下面是一个函数，它根据相关系数corr 和样本大小n 计算p值。 实际上，它基于您一直在使用的Scipy实现。

import pandas as pd
from scipy.special import betainc

def pvalue(corr, n=50):
    df = n - 2
    t_squared = corr**2 * (df / ((1.0 - corr) * (1.0 + corr)))
    prob = betainc(0.5*df, 0.5, df/(df+t_squared))
    return prob

你可以将这个函数应用于你已经拥有的相关性数值。

rolling_corr = df['x'].rolling(50).corr(df['y'])
pvalue(rolling_corr)

这可能不是完美的矢量化NumPy解决方案，但应该比一遍又一遍地计算相关性快数十倍。

- AlCorreia

---

谢谢。我理解从相关性中回溯p值的统计数据。不过，我真正想知道的答案是如何将一般的统计函数传递到滚动/扩展窗口中，绕过缓慢的“apply”方法。 - mch56

我理解这个问题，确实是一个非常好的问题。似乎你可以使用rolling.apply传递任何函数，但我找不到一种方法将其应用于两列。 - AlCorreia

是的，这是一个有趣的问题，我很感激你尝试解决它。不幸的是，apply 在底层只是一个循环，因此在 DataFrame 变得很大时会遇到性能问题。我希望能够找到一个 numpy 解决方案，或者可能对 .corr 本身进行一些修改。 - mch56

但是“rolling”函数背后也只是一个for循环，对吧？我认为减慢您当前解决方案的原因不是真正的“apply”，而是每次都需要创建一个新的数据框。此外，我更新了我的答案，并假设您拥有相关系数，您实际上不需要“apply”。 - AlCorreia

Answer 2

方法#1

corr2_coeff_rowwise介绍了如何在行之间进行逐元素相关性。我们可以将其简化为两列之间的逐元素相关操作。因此，我们得到一个使用corr2_coeff_rowwise的循环。接下来，我们会尝试将其向量化，并发现其中有一些可以向量化的部分：

1. 使用mean获取平均值。这可以通过使用统一过滤器进行向量化。
2. 接下来是获取这些平均值与输入数组的滑动元素之间的差异。为了将其转换为向量化版本，我们将利用broadcasting。

剩下的步骤与获取pearsonr的第一个输出相同。

要获得第二个输出，我们需要回到源代码。这应该很简单，考虑到第一个系数输出。

因此，考虑到这些，我们最终会得到类似于以下内容的代码 -

import scipy.special as special
from scipy.ndimage import uniform_filter

def sliding_corr1(a,b,W):
    # a,b are input arrays; W is window length

    am = uniform_filter(a.astype(float),W)
    bm = uniform_filter(b.astype(float),W)

    amc = am[W//2:-W//2+1]
    bmc = bm[W//2:-W//2+1]

    da = a[:,None]-amc
    db = b[:,None]-bmc

    # Get sliding mask of valid windows
    m,n = da.shape
    mask1 = np.arange(m)[:,None] >= np.arange(n)
    mask2 = np.arange(m)[:,None] < np.arange(n)+W
    mask = mask1 & mask2
    dam = (da*mask)
    dbm = (db*mask)

    ssAs = np.einsum('ij,ij->j',dam,dam)
    ssBs = np.einsum('ij,ij->j',dbm,dbm)
    D = np.einsum('ij,ij->j',dam,dbm)
    coeff = D/np.sqrt(ssAs*ssBs)

    n = W
    ab = n/2 - 1
    pval = 2*special.btdtr(ab, ab, 0.5*(1 - abs(np.float64(coeff))))
    return coeff,pval

因此，为了从Pandas系列的输入中获取最终输出 -

out = sliding_corr1(df['x'].to_numpy(copy=False),df['y'].to_numpy(copy=False),50)

第二种方法

与"第一种方法"非常相似，但我们将使用numba提高内存效率，以替换前一个方法中的步骤＃2。

from numba import njit
import math

@njit(parallel=True)
def sliding_corr2_coeff(a,b,amc,bmc):
    L = len(a)-W+1
    out00 = np.empty(L)
    for i in range(L):
        out_a = 0
        out_b = 0
        out_D = 0
        for j in range(W):
            d_a = a[i+j]-amc[i]
            d_b = b[i+j]-bmc[i]
            out_D += d_a*d_b
            out_a += d_a**2
            out_b += d_b**2
        out00[i] = out_D/math.sqrt(out_a*out_b)
    return out00

def sliding_corr2(a,b,W):
    am = uniform_filter(a.astype(float),W)
    bm = uniform_filter(b.astype(float),W)

    amc = am[W//2:-W//2+1]
    bmc = bm[W//2:-W//2+1]

    coeff = sliding_corr2_coeff(a,b,amc,bmc)

    ab = W/2 - 1
    pval = 2*special.btdtr(ab, ab, 0.5*(1 - abs(np.float64(coeff))))
    return coeff,pval

方法三

与前一个方法非常相似，只是我们将所有系数工作都推到了numba中-

@njit(parallel=True)
def sliding_corr3_coeff(a,b,W):
    L = len(a)-W+1
    out00 = np.empty(L)
    for i in range(L):
        a_mean = 0.0
        b_mean = 0.0
        for j in range(W):
            a_mean += a[i+j]
            b_mean += b[i+j]
        a_mean /= W
        b_mean /= W

        out_a = 0
        out_b = 0
        out_D = 0
        for j in range(W):
            d_a = a[i+j]-a_mean
            d_b = b[i+j]-b_mean
            out_D += d_a*d_b
            out_a += d_a*d_a
            out_b += d_b*d_b
        out00[i] = out_D/math.sqrt(out_a*out_b)
    return out00

def sliding_corr3(a,b,W):    
    coeff = sliding_corr3_coeff(a,b,W)
    ab = W/2 - 1
    pval = 2*special.btdtr(ab, ab, 0.5*(1 - np.abs(coeff)))
    return coeff,pval

时间 -

In [181]: df = pd.DataFrame({'x': np.random.rand(10000), 'y': np.random.rand(10000)})

In [182]: %timeit sliding_corr2(df['x'].to_numpy(copy=False),df['y'].to_numpy(copy=False),50)
100 loops, best of 3: 5.05 ms per loop

In [183]: %timeit sliding_corr3(df['x'].to_numpy(copy=False),df['y'].to_numpy(copy=False),50)
100 loops, best of 3: 5.51 ms per loop

注意:

• sliding_corr1在这个数据集上花费的时间似乎很长，很可能是由于它的第二步需要占用大量内存。

• 在使用numba函数之后，瓶颈转移到了使用special.btdtr进行p值计算。

- Divakar