正如@MB-F在评论中指出的，它是通过分析计算得出的。
在版本0.19.1的代码中，您可以隔离该部分代码，并根据r绘制p值

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import betainc
r = np.linspace(-1, 1, 1000)*(1-1e-10);

for n in [10, 100, 1000]:
    df = n - 2
    t_squared = r**2 * (df / ((1.0 - r) * (1.0 + r)))
    prob = betainc(0.5*df, 0.5, df/(df+t_squared))
    plt.semilogy(r, prob, label=f'n={n}')
plt.axvline(0.9672434106763087, ls='--', color='black', label='r value')
plt.legend()
plt.grid()

当前稳定版本1.9.3使用不同的公式

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import btdtr
r = np.linspace(-1, 1, 1000)*(1-1e-10);
for n in [10, 100, 1000]:
    ab = 0.5*n
    prob = btdtr(ab, ab, 0.5*(1-abs(r)))
    plt.semilogy(r, prob, label=f'n={n}')
plt.axvline(0.9672434106763087, ls='--', color='black', label='r value')
plt.legend()
plt.grid()

但产生相同的结果。
您可以看到，如果您有1000个点和您的相关性，p值将小于最小浮动值。

beta分布

Scipy提供了概率分布的集合，其中包括beta分布。
行

prob = btdtr(ab, ab, 0.5*(1-abs(r)))

可替换为

from scipy.stats import beta
prob = beta(ab, ab).cdf(0.5*(1-abs(r)))

在那里你可以得到更多的信息。