6.1 Normal

6.1.1 Normal univariada \(\cdot \; \mathcal{N}(\mu,\sigma)\)

Veja https://filipezabala.com/eb/distr-cont-esp.html#normal.

6.1.2 Normal bivariada \(\cdot \; \mathcal{NB}(\mu_1,\mu_2,\sigma_1,\sigma_2,\rho)\)

A distribuição normal bivariada é anotada por \(\mathcal{NB}(\mu_1,\mu_2,\sigma_1,\sigma_2,\rho)\) e dada pela expressão \[\begin{equation} f(x_1,x_2|\mu_1, \mu_2, \sigma_1, \sigma_2, \rho) = \\ \;\;\; \tfrac{1}{2\pi \sigma_1 \sigma_2 \sqrt{1-\rho^2}} \exp \Bigg\{ -\tfrac{1}{2(1-\rho^2)} \left[ \tfrac{(x_1-\mu_1)^2}{\sigma^2_1} + \tfrac{(x_2-\mu_2)^2}{\sigma^2_2} - \tfrac{2 \rho (x_1 - \mu_1) (x_2 - \mu_2)}{\sigma_1 \sigma_2} \right] \Bigg\} \tag{6.1} \end{equation}\]

para \(x_1, x_2 \in \rm I\!R\), \(\mu \in \rm I\!R^{2}\), \(\sigma_1, \sigma_2 \ge 0\), \(\sigma_{12} \in \rm I\!R\), \(-1 \le \rho \le +1\). O vetor \(\boldsymbol{\mu}' = \begin{bmatrix} \mu_{1} & \mu_{2} \end{bmatrix}\) é dado conforme Eq. (4.7), a matriz de covariâncias \(\Sigma = \begin{bmatrix} \sigma_{1}^2 & \sigma_{12} \\ \sigma_{12} & \sigma_{2}^2 \end{bmatrix}\) conforme Eq. (4.11) e a matriz de correlação \(\boldsymbol{\rho} = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}\) conforme Eq. (4.20).

Exercício 6.1 Verifique analiticamente (i.e., parta das definições e apresente o desenvolvimento das equações) que o produto de duas normais univariadas equivale à definição no caso bivariado quando \(\rho=0\). \(\\\)

Exemplo 6.1 Pode-se aplicar o resultado do Exercício 6.1.

# parâmetros
n <- 100
x1 <- seq(-5, 5, length = n)
x2 <- seq(-5, 5, length = n)
m1 <- 0; m2 <- 0  # médias
s1 <- 1; s2 <- 2  # desvios padrão

# Caso 1: produto de normais independentes, \rho = 0
z1 <- outer(x1, x2, function(x,y) dnorm(x,m1,s1) * dnorm(y,m2,s2))

# gráficos
rgl::persp3d(x1, x2, z1, col = 'gray')

6.1.3 Normal multivariada \(\cdot \; \mathcal{NM}(\boldsymbol{\mu},\Sigma)\)

A distribuição normal multivariada é anotada por \(\mathcal{NM}(\boldsymbol{\mu},\Sigma)\) e definida pela expressão \[\begin{equation} f(\boldsymbol{x}|\boldsymbol{\mu}, \Sigma) = \dfrac{1}{\sqrt{(2\pi)^k} |\Sigma|^{1/2}} \exp \bigg\{ -\frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu})' \Sigma^{-1} (\boldsymbol{x} - \boldsymbol{\mu}) \bigg\} \tag{6.2} \end{equation}\]

para \(\boldsymbol{x} \in \rm I\!R^{n}\), \(\boldsymbol{\mu} \in \rm I\!R^{n}\), \(\sigma_{ij} \ge 0\) se \(i = j\), \(\sigma_{ij} \in \rm I\!R\) se \(i \ne j\). O vetor de médias \(\boldsymbol{\mu}\) é dado conforme Eq. (4.7), e a matriz de covariâncias \(\Sigma\) conforme Eq. (4.11).

(Tong 1990) apresenta ‘algumas propriedades fundamentais’ da distribuição normal multivariada na Seção 1.1.1 (pp. 1-2).

Representa uma extensão natural da distribuição normal univariada e fornece um modelo adequado para muitos problemas da vida real
Pelo Teorema Central do Limite a distribuição do vetor de médias da amostra é assintoticamente normal, permitindo o uso da distribuição normal multivariada para aproximar distribuições de médias para amostras grandes
A função densidade de uma distribuição normal multivariada é determinada exclusivamente pelo vetor de médias e a matriz de covariância da variável aleatória \(X\)
Correlações zero implicam independência
A família de distribuições normais multivariadas é fechada sob transformações lineares e combinações lineares, i.e., as distribuições de transformações lineares ou combinações lineares de variáveis normais multivariadas são novamente normais multivariadas
A distribuição marginal de qualquer subconjunto de componentes de uma variável normal multivariada também é normal multivariada
A distribuição condicional em uma distribuição normal multivariada é multivariada normal
Para a distribuição normal multivariada, as propriedades de dependência positiva e negativa dos componentes de um vetor aleatório são completamente determinadas pelo sinal e pelo tamanho do coeficiente de correlação

6.1.3.1 Biblioteca `mvtnorm`

(Genz and Bretz 2009) revisam os métodos de computação numérica para probabilidades normais e \(t\) multivariadas, com foco nos métodos recentes de integração numérica. (Genz et al. 2021) apresentam a biblioteca mvtnorm para cálculo de probabilidades multivariadas normais e \(t\), quantis, variáveis aleatórias e densidades.

6.1.3.2 Biblioteca `scipy`

O SciPy (Virtanen et al. 2020) é desenvolvido abertamente no GitHub, por meio do consenso do SciPy e da comunidade científica Python. A função scipy.stats.multivariate_normal pode auxiliar no tratamento computacional da normal multivariada.

Exemplo 6.2 Podem-se calcular probabilidades com os dados do Exemplo 6.5. Considere que se deseja obter \[ Pr(X_1<0, X_2<0) = \dfrac{1}{\sqrt{(2\pi)^2 |\Sigma|}} \int_{-\infty}^{0} \int_{-\infty}^{0} e^{-\frac{1}{2} \boldsymbol{x}' \Sigma^{-1} \boldsymbol{x}} dx_2 dx_1 \]

library(mvtnorm)

# parâmetros
m <- c(0,0)
s <- diag(2)

# Pr(X1 < 0, X2 < 0)
lower <- c(-Inf, -Inf)
upper <- c(0, 0)
pmvnorm(lower, upper, m, s)

## [1] 0.25
## attr(,"error")
## [1] 1e-15
## attr(,"msg")
## [1] "Normal Completion"

Exemplo 6.3 (Genz 1992, 144–45) apresenta um exemplo para uma normal de 3 dimensões, também disponível na documentação da função mvtnorm::pmvnorm. Sendo \[\boldsymbol{\mu}' = \begin{bmatrix} 0 & 0 & 0 \end{bmatrix} \;\;\; \text{e} \;\;\; \boldsymbol{\rho} = \begin{bmatrix} 1 & \frac{3}{5} & \frac{1}{3} \\ \frac{3}{5} & 1 & \frac{11}{15} \\ \frac{1}{3} & \frac{11}{15} & 1 \end{bmatrix}\] calcula-se a probabilidade \[ Pr(X_1<1, X_2<4, X_3<2) = \dfrac{1}{\sqrt{(2\pi)^3 |\Sigma|}} \int_{-\infty}^{1} \int_{-\infty}^{4} \int_{-\infty}^{2} e^{-\frac{1}{2} \boldsymbol{x}' \Sigma^{-1} \boldsymbol{x}} dx_3 dx_2 dx_1 \]

Note que a função admite matrizes de correlação triangulares no argumento corr.

m <- 3
rho <- diag(m)
rho[2,1] <- 3/5
rho[3,1] <- 1/3
rho[3,2] <- 11/15
mvtnorm::pmvnorm(lower = rep(-Inf, m), upper = c(1,4,2), 
        mean = rep(0, m), corr = rho)

## [1] 0.8279846
## attr(,"error")
## [1] 4.098998e-07
## attr(,"msg")
## [1] "Normal Completion"

Exercício 6.2 Considere os dados do Exemplo 6.3. Obtenha a matriz \(\boldsymbol{\rho}\) completa e recalcule a probabilidade indicada.

Exemplo 6.4 Pode-se replicar o Exemplo 6.1 via mvtnorm::dmvnorm, fazendo \(\rho = 0\) (note que sigma é a matriz de covariâncias).

library(mvtnorm)
m <- c(m1,m2)
s <- diag(c(s1^2, s2^2))
z2 <- outer(x1, x2, function(x,y) dmvnorm(cbind(x,y), mean = m, sigma = s))

# verificando numericamente
all.equal(z1, z2)

## [1] TRUE

# gráficos
library(rgl)
persp3d(x1, x2, z2, col = 'lightblue')

Exemplo 6.5 Pelo fato de o parâmetro sigma da função dmvnorm ser a matriz de covariâncias (e não de correlações), sabe-se da Eq. (4.21) que a correlação \(\rho_{jk}\) é obtida dividindo-se a covariância \(\sigma_{jk}\) pelo produto dos desvios padrão \(\sigma_{j}\) e \(\sigma_{k}\). Assim, \[\begin{equation} \sigma_{jk} = \rho_{jk}\sigma_{j}\sigma_{k} \tag{6.3} \end{equation}\]

# parâmetros
n <- 100
x1 <- seq(-5, 5, length = n)
x2 <- seq(-5, 5, length = n)
m1 <- 0; m2 <- 0  # médias
s1 <- 1; s2 <- 2  # desvios padrão
rho <- 0.9        # correlação
s12 <- rho*s1*s2  # covariância

# via mvtnorm::dmvnorm, \rho = 0.9
library(mvtnorm)
m <- c(m1,m2)
s <- matrix(c(s1^2, s12, s12, s2^2), nrow = 2, byrow = T)
z3 <- outer(x1, x2, function(x,y) dmvnorm(cbind(x,y), mean = m, sigma = s))

# gráficos
library(rgl)
persp3d(x1, x2, z3, col = 'lightgreen')

Exemplo 6.6 Podem-se simular dados pseudo-aleatórios de uma distribuição normal multivariada a partir da função mvtnorm::rmvnorm.

library(mvtnorm)
set.seed(111) # para garantir reprodutibilidade
sim1 <- rmvnorm(1000, mean = c(100,20,0,-10), sigma = diag(4))
colMeans(sim1) # vetor de médias (amostrais)

## [1] 99.94026509 20.03499171  0.01309363 -9.97267381

cov(sim1)      # matriz de covariâncias (amostrais)

##            [,1]        [,2]        [,3]       [,4]
## [1,] 0.97980660  0.01473150  0.02560396 0.01487046
## [2,] 0.01473150  1.02316657 -0.05805611 0.02423877
## [3,] 0.02560396 -0.05805611  0.96812296 0.03730201
## [4,] 0.01487046  0.02423877  0.03730201 1.00238890

cor(sim1)      # matriz de correlação (amostral)

##            [,1]        [,2]        [,3]       [,4]
## [1,] 1.00000000  0.01471308  0.02628886 0.01500500
## [2,] 0.01471308  1.00000000 -0.05833235 0.02393422
## [3,] 0.02628886 -0.05833235  1.00000000 0.03786595
## [4,] 0.01500500  0.02393422  0.03786595 1.00000000

Exercício 6.3 Considere a distribuição normal multivariada definida na Eq. (6.2) e as funções de normal multivariada da sua ferramenta computacional de preferência. Sejam \(p=4\) dimensões tal que \(\mu_1 = 100\), \(\mu_2 = 0\), \(\mu_3 = -10\), \(\mu_4 = -200\), \(\rho_{12}=-0.86\), \(\rho_{13}=0.38\), \(\rho_{14}=-0.67\), \(\rho_{23}=-0.71\), \(\rho_{24}=0.25\), \(\rho_{34}=0.23\), \(\sigma_1=3\), \(\sigma_2=16\), \(\sigma_3=18\), \(\sigma_4=61\). Dica: considere a Eq. (6.3).

Leia a discussão https://stackoverflow.com/questions/47403003/why-is-the-pmvnorm-result-different-when-the-input-matrix-are-covariance-and-c.
Compare os algoritmos de GenzBretz, Miwa e TVPACK, detalhados na documentação da biblioteca mvtnorm (Genz et al. 2021) do R.
Calcule \(Pr(X_1<100, X_2<2, X_3<0, X_4<-50)\) comparando os três algoritmos descritos no item b.
Calcule \(Pr(X_1>100, X_2<2, X_3>0, X_4<-50)\) comparando os três algoritmos descritos no item b.
Calcule \(Pr(70<X_1<120, X_2<2, X_3>0, -125<X_4<-50)\) comparando os três algoritmos descritos no item b.
Simule \(n=2000\) valores da referida normal multivariada. Calcule o vetor de médias e as matrizes de covariância e correlação a partir dos valores simulados.
Compare com os valores utilizados na geração pseudo-aleatória. Dica: use o teste de Hotteling⁷ para comparar vetores de média. #e teste M de Box⁸ para comparar as matrizes de covariância.
Repita os itens c, d, e, agora utilizando os valores obtidos da simulação no item f. Utilize apenas um método (GenzBretz ou Miwa) utilizando a matriz de covariâncias amostrais. # i. Verifique que os cálculos podem ser realizados utilizando funções distribuição \(t\) de Student multivariadas com infinitos graus de liberdade.

Sugestão de solução: Capítulo 14

Exercício 6.4 Replique os exemplos e exercícios desta seção utilizando a biblioteca SciPy, descrita na Seção 6.1.3.2.

6.1.4 Para saber mais

(Patel and Read 1982) trazem uma coleção de resultados e propriedades relacionados à distribuição normal. (Tong 1990) fornece um tratamento abrangente de resultados relacionados à distribuição normal multivariada. Os temas principais são dependência, desigualdades de probabilidade e seus papéis na teoria e aplicações.

Referências

Genz, Alan. 1992. “Numerical Computation of Multivariate Normal Probabilities.” Journal of Computational and Graphical Statistics 1 (2): 141–49. https://doi.org/10.2307/1390838.

Genz, Alan, and Frank Bretz. 2009. Computation of Multivariate Normal and t Probabilities. Lecture Notes in Statistics. Heidelberg: Springer-Verlag.

Genz, Alan, Frank Bretz, Tetsuhisa Miwa, Xuefei Mi, Friedrich Leisch, Fabian Scheipl, and Torsten Hothorn. 2021. mvtnorm: Multivariate Normal and t Distributions. https://CRAN.R-project.org/package=mvtnorm.

Patel, Jagdish K., and Campbell B. Read. 1982. “Handbook of the Normal Distribution.” Marcel Dekker, New York.

Tong, Yung L. 1990. The Multivariate Normal Distribution. Springer-Verlag New York Inc.

Virtanen, Pauli, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, et al. 2020. “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.” Nature Methods 17: 261–72. https://doi.org/10.1038/s41592-019-0686-2.