如何使用SymPy完成皮尔逊相关系数求导与数值代入

皮尔逊相关系数的求和形式定义为：
$$\rho(X,Y) = \frac{n\sum_{i=1}^n x_i y_i - \sum_{i=1}^n x_i \sum_{i=1}^n y_i}{\sqrt{(n\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^{2)(n\sum_{i=1}}n y_i^2 - (\sum_{i=1}^n y_i)^2)}}$$

以下是完整可运行的实现步骤：

1. 导入SymPy库 #

import sympy as sp

2. 定义符号、求和项与相关系数表达式 #

# 定义基础符号：样本量n，求和索引i
n = sp.symbols('n', integer=True, positive=True)
i = sp.symbols('i', integer=True)
# 定义索引变量，存储x、y的样本序列
x = sp.IndexedBase('x')
y = sp.IndexedBase('y')

# 定义所有求和项
sum_x = sp.Sum(x[i], (i, 1, n))
sum_y = sp.Sum(y[i], (i, 1, n))
sum_xy = sp.Sum(x[i] * y[i], (i, 1, n))
sum_x_sq = sp.Sum(x[i] ** 2, (i, 1, n))
sum_y_sq = sp.Sum(y[i] ** 2, (i, 1, n))

# 构造皮尔逊相关系数的符号表达式
rho = (n * sum_xy - sum_x * sum_y) / sp.sqrt((n * sum_x_sq - sum_x ** 2) * (n * sum_y_sq - sum_y ** 2))

3. 对指定变量求导 #

比如要对第k个x样本值x[k]求导，实现代码如下：

# 定义目标变量的索引k
k = sp.symbols('k', integer=True, positive=True)
# 求导并化简表达式
d_rho_d_xk = sp.simplify(sp.diff(rho, x[k]))
# 可打印查看化简后的导数表达式
print("导数表达式：", d_rho_d_xk)

4. 代入实际数值计算 #

以3个样本为例，x = [1,2,3]，y = [2,4,5]，计算相关系数本身和对x_2的导数值：

# 样本数值配置
sample_n = 3
sample_x = [1, 2, 3]
sample_y = [2, 4, 5]

# 计算皮尔逊相关系数的实际值
rho_val = rho.subs(n, sample_n).doit()
for idx in range(sample_n):
    rho_val = rho_val.subs({x[idx+1]: sample_x[idx], y[idx+1]: sample_y[idx]})
print("皮尔逊相关系数值：", float(rho_val))

# 计算对x_2的导数值
d_rho_d_x2_val = d_rho_d_xk.subs({n: sample_n, k: 2}).doit()
for idx in range(sample_n):
    d_rho_d_x2_val = d_rho_d_x2_val.subs({x[idx+1]: sample_x[idx], y[idx+1]: sample_y[idx]})
print("对x_2的导数值：", float(d_rho_d_x2_val))

补充说明 #

• 如果需要批量处理大量样本，可以用sp.lambdify把符号表达式转换为数值函数，计算效率更高
• 若需要对y序列的变量求导，只需要修改diff方法的第二个参数为对应y的索引变量即可

内容的提问来源于stack exchange，提问作者Shyam