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

Я хочу рассчитать плотность многомерного нормального распределения вручную. В качестве входных данных моей функции у меня есть 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. Что не так с моим кодом?


r function probability-density normal-distribution
person Phd Student    schedule 09.04.2020    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Выражение 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
}
person user2554330    schedule 09.04.2020