如何计算NumPy二维数组中行X与列Y的协方差 可否直接调用np.cov()

import numpy as np
P = np.array(
    [
        [0.03607908, 0.03760034, 0.00503184, 0.0205082 , 0.01051408,
         0.03776221, 0.00131325, 0.03760817, 0.01770659],
        [0.03750162, 0.04317351, 0.03869997, 0.03069872, 0.02176718,
         0.04778769, 0.01021053, 0.00324185, 0.02475319],
        [0.03770951, 0.01053285, 0.01227089, 0.0339596 , 0.02296711,
         0.02187814, 0.01925662, 0.0196836 , 0.01996279],
        [0.02845139, 0.01209429, 0.02450163, 0.00874645, 0.03612603,
         0.02352593, 0.00300314, 0.00103487, 0.04071951],
        [0.00940187, 0.04633153, 0.01094094, 0.00172007, 0.00092633,
         0.02032679, 0.02536328, 0.03552956, 0.01107725]
    ]
)

不能直接调用np.cov()完成这个计算，直接调用得到的结果完全不符合需求。

• np.cov()的设计目标是计算样本协方差，默认参数rowvar=True，逻辑是把输入数组的每一行视为一个独立随机变量，每一列视为对应变量的一次独立观测，且默认所有观测权重相等，对应均匀采样得到的样本数据集，输出结果是所有输入变量的协方差矩阵。哪怕把P误当做样本集传入，得到的也是5个行变量的5×5协方差矩阵，根本不是单个X和Y的协方差值。
• 上述数组P不是样本观测集，是离散随机变量X、Y的联合概率质量矩阵：5行对应X的5个可能取值，9列对应Y的9个可能取值，矩阵元素P[i,j]是X取第i个值、Y取第j个值的联合概率，所有元素求和为1，和np.cov()要求的输入结构完全不匹配，直接调用算出的数值没有实际意义。

需要按照离散随机变量协方差的定义计算，步骤为：先计算X、Y的边缘期望，再代入协方差公式Cov(X,Y) = E[(X-E[X])(Y-E[Y])] = sum_i sum_j (x_i - E_X)*(y_j - E_Y)*P[i,j]即可，参考代码如下：

# 若X实际取值不是行索引0-4、Y实际取值不是列索引0-8，替换为实际取值数组即可
x = np.arange(P.shape[0])
y = np.arange(P.shape[1])
x_grid, y_grid = np.meshgrid(x, y, indexing='ij')

# 计算边缘期望
E_X = np.sum(x_grid * P)
E_Y = np.sum(y_grid * P)

# 计算协方差
cov_XY = np.sum((x_grid - E_X) * (y_grid - E_Y) * P)

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