Obliczanie korelacji i znaczenia Pearsona w Pythonie

Odpowiedzi:

202

Możesz spojrzeć na 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
źródło
2
A co ze współczynnikiem korelacji dwóch słowników ?!
user702846
2
@ user702846 Korelacja Pearsona jest zdefiniowana na macierzy 2xN. Nie ma powszechnie stosowanej metody, która konwertuje dwa słowniki na macierz 2xN, ale możesz użyć tablicy par wartości słownika odpowiadających kluczom przecięcia kluczy w słownikach.
winerd
115

Korelację Pearsona można obliczyć za pomocą liczb numpy'ego corrcoef.

import numpy
numpy.corrcoef(list1, list2)[0, 1]
winerd
źródło
wynik jest mylący, ale w rzeczywistości bardzo prosty. sprawdź to wyjaśnienie stackoverflow.com/a/3425548/1245622
linqu
56

Alternatywą może być natywna funkcja Scipy z linregress, która oblicza:

nachylenie: nachylenie linii regresji

intercept: intercept linii regresji

wartość r: współczynnik korelacji

wartość p: dwustronna wartość p dla testu hipotezy, którego hipotezą zerową jest to, że nachylenie wynosi zero

stderr: błąd standardowy oszacowania

A oto przykład:

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)

zwróci ci:

LinregressResult(slope=0.20833333333333337, intercept=13.375, rvalue=0.14499815458068521, pvalue=0.68940144811669501, stderr=0.50261704627083648)
Salvador Dali
źródło
2
Świetna odpowiedź - zdecydowanie najbardziej pouczająca. Działa również z dwurzędowymi pandami.DataFrame:lineregress(two_row_df)
dmeu
Genialna odpowiedź. Również bardzo intuicyjne, jeśli o tym pomyślisz
Raghuram
37

Jeśli nie masz ochoty instalować scipy, skorzystałem z tego szybkiego hacka, nieco zmodyfikowanego w stosunku do Programming Collective Intelligence :

(Edytowane pod kątem poprawności.)

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
źródło
2
Byłem zaskoczony, aby odkryć ten nie zgadza się z Excel NumPy i R. Patrz stackoverflow.com/questions/3949226/... .
dfrankow
2
Jak zauważył inny komentator, ma to błąd float / int. Myślę, że sum_y / n jest dzieleniem całkowitym dla liczb całkowitych. Jeśli użyjesz sum_x = float (sum (x)) i sum_y = float (sum (y)), to zadziała.
dfrankow
@dfrankow Myślę, że to dlatego, że imap nie radzi sobie z float. python podaje TypeError: unsupported operand type(s) for -: 'itertools.imap' and 'float'atnum = psum - (sum_x * sum_y/n)
alvas
4
Jako notatkę stylu Python marszczy brwi na temat tego niepotrzebnego użycia mapy (na korzyść list)
Maxim Khesin
14
Jako komentarz, weź pod uwagę, że biblioteki jako scipy i in. Są rozwijane przez ludzi znających wiele analiz numerycznych. Może to uniknąć wielu typowych pułapek (na przykład posiadanie bardzo dużych i bardzo małych liczb w X lub Y może spowodować katastrofalne anulowanie)
geekazoid
32

Poniższy kod jest prostą interpretacją definicji :

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)

Test:

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

zwroty

0.981980506062

To zgadza się z Excelem, tym kalkulatorem , SciPy (także NumPy Excelem ), które zwracają odpowiednio 0,981980506 i 0,9819805060619657 i 0,98198050606196574.

R :

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

EDYCJA : Naprawiono błąd wskazany przez komentatora.

dfrankow
źródło
4
Uważaj na rodzaj zmiennych! Wystąpił problem z int / float. W sum(x) / len(x)was dzielicie ints, a nie unosi się. Tak więc sum([1,5,7]) / len([1,5,7]) = 13 / 3 = 4, zgodnie z podziałem na liczby całkowite (podczas gdy chcesz 13. / 3. = 4.33...). Aby to naprawić, przepisz tę linię jako float(sum(x)) / float(len(x))(wystarczy jedna liczba zmiennoprzecinkowa, ponieważ Python konwertuje ją automatycznie).
Piotr Migdal
Twój kod nie będzie działać w przypadkach takich jak: [10,10,10], [0,0,0] lub [10,10], [10,0]. lub nawet [10,10], [10,10]
madCode
4
Współczynnik korelacji nie jest zdefiniowany dla żadnego z tych przypadków. Umieszczenie ich w R zwraca „NA” dla wszystkich trzech.
dfrankow
28

