Я работаю с некоторой статистикой и подумал, что напишу свою собственную функцию для выполнения корреляции Пирсона.По какой-то причине мой код дает неправильное значение p, которое не согласуется с scipy.Вот мой код:
x = [438.69461,
404.060132,
198.210118,
386.826542,
346.090443,
330.09382,
249.341107,
324.826081,
303.760307,
388.144244]
y = [115.343021,
89.971217,
174.421423,
155.054388,
22.921773,
352.884202,
254.11353,
176.763594,
265.541998,
222.213654]
N = len(list(zip(x, y)))
# Calculate covariance
cov = sum([(i - np.mean(x)) * (j - np.mean(y)) for i, j in zip(x, y)]) / (N-1)
# Calculate the standard devs of x and y
sx, sy = np.std(x, ddof=1), np.std(y, ddof=1)
r = cov / (sx * sy)
# r is not normally distributed, so convert to z
z_r = np.log((1+r)/(1-r)) / 2
SE_z_r = 1 / np.sqrt(N-3)
# See Field (2009) pp.172 for discussion on z score relative to particular non-zero value
z = z_r / (SE_z_r)
# convert to t value
t_r = r * np.sqrt(N-2) / np.sqrt(1 - (r**2))
# CIs for r as z
lower_z = z_r - (1.96 * SE_z_r)
upper_z = z_r + (1.96 * SE_z_r)
# Finally convert these values back to correlation coefficients...
lower = (np.exp(2*lower_z) - 1) / (np.exp(2*lower_z) + 1)
upper = (np.exp(2*upper_z) - 1) / (np.exp(2*upper_z) + 1)
# Compute the p value from the t-statistic
p = stats.norm.sf(z)*2
result = {'Pearson\'s r': r,
'95% CI (upper)':upper,
'95% CI (lower)':lower,
'Sig. (two-tailed)':p,
'N':N}
Это дает следующий результат:
Out[315]:
{'95% CI (lower)': -0.8080015100222782,
'95% CI (upper)': 0.3455450058720938,
'N': 10,
"Pearson's r": -0.3630848051556617,
'Sig. (two-tailed)': 1.6858418386527083}
Что отличается от scipy:
stats.pearsonr(x, y)
Out[316]: (-0.3630848051556618, 0.30243603632310984)
Я предполагаю, что я используюнеправильный расчет для значения р.Так что, если я использую p = stats.t.cdf(t_r, N-2)*2 для вычисления значения p, он выглядит хорошо с небольшой разницей:
Out[339]:
{'95% CI (lower)': -0.8080015100222782,
'95% CI (upper)': 0.3455450058720938,
'N': 10,
"Pearson's r": -0.3630848051556617,
'Sig. (two-tailed)': 0.3024360363231102}
stats.pearsonr(x, y)
Out[340]: (-0.3630848051556618, 0.30243603632310984)
Но затем изменение x и y на следующее:
x = [52, 41, 44, 68, 83]
y = [86, 92, 107, 135, 150]
даетстранные значения еще раз для р:
Out[343]:
{'95% CI (lower)': 0.07262140356975652,
'95% CI (upper)': 0.9932583155633771,
'N': 5,
"Pearson's r": 0.8973956836822184,
'Sig. (two-tailed)': 1.9611596807577028}
stats.pearsonr(x, y)
Out[344]: (0.8973956836822184, 0.03884031924229741)