рассчитать плотность многомерного нормального распределения вручную

Question

Я хочу рассчитать плотность многомерного нормального распределения вручную. В качестве входных данных моей функции у меня есть x, которая представляет собой n*p матрицу точек данных, вектор mu со средним значением n и ковариационную матрицу sigma dim p*p.

Я написал для этого следующую функцию:

`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.

В конце я хочу получить те же результаты, что и при использовании функции dmvnorm() из пакета mvtnorm. Что не так с моим кодом?

Answer 1

Выражение t(x-mu)%*%solve(sigma)%*%(x- mu) - это pxp, поэтому ваш результат именно этого размера. Вы хотите диагональ этой матрицы, которую вы можете получить, используя

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
}