Berechnung der Pearson-Korrelation und Signifikanz in Python

Antworten:

201

Sie können einen Blick darauf werfen scipy.stats:

from pydoc import help
from scipy.stats.stats import pearsonr
help(pearsonr)

>>>
Help on function pearsonr in module scipy.stats.stats:

pearsonr(x, y)
 Calculates a Pearson correlation coefficient and the p-value for testing
 non-correlation.

 The Pearson correlation coefficient measures the linear relationship
 between two datasets. Strictly speaking, Pearson's correlation requires
 that each dataset be normally distributed. Like other correlation
 coefficients, this one varies between -1 and +1 with 0 implying no
 correlation. Correlations of -1 or +1 imply an exact linear
 relationship. Positive correlations imply that as x increases, so does
 y. Negative correlations imply that as x increases, y decreases.

 The p-value roughly indicates the probability of an uncorrelated system
 producing datasets that have a Pearson correlation at least as extreme
 as the one computed from these datasets. The p-values are not entirely
 reliable but are probably reasonable for datasets larger than 500 or so.

 Parameters
 ----------
 x : 1D array
 y : 1D array the same length as x

 Returns
 -------
 (Pearson's correlation coefficient,
  2-tailed p-value)

 References
 ----------
 http://www.statsoft.com/textbook/glosp.html#Pearson%20Correlation
Sacha
quelle
2
Wie wäre es mit einem Korrelationskoeffizienten von zwei Wörterbüchern ?!
user702846
2
@ user702846 Die Pearson-Korrelation wird in einer 2xN-Matrix definiert. Es gibt keine allgemein anwendbare Methode, mit der zwei Wörterbücher in eine 2xN-Matrix konvertiert werden. Sie können jedoch auch das Array von Wörterbuchwertpaaren verwenden, die den Schlüsseln des Schnittpunkts der Schlüssel Ihrer Wörterbücher entsprechen.
Winerd
115

Die Pearson-Korrelation kann mit Numpys berechnet werden corrcoef.

import numpy
numpy.corrcoef(list1, list2)[0, 1]
winerd
quelle
Die Ausgabe ist verwirrend, aber eigentlich sehr einfach. Überprüfen Sie diese Erklärung stackoverflow.com/a/3425548/1245622
linqu
56

Eine Alternative kann eine native Scipy-Funktion von linregress sein, die Folgendes berechnet:

Steigung: Steigung der Regressionslinie

Achsenabschnitt: Achsenabschnitt der Regressionslinie

r-Wert: Korrelationskoeffizient

p-Wert: zweiseitiger p-Wert für einen Hypothesentest, dessen Nullhypothese lautet, dass die Steigung Null ist

stderr: Standardfehler der Schätzung

Und hier ist ein Beispiel:

a = [15, 12, 8, 8, 7, 7, 7, 6, 5, 3]
b = [10, 25, 17, 11, 13, 17, 20, 13, 9, 15]
from scipy.stats import linregress
linregress(a, b)

wird Sie zurückgeben:

LinregressResult(slope=0.20833333333333337, intercept=13.375, rvalue=0.14499815458068521, pvalue=0.68940144811669501, stderr=0.50261704627083648)
Salvador Dali
quelle
2
Tolle Antwort - bei weitem die informativste. Funktioniert auch mit einem zweireihigen Pandas.DataFrame:lineregress(two_row_df)
dmeu
Geniale Antwort. Sehr intuitiv, wenn Sie daran denken
Raghuram
37

Wenn Sie keine Lust haben, scipy zu installieren, habe ich diesen schnellen Hack verwendet, der leicht von Programming Collective Intelligence geändert wurde :

(Zur Richtigkeit bearbeitet.)

from itertools import imap

def pearsonr(x, y):
  # Assume len(x) == len(y)
  n = len(x)
  sum_x = float(sum(x))
  sum_y = float(sum(y))
  sum_x_sq = sum(map(lambda x: pow(x, 2), x))
  sum_y_sq = sum(map(lambda x: pow(x, 2), y))
  psum = sum(imap(lambda x, y: x * y, x, y))
  num = psum - (sum_x * sum_y/n)
  den = pow((sum_x_sq - pow(sum_x, 2) / n) * (sum_y_sq - pow(sum_y, 2) / n), 0.5)
  if den == 0: return 0
  return num / den
