相关性分析(基于Python)6—— 皮尔逊相关系数的核心性质与边界约束

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子保证可复现
np.random.seed(42)

# 生成具有线性关系的原始数据
n = 500
x = np.random.normal(50, 10, n)
y = 0.6 * x + np.random.normal(0, 8, n)

# 计算两种顺序的相关系数
r_xy, _ = stats.pearsonr(x, y)
r_yx, _ = stats.pearsonr(y, x)

# 创建可视化画布
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# 绘制X-Y散点图
ax1.scatter(x, y, alpha=0.6, color='#377eb8', s=30, edgecolors='white')
slope1, intercept1, _, _, _ = stats.linregress(x, y)
x_fit1 = np.linspace(min(x), max(x), 100)
ax1.plot(x_fit1, slope1 * x_fit1 + intercept1, color='black', linestyle='--', linewidth=2)
ax1.set_title(f'X为横轴，Y为纵轴\nPearson r_xy = {r_xy:.4f}', fontsize=14)
ax1.set_xlabel('变量 X', fontsize=12)
ax1.set_ylabel('变量 Y', fontsize=12)
ax1.grid(True, linestyle=':', alpha=0.6)

# 绘制Y-X散点图
ax2.scatter(y, x, alpha=0.6, color='#e41a1c', s=30, edgecolors='white')
slope2, intercept2, _, _, _ = stats.linregress(y, x)
x_fit2 = np.linspace(min(y), max(y), 100)
ax2.plot(x_fit2, slope2 * x_fit2 + intercept2, color='black', linestyle='--', linewidth=2)
ax2.set_title(f'Y为横轴，X为纵轴\nPearson r_yx = {r_yx:.4f}', fontsize=14)
ax2.set_xlabel('变量 Y', fontsize=12)
ax2.set_ylabel('变量 X', fontsize=12)
ax2.grid(True, linestyle=':', alpha=0.6)

plt.tight_layout()
plt.show()

# 输出验证结果
print(f"Y与X的相关系数r_yx: {r_yx:.6f}")
print(f"两个相关系数是否相等: {np.isclose(r_xy, r_yx)}")

执行结果如下：

性质2：边界性

皮尔逊相关系数的取值被严格限制在区间内，即：

当且仅当  与  存在完全线性函数关系时，等号成立。其中 对应完全正线性相关， 对应完全负线性相关。

证明： 该性质可由柯西-施瓦茨不等式导出。对于任意两个具有有限二阶矩的随机变量  和 ，柯西-施瓦茨不等式给出：

将不等式两侧同时开平方，得到：

将不等式两侧同时除以右侧的项，恰好就是皮尔逊相关系数的定义式，因此可以直接得出 。

等号成立的充要条件是柯西-施瓦茨不等式的等号成立条件，即存在常数和，使得以概率1成立。当时，；当时，。

样本相关系数的边界性推导完全一致，只需将期望替换为样本求和，将总体方差替换为样本平方和，最终可以得到的结论。

我们通过代码可视化不同相关系数下的散点分布，直观展示相关系数从到的变化过程，以及完全线性相关的边界情况。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)
n = 300
x = np.random.normal(0, 1, n)

# 定义不同的相关系数目标值
rho_list = [-1, -0.8, -0.4, 0, 0.4, 0.8, 1]
fig, axes = plt.subplots(1, 7, figsize=(21, 3), sharey=True)

for ax, rho in zip(axes, rho_list):
# 生成对应相关系数的Y变量
if rho == 1:
        y = x
elif rho == -1:
        y = -x
else:
# 根据线性关系生成带噪声的Y，保证理论相关系数为rho
        y = rho * x + np.random.normal(0, np.sqrt(1 - rho**2), n)

# 计算实际样本相关系数
    r, _ = stats.pearsonr(x, y)

# 绘制散点图
    ax.scatter(x, y, alpha=0.6, s=20, edgecolors='white')
# 绘制回归线
    slope, intercept, _, _, _ = stats.linregress(x, y)
    x_fit = np.linspace(-3, 3, 100)
    ax.plot(x_fit, slope * x_fit + intercept, color='black', linestyle='--', linewidth=1.5)

    ax.set_title(f'ρ = {rho}\nr = {r:.2f}', fontsize=12)
    ax.set_xlim(-3.5, 3.5)
    ax.set_ylim(-3.5, 3.5)
    ax.grid(True, linestyle=':', alpha=0.4)
    ax.set_aspect('equal')

axes[0].set_ylabel('变量 Y', fontsize=12)
fig.supxlabel('变量 X', fontsize=12)
plt.tight_layout()
plt.show()

执行结果如下：

性质3：标准化变量的协方差等价性

若对  和  进行Z-Score标准化，得到标准化变量  和 ：

则标准化变量的协方差，与原始变量的皮尔逊相关系数完全相等：

证明： 首先计算标准化变量的期望与方差。根据期望的线性性质，

同理 。根据方差的性质，

同理 。

根据协方差的定义，

将  和  的表达式代入：

样本层面的推导完全一致，标准化后的样本变量  和 ，其样本协方差与样本相关系数完全相等。

这个性质表明皮尔逊相关系数是剔除了变量自身尺度与位置影响后，纯粹的线性协同波动程度。

我们通过代码验证这个等价性，同时对比原始数据与标准化数据的分布变化。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)

# 生成带量纲差异的原始数据
n = 400
# X为身高数据，单位厘米，均值170，标准差8
x = np.random.normal(170, 8, n)
# Y为体重数据，单位千克，均值65，标准差10
y = 0.7 * x + np.random.normal(0, 7, n) - 54

# 计算原始数据的相关系数与协方差
r_original, _ = stats.pearsonr(x, y)
Cov_original = np.cov(x, y, ddof=1)[0, 1]

# 对数据进行Z-Score标准化
z_x = (x - np.mean(x)) / np.std(x, ddof=1)
z_y = (y - np.mean(y)) / np.std(y, ddof=1)

# 计算标准化数据的相关系数与协方差
r_standardized, _ = stats.pearsonr(z_x, z_y)
Cov_standardized = np.cov(z_x, z_y, ddof=1)[0, 1]

# 创建可视化画布
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# 绘制原始数据散点图
ax1.scatter(x, y, alpha=0.6, color='#377eb8', s=30, edgecolors='white')
ax1.set_title(f'原始数据分布\n协方差 = {Cov_original:.2f}\n相关系数 r = {r_original:.4f}', fontsize=14)
ax1.set_xlabel('身高 (cm)', fontsize=12)
ax1.set_ylabel('体重 (kg)', fontsize=12)
ax1.grid(True, linestyle=':', alpha=0.6)

