Simulação da Normal Multivariada
Mauro Campos. Departamento de Estatística, UFES
27 de agosto de 2018
- 1 A distribuição normal multivariada
- 2 Visualização da normal multivariada
- 3 Funções do R para simulação da normal multivariada
- 4 A decomposição espectral da matriz de covariância
- 5 Simulação da normal multivariada via decomposição espectral
- 6 Simulação de dados para o modelo de regressão linear simples
1 A distribuição normal multivariada
A distribuição normal multivariada é uma generalização da distribuição normal univariada com média \(-\infty<\mu<+\infty\) e variância \(\sigma^{2}>0\) cuja a densidade é dada por \[ f(x)=\frac{1}{\sqrt{2\pi}\sigma}\cdot\exp\left[-\frac{1}{2}\frac{(x-\mu)^{2}}{\sigma^{2}}\right]\hspace{1cm}-\infty<x<+\infty. \] Em \(d\) dimensões, um vetor aleatório \(\boldsymbol{X}=(X_{1},\dots,X_{d})^{T}\) possui distribuição normal multivariada se sua densidade é dada por \[ f(\boldsymbol{x})=\frac{1}{(2\pi)^{d/2}|\Sigma|^{1/2}}\cdot\exp\left[-\frac{1}{2} (\boldsymbol{x}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right]\hspace{1cm}\boldsymbol{x}\in\mathbb{R}^{d} \] onde \(\boldsymbol{\mu}\) é o vetor de médias e \(\boldsymbol{\Sigma}\) é a matriz de covariâncias (ou matriz de variâncias-covariâncias). Ou seja, \[ \boldsymbol{\mu}=\mbox{E}(\boldsymbol{X})=(\mu_{1},\dots,\mu_{d})^{T}_{d\times 1} \] onde \(\mu_{i}=\mbox{E}(X_{i})\) e \[ \boldsymbol{\Sigma}=\mbox{Cov}(\boldsymbol{X})=(\sigma_{ij})_{d\times d} \] onde \(\sigma_{ij}=\mbox{Cov}(X_{i},X_{j})\). É importante enfatizar que \(\boldsymbol{\Sigma}\) é uma matriz simétrica e positiva-definida, ou seja, \(\boldsymbol{x}^{T}\boldsymbol{\Sigma}\boldsymbol{x}>0\) para todo \(\boldsymbol{x}\) não-nulo.
Para denotar o fato que um vetor aleatório \(\boldsymbol{X}\) possui distribuição normal multivariada com média \(\boldsymbol{\mu}\) e covariância \(\boldsymbol{\Sigma}\), usamos a seguinte notação: \[ \boldsymbol{X}\sim\mbox{N}_{d}(\boldsymbol{\mu},\boldsymbol{\Sigma}). \]
Um vetor aleatório \(\boldsymbol{Z}=(Z_{1},\dots,Z_{d})^{T}\) possui distribuição normal multivariada padrão se \[ \boldsymbol{Z}\sim\mbox{N}_{d}(\boldsymbol{0},\boldsymbol{I}) \] onde \(\boldsymbol{0}=(0,\dots,0)^{T}\) é o vertor nulo \(d\)-dimensional e \(\boldsymbol{I}\) é a matriz identidade \(d\times d\).
2 Visualização da normal multivariada
A visualização direta da densidade da normal multivariada só é possível em \(d=1\) ou \(d=2\). Portanto vamos nos concentrar num caso específico bidimensional.
Considere uma distribuição normal multivariada com a seguinte média e covariância: \[ \mu=(1,1)^{T}\hspace{1cm}\mbox{e}\hspace{1cm} \Sigma= \begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix} . \]
No caso em questão, podemos usar a função dmvnorm() do pacote mvtnorm do R para plotar a densidade:
library(mvtnorm)
## Warning: package 'mvtnorm' was built under R version 3.4.4
nmax <- 50
x.points <- seq(-3,5,length.out=nmax)
y.points <- x.points
z <- matrix(0,nrow=nmax,ncol=nmax)
mu.vec <- c(1,1);sigma.mat <- matrix(c(2,1,1,1),nrow=2,ncol=2,byrow=TRUE)
mu.vec
## [1] 1 1
sigma.mat
## [,1] [,2]
## [1,] 2 1
## [2,] 1 1
for(i in 1:nmax){
for(j in 1:nmax){
z[i,j] <- dmvnorm(c(x.points[i],y.points[j]),mean=mu.vec,sigma=sigma.mat)
}
}
persp(x.points,y.points,z,xlab="x",ylab="y",theta=30,phi=45,col="aliceblue")
Figura: Densidade normal bivariada
Igualmente, podemos usar o R para plotar as curvas de nível da densidade:
nmax <- 100
x.points <- seq(-3,5,length.out=nmax)
y.points <- x.points
z <- matrix(0,nrow=nmax,ncol=nmax)
mu.vec <- c(1,1);sigma.mat <- matrix(c(2,1,1,1),nrow=2,ncol=2,byrow=TRUE)
for(i in 1:nmax){
for(j in 1:nmax){
z[i,j] <- dmvnorm(c(x.points[i],y.points[j]),mean=mu.vec,sigma=sigma.mat)
}
}
contour(x.points,y.points,z,xlab="x",ylab="y",col="blue")
abline(v=1,h=1)
Figura: Curvas de nível da densidade normal bivariada
É importante ressaltar que as funções do pacote mvtnorm do R NÃO estão limitadas a duas dimensões. A visualização direta da densidade multivariada é que tá limitada a no máximo duas dimensões.
3 Funções do R para simulação da normal multivariada
Podemos usar a função rmvnorm() do pacote mvtnorm do R para simular valores de uma distribuição normal multivariada. Vejamos dois exemplos:
Exemplo. Dez números aleatórios de uma normal trivariada:
set.seed(123)
mu.vec <- c(3,-1,1);mu.vec
## [1] 3 -1 1
sigma.mat <- matrix(c(3,-1,1,-1,1,0,1,0,2),nrow=3,ncol=3,byrow=TRUE);sigma.mat
## [,1] [,2] [,3]
## [1,] 3 -1 1
## [2,] -1 1 0
## [3,] 1 0 2
amostra <- rmvnorm(10,mean=mu.vec,sigma=sigma.mat)
colnames(amostra) <- c("X","Y","Z")
amostra
## X Y Z
## [1,] 2.693368 -0.8965388 2.9350820
## [2,] 3.645596 -0.8088004 3.3842808
## [3,] 4.031292 -2.3820455 0.1393931
## [4,] 1.899486 0.3192337 1.4147121
## [5,] 3.429788 -1.0905108 0.3796246
## [6,] 5.087757 -1.3694983 -1.0638211
## [7,] 3.984711 -1.7743711 -0.2552936
## [8,] 2.800873 -1.8958916 -0.1340146
## [9,] 2.921326 -2.2471782 1.8386234
## [10,] 4.129816 -2.0295690 2.7050676
pairs(amostra,pch=16,col="red")
Figura: Números aleatórios de uma normal trivariada
Exemplo. Números aleatórios de uma normal bivariada:
set.seed(123)
mu.vec <- c(1,1);mu.vec
## [1] 1 1
sigma.mat <- matrix(c(2,1,1,1),nrow=2,ncol=2,byrow=TRUE);sigma.mat
## [,1] [,2]
## [1,] 2 1
## [2,] 1 1
nmax <- 100
x.points <- seq(-3,5,length.out=nmax)
y.points <- x.points
z <- matrix(0,nrow=nmax,ncol=nmax)
for(i in 1:nmax){
for(j in 1:nmax){
z[i,j] <- dmvnorm(c(x.points[i],y.points[j]),mean=mu.vec,sigma=sigma.mat)
}
}
amostra <- rmvnorm(200,mean=mu.vec,sigma=sigma.mat)
contour(x.points,y.points,z,xlab="x",ylab="y",col="blue")
abline(v=1,h=1)
points(amostra,pch=16,col="red")
Figura: Números aleatórios de uma normal bivariada
4 A decomposição espectral da matriz de covariância
A decomposição espectral da matriz de covariâncias \(\boldsymbol{\Sigma}\) na sua forma canônica, significa escrever \(\boldsymbol{\Sigma}\) em termos dos seus autovalores e autovetores. Formalmente, a decomposição espectral de \(\boldsymbol{\Sigma}\) é dada por \[ \boldsymbol{\Sigma}=\boldsymbol{P}\boldsymbol{\Lambda}\boldsymbol{P}^{T} \] onde \(\boldsymbol{P}\) é uma matrix \(d\times d\) cujas colunas são os autovetores de \(\boldsymbol{\Sigma}\) e \(\boldsymbol{\Lambda}\) é uma matriz diagonal \(d\times d\) cujos elementos são os autovalores de \(\boldsymbol{\Sigma}\).
A raiz quadrada de \(\boldsymbol{\Sigma}\) é dada por \[ \boldsymbol{\Sigma}^{1/2}=\boldsymbol{P}\boldsymbol{\Lambda}^{1/2}\boldsymbol{P}^{T} \] e \(\boldsymbol{\Lambda}^{1/2}=\mbox{diag}(\mbox{autovalores}^{1/2})\).
Vejamos um exemplo:
sigma.mat <- matrix(c(2,1,1,1),nrow=2,ncol=2,byrow=TRUE);sigma.mat
## [,1] [,2]
## [1,] 2 1
## [2,] 1 1
de <- eigen(sigma.mat,symmetric=TRUE)
lambda <- de$values
lambda.mat <- diag(lambda); lambda.mat
## [,1] [,2]
## [1,] 2.618034 0.000000
## [2,] 0.000000 0.381966
P.mat <- de$vectors;P.mat
## [,1] [,2]
## [1,] -0.8506508 0.5257311
## [2,] -0.5257311 -0.8506508
P.mat%*%lambda.mat%*%t(P.mat)
## [,1] [,2]
## [1,] 2 1
## [2,] 1 1
lambda.mat[1,1]*(P.mat[,1]%*%t(P.mat[,1]))+lambda.mat[2,2]*(P.mat[,2]%*%t(P.mat[,2]))
## [,1] [,2]
## [1,] 2 1
## [2,] 1 15 Simulação da normal multivariada via decomposição espectral
O método da decomposição espectral afirma que \[ \boldsymbol{X}=\boldsymbol{R}\boldsymbol{Z}+\boldsymbol{\mu}\sim\mbox{N}_{d}(\boldsymbol{\mu},\boldsymbol{\Sigma}) \] se \(\boldsymbol{Z}\sim\mbox{N}_{d}(\boldsymbol{0},\boldsymbol{I})\), onde \(\boldsymbol{R}=\boldsymbol{\Sigma}^{1/2}\). Portanto o método da decomposição espectral é como segue:
Algoritmo: Método da decomposição espectral
Entrada: mu.vec, Sigma.mat
Saída: Um vetor gerado da distribuição normal multivariada
1 Gere Z.vec = (N(0,1),...,N(0,1))
2 Calcule a raiz R.mat de Sigma.mat via decomposição espetral
3 Calcule X.vec = R.mat%*%Z.vec + mu.vec
4 retorna XSe desejar uma amostra de tamanho \(n>1\) de uma distribuição normal multivariada, basta fazer: \[ \boldsymbol{A}=\boldsymbol{W}\boldsymbol{R}+\boldsymbol{M} \] onde \(\boldsymbol{W}\) é uma matrix \(n\times d\) onde \(W_{ij}\sim\mbox{N}(0,1)\) são iid’s e \(\boldsymbol{M}\) é uma matrix \(n\times d\) onde cada linhas é igual ao vetor \(\boldsymbol{\mu}\).
Exemplo. Números aleatórios de uma normal bivariada:
set.seed(321)
mu.vec <- c(1,1)
sigma.mat <- matrix(c(2,1,1,1),nrow=2,ncol=2,byrow=TRUE)
nmax <- 100
x.points <- seq(-3,5,length.out=nmax)
y.points <- x.points
z <- matrix(0,nrow=nmax,ncol=nmax)
for(i in 1:nmax){
for(j in 1:nmax){
z[i,j] <- dmvnorm(c(x.points[i],y.points[j]),mean=mu.vec,sigma=sigma.mat)
}
}
de <- eigen(sigma.mat,symmetric=TRUE)
lambda <- de$values
raizlambda.mat <- diag(sqrt(lambda));raizlambda.mat
## [,1] [,2]
## [1,] 1.618034 0.000000
## [2,] 0.000000 0.618034
P.mat <- de$vectors;P.mat
## [,1] [,2]
## [1,] -0.8506508 0.5257311
## [2,] -0.5257311 -0.8506508
R.mat <- P.mat%*%raizlambda.mat%*%t(P.mat)
n.num <- 200
W.mat <- matrix(rnorm(2*n.num),n.num,2,byrow=TRUE)
M.mat <- matrix(rep(mu.vec,n.num),n.num,2,byrow=TRUE)
amostra <- W.mat%*%R.mat+M.mat
contour(x.points,y.points,z,xlab="x",ylab="y",col="blue")
abline(v=1,h=1)
points(amostra,pch=16,col="red")
Figura: Números aleatórios de uma normal bivariada
6 Simulação de dados para o modelo de regressão linear simples
O modelo de regressão linear simples é definido como segue: \[ y_{i}=\beta_{0}+\beta_{1}x_{i}+\varepsilon_{i}\hspace{1 cm}i=1,\ldots,n \] onde \(\varepsilon_{i}\sim\mbox{N}(0,\sigma^{2})\) para todo \(i\). Note que este é um modelo estatístico paramétrico com três parâmetros desconhecidos, nomeadamente \(\beta_{0},\beta_{1}\) e \(\sigma^{2}\), mas que podem ser estimados a partir dos dados observados para as variáveis \(x\) e \(y\): \[ \mbox{Dados}\hspace{0.5cm}\Leftarrow\hspace{0.5cm}\left\{\hspace{0.3cm} \begin{array}{c|cc} \hline i & x & y \\ \hline 1 & x_{1} & y_{1} \\ 2 & x_{2} & y_{2} \\ \vdots & \vdots & \vdots \\ n & x_{n} & y_{n} \\ \hline \end{array} \right. \] Além disso, note que \(y_{i}\sim\mbox{N}(\beta_{0}+\beta_{1}x_{i},\sigma^{2})\) para todo \(i\).
O objetivo desta seção é simular dados (ou valores) para um modelo de regressão linear simples e também proporcionar o ajuste e a validação do modelo, tendo como base os dados simulados. Para simular os dados, precisamos fixar valores para os parâmetros. Sendo assim, os valores escolhidos são: \[ \beta_{0}=60\hspace{1cm}\beta_{1}=15\hspace{1cm}\sigma^{2}=2. \] A notação usada nesta seção é:
- \(\hat{\beta}_{0}\) estimativa de \(\beta_{0}\)
- \(\hat{\beta}_{1}\) estimativa de \(\beta_{1}\)
- \(\hat{\sigma}^{2}\) estimativa de \(\sigma^{2}\)
- \(\hat{\sigma}\) estimativa de \(\sigma\)
- \(\hat{y}_{i}=\hat{\beta}_{0}+\hat{\beta}_{1}x_{i}\) para todo \(i\)
- \(r_{i}=y_{i}-\hat{y}_{i}=y_{i}-(\hat{\beta}_{0}+\hat{\beta}_{1}x_{i})\) para todo \(i\) (\(i\)-ésimo resíduo)
- \(r_{i}/\hat{\sigma}\) para todo \(i\) (\(i\)-ésimo resíduo padronizado).
No R, a simulação dos dados pode ser feita da seguinte maneira:
beta0 <- 60; beta1 <- 15
sigma2 <- 2
sigma <- sqrt(sigma2)
#
# Simulação dos dados
#
set.seed(123)
x <- seq(0.1,1.9,0.08); n <- length(x)
y <- beta0 + beta1*x + rnorm(n,mean=0,sd=sigma)
#
# Os dados
#
dados <- data.frame(x,y)
dados
## x y
## 1 0.10 60.70737
## 2 0.18 62.37448
## 3 0.26 66.10435
## 4 0.34 65.19971
## 5 0.42 66.48284
## 6 0.50 69.92547
## 7 0.58 69.35183
## 8 0.66 68.11093
## 9 0.74 70.12864
## 10 0.82 71.66974
## 11 0.90 75.23111
## 12 0.98 75.20885
## 13 1.06 76.46678
## 14 1.14 77.25653
## 15 1.22 77.51392
## 16 1.30 82.02708
## 17 1.38 81.40407
## 18 1.46 79.11878
## 19 1.54 84.09187
## 20 1.62 83.63137
## 21 1.70 83.98987
## 22 1.78 86.39174
## 23 1.86 86.44901
cor(x,y,method="pearson") # Calcula a correlação linear de Pearson entre x e y
## [1] 0.9858966
cor.test(x,y) # Testa se a correlação linear de Pearson é nula
##
## Pearson's product-moment correlation
##
## data: x and y
## t = 26.996, df = 21, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9664481 0.9941055
## sample estimates:
## cor
## 0.9858966
#
# Visualização gráfica dos dados
# e a reta verdadeira y = 60 + 15x
#
plot(x,y,type="p",xlim=c(0,2),ylim=c(50,100),xlab="x",ylab="y",pch=16,main="")
abline(v=c(0,2),h=c(50,60,90))
abline(a=60,b=15,col="blue")
Figura: Visualização gráfica dos dados simulados
No R, o ajuste do modelo de regressão linear simples pode ser feito da seguinte maneira:
dados[1:5,]
## x y
## 1 0.10 60.70737
## 2 0.18 62.37448
## 3 0.26 66.10435
## 4 0.34 65.19971
## 5 0.42 66.48284
n # tamanho da amostra
## [1] 23
#
# Ajuste do modelo linear
#
modelo <- lm(y~x,data=dados)
modelo
##
## Call:
## lm(formula = y ~ x, data = dados)
##
## Coefficients:
## (Intercept) x
## 60.68 14.34
#
# Estimação pontual dos parâmetros
#
coefficients(modelo) # Estimativas pontuais para os beta's
## (Intercept) x
## 60.68243 14.33631
beta0.hat <- coefficients(modelo)[[1]]
beta0.hat
## [1] 60.68243
beta1.hat <- coefficients(modelo)[[2]]
beta1.hat
## [1] 14.33631
rss <- deviance(modelo) # Calcula RSS = Residual sum of squares
rss
## [1] 38.35762
sigma2.hat <- rss/(n-2)
sigma2.hat
## [1] 1.826553
sigma.hat <- sqrt(rss/(n-2)) # Residual standard error
sigma.hat # Estimativa pontual para o desvio-padrão (sigma) dos erros
## [1] 1.3515
#
# Estimação intervalar e Teste de hipótese
#
confint(modelo) # Estimativas intervalares para os beta's
## 2.5 % 97.5 %
## (Intercept) 59.45166 61.91321
## x 13.23193 15.44069
anova(modelo) # Testa se beta1 = 0
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x 1 1331.17 1331.17 728.79 < 2.2e-16 ***
## Residuals 21 38.36 1.83
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modelo)
##
## Call:
## lm(formula = y ~ x, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4947 -0.8937 -0.2208 0.7627 2.7074
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 60.6824 0.5918 102.5 <2e-16 ***
## x 14.3363 0.5311 27.0 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.352 on 21 degrees of freedom
## Multiple R-squared: 0.972, Adjusted R-squared: 0.9707
## F-statistic: 728.8 on 1 and 21 DF, p-value: < 2.2e-16
No R, a validação do modelo de regressão linear simples pode ser feita da seguinte maneira:
#
# Validação do modelo linear
#
y.hat <- fitted(modelo)
residuos <- residuals(modelo)
residuos.p <- residuos/(sigma.hat) # Residuos padronizados
data.frame(x,y,y.hat,residuos,residuos.p)
## x y y.hat residuos residuos.p
## 1 0.10 60.70737 62.11606 -1.4086960 -1.0423201
## 2 0.18 62.37448 63.26297 -0.8884885 -0.6574090
## 3 0.26 66.10435 64.40987 1.6944734 1.2537720
## 4 0.34 65.19971 65.55678 -0.3570638 -0.2641980
## 5 0.42 66.48284 66.70368 -0.2208419 -0.1634050
## 6 0.50 69.92547 67.85059 2.0748812 1.5352427
## 7 0.58 69.35183 68.99749 0.3543423 0.2621844
## 8 0.66 68.11093 70.14440 -2.0334631 -1.5045967
## 9 0.74 70.12864 71.29130 -1.1626576 -0.8602717
## 10 0.82 71.66974 72.43821 -0.7684668 -0.5686027
## 11 0.90 75.23111 73.58511 1.6460028 1.2179077
## 12 0.98 75.20885 74.73201 0.4768387 0.3528217
## 13 1.06 76.46678 75.87892 0.5878569 0.4349661
## 14 1.14 77.25653 77.02582 0.2307048 0.1707027
## 15 1.22 77.51392 78.17273 -0.6588069 -0.4874634
## 16 1.30 82.02708 79.31963 2.7074433 2.0032870
## 17 1.38 81.40407 80.46654 0.9375287 0.6936947
## 18 1.46 79.11878 81.61344 -2.4946595 -1.8458444
## 19 1.54 84.09187 82.76035 1.3315196 0.9852158
## 20 1.62 83.63137 83.90725 -0.2758801 -0.2041288
## 21 1.70 83.98987 85.05416 -1.0642875 -0.7874859
## 22 1.78 86.39174 86.20106 0.1906755 0.1410843
## 23 1.86 86.44901 87.34797 -0.8989555 -0.6651537
summary(residuos.p)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.8458 -0.6613 -0.1634 0.0000 0.5643 2.0033
qnorm(c(0.25,0.50,0.75),mean=0,sd=1) # Quantis da normal padrão para referência
## [1] -0.6744898 0.0000000 0.6744898
shapiro.test(residuos.p)
##
## Shapiro-Wilk normality test
##
## data: residuos.p
## W = 0.98482, p-value = 0.9698
#
# Visualização gráfica dos resíduos
#
par(mfrow=c(2,2))
plot(1:n,residuos.p,xlab="índice",ylab="r",ylim=c(-3.5,3.5),pch=16,main="Residuos padronizados")
abline(h=c(-3,0,3),col="red")
hist(residuos.p,freq=FALSE,col="orange",xlab="r",ylab="Densidade",xlim=c(-3.5,3.5),ylim=c(0,0.5),main="Residuos padronizados")
points(seq(-3.5,3.5,0.05),dnorm(seq(-3.5,3.5,0.05),0,1),type="l",col="blue")
boxplot(residuos.p,horizontal=TRUE,ylim=c(-3.5,3.5),col="yellow",xlab="r",main="Residuos padronizados")
qqnorm(residuos.p,xlim=c(-3.5,3.5),ylim=c(-3.5,3.5),pch=16,main = "Q-Q Plot Normal",xlab="Quantis teóricos",ylab ="Quantis amostrais",plot.it=T,datax=F)
abline(a=0,b=1,col="red")
Figura: Visualização gráfica dos resíduos
Por fim, podemos visualizar o modelo:
#
# Parâmetros estimados
#
beta0.hat
## [1] 60.68243
beta1.hat
## [1] 14.33631
sigma.hat
## [1] 1.3515
#
# Visualização gráfica
#
plot(x,y,type="p",xlim=c(0,2),ylim=c(50,100),xlab="x",ylab="y",pch=16,main="")
abline(v=c(0,2),h=c(50,60,90))
abline(a=60,b=15,col="blue")
abline(a=beta0.hat,b=beta1.hat,col="red")
legend(x=0.1,y=99,legend=c("Modelo real","Modelo estimado"),fill=c("blue","red"))
Figura: Regressão linear simples
Copyright © 2018 MCMC.