Python SciPy mstats.pearsonr用法及代码示例

本文简要介绍 python 语言中 scipy.stats.mstats.pearsonr 的用法。

用法:

scipy.stats.mstats.pearsonr(x, y)#

用于测试非相关性的 Pearson 相关系数和 p 值。

皮尔逊相关系数 [1] 衡量两个数据集之间的线性关系。 p 值的计算依赖于每个数据集呈正态分布的假设。 (有关输入非正态性对相关系数分布的影响的讨论，请参阅 Kowalski [3]。)与其他相关系数一样，该系数在 -1 和 +1 之间变化，其中 0 表示不相关。 -1 或 +1 的相关性意味着精确的线性关系。

参数 ：：

x： (N,) 数组

输入数组。

y： (N,) 数组

输入数组。

返回 ：：

r： 浮点数

皮尔逊相关系数。

p-value： 浮点数

Two-tailed p 值。

警告 ：：

ConstantInputWarning

如果输入是常量数组，则引发。这种情况下没有定义相关系数，所以返回np.nan。

NearConstantInputWarning

如果输入是“nearly” 常量，则引发。如果 norm(x - mean(x)) < 1e-13 * abs(mean(x)) ，则数组 x 被认为几乎是恒定的。在这种情况下，计算x - mean(x) 中的数值错误可能会导致 r 的计算不准确。

注意：

相关系数计算如下：

\[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 值的分布。该分布是区间 [-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 值。 p 值粗略地表示不相关系统生成的数据集的 Pearson 相关性至少与根据这些数据集计算出的相关性一样极端的概率。更准确地说，对于相关系数为 r 的给定样本，p 值是从相关性为零的总体中抽取的随机样本 x' 和 y' 的 abs(r') 大于或等于 abs(r') 的概率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。

参考：

[1] (1,2)

“Pearson correlation coefficient”，维基百科，https://en.wikipedia.org/wiki/Pearson_correlation_coefficient

[2]

学生，“相关系数的可能误差”，Biometrika，第 6 卷，第 2-3 期，1908 年 9 月 1 日，第 302-310 页。

[3]

C. J. Kowalski，“关于非正态性对样本 Product-Moment 相关系数分布的影响”皇家统计学会杂志。系列 C(应用统计)，卷。 21，第 1 期(1972 年)，第 1-12 页。

例子：

>>> import numpy as np
>>> from scipy import stats
>>> from scipy.stats import mstats
>>> mstats.pearsonr([1, 2, 3, 4, 5], [10, 9, 2.5, 6, 4])
(-0.7426106572325057, 0.1505558088534455)

如果 y = a + b*x + e，则 x 和 y 之间存在线性相关性，其中 a,b 是常数，e 是随机误差项，假定与 x 无关。为简单起见，假设 x 是标准正态，a=0，b=1，让 e 服从均值为零且标准差 s>0 的正态分布。

>>> s = 0.5
>>> x = stats.norm.rvs(size=500)
>>> e = stats.norm.rvs(scale=s, size=500)
>>> y = x + e
>>> mstats.pearsonr(x, y)
(0.9029601878969703, 8.428978827629898e-185) # may vary

这应该接近给出的确切值

>>> 1/np.sqrt(1 + s**2)
0.8944271909999159

对于 s=0.5，我们观察到高度相关。一般来说，较大的噪声方差会降低相关性，而当误差的方差变为零时，相关性会接近 1。

重要的是要记住，除非 (x, y) 是共同正态的，否则没有相关性并不意味着独立。当存在非常简单的依赖结构时，相关性甚至可以为零：如果 X 遵循标准正态分布，则令 y = abs(x)。请注意，x 和 y 之间的相关性为零。实际上，由于 x 的期望为零，所以 cov(x, y) = E[x*y]。根据定义，这等于 E[x*abs(x)]，对称性为零。以下代码行说明了这一观察结果：

>>> y = np.abs(x)
>>> mstats.pearsonr(x, y)
(-0.016172891856853524, 0.7182823678751942) # may vary

非零相关系数可能会产生误导。例如，如果 X 具有标准正态分布，则如果 x < 0 则定义 y = x，否则定义 y = 0。一个简单的计算表明 corr(x, y) = sqrt(2/Pi) = 0.797…，意味着高度相关：

>>> y = np.where(x < 0, x, 0)
>>> mstats.pearsonr(x, y)
(0.8537091583771509, 3.183461621422181e-143) # may vary

这是不直观的，因为如果 x 大于零，则 x 和 y 没有依赖性，如果我们对 x 和 y 进行采样，大约有一半的情况会发生这种情况。