Jeff Hammerbacher
quelle
2
Ich war überrascht, dass dies nicht mit Excel, NumPy und R übereinstimmt . Siehe stackoverflow.com/questions/3949226/… .
dfrankow
2
Wie ein anderer Kommentator betonte, hat dies einen float / int-Fehler. Ich denke, sum_y / n ist eine ganzzahlige Division für Ints. Wenn Sie sum_x = float (sum (x)) und sum_y = float (sum (y)) verwenden, funktioniert es.
dfrankow
@dfrankow Ich denke, das liegt daran, dass imap nicht mit float umgehen kann. Python gibt ein TypeError: unsupported operand type(s) for -: 'itertools.imap' and 'float'annum = psum - (sum_x * sum_y/n)
alvas
4
Als Stilnotiz runzelt Python die Stirn wegen dieser unnötigen Verwendung der Karte (zugunsten des Listenverständnisses)
Maxim Khesin
14
Bedenken Sie nur als Kommentar, dass Bibliotheken wie scipy et al. Von Personen entwickelt werden, die viel über numerische Analysen wissen. Dies kann Ihnen viele häufige Fallstricke ersparen (zum Beispiel kann eine sehr große und sehr kleine Anzahl in X oder Y zu einer katastrofischen Stornierung führen)
Geekazoid
32

Der folgende Code ist eine direkte Interpretation der Definition :

import math

def average(x):
    assert len(x) > 0
    return float(sum(x)) / len(x)

def pearson_def(x, y):
    assert len(x) == len(y)
    n = len(x)
    assert n > 0
    avg_x = average(x)
    avg_y = average(y)
    diffprod = 0
    xdiff2 = 0
    ydiff2 = 0
    for idx in range(n):
        xdiff = x[idx] - avg_x
        ydiff = y[idx] - avg_y
        diffprod += xdiff * ydiff
        xdiff2 += xdiff * xdiff
        ydiff2 += ydiff * ydiff

    return diffprod / math.sqrt(xdiff2 * ydiff2)

Prüfung:

print pearson_def([1,2,3], [1,5,7])

kehrt zurück

0.981980506062

Dies stimmt mit Excel überein, diesem Rechner , SciPy (auch NumPy ), der 0,981980506 und 0,9819805060619657 bzw. 0,98198050606196574 zurückgibt.

R :

> cor( c(1,2,3), c(1,5,7))
[1] 0.9819805

BEARBEITEN : Es wurde ein Fehler behoben, auf den ein Kommentator hinwies.

dfrankow
quelle
4
Vorsicht vor der Art der Variablen! Sie haben ein int / float-Problem festgestellt. In sum(x) / len(x)dir teile Ints, nicht Floats. Also sum([1,5,7]) / len([1,5,7]) = 13 / 3 = 4nach ganzzahliger Division (wohingegen Sie wollen 13. / 3. = 4.33...). Um dies zu beheben, schreiben Sie diese Zeile wie folgt um float(sum(x)) / float(len(x))(ein Gleitkomma reicht aus, da Python es automatisch konvertiert).
Piotr Migdal
Ihr Code funktioniert nicht in Fällen wie: [10,10,10], [0,0,0] oder [10,10], [10,0]. oder sogar [10,10], [10,10]
madCode
4
Der Korrelationskoeffizient ist für keinen dieser Fälle definiert. Wenn Sie sie in R setzen, wird für alle drei "NA" zurückgegeben.
dfrankow
28

Sie können dies auch tun mit pandas.DataFrame.corr:

import pandas as pd
a = [[1, 2, 3],
     [5, 6, 9],
     [5, 6, 11],
     [5, 6, 13],
     [5, 3, 13]]
df = pd.DataFrame(data=a)
df.corr()

Das gibt

          0         1         2
0  1.000000  0.745601  0.916579
1  0.745601  1.000000  0.544248
2  0.916579  0.544248  1.000000
Martin Thoma
quelle
5
Dies ist nur Korrelation ohne Bedeutung
Ivelin
12

Anstatt sich auf Numpy / Scipy zu verlassen, sollte meine Antwort am einfachsten zu codieren und die Schritte zur Berechnung des Pearson-Korrelationskoeffizienten (PCC) zu verstehen sein .

import math

# calculates the mean
def mean(x):
    sum = 0.0
    for i in x:
         sum += i
    return sum / len(x) 