# 绘制标准化数据散点图
ax2.scatter(z_x, z_y, alpha=0.6, color='#ff7f00', s=30, edgecolors='white')
ax2.axhline(0, color='black', lw=0.8)
ax2.axvline(0, color='black', lw=0.8)
ax2.set_title(f'标准化数据分布\n协方差 = {Cov_standardized:.4f}\n相关系数 r = {r_standardized:.4f}', fontsize=14)
ax2.set_xlabel('标准化身高 Z_X', fontsize=12)
ax2.set_ylabel('标准化体重 Z_Y', fontsize=12)
ax2.grid(True, linestyle=':', alpha=0.6)
ax2.set_aspect('equal')

plt.tight_layout()
plt.show()

# 输出验证结果
print(f"原始数据协方差: {Cov_original:.4f}")
print(f"原始数据相关系数: {r_original:.6f}")
print(f"标准化数据协方差: {Cov_standardized:.6f}")
print(f"标准化数据相关系数: {r_standardized:.6f}")
print(f"标准化协方差与原始相关系数是否相等: {np.isclose(Cov_standardized, r_original)}")

执行结果如下：

性质4：独立变量的相关系数为0

若随机变量  和  统计独立，则两者的总体皮尔逊相关系数。

证明： 根据统计独立性的定义，若  和  相互独立，则 。

将其代入协方差的展开式：

协方差为0，因此皮尔逊相关系数。

需要特别强调的是，这个性质的逆命题并不成立。相关系数为0仅能说明两个变量不存在线性依赖关系，无法证明它们相互独立，变量之间可能存在极强的非线性依赖关系，同时保持相关系数为0。

我们通过代码模拟两个独立的随机变量，验证其相关系数的分布特征，同时展示独立变量的联合分布形态。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)

# 生成两个完全独立的随机变量
n = 10000
# X服从标准正态分布
x = np.random.normal(0, 1, n)
# Y服从均匀分布，与X完全独立
y = np.random.uniform(-3, 3, n)

# 计算样本相关系数
r, p_value = stats.pearsonr(x, y)

# 重复抽样实验，展示相关系数的抽样分布
n_experiments = 5000
sample_size = 200
r_distribution = np.zeros(n_experiments)

for i in range(n_experiments):
    x_sample = np.random.normal(0, 1, sample_size)
    y_sample = np.random.uniform(-3, 3, sample_size)
    r_sample, _ = stats.pearsonr(x_sample, y_sample)
    r_distribution[i] = r_sample

# 创建可视化画布
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# 绘制独立变量的联合散点图
ax1.scatter(x, y, alpha=0.3, s=10, color='#4daf4a')
ax1.set_title(f'两个独立随机变量的联合分布\n总体相关系数 ρ=0\n样本相关系数 r={r:.4f}\np值={p_value:.4f}', fontsize=13)
ax1.set_xlabel('X (标准正态分布)', fontsize=12)
ax1.set_ylabel('Y (均匀分布)', fontsize=12)
ax1.grid(True, linestyle=':', alpha=0.6)

# 绘制相关系数的抽样分布直方图
ax2.hist(r_distribution, bins=50, density=True, alpha=0.7, color='#984ea3', edgecolor='white')
# 绘制正态分布拟合曲线
mean_r = np.mean(r_distribution)
std_r = np.std(r_distribution)
x_grid = np.linspace(-0.5, 0.5, 200)
pdf = stats.norm.pdf(x_grid, loc=mean_r, scale=std_r)
ax2.plot(x_grid, pdf, 'k--', linewidth=2, label='正态分布拟合')
# 标注总体相关系数0的位置
ax2.axvline(0, color='red', linestyle='-', linewidth=2, label='总体相关系数 ρ=0')
ax2.axvline(mean_r, color='black', linestyle='--', linewidth=2, label=f'样本均值 E[r]={mean_r:.4f}')

ax2.set_title(f'独立变量样本相关系数的抽样分布\n(样本量n={sample_size}, 重复实验{n_experiments}次)', fontsize=13)
ax2.set_xlabel('样本相关系数 r', fontsize=12)
ax2.set_ylabel('概率密度', fontsize=12)
ax2.legend()
ax2.grid(True, linestyle=':', alpha=0.6)

plt.tight_layout()
plt.show()

执行结果如下：

2线性变换下的不变性与符号规则

皮尔逊相关系数最具实用价值的性质之一，是它对变量的线性变换具有严格的不变性。

我们先定义变量的线性变换。对于随机变量  和 ，我们构造两个新的随机变量  和 ：

其中  为任意常数，且 （若  或  为，变换后的变量为常数，方差为，无法计算相关系数）。

核心结论是：对变量进行任意线性变换，皮尔逊相关系数的绝对值保持不变，符号仅由两个变量的尺度变换系数的乘积符号决定。平移变换对相关系数无任何影响。

我们分两步完成这个结论的严谨推导。

性质1：平移变换的绝对不变性

当尺度系数 时，线性变换退化为纯平移变换：

此时可以证明，变换后的相关系数与原始相关系数完全相等：

证明： 首先计算  和  的期望：

去中心化后的偏差项为：

可以看到，平移变换完全不改变变量去中心化后的偏差项。因此协方差保持不变：

同时，平移变换也不改变变量的方差：

协方差与方差均未发生变化，因此相关系数。

这个性质的实践意义极强。它说明，对变量进行任何平移操作，比如将摄氏度转换为开尔文温度（加273.15）、将日期转换为距离某一天的天数、对收入数据进行常数项调整，都完全不会改变变量之间的皮尔逊相关系数。

性质2：尺度变换的符号规则与绝对值不变性

当平移系数  时，线性变换退化为纯尺度变换：

此时变换后的相关系数为：

其中为符号函数，当括号内的数值大于时返回，小于时返回。

证明： 首先计算  和  的期望与去中心化偏差：

计算协方差：

计算方差：

将协方差与标准差代入相关系数定义式：

上述结果表明：尺度变换不会改变相关系数的绝对值，即无论我们将变量放大或缩小多少倍，相关系数的绝对值始终保持不变；尺度变换仅会改变相关系数的符号，即当两个变量的尺度系数乘积为正（同号）时，相关系数符号不变；当乘积为负（异号）时，相关系数符号反转。

完整线性变换的综合性质

结合平移变换与尺度变换的性质，我们可以得到完整线性变换下的最终结论。对于任意线性变换 ，变换后的相关系数为：

平移系数  和  完全不影响最终的相关系数，只有尺度系数 和  会影响相关系数的符号，且不会改变其绝对值。

这个性质直接解释了为什么相关系数是无量纲的统计量，因此我们不需要考虑变量的量纲差异，就可以直接比较不同变量对之间的相关强度。

我们通过代码模拟三种典型的线性变换场景，可视化变换前后的散点分布与相关系数变化，验证上述性质。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)

# 生成原始数据
n = 400
x = np.random.normal(10, 2, n)
y = 0.7 * x + np.random.normal(0, 1.5, n)
r_original, _ = stats.pearsonr(x, y)

# 场景1：纯平移变换
u1 = x + 50# X加50
v1 = y - 20# Y减20
r1, _ = stats.pearsonr(u1, v1)

