Вычисление корреляции и значимости Пирсона в Python

Я ищу функцию, которая принимает в качестве входных данных два списка и возвращает корреляцию Пирсона и значение корреляции.

4b9b3361

Вы можете посмотреть 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 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 correlation coefficient,
  2-tailed p-value)

 References
 ----------
 http://www.statsoft.com/textbook/glosp.html#Pearson%20Correlation
144
ответ дан 16 окт. '10 в 17:29
источник

Корреляция Пирсона может быть рассчитана с помощью numpy corrcoef.

import numpy
numpy.corrcoef(list1, list2)[0, 1]
68
ответ дан 16 апр. '13 в 3:17
источник

Если вам не хочется устанавливать scipy, я использовал этот быстрый хак, слегка модифицированный Программирование коллективного интеллекта:

(Отредактировано для правильности.)

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
30
ответ дан 19 апр. '11 в 11:52
источник

Альтернативой может быть встроенная функция scipy из linregress, которая вычисляет:

наклон: наклон линии регрессии

перехват: перехват линии регрессии

r-значение: коэффициент корреляции

p-значение: двухстороннее p-значение для теста гипотезы, нулевой гипотезой которого является то, что наклон равен нулю

stderr: стандартная ошибка оценки

И вот пример:

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)

вернет вас:

LinregressResult(slope=0.20833333333333337, intercept=13.375, rvalue=0.14499815458068521, pvalue=0.68940144811669501, stderr=0.50261704627083648)
25
ответ дан 21 нояб. '15 в 16:10
источник

Следующий код представляет собой прямолинейную интерпретацию определения:

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)

Тест:

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

возвращает

0.981980506062

Это согласуется с Excel, этот калькулятор, SciPy (также NumPy), которые возвращают 0.981980506 и 0.9819805060619657, и 0.98198050606196574 соответственно.

R:

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

EDIT: исправлена ​​ошибка, отмеченная комментатором.

24
ответ дан 29 окт. '11 в 16:42
источник

Вместо того, чтобы полагаться на numpy/scipy, я думаю, что мой ответ должен быть самым простым для кода и понять шаги при вычислении коэффициента корреляции Пирсона (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)

Значение PCC в основном показывает вам, как сильно коррелирует две переменные/списки. Важно отметить, что значение PCC составляет от -1 до 1. Значение от 0 до 1 означает положительную корреляцию. Значение 0 = наибольшая вариация (без корреляции). Значение от -1 до 0 означает отрицательную корреляцию.

11
ответ дан 30 июня '13 в 14:39
источник

Хм, многие из этих ответов долго и трудно читать код...

Я предлагаю использовать numpy с его отличными функциями при работе с массивами:

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)
6
ответ дан 31 окт. '13 в 18:28
источник

Вы можете сделать это с помощью 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()

Это дает

          0         1         2
0  1.000000  0.745601  0.916579
1  0.745601  1.000000  0.544248
2  0.916579  0.544248  1.000000
6
ответ дан 29 янв. '17 в 14:28
источник

Здесь вариант ответа mkh, который работает намного быстрее, чем этот, и scipy.stats.pearsonr, используя 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)
3
ответ дан 22 марта '15 в 15:55
источник

Вот реализация корреляции pearson, основанная на разреженном векторе. Здесь векторы выражаются в виде списка кортежей, выраженных как (индекс, значение). Два разреженных вектора могут иметь разную длину, но по всему размеру вектора должны быть одинаковыми. Это полезно для приложений для интеллектуального анализа текста, где размер вектора чрезвычайно велик из-за того, что большинство функций являются мешками слов, и, следовательно, вычисления обычно выполняются с использованием разреженных векторов.

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)

Единичные тесты:

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)
3
ответ дан 15 сент. '13 в 23:19
источник

Это реализация функции корреляции Пирсона с использованием 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

2
ответ дан 03 марта '15 в 0:13
источник

Вы можете взглянуть на эту статью. Это хорошо документированный пример расчета корреляции на основе исторических валютных парных данных из нескольких файлов с использованием библиотеки pandas (для Python), а затем создания графика тепловой карты с использованием библиотеки морского судна.

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

1
ответ дан 07 авг. '15 в 0:00
источник

Вы можете задаться вопросом, как интерпретировать вашу вероятность в контексте поиска корреляции в определенном направлении (отрицательная или положительная корреляция). Вот функция, которую я написал, чтобы помочь с этим. Это может быть даже правильно!

На основе информации, которую я почерпнул из http://www.vassarstats.net/rsig.html и http://en.wikipedia.org/wiki/Student%27s_t_distribution, благодаря другим ответам, размещенным здесь.

# 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)
0
ответ дан 09 янв. '13 в 5:18
источник