带有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值的方向吗？

Answer 1

我想不出一种在熊猫中rolling直接使用的聪明方法，但是请注意，给定相关系数，您可以计算p值。

Pearson的相关系数服从Student的t分布，您可以通过将p值插入不完全beta函数定义的cdf来获得p值scipy.special.betainc。听起来很复杂，但是可以用几行代码完成。以下是在给定相关系数corr和样本大小的情况下计算p值的函数n。它实际上是基于您一直在使用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解决方案，但比一遍又一遍地计算相关性要快几十倍。