# calculates the sample standard deviation
def sampleStandardDeviation(x):
    sumv = 0.0
    for i in x:
         sumv += (i - mean(x))**2
    return math.sqrt(sumv/(len(x)-1))

# calculates the PCC using both the 2 functions above
def pearson(x,y):
    scorex = []
    scorey = []

    for i in x: 
        scorex.append((i - mean(x))/sampleStandardDeviation(x)) 

    for j in y:
        scorey.append((j - mean(y))/sampleStandardDeviation(y))

# multiplies both lists together into 1 list (hence zip) and sums the whole list   
    return (sum([i*j for i,j in zip(scorex,scorey)]))/(len(x)-1)

Die Bedeutung von PCC besteht im Wesentlichen darin, Ihnen zu zeigen, wie stark die beiden Variablen / Listen korreliert sind. Es ist wichtig zu beachten, dass der PCC-Wert zwischen -1 und 1 liegt . Ein Wert zwischen 0 und 1 bedeutet eine positive Korrelation. Wert 0 = höchste Variation (keinerlei Korrelation). Ein Wert zwischen -1 und 0 bedeutet eine negative Korrelation.

compski
quelle
2
Beachten Sie, dass Python eine integrierte sumFunktion hat.
Bfontaine
5
Es hat eine erstaunliche Komplexität und langsame Leistung auf 2 Listen mit mehr als 500 Werten.
Nikolay Fominyh
9

Pearson-Koeffizientenberechnung mit Pandas in Python: Ich würde empfehlen, diesen Ansatz zu versuchen, da Ihre Daten Listen enthalten. Es ist einfach, mit Ihren Daten zu interagieren und sie über die Konsole zu bearbeiten, da Sie Ihre Datenstruktur visualisieren und nach Ihren Wünschen aktualisieren können. Sie können den Datensatz auch exportieren und speichern und neue Daten zur späteren Analyse aus der Python-Konsole hinzufügen. Dieser Code ist einfacher und enthält weniger Codezeilen. Ich gehe davon aus, dass Sie einige kurze Codezeilen benötigen, um Ihre Daten für die weitere Analyse zu überprüfen

Beispiel:

data = {'list 1':[2,4,6,8],'list 2':[4,16,36,64]}

import pandas as pd #To Convert your lists to pandas data frames convert your lists into pandas dataframes

df = pd.DataFrame(data, columns = ['list 1','list 2'])

from scipy import stats # For in-built method to get PCC

pearson_coef, p_value = stats.pearsonr(df["list 1"], df["list 2"]) #define the columns to perform calculations on
print("Pearson Correlation Coefficient: ", pearson_coef, "and a P-value of:", p_value) # Results 

Sie haben Ihre Daten jedoch nicht veröffentlicht, damit ich die Größe des Datensatzes oder die Transformationen sehen kann, die möglicherweise vor der Analyse erforderlich sind.

Web Ster
quelle
Hallo, willkommen bei StackOverflow! Fügen Sie zu Beginn Ihrer Antwort eine kurze Beschreibung hinzu, warum Sie diesen Code gewählt haben und wie er in diesem Fall gilt!
Tristo
8

Hmm, viele dieser Antworten haben langen und schwer lesbaren Code ...

Ich würde vorschlagen, numpy mit seinen raffinierten Funktionen zu verwenden, wenn Sie mit Arrays arbeiten:

import numpy as np
def pcc(X, Y):
   ''' Compute Pearson Correlation Coefficient. '''
   # Normalise X and Y
   X -= X.mean(0)
   Y -= Y.mean(0)
   # Standardise X and Y
   X /= X.std(0)
   Y /= Y.std(0)
   # Compute mean product
   return np.mean(X*Y)

# Using it on a random example
from random import random
X = np.array([random() for x in xrange(100)])
Y = np.array([random() for x in xrange(100)])
pcc(X, Y)
Martin F. Thomsen
quelle
Obwohl mir diese Antwort sehr gefällt, würde ich empfehlen, sowohl X als auch Y innerhalb der Funktion zu kopieren / klonen. Andernfalls werden beide geändert, was möglicherweise kein gewünschtes Verhalten ist.
Antonimmo
7

Dies ist eine Implementierung der Pearson-Korrelationsfunktion unter Verwendung von numpy:


def corr(data1, data2):
    "data1 & data2 should be numpy arrays."
    mean1 = data1.mean() 
    mean2 = data2.mean()
    std1 = data1.std()
    std2 = data2.std()

