Cálculo de la correlación y la importancia de Pearson en Python

Respuestas:

202

Puedes echar un vistazo a 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
fuente
2
¿Qué tal el coeficiente de correlación de dos diccionarios?
user702846
2
@ user702846 La correlación de Pearson se define en una matriz 2xN. No existe un método generalmente aplicable que convierta dos diccionarios en una matriz de 2xN, pero puede usar la matriz de pares de valores de diccionario correspondientes a las claves de la intersección de las claves de sus diccionarios.
Winerd
115

La correlación de Pearson se puede calcular con numpy's corrcoef.

import numpy
numpy.corrcoef(list1, list2)[0, 1]
Winerd
fuente
La salida es confusa, pero en realidad muy simple. consulte esta explicación stackoverflow.com/a/3425548/1245622
linqu
56

Una alternativa puede ser una función scipy nativa de linregress que calcula:

pendiente: pendiente de la línea de regresión

intercepción: intercepción de la línea de regresión

valor r: coeficiente de correlación

valor p: valor p de dos lados para una prueba de hipótesis cuya hipótesis nula es que la pendiente es cero

stderr: error estándar de la estimación

Y aquí hay un ejemplo:

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)

te devolveré:

LinregressResult(slope=0.20833333333333337, intercept=13.375, rvalue=0.14499815458068521, pvalue=0.68940144811669501, stderr=0.50261704627083648)
Salvador Dalí
fuente
2
Gran respuesta, con mucho, la más informativa. También funciona con pandas de dos filas. Marco de datos:lineregress(two_row_df)
dmeu
Brillante respuesta. Muy intuitivo también, si lo piensas
Raghuram
37

Si no tiene ganas de instalar scipy, he usado este truco rápido, ligeramente modificado de Programming Collective Intelligence :

(Editado para su corrección).

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
fuente
2
Me sorprendió descubrir que esto no está de acuerdo con Excel, NumPy y R. Ver stackoverflow.com/questions/3949226/… .
dfrankow
2
Como señaló otro comentarista, esto tiene un error float / int. Creo que sum_y / n es una división entera para ints. Si usa sum_x = float (sum (x)) y sum_y = float (sum (y)), funciona.
dfrankow
@dfrankow Creo que es porque imap no puede manejar flotante. python da un TypeError: unsupported operand type(s) for -: 'itertools.imap' and 'float'atnum = psum - (sum_x * sum_y/n)
alvas el
44
Como nota de estilo, Python frunce el ceño ante este uso innecesario del mapa (a favor de la comprensión de la lista)
Maxim Khesin
14
Solo como comentario, considere que las bibliotecas como scipy et al son desarrolladas por personas que conocen muchos análisis numéricos. Esto puede evitar muchos problemas comunes (por ejemplo, tener números muy grandes y muy pequeños en X o Y puede resultar en una cancelación catastrófica)
geekazoid
32

El siguiente código es una interpretación directa de la definición :

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)

Prueba:

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

devoluciones

0.981980506062

Esto concuerda con Excel, esta calculadora , SciPy (también NumPy ), que devuelve 0.981980506 y 0.9819805060619657 y 0.98198050606196574, respectivamente.

R :

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

EDITAR : se corrigió un error señalado por un comentarista.

dfrankow
fuente
44
¡Cuidado con el tipo de variables! Ha encontrado un problema int / float. En sum(x) / len(x)ustedes dividen ints, no flotadores. Entonces sum([1,5,7]) / len([1,5,7]) = 13 / 3 = 4, de acuerdo con la división de enteros (mientras lo desee 13. / 3. = 4.33...). Para solucionarlo, vuelva a escribir esta línea como float(sum(x)) / float(len(x))(un flotador es suficiente, ya que Python lo convierte automáticamente).
Piotr Migdal
Su código no funcionará para casos como: [10,10,10], [0,0,0] o [10,10], [10,0]. o incluso [10,10], [10,10]
madCode
44
El coeficiente de correlación no está definido para ninguno de esos casos. Ponerlos en R devuelve "NA" para los tres.
dfrankow
28

También puedes hacer esto con 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()

Esto da

          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
