在R中模拟高维多元正态数据遇正定矩阵错误的求助

原代码与报错 #

用户尝试模拟n=100、p=400的高维多元正态数据，变量分两组且存在部分相关性，原代码如下：

## load library MASS
library(MASS)

## sample size set to n = 100
sample_size <- 100      

## 模拟两组各200个变量的均值向量                  
sample_meanvector <- c(runif(200,0,1), runif(200,6,8))

## 构造存在部分相关性的协方差矩阵
sample_covariance_matrix <- matrix(NA, nrow = 400, ncol = 400)
diag(sample_covariance_matrix) <- 1
set.seed(666)
sample_covariance_matrix[lower.tri(sample_covariance_matrix)] <- runif(79800, 0.00001, 0.2)
sample_covariance_matrix[lower.tri(sample_covariance_matrix)][sample(1:79800, 10000)] <- runif(10000, 0.6, 0.9)

## 转换为对称矩阵
sample_covariance_matrix[upper.tri(sample_covariance_matrix)]<-t(sample_covariance_matrix)[upper.tri(sample_covariance_matrix)]

## 生成多元正态分布数据
sample_distribution <- mvrnorm(n = sample_size,
                               mu = sample_meanvector,
                               Sigma = sample_covariance_matrix)

运行时出现报错：

Error in mvrnorm(n = sample_size, mu = sample_meanvector, Sigma = sample_covariance_matrix) : 'Sigma' is not positive definite

问题1：报错原因 #

多元正态分布要求协方差矩阵必须是正定矩阵，你的构造方式存在两个核心问题：

• 随机填充下三角元素后强行对称，极易导致矩阵出现负特征值，破坏正定性。尤其是大量设置0.6-0.9的高相关元素，容易让某些变量的线性组合方差为负或接近0。
• 高维场景下（p=400远大于n=100），随机构造的协方差矩阵天生容易违反正定的数学要求，变量间的相关性组合更容易出现矛盾。

问题2：修正方案 #

方法1：低秩矩阵乘法生成正定矩阵 #

通过低秩矩阵相乘的方式，既能控制分组相关性，又能从数学上保证矩阵正定：

library(MASS)

sample_size <- 100
p <- 400
group1_size <- 200
group2_size <- 200

# 构造分组均值向量
sample_meanvector <- c(runif(group1_size, 0, 1), runif(group2_size, 6, 8))

set.seed(666)
# 构造低秩矩阵，控制分组相关性强度（秩为50，可根据需求调整）
A <- matrix(rnorm(p * 50), nrow = p)
A[1:group1_size, ] <- A[1:group1_size, ] * 0.2  # 第一组弱相关
A[(group1_size+1):p, ] <- A[(group1_size+1):p, ] * 0.8  # 第二组强相关

# 生成正定协方差矩阵，再标准化对角线为1
sample_covariance_matrix <- A %*% t(A)
diag(sample_covariance_matrix) <- 1

# 生成高维多元正态数据
sample_distribution <- mvrnorm(n = sample_size,
                               mu = sample_meanvector,
                               Sigma = sample_covariance_matrix)

方法2：用专用函数生成正定矩阵 #

借助clusterGeneration包的genPositiveDefMat函数直接生成正定矩阵，再调整分组相关性：

library(MASS)
library(clusterGeneration)

sample_size <- 100
p <- 400
group1_size <- 200
group2_size <- 200

sample_meanvector <- c(runif(group1_size, 0, 1), runif(group2_size, 6, 8))

set.seed(666)
# 生成基础正定矩阵
base_cov <- genPositiveDefMat(p, covMethod = "unifcorrmat")$Sigma

# 调整第一组内部弱相关
base_cov[1:group1_size, 1:group1_size] <- apply(base_cov[1:group1_size, 1:group1_size], 
                                                c(1,2), 
                                                function(x) runif(1, 0.00001, 0.2))
# 调整第二组内部强相关
base_cov[(group1_size+1):p, (group1_size+1):p] <- apply(base_cov[(group1_size+1):p, (group1_size+1):p], 
                                                      c(1,2), 
                                                      function(x) runif(1, 0.6, 0.9))

# 保证矩阵对称并修正数值误差，强制对角线为1
sample_covariance_matrix <- (base_cov + t(base_cov))/2
diag(sample_covariance_matrix) <- 1
# 最终确保正定（处理轻微数值问题）
sample_covariance_matrix <- nearPD(sample_covariance_matrix)$mat

# 生成数据
sample_distribution <- mvrnorm(n = sample_size,
                               mu = sample_meanvector,
                               Sigma = sample_covariance_matrix)

进阶：块对角协方差矩阵（组间独立） #

如果需要严格的分组独立（组内相关、组间无相关），可以构造块对角正定矩阵：

library(MASS)
library(clusterGeneration)
library(Matrix)

sample_size <- 100
group1_size <- 200
group2_size <- 200

sample_meanvector <- c(runif(group1_size, 0, 1), runif(group2_size, 6, 8))

set.seed(666)
# 生成第一组正定协方差矩阵（弱相关）
cov_group1 <- genPositiveDefMat(group1_size)$Sigma
diag(cov_group1) <- 1
cov_group1[lower.tri(cov_group1)] <- runif(group1_size*(group1_size-1)/2, 0.00001, 0.2)
cov_group1 <- (cov_group1 + t(cov_group1))/2

# 生成第二组正定协方差矩阵（强相关）
cov_group2 <- genPositiveDefMat(group2_size)$Sigma
diag(cov_group2) <- 1
cov_group2[lower.tri(cov_group2)] <- runif(group2_size*(group2_size-1)/2, 0.6, 0.9)
cov_group2 <- (cov_group2 + t(cov_group2))/2

# 组合成块对角矩阵，组间相关性为0
sample_covariance_matrix <- as.matrix(bdiag(cov_group1, cov_group2))

# 生成数据
sample_distribution <- mvrnorm(n = sample_size,
                               mu = sample_meanvector,
                               Sigma = sample_covariance_matrix)

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