# 场景2：同号尺度变换
u2 = 10 * x  # X放大10倍
v2 = 2 * y   # Y放大2倍
r2, _ = stats.pearsonr(u2, v2)

# 场景3：异号尺度变换
u3 = 5 * x   # X放大5倍
v3 = -3 * y  # Y放大3倍并反转符号
r3, _ = stats.pearsonr(u3, v3)

# 创建可视化画布
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
axes = axes.flatten()

# 绘制原始数据
axes[0].scatter(x, y, alpha=0.6, color='#377eb8', s=30, edgecolors='white')
slope, intercept, _, _, _ = stats.linregress(x, y)
x_fit = np.linspace(min(x), max(x), 100)
axes[0].plot(x_fit, slope * x_fit + intercept, color='black', linestyle='--', linewidth=2)
axes[0].set_title(f'原始数据\nPearson r = {r_original:.4f}', fontsize=14)
axes[0].set_xlabel('原始 X', fontsize=12)
axes[0].set_ylabel('原始 Y', fontsize=12)
axes[0].grid(True, linestyle=':', alpha=0.6)

# 绘制纯平移变换数据
axes[1].scatter(u1, v1, alpha=0.6, color='#4daf4a', s=30, edgecolors='white')
slope1, intercept1, _, _, _ = stats.linregress(u1, v1)
x_fit1 = np.linspace(min(u1), max(u1), 100)
axes[1].plot(x_fit1, slope1 * x_fit1 + intercept1, color='black', linestyle='--', linewidth=2)
axes[1].set_title(f'纯平移变换 U=X+50, V=Y-20\nPearson r = {r1:.4f}', fontsize=14)
axes[1].set_xlabel('平移后的 U', fontsize=12)
axes[1].set_ylabel('平移后的 V', fontsize=12)
axes[1].grid(True, linestyle=':', alpha=0.6)

# 绘制同号尺度变换数据
axes[2].scatter(u2, v2, alpha=0.6, color='#984ea3', s=30, edgecolors='white')
slope2, intercept2, _, _, _ = stats.linregress(u2, v2)
x_fit2 = np.linspace(min(u2), max(u2), 100)
axes[2].plot(x_fit2, slope2 * x_fit2 + intercept2, color='black', linestyle='--', linewidth=2)
axes[2].set_title(f'同号尺度变换 U=10X, V=2Y\nPearson r = {r2:.4f}', fontsize=14)
axes[2].set_xlabel('尺度变换后的 U', fontsize=12)
axes[2].set_ylabel('尺度变换后的 V', fontsize=12)
axes[2].grid(True, linestyle=':', alpha=0.6)

# 绘制异号尺度变换数据
axes[3].scatter(u3, v3, alpha=0.6, color='#e41a1c', s=30, edgecolors='white')
slope3, intercept3, _, _, _ = stats.linregress(u3, v3)
x_fit3 = np.linspace(min(u3), max(u3), 100)
axes[3].plot(x_fit3, slope3 * x_fit3 + intercept3, color='black', linestyle='--', linewidth=2)
axes[3].set_title(f'异号尺度变换 U=5X, V=-3Y\nPearson r = {r3:.4f}', fontsize=14)
axes[3].set_xlabel('尺度变换后的 U', fontsize=12)
axes[3].set_ylabel('尺度变换后的 V', fontsize=12)
axes[3].grid(True, linestyle=':', alpha=0.6)

plt.tight_layout()
plt.show()

# 输出验证结果
print(f"原始数据相关系数: {r_original:.6f}")
print(f"纯平移变换后相关系数: {r1:.6f} (与原始相等: {np.isclose(r1, r_original)})")
print(f"同号尺度变换后相关系数: {r2:.6f} (与原始相等: {np.isclose(r2, r_original)})")
print(f"异号尺度变换后相关系数: {r3:.6f} (与原始符号相反、绝对值相等: {np.isclose(np.abs(r3), np.abs(r_original))})")

执行结果如下：

3几何性质与线性回归的内在关联

3.1 相关系数与向量夹角的严格对应

在容量为  的样本空间中，我们可以将两个变量的观测值映射为n维欧几里得空间中的两个向量：

这两个向量是原始数据去中心化后的结果，向量的起点位于样本空间的原点，每个分量对应一个样本的观测偏差。

根据线性代数的定义，两个向量的内积等于向量模长的乘积乘以向量夹角的余弦值：

其中  为两个向量在  维空间中的夹角，取值范围为。

我们将内积与模长的表达式展开：

将上述表达式代入内积公式，两侧同时除以，可以得到：

上述结果表明：样本皮尔逊相关系数，等于n维样本空间中两个去中心化数据向量夹角的余弦值。

当 时，，两个向量完全共线且同向，对应完全正线性相关；当 时，，两个向量在  维空间中相互正交，对应无线性相关关系；当 时，，两个向量完全共线且反向，对应完全负线性相关。

3.2 相关系数与线性回归系数的代数关系

一元线性回归中，我们试图用X的线性函数来拟合 ，拟合方程为：

其中  为  对  的回归系数，代表  每变化一个单位， 的平均变化量； 为截距项。

根据最小二乘法，回归系数   的最优估计值为：

我们将这个公式与皮尔逊相关系数的定义式结合，可以推导出回归系数与相关系数的内在关联：

其中  和  分别为  和  的样本标准差。

同理，若我们用  来拟合 ，得到  对  的回归系数：

从这两个公式可以得出以下结论：一是回归系数的符号与相关系数的符号完全一致；二是回归系数的大小由相关系数和两个变量的标准差之比共同决定，相关系数反映了线性关系的紧密程度，标准差之比则反映了两个变量的尺度差异；三是两个方向的回归系数的乘积，等于相关系数的平方：

这个结论表明相关系数它是两个方向回归系数的几何均值。

3.3 相关系数的平方与决定系数

在线性回归中，决定系数  是评估模型拟合效果的核心指标，它代表了因变量的总方差中，能够被自变量的线性关系解释的比例。

我们可以将Y的总方差拆解为两部分：

等式左侧为总平方和SST，代表  的总波动；右侧第一项为回归平方和SSR，代表模型能够解释的波动；第二项为残差平方和SSE，代表模型无法解释的随机波动。

决定系数  的定义为：

我们可以通过严格的代数推导证明，在一元线性回归中，决定系数  恰好等于皮尔逊相关系数的平方：

证明： 根据回归方程的拟合值，可以得到。 将其代入回归平方和SSR：

将代入，同时，即：

而总平方和SST = 。 因此：

这个推导表明：相关系数的平方，代表了两个变量之间线性关系能够解释的方差比例。

当  时，，说明  的方差中有  可以被  的线性关系解释；当  时，，说明仅有  的方差可以被解释，即使相关系数在统计上显著，其实际的解释能力也非常有限。

我们通过代码可视化相关系数、回归系数、决定系数三者之间的关系，直观展示不同相关系数下的回归拟合效果与方差解释率。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)
n = 300
x = np.random.normal(0, 1, n)

