Pandas expanding/rolling window Korrelationsberechnung mit p-value

Wir empfehlen Ihnen, diese Lösung in einer kontrollierten Umgebung zu testen, bevor Sie sie in die Produktion schicken, Grüße.

Lösung:

Mir ist kein gescheiter Weg eingefallen, dies in Pandas mit rolling direkt, aber beachten Sie, dass Sie den p-Wert anhand des Korrelationskoeffizienten berechnen können.

Der Pearson-Korrelationskoeffizient folgt der Student'schen t-Verteilung und man kann den p-Wert erhalten, indem man ihn in die cdf einfügt, die durch die unvollständige Beta-Funktion definiert ist, scipy.special.betainc. Das hört sich kompliziert an, lässt sich aber mit ein paar Zeilen Code bewerkstelligen. Im Folgenden finden Sie eine Funktion, die den p-Wert anhand des Korrelationskoeffizienten berechnet corr und des Stichprobenumfangs n. Sie basiert eigentlich auf der Implementierung von scipy, die du benutzt hast.

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

Sie können diese Funktion dann auf die Korrelationswerte anwenden, die Sie bereits haben.

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

Es ist vielleicht nicht die perfekte vektorisierte Numpy-Lösung, aber es sollte zehnmal schneller sein, als die Korrelationen immer wieder zu berechnen.

Ansatz 1

corr2_coeff_rowwise listet auf, wie man die elementweise Korrelation zwischen Zeilen durchführt. Wir könnten es auf einen Fall für elementweise Korrelation zwischen zwei Spalten herunterbrechen. Wir würden also mit einer Schleife enden, die Folgendes verwendet corr2_coeff_rowwise. Dann würden wir versuchen, sie zu vektorisieren und sehen, dass es Teile in ihr gibt, die vektorisiert werden können:

1. Durchschnittswerte erhalten mit mean. Dies könnte mit Hilfe eines einheitlichen Filters vektorisiert werden.
2. Als nächstes wurden die Differenzen zwischen diesen Durchschnittswerten und den gleitenden Elementen der Eingabefelder ermittelt. Zur Portierung auf ein vektorisiertes Array würden wir folgendes verwenden broadcasting.

Der Rest bleibt gleich, um die erste der beiden Ausgaben von pearsonr.

Um die zweite Ausgabe zu erhalten, gehen wir zurück zu source code. Dies sollte angesichts des ersten Koeffizientenausgangs einfach sein.

Mit diesen Überlegungen würden wir also etwa so enden.

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

Um also die endgültige Ausgabe aus den Eingaben der Pandas-Reihe zu erhalten -

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

Ansatz #2

Sehr ähnlich wie Approach #1, aber wir verwenden numba um die Speichereffizienz zu verbessern und Schritt 2 des vorherigen Ansatzes zu ersetzen.

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

Ansatz Nr. 3

Sehr ähnlich wie der vorherige Ansatz, außer dass wir die gesamte Koeffizientenarbeit nach 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

---

Zeitangaben.

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

Anmerkung:

• sliding_corr1 scheint bei diesem Datensatz viel Zeit in Anspruch zu nehmen, was höchstwahrscheinlich an der Speicheranforderung von Schritt 2 liegt.

• Der Engpass nach der Verwendung der numba-Funktionen geht dann auf die p-val-Berechnung mit special.btdtr.

Hier sind die Kommentare und Bewertungen

Wir wissen es zu schätzen, dass Sie unsere Forschung unterstützen möchten, indem Sie einen Kommentar hinzufügen oder sie bewerten. Wir heißen Sie willkommen.

---

---

Nutzen Sie unsere Suchmaschine