1 Introdução

A distribuição Normal multivariada tem uma ampla gama de aplicações em estatística. Exemplos incluem análise de componentes principais, regressão multivariada, análise de variância multivariada, análise fatorial entre outras. Dizemos que um vetor aleatório \(\mathbf{Y}\) de dimensão \(n \times 1\) tem distribuição normal multivariada se sua função densidade probabilidade é dada por

\[\begin{equation} f(\mathbf{y}|\mathbf{\mu}, \Sigma)=\left(\frac{1}{2\pi}\right)^{n/2}|\Sigma|^{-1/2}\exp\{-\frac{1}{2}(\mathbf{y}-\mathbf{\mu})^{\top}\Sigma^{-1}(\mathbf{y}-\mathbf{\mu})\}, \end{equation}\]

onde \(\mathbf{\mu}\) é um vetor \(n \times 1\) de valores esperados e \(\Sigma\) é uma matriz \((n \times n)\) simétrica e positiva definida. É usual representar esta situação pela notação \(\mathbf{Y} \sim N(\mathbf{\mu}, \Sigma)\). E neste caso temos que \(E(\mathbf{Y}) = \mathbf{\mu}\) e \(V(\mathbf{Y}) = \Sigma\).

O objetivo deste tutorial é discutir diferentes formas de implementar a distribuição normal multivariada em R.

Note que a implementação da distribuição normal multivariada envolve ao menos três tarefas computacionalmente caras são elas:

Cálculo do determinante de \(\Sigma\);
Multiplicação de matrizes;
Inversa da matriz \(\Sigma\).

O pacote mvtnorm do software R fornece uma implementação de tal distribuição através da função dmvnorm() e usaremos tal função como comparação para mensurar a eficiência das nossas implementações.

O uso de tal função é muito simples e segue o padrão do R de uso de distribuições de probabilidade. O código abaixo ilustra o uso da função dmvnorm() para o cálculo da densidade normal multivariada.

require(mvtnorm)
# Loading required package: mvtnorm
mu <- c(0,0)
Sigma <- matrix(c(1,0.8,0.8,1),2,2)
dmvnorm(x = c(0,0), mean = mu, sigma = Sigma)
# [1] 0.2652582

2 Implementação ingênua

Como primeira tentativa de implementação da distribuição normal vamos usar a função solve() do software R para explicitamente calcular \(\Sigma^{-1}\). Em geral estamos interessados no log da densidade, por isso vamos calcular a densidade sempre em escala log.

my_dmvnorm1 <- function(x, mean, sigma) {
  n <- length(x)
  c1 <- (1/(2*pi))^(n/2)
  c2 <- determinant(sigma, logarithm = FALSE)$modulus^-0.5
  res <- x - mu
  c3 <- t(res)%*%solve(sigma)%*%res
  out <- c1*c2*exp(-0.5*c3)
  return(log(out))
}

Podemos comparar nossa primeira implementação com a implementação padrão do R para verificar se está correta.

my_dmvnorm1(x = c(0,0), mean = mu, sigma = Sigma)
#           [,1]
# [1,] -1.327051
# attr(,"logarithm")
# [1] FALSE
dmvnorm(c(0,0), mean = mu, sigma = Sigma, log = TRUE)
# [1] -1.327051

3 Implementação usando solve()

Uma primeira forma de melhorar a implementação anterior é calcular todas as expressões em escala log para dar maior estabilidade numérica a nossa implementação. Para isto note que

\[\begin{equation} \log(f(\mathbf{y}|\mathbf{\mu},\Sigma)) = -\frac{n}{2} \log(2\pi) - \frac{1}{2} \log |\Sigma| - \frac{1}{2} (\mathbf{y} - \mathbf{\mu})^{\top} \Sigma^{-1}(\mathbf{y} - \mathbf{\mu}) \end{equation}\]

A implementação fica então

my_dmvnorm2 <- function(x, mean, sigma) {
  n <- length(x)
  c1 <- -(n/2)*log(2*pi) 
  c2 <- -0.5*determinant(sigma)$modulus
  res <- x - mean
  c3 <- -0.5*t(res)%*%solve(sigma)%*%res
  return(c1+c2+c3)
}
# Verificando
my_dmvnorm2(x = c(0,0), mean = mu, sigma = Sigma)
#           [,1]
# [1,] -1.327051
# attr(,"logarithm")
# [1] TRUE

4 Implementação via decomposição LU

Note que nas implementações anteriores estamos explicitamente calculando a inversa da matriz \(\Sigma\). Nos vimos que em geral isso não é necessário e podemos economizar os cálculos resolvendo apenas o sistema linear \(\Sigma^{-1} (\mathbf{x} - \mathbf{\mu})\) ao invês de encontrar \(\Sigma^{-1}\) usando a decomposição LU e depois multiplicar pelo vetor \((\mathbf{x} - \mathbf{\mu})\).

my_dmvnorm3 <- function(x, mean, sigma) {
  n <- length(x)
  c1 <- -(n/2)*log(2*pi) 
  c2 <- -0.5*determinant(sigma)$modulus 
  res <- x - mean
  c3 <- -0.5* t(res)%*%solve(sigma,res)
  return(c1+c2+c3)
}
## Verificando
my_dmvnorm3(x = c(0,0), mean = mu, sigma = Sigma)
#           [,1]
# [1,] -1.327051
# attr(,"logarithm")
# [1] TRUE

5 Implementação via decomposição espectral