# 定义不同的相关系数目标值
rho_list = [0.3, 0.5, 0.7, 0.9]
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
axes = axes.flatten()

for ax, rho in zip(axes, rho_list):
# 生成对应相关系数的Y变量
    y = rho * x + np.random.normal(0, np.sqrt(1 - rho**2), n)

# 计算相关系数、回归系数、决定系数
    r, _ = stats.pearsonr(x, y)
    slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
    r_squared = r_value ** 2

# 生成拟合线与置信区间
    x_fit = np.linspace(-3, 3, 100)
    y_fit = slope * x_fit + intercept

# 绘制散点图
    ax.scatter(x, y, alpha=0.5, s=30, edgecolors='white', label='观测数据')
# 绘制回归线
    ax.plot(x_fit, y_fit, color='red', linestyle='-', linewidth=3, label=f'回归直线 y={slope:.2f}x+{intercept:.2f}')
# 绘制均值基准线
    ax.axhline(np.mean(y), color='black', linestyle=':', linewidth=2, label='Y的均值')
    ax.axvline(0, color='black', linestyle=':', linewidth=1)

# 添加统计信息文本
    info_text = (
f"皮尔逊相关系数 r = {r:.4f}\n"
f"决定系数 R² = {r_squared:.4f}\n"
f"可解释方差比例 = {r_squared*100:.2f}%"
    )
    ax.text(0.05, 0.95, info_text, transform=ax.transAxes, fontsize=12,
            verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))

    ax.set_title(f'相关系数 r = {rho}', fontsize=14)
    ax.set_xlabel('自变量 X', fontsize=12)
    ax.set_ylabel('因变量 Y', fontsize=12)
    ax.set_xlim(-3.5, 3.5)
    ax.set_ylim(-3.5, 3.5)
    ax.grid(True, linestyle=':', alpha=0.6)
    ax.legend(loc='lower right')

plt.tight_layout()
plt.show()

执行结果如下：

4样本相关系数的抽样性质与偏差校正

性质1：一致性

样本相关系数  是总体相关系数  的一致估计量。即当样本量  趋向于无穷大时， 依概率收敛于：

证明： 根据大数定律，样本协方差依概率收敛于总体协方差；样本标准差、依概率收敛于总体标准差、。

根据连续映射定理，随机变量序列的连续函数的极限，等于随机变量极限的连续函数。皮尔逊相关系数是协方差与标准差的连续函数，因此：

这个性质表明，只要样本量足够大，样本相关系数会无限接近总体的真实相关系数。无论总体相关系数的真实值是多少，随着样本量的增加，抽样误差会逐渐缩小，估计结果会越来越稳定。

性质2：有偏性

样本相关系数  是总体相关系数  的有偏估计量。即样本相关系数的期望不等于总体真实相关系数：

偏差的来源与样本量和总体相关系数的大小密切相关。当样本量较小时，偏差会非常明显，尤其是当  时，偏差会急剧增大；当 时，样本相关系数  的期望 ，此时估计是无偏的；当  时，样本相关系数的绝对值会系统性地小于总体相关系数的绝对值，即 ；随着样本量  的增大，偏差会逐渐趋近于，这与一致性的结论一致。

我们可以通过模拟实验直观展示这种偏差的存在，以及样本量和总体相关系数对偏差大小的影响。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import time

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)

# 实验参数设置
n_experiments = 5000# 减少到5000次实验（仍具有统计意义）
rho_list = [0, 0.3, 0.6, 0.9]  # 总体相关系数
sample_sizes = [5, 10, 20, 50, 100, 500]  # 不同样本量

print("开始模拟...")
start_time = time.time()

# 存储结果
bias_results = np.zeros((len(rho_list), len(sample_sizes)))

# 进行模拟实验 - 内存友好版本
for i, rho in enumerate(rho_list):
    print(f"正在处理 ρ={rho}...")
for j, n in enumerate(sample_sizes):
        r_values = np.zeros(n_experiments)

# 方式1：直接循环，每次只生成当前实验的数据（内存友好）
for k in range(n_experiments):
# 生成二元正态分布数据（使用更简单的方法）
            x = np.random.randn(n)
# 使用条件分布生成相关数据：y = ρ*x + √(1-ρ²)*ε
            epsilon = np.random.randn(n)
            y = rho * x + np.sqrt(1 - rho ** 2) * epsilon

# 计算相关系数
            r_values[k] = np.corrcoef(x, y)[0, 1]

# 计算偏差
        mean_r = np.mean(r_values)
        bias = mean_r - rho
        bias_results[i, j] = bias
        print(f"  样本量 n={n}: 偏差={bias:.6f}")

end_time = time.time()
print(f"\n模拟完成！总耗时: {end_time - start_time:.2f} 秒")

# 可视化偏差结果
fig, ax = plt.subplots(figsize=(12, 7))

for i, rho in enumerate(rho_list):
    ax.plot(sample_sizes, bias_results[i, :], marker='o', linewidth=2,
            markersize=8, label=f'总体ρ={rho}')

# 绘制零偏差基准线
ax.axhline(0, color='black', linestyle='--', linewidth=1.5)

ax.set_title('样本相关系数的偏差与样本量、总体ρ的关系', fontsize=16)
ax.set_xlabel('样本量 n', fontsize=14)
ax.set_ylabel('偏差 E[r] - ρ', fontsize=14)
ax.set_xscale('log')
ax.grid(True, linestyle=':', alpha=0.6)
ax.legend(fontsize=12)
plt.tight_layout()
plt.show()

# 输出偏差数值表
print("\n样本相关系数的偏差表 (行：总体ρ，列：样本量n)")
print("样本量:", sample_sizes)
print(np.round(bias_results, 6))

执行结果如下：

偏差校正方法

在小样本场景下，样本相关系数的偏差会对分析结果产生显著影响，因此需要对r进行偏差校正，得到更接近总体真实ρ的估计值。

目前应用最广泛的是Olkin & Pratt 提出的近似无偏校正公式：

这个公式适用于  的场景，能够有效减小小样本下的偏差。

另一种常用的校正方法是基于Fisher Z变换的校正。我们先对  进行Fisher Z变换得到 ，对  进行偏差校正后，再通过反变换得到校正后的相关系数。校正公式为：

我们通过模拟实验对比两种校正方法的效果，验证它们对偏差的改善程度。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)

# 实验参数
rho_true = 0.7
sample_size = 10
n_experiments = 10000

# 存储结果
r_raw = np.zeros(n_experiments)
r_olkin = np.zeros(n_experiments)
r_fisher = np.zeros(n_experiments)

# 模拟实验
for i in range(n_experiments):
# 生成数据
    cov_matrix = [[1, rho_true], [rho_true, 1]]
    data = np.random.multivariate_normal([0, 0], cov_matrix, size = sample_size)
    r, _ = stats.pearsonr(data[:, 0], data[:, 1])

