R语言mvtnorm包dmvnorm()输出异常：仅差一行的矩阵结果差异显著

问题现象 #

使用mvtnorm包的dmvnorm()函数时，仅相差一行数据的两个矩阵，计算得到的多元正态密度结果差异极大：包含21行数据的train_mat2输出全为0，而取前20行的train_mat3输出为正常的概率密度值。

复现代码 #

library(mvtnorm)
library(magrittr)

train_mat2 <- 
  data.frame(c1 = c(13.29, 13.68, 14.06, 14.10, 13.86, 13.05, 14.20, 13.77, 
                    12.64, 13.67, 12.17, 12.37, 12.37, 12.29, 13.86, 13.49,
                    11.66, 11.84, 12.33, 12.70, 12.00),
             c2 = c(2.68, 2.36, 2.61, 2.40, 2.27, 2.55, 2.45, 2.68, 2.02, 1.92, 
                    2.53, 2.56, 1.92, 2.21, 2.67, 2.24, 1.92, 2.23, 1.95, 2.40, 2.00)) %>%
  as.matrix()

train_mat3 <- train_mat2[1:20,]

# 第一个调用输出全0
dmvnorm(train_mat2,
        mean=colMeans(train_mat2),
        sigma=cov(train_mat2))

# 第二个调用输出正常密度值
dmvnorm(train_mat3,
        mean=colMeans(train_mat3),
        sigma=cov(train_mat3))

输出结果 #

• train_mat2的计算输出：

[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

• train_mat3的计算输出：

[1] 0.32385982 0.59184100 0.33153870 0.34550971 0.38847115 0.52154257 0.30568408 0.32928314 0.41859960 0.07385217 0.16868584 0.21340448 0.25124098
[14] 0.51129759 0.32758244 0.58043723 0.12999843 0.24604814 0.28729287 0.63333639

原因解释 #

核心原因是数值精度下溢：

1. 加入第21行数据后，train_mat2的协方差矩阵行列式变得极小（可通过det(cov(train_mat2))验证），协方差矩阵接近奇异状态。
2. 多元正态密度计算公式中包含exp(-0.5*(x-μ)^T Σ^{-1} (x-μ))项，当Σ接近奇异时，Σ的逆矩阵数值会变得极大，导致指数部分成为绝对值极大的负数，最终exp(极大负数)的结果小到超出R的浮点数精度范围，被显示为0。
3. train_mat3的协方差矩阵行列式更大，逆矩阵数值更稳定，计算出的密度值处于正常浮点数范围内，因此能显示正常结果。

验证与解决方法 #

• 验证：使用log=TRUE参数计算对数密度，避免下溢问题，可查看实际的对数值：

dmvnorm(train_mat2, mean=colMeans(train_mat2), sigma=cov(train_mat2), log=TRUE)

• 解决：若需要非对数的密度值，可将对数结果取指数；或检查数据是否存在线性相关、异常值，对数据进行预处理（如移除异常值、正则化协方差矩阵）。

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