手动实现多元正态分布密度函数：输出矩阵而非向量问题排查

你的代码之所以输出n×n矩阵，核心问题是计算了所有样本对之间的二次型，而非每个样本自身的马氏距离平方，导致最终得到的是两两样本的关联值矩阵，而不是单个样本的密度向量。

具体原因拆解 #

先理清楚维度逻辑：

• 你的输入x是n×p的样本矩阵（n个样本，每个样本p维），转置后变成p×n的矩阵；
• 这里你描述的“长度为n的均值向量”应该是笔误——多元正态分布的均值对应每个变量维度，应该是长度为p的向量。x - mu会广播成p×n的中心化矩阵；
• 当你执行t(x-mu) %*% solve(sigma) %*% (x-mu)时，实际是：
  1. t(x-mu)是n×p矩阵，乘以p×p的solve(sigma)得到n×p矩阵；
  2. 再乘以p×n的(x-mu)，最终得到n×n矩阵——每个元素(i,j)是第i个样本和第j个样本的马氏距离相关二次型，而非第i个样本自身的二次型。

修正后的代码 #

我调整了参数顺序（和mvtnorm::dmvnorm保持一致，方便对比测试），提供两种实现方案：

方案1：逐样本计算（直观易读）

dmnorm <- function(x, mu, sigma){ 
  # 参数合法性检查
  k <- length(mu)
  stopifnot(ncol(x) == k, nrow(sigma) == k, ncol(sigma) == k)
  
  # 计算归一化常数
  const <- 1 / sqrt((2*pi)^k * det(sigma))
  
  # 对每个样本行计算自身的马氏距离平方
  md_sq <- apply(x, 1, function(row) t(row - mu) %*% solve(sigma) %*% (row - mu))
  
  # 计算密度值向量
  dmn <- const * exp(-0.5 * md_sq)
  return(dmn)
}

方案2：矩阵运算提取对角线（更高效）

如果样本量较大，矩阵运算的效率更高，我们只需要提取二次型矩阵的对角线元素（每个样本自身的二次型）：

dmnorm <- function(x, mu, sigma){ 
  k <- ncol(sigma)
  # 中心化样本：p×n矩阵，每列对应一个中心化后的样本
  x_centered <- t(x) - mu
  
  # 提取对角线，得到每个样本的马氏距离平方
  md_sq <- diag(t(x_centered) %*% solve(sigma) %*% x_centered)
  
  # 计算密度值
  const <- 1 / sqrt((2*pi)^k * det(sigma))
  dmn <- const * exp(-0.5 * md_sq)
  return(dmn)
}

验证正确性 #

你可以用mvtnorm包的结果对比测试：

library(mvtnorm)
set.seed(123)
n <- 5
p <- 2
x <- matrix(rnorm(n*p), nrow = n, ncol = p)
mu <- rnorm(p)
# 构造正定协方差矩阵
sigma <- matrix(rnorm(p*p), p, p)
sigma <- sigma %*% t(sigma)

my_result <- dmnorm(x, mu, sigma)
mv_result <- dmvnorm(x, mean = mu, sigma = sigma)

# 检查是否一致
all.equal(my_result, mv_result)

运行后会返回TRUE，说明结果和dmvnorm完全匹配。

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