# 原始相关系数
    r_raw[i] = r

# Olkin & Pratt校正
if sample_size > 4:
        r_olkin[i] = r * (1 + (1 - r ** 2) / (2 * (sample_size - 4)))
else:
        r_olkin[i] = r

# Fisher Z变换校正
    z_r = np.arctanh(np.clip(r, -0.99999, 0.99999))
    z_corrected = z_r - r / (2 * (sample_size - 1))
    r_fisher[i] = np.tanh(z_corrected)

# 计算偏差与均方误差
bias_raw = np.mean(r_raw) - rho_true
bias_olkin = np.mean(r_olkin) - rho_true
bias_fisher = np.mean(r_fisher) - rho_true

mse_raw = np.mean((r_raw - rho_true) ** 2)
mse_olkin = np.mean((r_olkin - rho_true) ** 2)
mse_fisher = np.mean((r_fisher - rho_true) ** 2)

# 可视化分布对比
fig, ax = plt.subplots(figsize=(12, 7))

# 绘制直方图
ax.hist(r_raw, bins=50, density=True, alpha=0.4, color='#377eb8',
        label=f'原始r\n偏差={bias_raw:.4f}, MSE={mse_raw:.4f}')
ax.hist(r_olkin, bins=50, density=True, alpha=0.4, color='#4daf4a',
        label=f'Olkin校正\n偏差={bias_olkin:.4f}, MSE={mse_olkin:.4f}')
ax.hist(r_fisher, bins=50, density=True, alpha=0.4, color='#984ea3',
        label=f'Fisher校正\n偏差={bias_fisher:.4f}, MSE={mse_fisher:.4f}')

# 绘制真实ρ的基准线
ax.axvline(rho_true, color='red', linestyle='-', linewidth=2, label=f'真实总体ρ={rho_true}')
# 绘制各方法的均值
ax.axvline(np.mean(r_raw), color='#377eb8', linestyle='--', linewidth=1.5)
ax.axvline(np.mean(r_olkin), color='#4daf4a', linestyle='--', linewidth=1.5)
ax.axvline(np.mean(r_fisher), color='#984ea3', linestyle='--', linewidth=1.5)

ax.set_title(f'小样本下相关系数校正效果对比 (n={sample_size}, 真实ρ={rho_true})', fontsize=16)
ax.set_xlabel('相关系数估计值', fontsize=14)
ax.set_ylabel('概率密度', fontsize=14)
ax.legend(fontsize=12)
ax.grid(True, linestyle=':', alpha=0.6)
plt.tight_layout()
plt.show()

执行结果如下：

性质3：渐近正态性

当样本量  足够大时，样本相关系数  的抽样分布近似服从正态分布：

其中表示依分布收敛。

需要特别强调的是，这个渐近正态性仅在大样本下成立。当样本量较小，或者  时， 的抽样分布会呈现出强烈的偏态，完全不符合正态分布，此时须使用Fisher Z变换进行分布正态化。

我们通过代码对比不同样本量下， 的抽样分布与渐近正态分布的拟合效果，验证渐近正态性的适用条件。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)

# 实验参数
rho_true = 0.6
sample_sizes = [10, 50, 200, 1000]
n_experiments = 10000

# 创建可视化画布
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
axes = axes.flatten()

for ax, n in zip(axes, sample_sizes):
# 模拟生成样本相关系数
    r_values = np.zeros(n_experiments)
for i in range(n_experiments):
        cov_matrix = [[1, rho_true], [rho_true, 1]]
        data = np.random.multivariate_normal([0, 0], cov_matrix, size = n)
        r, _ = stats.pearsonr(data[:, 0], data[:, 1])
        r_values[i] = r

# 计算渐近正态分布的参数
    mean_asym = rho_true
    std_asym = (1 - rho_true ** 2) / np.sqrt(n - 1)

# 绘制直方图
    ax.hist(r_values, bins=50, density=True, alpha=0.7, color='#377eb8', edgecolor='white', label='r的经验分布')
# 绘制渐近正态分布曲线
    x_grid = np.linspace(-0.2, 1.0, 200)
    pdf_asym = stats.norm.pdf(x_grid, loc=mean_asym, scale=std_asym)
    ax.plot(x_grid, pdf_asym, 'k--', linewidth=2, label='渐近正态分布')
# 绘制真实ρ的基准线
    ax.axvline(rho_true, color='red', linestyle='-', linewidth=2, label=f'真实ρ={rho_true}')

# 计算Q-Q图的R²，评估正态拟合度
    _, (_, _, r_squared) = stats.probplot(r_values, dist="norm", fit=True)

# 添加统计信息
    info_text = (
f"样本量 n = {n}\n"
f"r的均值 = {np.mean(r_values):.4f}\n"
f"r的标准差 = {np.std(r_values):.4f}\n"
f"渐近标准差 = {std_asym:.4f}\n"
f"Q-Q图 R² = {r_squared:.4f}"
    )
    ax.text(0.05, 0.95, info_text, transform=ax.transAxes, fontsize=11,
            verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))

    ax.set_title(f'样本量n={n}时r的抽样分布', fontsize=14)
    ax.set_xlabel('样本相关系数 r', fontsize=12)
    ax.set_ylabel('概率密度', fontsize=12)
    ax.legend()
    ax.grid(True, linestyle=':', alpha=0.6)

plt.tight_layout()
plt.show()

执行结果如下：

5异常值与杠杆点的扰动性质

异常值与杠杆点的定义

在相关性分析中，我们将对相关系数产生强扰动的点分为两类：一是异常值，即指在因变量  方向上远离主体数据的点，即该点的  值与基于  的预测值之间存在极大的残差；二是杠杆点，即指在自变量  方向上远离主体数据的点，即该点的  值远大于或远小于其他样本的  值。高杠杆点即使残差很小，也会对相关系数和回归结果产生极强的影响。

单个数据点的影响力量化

我们可以通过影响函数来量化单个数据点对样本相关系数的影响力。具体来说，我们计算剔除第  个数据点后，样本相关系数的变化量，以此衡量该点的影响力：

其中表示剔除第  个数据点后，剩余样本计算得到的相关系数。

影响力的绝对值越大，说明该点对相关系数的扰动越强。当影响力的符号与原相关系数相反时，说明该点会削弱原有的相关关系；当符号相同时，说明该点会强化原有的相关关系。

在极端情况下，单个高杠杆点甚至可以完全决定相关系数的大小和符号，让原本无线性相关的数据呈现出极强的相关关系，或者让原本正相关的数据变成负相关。

我们通过代码模拟四种典型场景，可视化异常值与杠杆点对相关系数的扰动效果，同时量化每个点的影响力。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)
n_base = 30

# 生成基础数据：无相关的随机数据
x_base = np.random.normal(0, 1, n_base)
y_base = np.random.normal(0, 1, n_base)
r_base, _ = stats.pearsonr(x_base, y_base)