#     corr = ((data1-mean1)*(data2-mean2)).mean()/(std1*std2)
    corr = ((data1*data2).mean()-mean1*mean2)/(std1*std2)
    return corr

Mojtaba Khodadadi
quelle
7

Hier ist eine Variante von mkhs Antwort, die viel schneller läuft als sie, und scipy.stats.pearsonr mit numba.

import numba

@numba.jit
def corr(data1, data2):
    M = data1.size

    sum1 = 0.
    sum2 = 0.
    for i in range(M):
        sum1 += data1[i]
        sum2 += data2[i]
    mean1 = sum1 / M
    mean2 = sum2 / M

    var_sum1 = 0.
    var_sum2 = 0.
    cross_sum = 0.
    for i in range(M):
        var_sum1 += (data1[i] - mean1) ** 2
        var_sum2 += (data2[i] - mean2) ** 2
        cross_sum += (data1[i] * data2[i])

    std1 = (var_sum1 / M) ** .5
    std2 = (var_sum2 / M) ** .5
    cross_mean = cross_sum / M

    return (cross_mean - mean1 * mean2) / (std1 * std2)
Schildkröten sind süß
quelle
5

Hier ist eine Implementierung für die Pearson-Korrelation basierend auf einem spärlichen Vektor. Die Vektoren hier werden als Liste von Tupeln ausgedrückt als (Index, Wert). Die zwei spärlichen Vektoren können unterschiedlich lang sein, aber insgesamt muss die Vektorgröße gleich sein. Dies ist nützlich für Text Mining-Anwendungen, bei denen die Vektorgröße extrem groß ist, da die meisten Merkmale eine Wortsammlung sind und daher Berechnungen normalerweise mit spärlichen Vektoren durchgeführt werden.

def get_pearson_corelation(self, first_feature_vector=[], second_feature_vector=[], length_of_featureset=0):
    indexed_feature_dict = {}
    if first_feature_vector == [] or second_feature_vector == [] or length_of_featureset == 0:
        raise ValueError("Empty feature vectors or zero length of featureset in get_pearson_corelation")

    sum_a = sum(value for index, value in first_feature_vector)
    sum_b = sum(value for index, value in second_feature_vector)

    avg_a = float(sum_a) / length_of_featureset
    avg_b = float(sum_b) / length_of_featureset

    mean_sq_error_a = sqrt((sum((value - avg_a) ** 2 for index, value in first_feature_vector)) + ((
        length_of_featureset - len(first_feature_vector)) * ((0 - avg_a) ** 2)))
    mean_sq_error_b = sqrt((sum((value - avg_b) ** 2 for index, value in second_feature_vector)) + ((
        length_of_featureset - len(second_feature_vector)) * ((0 - avg_b) ** 2)))

    covariance_a_b = 0

    #calculate covariance for the sparse vectors
    for tuple in first_feature_vector:
        if len(tuple) != 2:
            raise ValueError("Invalid feature frequency tuple in featureVector: %s") % (tuple,)
        indexed_feature_dict[tuple[0]] = tuple[1]
    count_of_features = 0
    for tuple in second_feature_vector:
        count_of_features += 1
        if len(tuple) != 2:
            raise ValueError("Invalid feature frequency tuple in featureVector: %s") % (tuple,)
        if tuple[0] in indexed_feature_dict:
            covariance_a_b += ((indexed_feature_dict[tuple[0]] - avg_a) * (tuple[1] - avg_b))
            del (indexed_feature_dict[tuple[0]])
        else:
            covariance_a_b += (0 - avg_a) * (tuple[1] - avg_b)

    for index in indexed_feature_dict:
        count_of_features += 1
        covariance_a_b += (indexed_feature_dict[index] - avg_a) * (0 - avg_b)

    #adjust covariance with rest of vector with 0 value
    covariance_a_b += (length_of_featureset - count_of_features) * -avg_a * -avg_b

    if mean_sq_error_a == 0 or mean_sq_error_b == 0:
        return -1
    else:
        return float(covariance_a_b) / (mean_sq_error_a * mean_sq_error_b)

Unit Tests:

