使用scipy.stats.multivariate_normal生成9维样本后CDF恒为零的问题

为什么9维样本的CDF全为0？ #

这是高维空间下的数值精度问题，再加上你的协方差里存在退化维度（第4、8维的标准差为0），双重因素导致计算出的累积概率太小，低于Python浮点数的最小可表示非零值（约1e-308），最终被显示为0。

具体来说：

• 你的9维分布是独立单变量正态的乘积（协方差为对角矩阵），但高维下，即使每个单变量的CDF是中等大小，多个小数相乘后结果会指数级缩小。比如9个0.5相乘是0.00195，已经很小；如果有几个维度的样本落在分布尾部，乘积会直接降到浮点数精度以下。
• 另外，第4、8维的标准差为0，这两个维度是固定等于均值的退化分布。虽然allow_singular=True允许这种情况，但高维数值积分算法在处理奇异协方差时，误差会被放大，进一步压低计算结果。

而2维时，即使方差大，两个小数相乘的结果还能保持在浮点数可表示范围内，所以能得到非零值。

解决办法 #

1. 分解独立变量的CDF计算（最有效）

因为你的协方差是对角矩阵，多元正态的CDF等于各个独立单变量正态CDF的乘积。手动逐个计算单变量CDF再相乘，能避免高维积分的数值误差：

import numpy as np
from scipy.stats import norm, multivariate_normal

mean = np.array([1.0,1.0,-1.0,1.0,1.0,-1.0,-1.0,-1.0,-1.0])
std = np.array([35.98047,63.84033,49.01465,0.0,59.04541,33.57812,59.36084,0.0,11.00098])
normal = multivariate_normal(mean=mean, cov=np.diag(std**2), allow_singular=True)

x_samples = [normal.rvs() for _ in range(30)]
cdf = []

for x_i in x_samples:
    # 逐个计算单变量CDF，方差为0时样本等于均值，对应维度CDF为1
    univariate_cdfs = [norm.cdf(x, m, s) if s !=0 else 1.0 for x, m, s in zip(x_i, mean, std)]
    joint_cdf = np.prod(univariate_cdfs)
    cdf.append(joint_cdf)

print(cdf)

这样计算出来的结果就不会全为0，能保留有效数值。

2. 移除退化维度

如果那些方差为0的维度没有实际意义，可以直接从分布中移除，减少维度后再计算CDF：

# 筛选出标准差不为0的维度
valid_idx = std != 0
mean_valid = mean[valid_idx]
std_valid = std[valid_idx]

normal_valid = multivariate_normal(mean=mean_valid, cov=np.diag(std_valid**2))
x_samples_valid = [normal_valid.rvs() for _ in range(30)]
cdf = [normal_valid.cdf(x) for x in x_samples_valid]
print(cdf)

去掉2个退化维度后，剩下7维的CDF计算精度会提升，能得到非零值。

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