皮尔逊相关系数#
- scipy.stats.pearsonr(x, y, *, alternative='two-sided', method=None, axis=0)[源代码][源代码]#
皮尔逊相关系数和用于检验非相关的p值。
皮尔逊相关系数 [1] 衡量两个数据集之间的线性关系。与其他相关系数一样,其值在 -1 和 +1 之间变化,0 表示没有相关性。相关系数为 -1 或 +1 表示存在完全线性关系。正相关表示随着 x 增加,y 也增加。负相关表示随着 x 增加,y 减少。
此函数还对零假设进行检验,即样本下的分布是不相关的且正态分布的。(参见 Kowalski [3] 关于输入非正态性对相关系数分布影响的讨论。)p 值大致表示不相关的系统产生数据集的概率,这些数据集的皮尔逊相关性至少与从这些数据集计算出的相关性一样极端。
- 参数:
- xarray_like
输入数组。
- yarray_like
输入数组。
- 轴int 或 None, 默认
执行计算的轴。默认值为0。如果为None,则在执行计算之前将两个数组展平。
Added in version 1.13.0.
- 替代方案{‘双侧’, ‘大于’, ‘小于’}, 可选
定义备择假设。默认是’双侧’。以下选项可用:
‘双边’: 相关性非零
‘less’: 相关性为负(小于零)
‘greater’: 相关性为正(大于零)
Added in version 1.9.0.
- 方法重采样方法, 可选
定义用于计算p值的方法。如果 method 是
PermutationMethod
/MonteCarloMethod
的实例,则使用scipy.stats.permutation_test
/scipy.stats.monte_carlo_test
并结合提供的配置选项和其他适当的设置来计算p值。否则,p值将按照注释中的说明进行计算。Added in version 1.11.0.
- 返回:
- 结果 :
PearsonRResult
PearsonRResult
一个具有以下属性的对象:
- 统计浮动
皮尔逊积矩相关系数。
- p值浮动
与所选替代方案相关的p值。
该对象具有以下方法:
- confidence_interval(置信水平, 方法)
这将计算给定置信水平下相关系数 statistic 的置信区间。置信区间以
namedtuple
的形式返回,字段为 low 和 high。如果未提供 method,则使用 Fisher 变换 [1] 计算置信区间。如果 method 是BootstrapMethod
的实例,则使用scipy.stats.bootstrap
根据提供的配置选项和其他适当设置计算置信区间。在某些情况下,由于退化重采样,置信限可能为 NaN,这在非常小的样本(约6个观察值)中是典型的。
- 结果 :
- 警告:
ConstantInputWarning
如果输入是一个常数数组,则会引发此错误。在这种情况下,相关系数未定义,因此返回
np.nan
。NearConstantInputWarning
如果输入是“几乎”恒定的,则会引发此错误。如果
norm(x - mean(x)) < 1e-13 * abs(mean(x))
,则认为数组x
几乎是恒定的。在这种情况下,计算x - mean(x)
时的数值误差可能导致 r 的计算不准确。
参见
spearmanr
斯皮尔曼等级相关系数。
kendalltau
Kendall’s tau,一种用于有序数据的关联度量。
注释
相关系数计算如下:
\[r = \frac{\sum (x - m_x) (y - m_y)} {\sqrt{\sum (x - m_x)^2 \sum (y - m_y)^2}}\]其中 \(m_x\) 是向量 x 的均值,\(m_y\) 是向量 y 的均值。
在假设 x 和 y 来自独立正态分布(因此总体相关系数为 0)的情况下,样本相关系数 r 的概率密度函数为([1],[2]):
\[f(r) = \frac{{(1-r^2)}^{n/2-2}}{\mathrm{B}(\frac{1}{2},\frac{n}{2}-1)}\]其中 n 是样本数量,B 是 beta 函数。这有时被称为 r 的精确分布。这是在
pearsonr
中用于计算 p 值的分布,当 method 参数保留其默认值(None)时。该分布是区间 [-1, 1] 上的 beta 分布,形状参数 a = b = n/2 - 1。在 SciPy 的 beta 分布实现中,r 的分布为:dist = scipy.stats.beta(n/2 - 1, n/2 - 1, loc=-1, scale=2)
pearsonr
返回的默认 p 值是双侧 p 值。对于具有相关系数 r 的给定样本,p 值是随机抽取的样本 x’ 和 y’ 的 abs(r’) 大于或等于 abs(r) 的概率,这些样本来自零相关的总体。在上述dist
对象的术语中,给定 r 和长度 n 的 p 值可以计算为:p = 2*dist.cdf(-abs(r))
当 n 为 2 时,上述连续分布未被明确定义。可以将 beta 分布的极限解释为形状参数 a 和 b 趋近于 a = b = 0 时的离散分布,其在 r = 1 和 r = -1 处具有相等的概率质量。更直接地,可以观察到,给定数据 x = [x1, x2] 和 y = [y1, y2],并且假设 x1 != x2 且 y1 != y2,r 的唯一可能值为 1 和 -1。因为对于任何长度为 2 的样本 x’ 和 y’,abs(r’) 都将为 1,所以长度为 2 的样本的双侧 p 值总是 1。
为了向后兼容,返回的对象也表现得像一个包含统计量和p值的长度为二的元组。
参考文献
[2]Student, “相关系数的可能误差”, Biometrika, 第6卷, 第2-3期, 1908年9月1日, 第302-310页。
[3]C. J. Kowalski, “On the Effects of Non-Normality on the Distribution of the Sample Product-Moment Correlation Coefficient” Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 21, No. 1 (1972), pp. 1-12.
示例
>>> import numpy as np >>> from scipy import stats >>> x, y = [1, 2, 3, 4, 5, 6, 7], [10, 9, 2.5, 6, 4, 3, 2] >>> res = stats.pearsonr(x, y) >>> res PearsonRResult(statistic=-0.828503883588428, pvalue=0.021280260007523286)
要执行测试的精确排列版本:
>>> rng = np.random.default_rng() >>> method = stats.PermutationMethod(n_resamples=np.inf, random_state=rng) >>> stats.pearsonr(x, y, method=method) PearsonRResult(statistic=-0.828503883588428, pvalue=0.028174603174603175)
在零假设下进行测试,假设数据是从 均匀 分布中抽取的:
>>> method = stats.MonteCarloMethod(rvs=(rng.uniform, rng.uniform)) >>> stats.pearsonr(x, y, method=method) PearsonRResult(statistic=-0.828503883588428, pvalue=0.0188)
要生成一个渐近的90%置信区间:
>>> res.confidence_interval(confidence_level=0.9) ConfidenceInterval(low=-0.9644331982722841, high=-0.3460237473272273)
而对于一个引导置信区间:
>>> method = stats.BootstrapMethod(method='BCa', random_state=rng) >>> res.confidence_interval(confidence_level=0.9, method=method) ConfidenceInterval(low=-0.9983163756488651, high=-0.22771001702132443) # may vary
如果提供了N维数组,则根据与大多数
scipy.stats
函数相同的约定,在一次调用中执行多个测试:>>> rng = np.random.default_rng() >>> x = rng.standard_normal((8, 15)) >>> y = rng.standard_normal((8, 15)) >>> stats.pearsonr(x, y, axis=0).statistic.shape # between corresponding columns (15,) >>> stats.pearsonr(x, y, axis=1).statistic.shape # between corresponding rows (8,)
要执行数组切片之间的所有成对比较,请使用标准的 NumPy 广播技术。例如,要计算所有行对之间的相关性:
>>> stats.pearsonr(x[:, np.newaxis, :], y, axis=-1).statistic.shape (8, 8)
如果 y = a + b*x + e,其中 a,b 是常数,e 是随机误差项,假设与 x 无关,则 x 和 y 之间存在线性依赖关系。为简单起见,假设 x 是标准正态分布,a=0,b=1,并让 e 服从均值为零且标准差 s>0 的正态分布。
>>> rng = np.random.default_rng() >>> s = 0.5 >>> x = stats.norm.rvs(size=500, random_state=rng) >>> e = stats.norm.rvs(scale=s, size=500, random_state=rng) >>> y = x + e >>> stats.pearsonr(x, y).statistic 0.9001942438244763
这应该接近于给定的确切值。
>>> 1/np.sqrt(1 + s**2) 0.8944271909999159
对于 s=0.5,我们观察到高度的相关性。通常,噪声的方差越大,相关性越低,而当误差的方差趋近于零时,相关性趋近于一。
重要的是要记住,如果没有相关性并不意味着独立,除非 (x, y) 是联合正态分布。即使在存在非常简单的依赖结构时,相关性也可能是零:如果 X 服从标准正态分布,令 y = abs(x)。注意 x 和 y 之间的相关性是零。实际上,由于 x 的期望值为零,cov(x, y) = E[x*y]。根据定义,这等于 E[x*abs(x)],由于对称性,这为零。以下代码行说明了这一观察结果:
>>> y = np.abs(x) >>> stats.pearsonr(x, y) PearsonRResult(statistic=-0.05444919272687482, pvalue=0.22422294836207743)
非零的相关系数可能会产生误导。例如,如果 X 服从标准正态分布,定义 y = x 如果 x < 0,否则 y = 0。一个简单的计算显示 corr(x, y) = sqrt(2/Pi) = 0.797…,这意味着高度的相关性:
>>> y = np.where(x < 0, x, 0) >>> stats.pearsonr(x, y) PearsonRResult(statistic=0.861985781588, pvalue=4.813432002751103e-149)
这并不直观,因为如果 x 大于零,x 和 y 之间就没有依赖关系,如果我们对 x 和 y 进行采样,这种情况大约会发生在一半的案例中。