simple_correlation.sagews

Пример вычисления коэффициента корреляции Пирсона

Чтение данных из файла

f = open('longley.csv', 'r')
firstline = f.readline()
data = [firstline.split(',')]
for line in f :
    row = line[:-1].split(',')
    data.append([row[0]]+[float(x) for x in row[1:]])
f.close()
# Выделяем из данных только числовой массив
num_data = [row[1:] for row in data[1:]]
num_data = matrix(num_data)
# Проверяем, что всё считано корректно
show(num_data)

$\displaystyle \left(\begin{array}{rrrrrrr} 83.0 & 234.289 & 235.6 & 159.0 & 107.608 & 1947.0 & 60.323 \\ 88.5 & 259.426 & 232.5 & 145.6 & 108.632 & 1948.0 & 61.122 \\ 88.2 & 258.054 & 368.2 & 161.6 & 109.773 & 1949.0 & 60.171 \\ 89.5 & 284.599 & 335.1 & 165.0 & 110.929 & 1950.0 & 61.187 \\ 96.2 & 328.975 & 209.9 & 309.9 & 112.075 & 1951.0 & 63.221 \\ 98.1 & 346.999 & 193.2 & 359.4 & 113.27 & 1952.0 & 63.639 \\ 99.0 & 365.385 & 187.0 & 354.7 & 115.094 & 1953.0 & 64.989 \\ 100.0 & 363.112 & 357.8 & 335.0 & 116.219 & 1954.0 & 63.761 \\ 101.2 & 397.469 & 290.4 & 304.8 & 117.388 & 1955.0 & 66.019 \\ 104.6 & 419.18 & 282.2 & 285.7 & 118.734 & 1956.0 & 67.857 \\ 108.4 & 442.769 & 293.6 & 279.8 & 120.445 & 1957.0 & 68.169 \\ 110.8 & 444.546 & 468.1 & 263.7 & 121.95 & 1958.0 & 66.513 \\ 112.6 & 482.704 & 381.3 & 255.2 & 123.366 & 1959.0 & 68.655 \\ 114.2 & 502.601 & 393.1 & 251.4 & 125.368 & 1960.0 & 69.564 \\ 115.7 & 518.173 & 480.6 & 257.2 & 127.852 & 1961.0 & 69.331 \\ 116.9 & 554.894 & 400.7 & 282.7 & 130.081 & 1962.0 & 70.551 \end{array}\right)$

По формулам из методички

$$\widehat{r}_{XY} =\frac{n\, \sum\limits _{i=1}^{n}x_{i} y_{i} -\sum\limits _{i=1}^{n}x_{i} \cdot \sum\limits _{i=1}^{n}y_{i} }{\sqrt{\left[n\, \sum\limits _{i=1}^{n}x_{i}^{2} -\left(\sum\limits _{i=1}^{n}x_{i} \right)^{2} \right]\cdot \left[n\, \sum\limits _{i=1}^{n}y_{i}^{2} -\left(\sum\limits _{i=1}^{n}y_{i} \right)^{2} \right]} };$$$$t_{\text{СТ}}=\left|\widehat{r}\right|\sqrt{\dfrac{n-2}{1-\widehat{r}^2}}.$$

# Корреляция между Population (столбец 5) и Employed (столбец 7)
x = list(num_data.column(4))
y = list(num_data.column(6))
sp = scatter_plot(zip(x, y), axes_labels=['Population', 'Employed'])
sp.show()
print
n = len(x)
x_sq = [elem^2 for elem in x]
y_sq = [elem^2 for elem in y]
xy = [x[i]*y[i] for i in range(n)]
r = (n*sum(xy)-sum(x)*sum(y))/sqrt((n*sum(x_sq)-(sum(x))^2)*(n*sum(y_sq)-(sum(y))^2))
print 'Корреляция между', data[0][5], 'и', data[0][7], ':', r
t_stat = abs(r)*sqrt((n-2)/1-r^2)
t_crit = 2.1448 # нашли по таблице при alpha=0.1, n=16 --> t(0.05; 14)
if t_stat > t_crit :
    print 'Корреляция значима'
else :
    print 'Корреляция не выявлена'

Корреляция между "Population" и "Employed" : 0.960390571594 Корреляция значима

Использование встроенных функций

from scipy.stats import pearsonr
# Корреляция между Population (столбец 5) и Employed (столбец 7)
Pop = num_data.column(4)
Empl = num_data.column(6)
r = pearsonr(Pop, Empl)[0]
p_value = pearsonr(Pop, Empl)[1]
print 'Корреляция между', data[0][5], 'и', data[0][7], ':', r, ', p-значение равно', p_value
if p_value < 0.05 :
    print 'Корреляция значима'
else :
    print 'Корреляция не выявлена'

Корреляция между "Population" и "Employed" : 0.960390571594 , p-значение равно 3.69324508517e-09 Корреляция значима

Использование R

%r
with(longley, cor.test(Population, Employed))

Pearson's product-moment correlation data: Population and Employed t = 12.896, df = 14, p-value = 3.693e-09 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.8869236 0.9864676 sample estimates: cor 0.9603906