# 场景1：原始无相关数据
x1 = x_base.copy()
y1 = y_base.copy()
r1, _ = stats.pearsonr(x1, y1)

# 场景2：加入Y方向的单个异常值
x2 = np.append(x_base, 0)
y2 = np.append(y_base, 8)
r2, _ = stats.pearsonr(x2, y2)

# 场景3：加入X方向的单个高杠杆点（正相关）
x3 = np.append(x_base, 10)
y3 = np.append(y_base, 10)
r3, _ = stats.pearsonr(x3, y3)

# 场景4：加入X方向的单个高杠杆点（负相关）
x4 = np.append(x_base, 10)
y4 = np.append(y_base, -10)
r4, _ = stats.pearsonr(x4, y4)

# 计算每个场景中异常点的影响力
defcalculate_influence(x, y):
    n = len(x)
    r_full, _ = stats.pearsonr(x, y)
    influence = np.zeros(n)
for i in range(n):
        x_drop = np.delete(x, i)
        y_drop = np.delete(y, i)
        r_drop, _ = stats.pearsonr(x_drop, y_drop)
        influence[i] = r_full - r_drop
return influence

inf1 = calculate_influence(x1, y1)
inf2 = calculate_influence(x2, y2)
inf3 = calculate_influence(x3, y3)
inf4 = calculate_influence(x4, y4)

# 创建可视化画布
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
axes = axes.flatten()
scenarios = [
    (axes[0], x1, y1, r1, inf1, '场景1：原始无相关数据', '#377eb8'),
    (axes[1], x2, y2, r2, inf2, '场景2：加入Y方向异常值', '#e41a1c'),
    (axes[2], x3, y3, r3, inf3, '场景3：加入正相关高杠杆点', '#4daf4a'),
    (axes[3], x4, y4, r4, inf4, '场景4：加入负相关高杠杆点', '#984ea3')
]

for ax, x, y, r, inf, title, color in scenarios:
# 绘制主体数据点
    ax.scatter(x[:-1], y[:-1], alpha=0.7, s=50, edgecolors='white', label='主体数据')
# 绘制异常/杠杆点
if len(x) > n_base:
        ax.scatter(x[-1], y[-1], color=color, s=150, marker='*', edgecolors='black', zorder=5, label='异常/杠杆点')
# 绘制回归线
    slope, intercept, _, _, _ = stats.linregress(x, y)
    x_fit = np.linspace(min(x)-1, max(x)+1, 100)
    ax.plot(x_fit, slope * x_fit + intercept, color='black', linestyle='--', linewidth=2, label='整体回归线')

# 添加统计信息
    info_text = (
f"整体相关系数 r = {r:.4f}\n"
f"剔除异常点后 r = {r - inf[-1]:.4f}\n"
f"异常点影响力 = {inf[-1]:.4f}"
    )
    ax.text(0.05, 0.95, info_text, transform=ax.transAxes, fontsize=11,
            verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))

    ax.set_title(title, fontsize=14)
    ax.set_xlabel('X', fontsize=12)
    ax.set_ylabel('Y', fontsize=12)
    ax.grid(True, linestyle=':', alpha=0.6)
    ax.legend()

plt.tight_layout()
plt.show()

执行结果如下：

从模拟结果可以清晰地看到，单个高杠杆点可以让原本相关系数接近0的无相关数据，瞬间变成相关系数超过  的强相关数据，甚至可以反转相关系数的符号。这种极强的扰动性，是皮尔逊相关系数最致命的缺陷之一。

在实际数据分析中，我们必须在计算相关系数之前，通过散点图检查数据中是否存在异常值与高杠杆点。如果存在，需要根据业务逻辑判断该点是否为有效数据，若为无效数据则需要剔除，若为有效数据则需要使用更稳健的相关性度量方法。

6分组数据与聚合效应下的相关性质

全相关定律与辛普森悖论的数学条件

我们可以将总体相关系数拆解为组内效应与组间效应的组合，这个拆解基于全协方差定律与全方差定律。

假设我们根据离散分组变量 ，将总体分为  个组，每个组的权重为，即第  组的样本占比。

根据全协方差定律，总体协方差可以拆解为：

等式右侧第一项为组内协方差的加权平均，代表控制了分组变量Z后， 和  之间真实的组内线性关系；第二项为组间均值的协方差，代表不同分组之间， 的组均值与  的组均值的协同波动程度。

同理，根据全方差定律， 和  的总体方差可以拆解为：

总体皮尔逊相关系数是总体协方差与总体标准差乘积的比值，因此它由组内协方差、组间协方差、组内方差、组间方差共同决定。

辛普森悖论发生的数学条件： 当组内协方差的加权平均的符号，与组间均值的协方差的符号相反，且组间协方差的绝对值足够大，足以压倒组内协方差的加权平均时，就会发生辛普森悖论。

具体来说，若每个组内X和Y都呈负相关，即，但组间均值的协方差为正且绝对值足够大，导致总体协方差 ，最终总体相关系数为正，与组内相关系数的符号完全相反。

我们通过代码模拟经典的辛普森悖论场景，可视化分组数据与全量数据的相关关系差异，验证上述数学条件。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)
n_per_group = 100

# 生成3个分组的数据，组内均为负相关
# 组1：低分组
x1 = np.random.normal(2, 0.8, n_per_group)
y1 = -0.8 * x1 + np.random.normal(12, 0.8, n_per_group)
r1, _ = stats.pearsonr(x1, y1)

# 组2：中分组
x2 = np.random.normal(5, 0.8, n_per_group)
y2 = -0.8 * x2 + np.random.normal(20, 0.8, n_per_group)
r2, _ = stats.pearsonr(x2, y2)

# 组3：高分组
x3 = np.random.normal(8, 0.8, n_per_group)
y3 = -0.8 * x3 + np.random.normal(28, 0.8, n_per_group)
r3, _ = stats.pearsonr(x3, y3)

# 聚合全量数据
x_all = np.concatenate([x1, x2, x3])
y_all = np.concatenate([y1, y2, y3])
r_all, _ = stats.pearsonr(x_all, y_all)

# 计算全协方差拆解
# 组内协方差的加权平均
cov_within = (np.cov(x1, y1)[0,1] + np.cov(x2, y2)[0,1] + np.cov(x3, y3)[0,1]) / 3
# 组间均值的协方差
group_means_x = [np.mean(x1), np.mean(x2), np.mean(x3)]
group_means_y = [np.mean(y1), np.mean(y2), np.mean(y3)]
cov_between = np.cov(group_means_x, group_means_y)[0,1]
# 总体协方差
cov_total = np.cov(x_all, y_all)[0,1]

# 可视化
fig, ax = plt.subplots(figsize=(12, 8))

# 绘制各组散点与回归线
groups = [
    ('低分组', x1, y1, r1, '#377eb8'),
    ('中分组', x2, y2, r2, '#4daf4a'),
    ('高分组', x3, y3, r3, '#984ea3')
]