fuente
55
Esto es solo correlación sin importancia
Ivelin
12

En lugar de confiar en numpy / scipy, creo que mi respuesta debería ser la más fácil de codificar y comprender los pasos para calcular el coeficiente de correlación de Pearson (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)

La importancia de PCC es básicamente mostrarle cuán fuertemente correlacionadas están las dos variables / listas. Es importante tener en cuenta que el valor de PCC varía de -1 a 1 . Un valor entre 0 y 1 denota una correlación positiva. Valor de 0 = variación más alta (sin correlación alguna). Un valor entre -1 y 0 denota una correlación negativa.

compski
fuente
2
Tenga en cuenta que Python tiene una sumfunción integrada.
bfontaine
55
Tiene una complejidad sorprendente y un rendimiento lento en 2 listas con más de 500 valores.
Nikolay Fominyh
9

Cálculo del coeficiente de Pearson usando pandas en python: sugeriría probar este enfoque ya que sus datos contienen listas. Será fácil interactuar con sus datos y manipularlos desde la consola, ya que puede visualizar su estructura de datos y actualizarla como desee. También puede exportar el conjunto de datos y guardarlo y agregar nuevos datos de la consola de Python para su posterior análisis. Este código es más simple y contiene menos líneas de código. Supongo que necesita algunas líneas rápidas de código para filtrar sus datos para un análisis posterior.

Ejemplo:

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 

Sin embargo, no publicó sus datos para que yo vea el tamaño del conjunto de datos o las transformaciones que podrían ser necesarias antes del análisis.

Web Ster
fuente
Hola, bienvenido a StackOverflow! ¡Intente agregar una breve descripción de por qué eligió este código y cómo se aplica en este caso al comienzo de su respuesta!
Tristo
8

Hmm, muchas de estas respuestas tienen un código largo y difícil de leer ...

Sugeriría usar numpy con sus ingeniosas características cuando trabaje con matrices:

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
fuente
Aunque me gusta mucho esta respuesta, recomendaría copiar / clonar tanto X como Y dentro de la función. De lo contrario, ambos se alteran, lo que podría no ser un comportamiento deseado.
antonimmo
7

Esta es una implementación de la función de correlación de Pearson usando 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
fuente
7

Aquí hay una variante en la respuesta de mkh que se ejecuta mucho más rápido que scipy.stats.pearsonr, usando 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)
Las tortugas son lindas
fuente
5

Aquí hay una implementación para la correlación de Pearson basada en un vector disperso. Los vectores aquí se expresan como una lista de tuplas expresadas como (índice, valor). Los dos vectores dispersos pueden tener una longitud diferente, pero sobre todo el tamaño del vector tendrá que ser el mismo. Esto es útil para aplicaciones de minería de texto donde el tamaño del vector es extremadamente grande debido a que la mayoría de las características son bolsas de palabras y, por lo tanto, los cálculos generalmente se realizan usando vectores dispersos.

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)

Pruebas unitarias:

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
fuente
3

Tengo una solución muy simple y fácil de entender para esto. Para dos matrices de igual longitud, el coeficiente de Pearson se puede calcular fácilmente de la siguiente manera:

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
fuente
1

Quizás se pregunte cómo interpretar su probabilidad en el contexto de buscar una correlación en una dirección particular (correlación negativa o positiva). Aquí hay una función que escribí para ayudar con eso. ¡Incluso podría ser correcto!

Se basa en la información que obtuve de http://www.vassarstats.net/rsig.html y http://en.wikipedia.org/wiki/Student%27s_t_distribution , gracias a otras respuestas publicadas aquí.

# 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
fuente
1

Puedes echar un vistazo a este artículo. Este es un ejemplo bien documentado para calcular la correlación basada en datos históricos de pares de divisas forex de múltiples archivos usando la biblioteca de pandas (para Python), y luego generar un diagrama de mapa de calor usando la biblioteca naciente.

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

Oleg K
fuente
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
Fuente
fuente
Las respuestas de solo código no se consideran una buena práctica. Considere agregar algunas palabras para explicar cómo su código aborda la pregunta. (lea la página de ayuda sobre cómo responder una pregunta sobre SO)
Yannis