Я пытаюсь получить под капотом stats.pearsonr
для моего собственного понимания.Я хотел бы написать r, p = stats.pearsonr(x, y)
как отдельный код (см. Ниже).Я могу воспроизвести r
, однако для вычисления p-значения используется scipy.special.betainc
.Я порылся в https://github.com/scipy/scipy/tree/master/scipy/special, но не смог найти betainc
(полагаю, это может быть какой-то основной код C).Любая идея, где я могу найти исходный код для этого?
import numpy as np
x = np.random.rand(100)
y = np.random.rand(100)
n = len(x)
mx = x.mean()
my = y.mean()
xm, ym = x - mx, y - my
r_num = np.add.reduce(xm * ym)
r_den = np.sqrt(np.sum(xm*xm, 0) * np.sum(ym*ym, 0))
r = r_num / r_den
r = max(min(r, 1.0), -1.0)
df = n - 2
if abs(r) == 1.0:
p = 0.0
else:
t_squared = r**2 * (df / ((1.0 - r) * (1.0 + r)))
# Now code https://github.com/scipy/scipy/blob/14142ff70d84a6ce74044a27919c850e893648f7/scipy/stats/stats.py#L3020
_x = df/(df+t_squared)
_x = np.asarray(_x)
_x = np.where(_x < 1.0, _x, 1.0)
_a = 0.5*df
_b = 0.5
# Code special.betainc(_a, _b, _x)