for label, x, y, r, color in groups:
    ax.scatter(x, y, alpha=0.7, color=color, s=50, edgecolors='white', label=f'{label} (组内r={r:.2f})')
    slope, intercept, _, _, _ = stats.linregress(x, y)
    x_fit = np.linspace(min(x), max(x), 100)
    ax.plot(x_fit, slope * x_fit + intercept, color=color, linewidth=3)

# 绘制整体回归线
slope_all, intercept_all, _, _, _ = stats.linregress(x_all, y_all)
x_fit_all = np.linspace(min(x_all), max(x_all), 100)
ax.plot(x_fit_all, slope_all * x_fit_all + intercept_all, color='black', linestyle='--', linewidth=3, label=f'整体回归线 (整体r={r_all:.2f})')

# 绘制各组均值点
ax.scatter(group_means_x, group_means_y, color='red', marker='X', s=200, edgecolors='black', zorder=5, label='各组均值点')

# 添加协方差拆解信息
info_text = (
f"组内协方差加权平均 = {cov_within:.4f}\n"
f"组间均值协方差 = {cov_between:.4f}\n"
f"总体协方差 = {cov_total:.4f}\n"
f"组内相关系数均为负，整体相关系数为正，辛普森悖论成立"
)
ax.text(0.05, 0.95, info_text, transform=ax.transAxes, fontsize=12,
        verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))

ax.set_title('辛普森悖论：分组负相关 vs 整体正相关', fontsize=16)
ax.set_xlabel('X', fontsize=14)
ax.set_ylabel('Y', fontsize=14)
ax.grid(True, linestyle=':', alpha=0.6)
ax.legend(fontsize=12)
plt.tight_layout()
plt.show()

执行结果如下：

生态谬误

生态谬误是辛普森悖论的延伸，指的是用聚合层面的相关系数推断个体层面的相关关系，会产生严重的偏差甚至完全错误的结论。

在社会学、经济学、公共卫生等领域，我们经常会遇到聚合层面的数据，比如城市、省份、国家层面的统计指标。如果我们直接计算这些聚合指标的相关系数，并将其结论推广到个体层面，就会陷入生态谬误。

最经典的例子是，在国家层面，人均巧克力消费量与诺贝尔奖获得数量呈极强的正相关，但这并不意味着个体吃巧克力越多，获得诺贝尔奖的概率越高。国家层面的相关关系，完全无法推广到个体层面。

我们通过代码模拟生态谬误的场景，对比个体层面与聚合层面的相关系数差异，直观展示这种偏差。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)

# 模拟10个城市，每个城市200个个体
n_cities = 10
n_individuals_per_city = 200

# 每个城市的收入水平均值（城市层面的潜在变量）
city_income_mean = np.linspace(3000, 15000, n_cities)

# 存储个体层面数据
individual_income = []
individual_spending = []
city_labels = []

# 生成个体数据
for city_idx, income_mean in enumerate(city_income_mean):
# 个体收入围绕城市均值波动
    income = np.random.normal(income_mean, 2000, n_individuals_per_city)
# 个体消费与个体收入的真实关系：弱正相关，相关系数约0.2
    spending = 0.2 * income + np.random.normal(income_mean * 0.3, 1500, n_individuals_per_city)

    individual_income.extend(income)
    individual_spending.extend(spending)
    city_labels.extend([city_idx] * n_individuals_per_city)

# 转换为数组
individual_income = np.array(individual_income)
individual_spending = np.array(individual_spending)
city_labels = np.array(city_labels)

# 计算个体层面的相关系数
r_individual, _ = stats.pearsonr(individual_income, individual_spending)

# 计算城市层面的聚合数据（均值）
city_income_agg = []
city_spending_agg = []
for city_idx in range(n_cities):
    mask = city_labels == city_idx
    city_income_agg.append(np.mean(individual_income[mask]))
    city_spending_agg.append(np.mean(individual_spending[mask]))

city_income_agg = np.array(city_income_agg)
city_spending_agg = np.array(city_spending_agg)

# 计算城市层面的相关系数
r_city, _ = stats.pearsonr(city_income_agg, city_spending_agg)

# 可视化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 7))

# 个体层面散点图
ax1.scatter(individual_income, individual_spending, alpha=0.3, s=10, color='#377eb8')
slope1, intercept1, _, _, _ = stats.linregress(individual_income, individual_spending)
x_fit1 = np.linspace(min(individual_income), max(individual_income), 100)
ax1.plot(x_fit1, slope1 * x_fit1 + intercept1, color='black', linestyle='--', linewidth=3)
ax1.set_title(f'个体层面：收入 vs 消费\nPearson r = {r_individual:.4f}', fontsize=14)
ax1.set_xlabel('个体月收入', fontsize=12)
ax1.set_ylabel('个体月消费', fontsize=12)
ax1.grid(True, linestyle=':', alpha=0.6)

# 城市层面散点图
ax2.scatter(city_income_agg, city_spending_agg, color='#e41a1c', s=150, marker='o', edgecolors='black')
slope2, intercept2, _, _, _ = stats.linregress(city_income_agg, city_spending_agg)
x_fit2 = np.linspace(min(city_income_agg), max(city_income_agg), 100)
ax2.plot(x_fit2, slope2 * x_fit2 + intercept2, color='black', linestyle='--', linewidth=3)
ax2.set_title(f'城市聚合层面：平均收入 vs 平均消费\nPearson r = {r_city:.4f}', fontsize=14)
ax2.set_xlabel('城市月均收入', fontsize=12)
ax2.set_ylabel('城市月均消费', fontsize=12)
ax2.grid(True, linestyle=':', alpha=0.6)

plt.tight_layout()
plt.show()

# 输出结果
print(f"个体层面相关系数: {r_individual:.4f}")
print(f"城市聚合层面相关系数: {r_city:.4f}")

执行结果如下：

从模拟结果可以看到，个体层面收入与消费的相关系数仅为  左右，属于弱相关；但聚合到城市层面后，相关系数接近，呈现出完美的正相关。如果我们直接用城市层面的强相关，推断个体收入越高消费越高，就会严重高估个体层面的相关强度，陷入生态谬误。

7线性性本质与适用边界

通过前面的推导与模拟，我们已经多次强调，皮尔逊相关系数只能捕捉变量之间的线性依赖关系。本节我们将从数学上严格证明这一点，通过经典的非线性案例展示相关系数的失效场景，最终明确皮尔逊相关系数的核心适用边界。

线性性的数学本质

皮尔逊相关系数为的充要条件是，两个变量之间不存在线性依赖关系，而非不存在任何依赖关系。

我们可以从线性回归的角度来理解这一点。当皮尔逊相关系数  时， 对  的回归系数，这意味着最优的线性拟合线是一条水平直线， 的线性变化完全无法解释  的变化。但这并不代表  和  之间没有关系，它们之间可能存在极强的非线性函数关系，只是这种关系无法用一条直线来拟合。