Możesz to również zrobić pandas.DataFrame.corrza pomocą:

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()

To daje

          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
źródło
5
To tylko korelacja bez znaczenia
Ivelin
12

Uważam, że zamiast polegać na numpy / scipy, moja odpowiedź powinna być najłatwiejsza do zakodowania i zrozumienia kroków w obliczaniu współczynnika korelacji Pearsona (PCC).

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)

Znaczenie PCC jest po prostu pokazać, jak silnie skorelowane dwie zmienne / listy są. Należy zauważyć, że wartość PCC wynosi od -1 do 1 . Wartość od 0 do 1 oznacza korelację dodatnią. Wartość 0 = najwyższa zmienność (bez żadnej korelacji). Wartość od -1 do 0 oznacza ujemną korelację.

compski
źródło
2
Zauważ, że Python ma wbudowaną sumfunkcję.
bfontaine,
5
Ma niesamowitą złożoność i niską wydajność na 2 listach z ponad 500 wartościami.
Nikolay Fominyh
9

Obliczanie współczynnika Pearsona za pomocą pand w pythonie: Sugeruję wypróbowanie tego podejścia, ponieważ dane zawierają listy. Łatwo będzie wchodzić w interakcję z danymi i manipulować nimi z poziomu konsoli, ponieważ możesz wizualizować strukturę danych i aktualizować ją według własnego uznania. Możesz także wyeksportować zestaw danych i zapisać go oraz dodać nowe dane z konsoli Pythona do późniejszej analizy. Ten kod jest prostszy i zawiera mniej wierszy kodu. Zakładam, że potrzebujesz kilku szybkich linii kodu, aby przesłać dane do dalszej analizy

Przykład:

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 

Jednak nie przesłałeś mi swoich danych, aby zobaczyć rozmiar zestawu danych lub transformacje, które mogą być potrzebne przed analizą.

Web Ster
źródło
Witaj, witaj w StackOverflow! Spróbuj dodać krótki opis, dlaczego wybrałeś ten kod i jak ma on zastosowanie w tym przypadku na początku swojej odpowiedzi!
Tristo
8

Hmm, wiele z tych odpowiedzi ma długi i trudny do odczytania kod ...

Podczas pracy z tablicami sugerowałbym używanie numpy z jego ciekawymi funkcjami:

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
źródło
Chociaż bardzo podoba mi się ta odpowiedź, radziłbym skopiować / sklonować X i Y wewnątrz funkcji. W przeciwnym razie oba zostaną zmienione, co może nie być pożądanym zachowaniem.
antonimmo
7

Jest to implementacja funkcji korelacji Pearsona za pomocą 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
źródło
7

Oto wariant odpowiedzi mkh, który działa znacznie szybciej od niego, i scipy.stats.pearsonr, używając 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)
Żółwie są słodkie
źródło
5

Oto implementacja korelacji Pearsona na podstawie rzadkiego wektora. Wektory tutaj są wyrażone jako lista krotek wyrażona jako (indeks, wartość). Dwa rzadkie wektory mogą mieć różną długość, ale dla całego rozmiaru wektora będą musiały być takie same. Jest to przydatne w aplikacjach do eksploracji tekstu, w których rozmiar wektora jest niezwykle duży ze względu na to, że większość funkcji to zbiór słów, a zatem obliczenia są zwykle wykonywane przy użyciu rzadkich wektorów.

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)

Testy jednostkowe:

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
źródło
3

Mam na to bardzo proste i łatwe do zrozumienia rozwiązanie. W przypadku dwóch tablic o równej długości współczynnik Pearsona można łatwo obliczyć w następujący sposób:

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
źródło
1

Możesz się zastanawiać, jak zinterpretować swoje prawdopodobieństwo w kontekście poszukiwania korelacji w określonym kierunku (korelacja ujemna lub dodatnia). Oto funkcja, którą napisałem, aby w tym pomóc. To może nawet mieć rację!

Opiera się na informacjach zebranych z http://www.vassarstats.net/rsig.html i http://en.wikipedia.org/wiki/Student%27s_t_distribution , dzięki innym odpowiedziom zamieszczonym tutaj.

# 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
źródło
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
Źródło
źródło
Odpowiedzi zawierające tylko kod nie są uważane za dobrą praktykę. Rozważ dodanie kilku słów, aby wyjaśnić, w jaki sposób Twój kod rozwiązuje pytanie. (przeczytaj stronę pomocy, w jaki sposób odpowiedzieć na pytanie dotyczące SO)
Yannis,