r-手动计算多元正态分布的密度

Question

我想手动计算多元正态分布的密度。作为我的函数的输入，我有作为数据点的n*p矩阵的x、具有n均值的向量mu和dimp*p的协方差矩阵sigma。

我为此编写了以下函数：

`dmnorm <- function(mu, sigma, x){
k <- ncol(sigma)
x <- t(x)
dmn <- exp((-1/2)*t(x-mu)%*%solve(sigma)%*%(x- 
mu))/sqrt(((2*pi)^k)*det(sigma))  
return(dmn)
}`

我自己的函数给了我一个n*n的矩阵。然而，我应该得到一个长度为n的向量。

最后，我希望得到与使用mvtnorm包中的dmvnorm()函数相同的结果。我的代码出了什么问题？

Answer 1

表达式t(x-mu)%*%solve(sigma)%*%(x- mu)是p x p，所以这就是为什么结果是这样大的原因。你想要矩阵的对角线，你可以用得到

diag(t(x-mu)%*%solve(sigma)%*%(x-mu))

所以完整的功能应该是

dmnorm <- function(mu, sigma, x){
k <- ncol(sigma)
x <- t(x)
dmn <- exp((-1/2)*diag(t(x-mu)%*%solve(sigma)%*%(x- 
mu)))/sqrt(((2*pi)^k)*det(sigma))  
dmn
}