Na implementação 3 usamos a decomposição LU para resolver o sistema, porém a decomposição LU não nos ajudou a calcular o determinante de \(\Sigma\). Assim, podemos procurar decomposições que nos ajudem a resolver o sistema linear, bem como, a calcular o determinante de \(\Sigma\). Uma decomposição popular é a decomposição em autovalores e autovetores, também chamada de decomposição espectral.

Sendo \(\Sigma\) uma matrix quadrada simétrica a decomposição espectral diz que \(\Sigma = V \Lambda V^{\top}\), onde \(\Lambda\) é uma matriz diagonal com os autovetores. Um resultado bem conhecido é que o produto dos autovetores é igual ao determinante de \(\Sigma\). Além disso, a inversa é facilmente obtida como \(\Sigma = V \Lambda^{-1} V^{\top}\).

Sendo assim,

my_dmvnorm4 <- function(x, mean, sigma) {
  n <- length(x)
  c1 <- -(n/2)*log(2*pi) 
  dd <- eigen(sigma)
  c2 <- -0.5*log(prod(dd$values)) 
  res <- x - mean
  Sigma_inv <- dd$vectors%*%diag(1/dd$values)%*%t(dd$vectors)
  c3 <- -05*t(res)%*%Sigma_inv%*%res
  return(c1+c2+c3)
}
## Verificando
my_dmvnorm4(x = c(0,0), mean = mu, sigma = Sigma)
#           [,1]
# [1,] -1.327051

6 Implementação via decomposição de Cholesky

Uma outra opção é usar a decomposição de Cholesky que é dada por \(\Sigma = LL^{\top}\), onde \(L\) é uma matriz triangular inferior. Uma vez que \(L\) é triangular inferior o determinante de \(\Sigma\) fica sendo o quadrado do produto dos elementos da diagonal de \(L\).

my_dmvnorm5 <- function(x, mean, sigma) {
  n <- length(x)
  c1 <- -(n/2)*log(2*pi) 
  L <- chol(sigma)
  c2 <- -0.5*sum(log(diag(L)^2))
  res <- x - mean
  tt <- solve(L, res)
  c3 <- -0.5*t(tt)%*%tt
  return(c1+c2+c3)
}
# Verificando
my_dmvnorm5(x = c(0,0), mean = mu, sigma = Sigma)
#           [,1]
# [1,] -1.327051

7 Comparando as implementações

Podemos comparar as implementações anteriores para ver qual é mais eficiente computacionalmente. Para isto vamos criar uma matriz de covariancia \(900 \times 900\). E usar o pacote rbenchmark para comparar o tempo computacional de cada implementação.

x <- 1:30
grid <- expand.grid(x,x)
DD <- as.matrix(dist(grid, diag = T, upper = TRUE))
Sigma <- exp(-DD)
mean <- rep(0, dim(Sigma)[1])
x <- rep(0, dim(Sigma)[1])

Avaliando as funções para ter certeza que todas retornam o mesmo valor.

dmvnorm(x = x, mean = mean, sigma = Sigma, log = TRUE)
# [1] -711.7465
my_dmvnorm1(x = x, mean = mean, sigma = Sigma) # Não funciona
#      [,1]
# [1,] -Inf
# attr(,"logarithm")
# [1] FALSE
my_dmvnorm2(x = x, mean = mean, sigma = Sigma)
#           [,1]
# [1,] -711.7465
# attr(,"logarithm")
# [1] TRUE
my_dmvnorm3(x = x, mean = mean, sigma = Sigma)
#           [,1]
# [1,] -711.7465
# attr(,"logarithm")
# [1] TRUE
my_dmvnorm4(x = x, mean = mean, sigma = Sigma)
#           [,1]
# [1,] -711.7465
my_dmvnorm5(x = x, mean = mean, sigma = Sigma)
#           [,1]
# [1,] -711.7465

Claramente a implementação mais ingênua não funciona com matrizes não triviais. Assim, não faz sentido inclui-lá na comparação.

require(rbenchmark)
# Loading required package: rbenchmark
benchmark("dmvnorm" = {dmvnorm(x = x, mean = mean, sigma = Sigma, log = TRUE)},
          "my2" = {my_dmvnorm2(x = x, mean = mean, sigma = Sigma)},
          "my3" = {my_dmvnorm3(x = x, mean = mean, sigma = Sigma)},
          "my4" = {my_dmvnorm4(x = x, mean = mean, sigma = Sigma)},
          "my5" = {my_dmvnorm5(x = x, mean = mean, sigma = Sigma)},
          replications = 10,
          columns = c("test", "replications", "elapsed",
                      "relative", "user.self", "sys.self"))
#      test replications elapsed relative user.self sys.self
# 1 dmvnorm           10   0.670    2.154     1.061    1.235
# 2     my2           10   0.771    2.479     1.874    0.907
# 3     my3           10   0.408    1.312     0.956    0.538
# 4     my4           10   5.483   17.630     9.786   10.372
# 5     my5           10   0.311    1.000     0.715    0.429

Claramente o uso da decomposição de Cholesky é o que permite fazer a melhor implementação da distribuição normal multivariada entre as opções que consideramos. Note que apenas com algum conhecimento de álgebra linear e e propriedades da decomposição de Cholesky nos pudemos facilmente fazer uma função em R mais eficiente que a implementada no pacote mvtnorm(). Deixo como exercício para você identificar o porque que nossa implementação 5 é mais rápido do que a função dmvnorm(), uma vez que ambas estão usando a decomposição de Cholesky para calcular a densidade.

Métodos Computacionais para Inferência Estatística

Implementação: Distribuição Normal Multivariada

Prof. Wagner Hugo Bonat LEG/UFPR

2019