经典非线性失效案例

我们通过四个经典的非线性模型，从数学上推导相关系数为0的过程，同时通过代码可视化这些强非线性关系下，皮尔逊相关系数失效的场景。

案例1：二次函数关系 

假设  服从对称分布，比如区间  上的均匀分布，或者均值为的正态分布，此时 ，且  的奇数阶矩为0。

我们计算  和  的协方差：

协方差为0，因此皮尔逊相关系数 。

尽管  完全由  确定性决定，两者之间存在完美的二次函数关系，但皮尔逊相关系数却为0，完全无法捕捉这种非线性依赖。

案例2：绝对值关系 

同样假设X服从关于0对称的分布，此时  是一个奇函数，在对称区间上的期望为0：

因此协方差 ，相关系数 。

案例3：正弦函数关系 

假设  服从区间  上的均匀分布，此时  是一个偶函数，而我们计算协方差时：

这里需要注意，当  的区间为  时，相关系数并不为0；但当我们将区间扩展为 时， 在对称区间上的积分会相互抵消，，此时相关系数 。

案例4：圆环分布 

假设  和  服从二维圆环分布，即 ，其中  服从区间 上的均匀分布， 为固定半径。

计算协方差：

因此，相关系数 。

尽管  和  之间存在完美的确定性函数关系，给定  的值可以完全确定  的取值范围，但皮尔逊相关系数依然为0。

我们通过代码可视化这四个案例，验证相关系数的计算结果，直观展示皮尔逊相关系数对非线性关系的盲区。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 设定随机种子
np.random.seed(42)
n = 5000

# 创建可视化画布
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
axes = axes.flatten()

# 案例1：Y = X²
x1 = np.random.normal(0, 2, n)
y1 = x1 ** 2
r1, _ = stats.pearsonr(x1, y1)
axes[0].scatter(x1, y1, alpha=0.4, s=20, edgecolors='white')
axes[0].set_title(f'案例1：Y = X²\nPearson r = {r1:.4f}', fontsize=14)
axes[0].set_xlabel('X', fontsize=12)
axes[0].set_ylabel('Y', fontsize=12)
axes[0].grid(True, linestyle=':', alpha=0.6)

# 案例2：Y = |X|
x2 = np.random.normal(0, 2, n)
y2 = np.abs(x2)
r2, _ = stats.pearsonr(x2, y2)
axes[1].scatter(x2, y2, alpha=0.4, s=20, edgecolors='white')
axes[1].set_title(f'案例2：Y = |X|\nPearson r = {r2:.4f}', fontsize=14)
axes[1].set_xlabel('X', fontsize=12)
axes[1].set_ylabel('Y', fontsize=12)
axes[1].grid(True, linestyle=':', alpha=0.6)

# 案例3：Y = sin(X)，X~U(-2π, 2π)
x3 = np.random.uniform(-2 * np.pi, 2 * np.pi, n)
y3 = np.sin(x3)
r3, _ = stats.pearsonr(x3, y3)
axes[2].scatter(x3, y3, alpha=0.4, s=20, edgecolors='white')
axes[2].set_title(f'案例3：Y = sin(X)\nPearson r = {r3:.4f}', fontsize=14)
axes[2].set_xlabel('X', fontsize=12)
axes[2].set_ylabel('Y', fontsize=12)
axes[2].grid(True, linestyle=':', alpha=0.6)

# 案例4：圆环分布
theta = np.random.uniform(0, 2 * np.pi, n)
r = np.random.normal(5, 0.2, n)
x4 = r * np.cos(theta)
y4 = r * np.sin(theta)
r4, _ = stats.pearsonr(x4, y4)
axes[3].scatter(x4, y4, alpha=0.4, s=20, edgecolors='white')
axes[3].set_title(f'案例4：圆环分布 X²+Y²=R²\nPearson r = {r4:.4f}', fontsize=14)
axes[3].set_xlabel('X', fontsize=12)
axes[3].set_ylabel('Y', fontsize=12)
axes[3].set_aspect('equal')
axes[3].grid(True, linestyle=':', alpha=0.6)

plt.tight_layout()
plt.show()

执行结果如下：

皮尔逊相关系数的核心适用边界

基于前面的所有性质推导与模拟验证，我们可以明确皮尔逊相关系数的核心适用前提与边界：

一是变量之间存在线性关系。这是最核心的前提，皮尔逊相关系数仅能度量线性依赖关系，对任何形式的非线性依赖都会失效，甚至给出完全错误的结论。

二是变量为连续型随机变量。皮尔逊相关系数的定义基于方差与协方差，仅适用于连续型变量，不适用于分类变量、有序等级变量。

三是变量服从二元正态分布。这是相关系数统计推断的核心前提。只有当两个变量服从二元正态分布时，不相关才等价于独立，样本相关系数的抽样分布才能通过Fisher Z变换正态化，相关的假设检验与置信区间计算才有效。

四十数据中无强异常值与高杠杆点。皮尔逊相关系数对异常值与高杠杆点极其敏感，单个极端点可以完全改变相关系数的大小与符号，导致分析结果失真。

五是*数据具有同质性，无潜在分组结构。如果数据中存在潜在的分组变量，可能会出现辛普森悖论，导致整体相关系数与分组内的真实相关关系完全相反，产生误导性结论。

本节总结

在本节中，我们从代数、几何、统计推断三个维度，系统拆解了皮尔逊相关系数的核心性质与边界约束，完成了对这个最常用统计量的全面解构。

在代数性质层面，我们证明了皮尔逊相关系数的对称性、绝对边界性、标准化变量的协方差等价性，以及独立变量的相关系数为0的性质，明确了它的基础代数约束。

在线性变换性质层面，我们推导了平移变换的绝对不变性与尺度变换的符号规则，证明了任意线性变换都不会改变相关系数的绝对值，仅可能反转其符号，解释了相关系数的无量纲特性。

在几何性质层面，我们将相关系数映射为  维样本空间中去中心化向量的夹角余弦值，揭示了它与线性回归系数的内在关联，证明了相关系数的平方等于线性回归的决定系数 ，为相关系数提供了直观的方差解释率解读。

在抽样性质层面，我们明确了样本相关系数是总体相关系数的一致但有偏估计量，推导了它的渐近正态性，提供了小样本下的偏差校正方法，解决了小样本场景下的估计偏差问题。

在扰动性质层面，我们量化了异常值与高杠杆点对相关系数的影响力，通过模拟实验展示了单个极端点对相关系数的强扰动性，明确了皮尔逊相关系数的非稳健性。

在分组数据性质层面，我们通过全相关定律拆解了总体相关系数的组内与组间效应，明确了辛普森悖论的数学条件，揭示了聚合数据带来的生态谬误，强调了分层分析的重要性。

最后，我们通过数学推导与经典案例，明确了皮尔逊相关系数的线性性本质，展示了它在非线性关系下的完全失效，最终划定了皮尔逊相关系数的核心适用边界。