使用faux包rnorm_multi()时自定义相关矩阵非正定的问题咨询

问题描述 #

尝试使用faux包的rnorm_multi()函数生成数据框，希望传入的相关矩阵包含0到1、步长为0.01的所有相关值，但运行后始终报错：

Error in cormat(r, vars) : correlation matrix not positive definite

构造矩阵的代码如下：

library(faux)
# 创建空的100x100矩阵，填充下三角为目标值，转置后重复操作使其对称
matrix <- matrix(rep(0, 100*100), nrow = 100) 
matrix[lower.tri(matrix)] <- seq(0, .99, .01) 
matrix <- t(matrix) 
matrix[lower.tri(matrix)] <- seq(0, .99, .01) 
diag(matrix) <- 1 # 相关矩阵对角线应为1

# 基于上述相关矩阵生成数据框
df <- rnorm_multi(100, 100, 0, 1, r = matrix, 
                  varnames = paste("x", seq(0, 1, .01), sep = "."))

核心原因 #

你构造的矩阵不满足相关矩阵的正定（或半正定）要求。相关矩阵必须是正定矩阵（所有特征值严格大于0）或半正定矩阵（所有特征值非负），这是多元正态分布存在的必要条件。

你的构造方式是按递增顺序填充下三角，这种规则化的填充会导致某些变量组合的相关性违反正定约束——比如多个高相关变量的组合会让矩阵的特征值出现负数，从而触发报错。

可行解决方案 #

直接使用原矩阵不可行，但可以通过以下方式修正：

1. 修正为最近的正定矩阵 #

使用Matrix包的nearPD()函数，找到与原矩阵最接近的正定矩阵，微小调整后即可满足rnorm_multi()的要求：

library(faux)
library(Matrix)

# 构造初始矩阵
matrix <- matrix(rep(0, 100*100), nrow = 100) 
matrix[lower.tri(matrix)] <- seq(0, .99, .01) 
matrix <- t(matrix) 
matrix[lower.tri(matrix)] <- seq(0, .99, .01) 
diag(matrix) <- 1 

# 修正为正定矩阵（保留原矩阵的大部分特征）
pd_matrix <- as.matrix(nearPD(matrix)$mat)

# 生成数据框
df <- rnorm_multi(100, 100, 0, 1, r = pd_matrix, 
                  varnames = paste("x", seq(0, 1, .01), sep = "."))

2. 随机分配相关值 #

如果允许打乱相关值的位置，也可以将0到0.99的步长值随机分配到下三角，再构造对称矩阵，这种方式更大概率生成正定矩阵（但仍需检查）：

library(faux)

# 随机打乱相关值
cor_vals <- sample(seq(0, .99, .01))

# 构造矩阵
matrix <- matrix(rep(0, 100*100), nrow = 100) 
matrix[lower.tri(matrix)] <- cor_vals
matrix <- t(matrix) 
matrix[lower.tri(matrix)] <- cor_vals
diag(matrix) <- 1 

# 检查矩阵是否正定（允许微小数值误差）
all(eigen(matrix)$values >= -1e-8) 

# 若正定则生成数据
if(all(eigen(matrix)$values >= -1e-8)){
  df <- rnorm_multi(100, 100, 0, 1, r = matrix, 
                    varnames = paste("x", seq(0, 1, .01), sep = "."))
}

内容的提问来源于stack exchange，提问作者Sören Michallek