如何在R中手动编写多元正态分布核似然？代码正确性验证

你的代码核心错误 #

你写的这行代码是不正确的：

(colSums(dat-mu_mat)%*%solve(sigma)%*%colSums(dat-mu_mat))

问题出在先对所有样本与均值的差值求和，再计算二次型，这和多元正态对数似然要求的「每个样本差值的二次型之和」完全不是一回事。正确逻辑应该是对每个样本计算$(x_i - \mu)^T \Sigma^{-1} (x_i - \mu)$，再把所有结果加起来，而不是先把差值向量加总再算二次型。

正确的二次型计算方式 #

对于n个d维样本（存储为n行d列的矩阵dat），均值向量mu，协方差矩阵sigma，可以用以下几种高效方式计算二次型总和：

• 方式1：用R内置的mahalanobis函数，直接计算每个样本的马氏距离（也就是我们需要的二次型），再求和：

  centered <- dat - mu  # R自动广播，无需手动构造mu_mat
  sum_quad <- sum(mahalanobis(dat, center = mu, cov = sigma))
  

• 方式2：用矩阵运算计算迹（二次型之和等于t(centered) %*% solve(sigma) %*% centered的对角线元素之和）：

  centered <- dat - mu
  sum_quad <- sum(diag(t(centered) %*% solve(sigma) %*% centered))
  

修正后的完整代码 #

library(mvtnorm)

# 生成模拟数据
dat <- rmvnorm(100, mean = c(200, 0.1), sigma = matrix(c(5, 0, 0, 0.02), nrow = 2))
n <- nrow(dat)
d <- ncol(dat)
mu <- c(200, 0.1)
sigma <- matrix(c(5, 0, 0, 0.02), nrow = 2)

# 手动实现多元正态对数似然
loglik_mvn <- function(n, d, x, mu, sigma) {
  # 计算所有样本的二次型总和
  sum_quad <- sum(mahalanobis(x, center = mu, cov = sigma))
  
  # 对数似然公式计算
  (-n*d/2)*log(2*pi) - (n/2)*determinant(sigma, logarithm = TRUE)$modulus - 0.5*sum_quad
}

# 和内置函数结果对比
manual_result <- loglik_mvn(n, d, dat, mu, sigma)
builtin_result <- sum(dmvnorm(dat, mean = mu, sigma = sigma, log = TRUE))

cat("手动计算对数似然：", manual_result, "\n")
cat("内置函数计算对数似然：", builtin_result, "\n")

额外注意点 #

• 无需手动构造mu_mat：R的广播机制会自动将长度为d的均值向量mu与n行d列的dat逐行相减，得到每个样本的差值矩阵。
• 行列式计算：用determinant(sigma, logarithm = TRUE)$modulus比det(sigma, log=TRUE)兼容性更好，两者都能得到对数行列式值。
• 结果一致性：运行修正后的代码，手动计算结果和dmvnorm求和的结果应该完全一致（浮点误差范围内）。

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