def test_get_get_pearson_corelation(self):
    vector_a = [(1, 1), (2, 2), (3, 3)]
    vector_b = [(1, 1), (2, 5), (3, 7)]
    self.assertAlmostEquals(self.sim_calculator.get_pearson_corelation(vector_a, vector_b, 3), 0.981980506062, 3, None, None)

    vector_a = [(1, 1), (2, 2), (3, 3)]
    vector_b = [(1, 1), (2, 5), (3, 7), (4, 14)]
    self.assertAlmostEquals(self.sim_calculator.get_pearson_corelation(vector_a, vector_b, 5), -0.0137089240555, 3, None, None)
Vipul Sharma
quelle
3

Ich habe eine sehr einfache und leicht verständliche Lösung dafür. Für zwei Arrays gleicher Länge kann der Pearson-Koeffizient leicht wie folgt berechnet werden:

def manual_pearson(a,b):
"""
Accepts two arrays of equal length, and computes correlation coefficient. 
Numerator is the sum of product of (a - a_avg) and (b - b_avg), 
while denominator is the product of a_std and b_std multiplied by 
length of array. 
"""
  a_avg, b_avg = np.average(a), np.average(b)
  a_stdev, b_stdev = np.std(a), np.std(b)
  n = len(a)
  denominator = a_stdev * b_stdev * n
  numerator = np.sum(np.multiply(a-a_avg, b-b_avg))
  p_coef = numerator/denominator
  return p_coef
aumpen
quelle
1

Sie fragen sich vielleicht, wie Sie Ihre Wahrscheinlichkeit im Zusammenhang mit der Suche nach einer Korrelation in eine bestimmte Richtung (negative oder positive Korrelation) interpretieren können. Hier ist eine Funktion, die ich geschrieben habe, um Ihnen dabei zu helfen. Es könnte sogar richtig sein!

Es basiert auf Informationen, die ich von http://www.vassarstats.net/rsig.html und http://en.wikipedia.org/wiki/Student%27s_t_distribution erhalten habe , dank anderer Antworten, die hier veröffentlicht wurden.

# Given (possibly random) variables, X and Y, and a correlation direction,
# returns:
#  (r, p),
# where r is the Pearson correlation coefficient, and p is the probability
# that there is no correlation in the given direction.
#
# direction:
#  if positive, p is the probability that there is no positive correlation in
#    the population sampled by X and Y
#  if negative, p is the probability that there is no negative correlation
#  if 0, p is the probability that there is no correlation in either direction
def probabilityNotCorrelated(X, Y, direction=0):
    x = len(X)
    if x != len(Y):
        raise ValueError("variables not same len: " + str(x) + ", and " + \
                         str(len(Y)))
    if x < 6:
        raise ValueError("must have at least 6 samples, but have " + str(x))
    (corr, prb_2_tail) = stats.pearsonr(X, Y)

    if not direction:
        return (corr, prb_2_tail)

    prb_1_tail = prb_2_tail / 2
    if corr * direction > 0:
        return (corr, prb_1_tail)

    return (corr, 1 - prb_1_tail)
Joshua Richardson
quelle
1

Sie können sich diesen Artikel ansehen. Dies ist ein gut dokumentiertes Beispiel für die Berechnung der Korrelation basierend auf historischen Forex-Währungspaardaten aus mehreren Dateien mithilfe der Pandas-Bibliothek (für Python) und die anschließende Erstellung eines Heatmap-Diagramms mithilfe der Seaborn-Bibliothek.

http://www.tradinggeeks.net/2015/08/calculating-correlation-in-python/

Oleg K.
quelle
0
def pearson(x,y):
  n=len(x)
  vals=range(n)

  sumx=sum([float(x[i]) for i in vals])
  sumy=sum([float(y[i]) for i in vals])

  sumxSq=sum([x[i]**2.0 for i in vals])
  sumySq=sum([y[i]**2.0 for i in vals])

  pSum=sum([x[i]*y[i] for i in vals])
  # Calculating Pearson correlation
  num=pSum-(sumx*sumy/n)
  den=((sumxSq-pow(sumx,2)/n)*(sumySq-pow(sumy,2)/n))**.5
  if den==0: return 0
  r=num/den
  return r
Quelle
quelle
Nur-Code-Antworten gelten nicht als bewährte Methode. Bitte fügen Sie ein paar Wörter hinzu, um zu erklären, wie Ihr Code die Frage beantwortet. (Lesen Sie die Hilfeseite zur Beantwortung einer Frage zu SO)
Yannis