5.1 Normal

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

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

5.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{5.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.9), a matriz de covariâncias \(\Sigma = \begin{bmatrix} \sigma_{1}^2 & \sigma_{12} \\ \sigma_{12} & \sigma_{2}^2 \end{bmatrix}\) conforme Eq. (4.13) e a matriz de correlação \(\boldsymbol{\rho} = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}\) conforme Eq. (4.22).

Exercício 5.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 5.1 Pode-se aplicar o resultado do Exercício 5.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
if(type_book == 'bookdown::gitbook'){
  rgl::persp3d(x1, x2, z1, col = 'gray')
}

5.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{5.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.9), e a matriz de covariâncias \(\Sigma\) conforme Eq. (4.13).

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

  1. Representa uma extensão natural da distribuição normal univariada e fornece um modelo adequado para muitos problemas da vida real
  2. 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
  3. 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\)
  4. Correlações zero implicam independência
  5. 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
  6. A distribuição marginal de qualquer subconjunto de componentes de uma variável normal multivariada também é normal multivariada
  7. A distribuição condicional em uma distribuição normal multivariada é multivariada normal
  8. 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

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

5.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 5.2 Podem-se calcular probabilidades com os dados do Exemplo 5.1. Considere que se deseja obter \[ Pr(X_1<0, X_2<0) = \dfrac{1}{\sqrt{(2\pi)^2 \left|\begin{bmatrix} 1 & 0 \\ 0 & 4 \end{bmatrix} \right|}} \int_{-\infty}^{0} \int_{-\infty}^{0} e^{-\frac{1}{2} \boldsymbol{x}' \begin{bmatrix} 1 & 0 \\ 0 & 4 \end{bmatrix}^{-1} \boldsymbol{x}} dx_2 dx_1 \]

library(mvtnorm)

# parâmetros
m <- c(0,0)   # vetor de médias
s <- diag(2)  # matriz de covariâncias
s[2,2] <- 2^2 # variância de X2

# Pr(X1 < 0, X2 < 0)
lower <- c(-Inf, -Inf)
upper <- c(0, 0)
pmvnorm(lower, upper, mean = m, sigma = s, seed = 42)
## [1] 0.25
## attr(,"error")
## [1] 1e-15
## attr(,"msg")
## [1] "Normal Completion"

Exemplo 5.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, seed = 42)
## [1] 0.8279853
## attr(,"error")
## [1] 1.50763e-06
## attr(,"msg")
## [1] "Normal Completion"

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

Exemplo 5.4 Pode-se replicar o Exemplo 5.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)
if(type_book == 'bookdown::gitbook'){
  persp3d(x1, x2, z2, col = 'lightblue')
}

Exemplo 5.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.23) 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{5.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)
if(type_book == 'bookdown::gitbook'){
  persp3d(x1, x2, z3, col = 'lightgreen')
}

Exemplo 5.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))
colnames(sim1) <- paste0('X',1:4)
colMeans(sim1) # vetor de médias (amostrais)
##           X1           X2           X3           X4 
## 100.01080923  19.93781889   0.02050894  -9.95346044
cov(sim1)      # matriz de covariâncias (amostrais)
##              X1         X2           X3           X4
## X1  0.976004830 0.04433941 -0.008425825  0.006718199
## X2  0.044339406 0.96758992  0.004532920  0.018780622
## X3 -0.008425825 0.00453292  1.046619468 -0.046749551
## X4  0.006718199 0.01878062 -0.046749551  0.982410919
cor(sim1)      # matriz de correlação (amostral)
##              X1          X2           X3           X4
## X1  1.000000000 0.045626612 -0.008336659  0.006860888
## X2  0.045626612 1.000000000  0.004504411  0.019262724
## X3 -0.008336659 0.004504411  1.000000000 -0.046103771
## X4  0.006860888 0.019262724 -0.046103771  1.000000000

Exercício 5.3 Considere a distribuição normal multivariada definida na Eq. (5.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. (5.3).

  1. Leia a discussão https://stackoverflow.com/questions/47403003/why-is-the-pmvnorm-result-different-when-the-input-matrix-are-covariance-and-c.
  2. Compare os algoritmos de GenzBretz, Miwa e TVPACK, detalhados na documentação da biblioteca mvtnorm (Genz et al. 2021) do R.
  3. Calcule \(Pr(X_1<100, X_2<2, X_3<0, X_4<-50)\) comparando os três algoritmos descritos no item b.
  4. Calcule \(Pr(X_1>100, X_2<2, X_3>0, X_4<-50)\) comparando os três algoritmos descritos no item b.
  5. 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.
  6. 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.
  7. Compare com os valores utilizados na geração pseudo-aleatória. Dica: use o teste de Hotteling5 para comparar vetores de média.
  8. 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.

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

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

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

References

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.

  1. https://search.r-project.org/CRAN/refmans/ICSNP/html/HotellingsT.html, https://filipezabala.com/ea/teste-vetor-media.html↩︎