Distribución Normal Multivariada
Jairo A. Ayala Godoy
Introducción
La distribución normal multivariada, también conocida como la distribución gaussiana multivariada, es un concepto fundamental en estadística y probabilidad. Se utiliza para modelar conjuntos de variables aleatorias que están correlacionadas entre sí. Esta distribución es una extensión de la distribución normal univariada, que describe una sola variable aleatoria.
En el contexto de la distribución normal multivariada, estamos interesados en comprender cómo se distribuyen conjuntamente varias variables aleatorias continuas. Cada conjunto de variables aleatorias multivariadas se caracteriza por un vector de medias y una matriz de covarianza, que describen la tendencia central y la relación entre las variables, respectivamente.
La distribución normal multivariada es ampliamente utilizada en diversos campos, incluyendo la estadística, la econometría, la ingeniería y la ciencia de datos. Se utiliza para modelar datos en los que las variables están correlacionadas de alguna manera, lo que la hace especialmente valiosa para el análisis de datos multivariados y la predicción.
Función de densidad
Supongamos que tenemos \(p\) variables aleatorias normales independientes, entonces la distribución conjunta de estas variables se conoce como distribución normal multivariada que tiene como función de densidad: \[ f({\bf x}|\mu,\Sigma) = \dfrac{1}{\sqrt{(2\pi)^p |\Sigma|}} e^{-\frac{1}{2} ({\bf x}-\mu)^t \Sigma^{-1} ({\bf x}-\mu)} \] donde \(\mu\) es el vector de medias de dimensión \(p\times 1\) y \(\Sigma\) la matriz de covarianzas de dimensión \(p\times p\). El símbolo \(|\Sigma|\) se refiere al determinante de la matriz \(\Sigma\). Esta función asume que \(\Sigma\)se puede invertir, y una condición suficiente para la existencia de una inversa es que \(|\Sigma| \not = 0\). Además, esta matriz debe ser semidefinida positiva \((\Sigma_{ij}\geq 0)\), para asegurar que el punto más probable es \(\mu\) y que, a medida que \({\bf x}\) se aleja de \(\mu\) en cualquier dirección, entonces la probabilidad de observar \({\bf x}\) disminuye.
Se dice que un vector \({\bf x}\) tiene distribución normal \(p\)-variada de parámetros \(\mu\) y \(S\), donde los parámetros tendrán la dimensión adecuadas. Se denota como: \({\bf x} \sim N_p (\mu,S)\). Particularmente si tenemos solo dos variables normales univariadas \((p=2)\), la distribución se conoce como normal bivariada \({\bf x} \sim N_2 (\mu,S)\).
Ingresa a: Densidad de Normal Bivariada y explora la aplicación.
Ejemplo 1:
Veamos a modo de ejemplo como se puede visualizar la densidad de una
distribución normal bivariada. Supongamos que tenemos las variables
Peso y
Estatura con vector de medias y matriz de
covarianzas respectivamente: \[
\mu= \begin{pmatrix} 150\\ 165\end{pmatrix} \quad \text{y} \quad\Sigma=
\begin{pmatrix} 25 & 10 \\ 10 & 20 \end{pmatrix}
\] Para crear una superficie de una normal bivariada en
R se puede usar el siguiente código. La función dmvnorm
del paquete mvtnorm sirve para obtener la densidad de
cada punto.
library(mvtnorm)
mu <- c(150, 165)
Sigma <- matrix(c(25,10,10,20), nrow=2 , byrow = TRUE)
n <- 50
Peso <- seq(from=130, to=170, length.out=n)
Estatura <- seq(from=145, to=185, length.out=n)
densidad <- function(x1, x2) dmvnorm(cbind(x1, x2),mu, Sigma)
Z <- outer(Peso, Estatura, FUN="densidad")
persp(Peso, Estatura, Z, theta=45, phi=30,xlab=" ",ylab=" ", zlab=" ", ticktype="detailed", nticks=4, col="lightblue")
title(main = "Distribucion Normal Bivariada")
Ejemplo 2:
Ahora, a modo de ejemplo simularemos una normal multivariada.
Supongamos que queremos simular 1000 observaciones de las variables
Peso y
Estatura con vector de medias y matriz de
covarianzas respectivamente: \[
\mu= \begin{pmatrix} 150\\ 165\end{pmatrix} \quad \text{y} \quad \Sigma=
\begin{pmatrix} 25 & 10 \\ 10 & 20 \end{pmatrix}
\]
Para simular una normal bivariada se puede usar la función rmvnorm
del paquete mvtnorm.
Para hacer lo solicitado se puede usar el siguiente código:
library(mvtnorm)
# Datos simulados
n <- 1000
set.seed(2023)
datos <- data.frame(rmvnorm(n,mu,Sigma))
colnames(datos) <- c("Peso","Estatura")
plot(datos, xlab="Peso",ylab="Estatura", pch=19, las=1)
points(x=mu[1], y=mu[2], lwd=5, col="red", pch=19)
Ahora, usando estos datos simulados podemos construir la densidad. Usaremos el paquete plotly.
library(plotly)
densidad <- dmvnorm(datos,mu,Sigma)
plot_ly(x=~datos$Peso, y=~datos$Estatura, z=~densidad,type = "scatter3d", color=densidad, xlab=" ",ylab=" ", zlab=" ") También podemos construir el diagrama de dispersión usando el paquete scatterplot3d.
library(scatterplot3d)
scatterplot3d(x=datos[,1], y=datos[,2], z=densidad, pch=19, cex.lab=1, highlight.3d=TRUE, type="h", xlab=" ",ylab=" ", zlab=" ")
Propiedades principales
La distribución normal multivariada cumple las siguientes propiedades:
La distribución es simétrica alrededor de \(\mu\). La simetría se comprueba sustituyendo en la densidad \({\bf x}\) por \(\mu \pm a\) y observando que: \(f(\mu+ a) =f(\mu − a)\).
La distribución tiene un único máximo en \(\mu\). Al \(\Sigma\) ser semidefinida positiva el término del exponente \(({\bf x} − \mu)\Sigma^{−1}({\bf x} − \mu)\) es siempre positivo, y la densidad \(f({\bf x})\) será máxima cuando dicho término sea cero, lo que ocurre para \({\bf x} = {\mu}\).
El vector de medias es \(\mu\) y su matriz de covarianzas es \(\Sigma\). Esta propiedad puede demostrarse rigurosamente, se deduce de la comparación de las densidades univariada y multivariada.
Cualquier vector \({\bf x}\) con distribución normal multivariada puede convertirse mediante una transformación lineal en un vector \({\bf z}\) normal multivariado con vector de medias \({\bf 0}\) y matriz de covarianza \(I\). Comúnmente se conoce como normal estándar multivariada . La prueba de esta propiedad se puede visualizar así: Como\(\Sigma\)es invertible, podemos definimos una nueva variable: \[ {\bf z} = \sqrt{\Sigma^{-1}}({\bf x}-\mu) \] Entonces \({\bf x}= \mu + \sqrt{\Sigma}\,{\bf z}\). Donde es fácil ver que: \[ E[{\bf z}]= {\bf 0} \quad \text{y} \quad \sqrt{\Sigma}\cdot \Sigma^{-1}\cdot \sqrt{\Sigma}=I \] Por lo tanto, cualquier vector \({\bf x}\) de variables normales puede transformarse en otro vector \({z}\) variables normales independientes con vector de medias \({\bf 0}\) y matriz de covarianza \(I\).
Las distribuciones marginales son normales. Como las \(p\) variables son independientes la comprobación de esta propiedad es inmediata.
Cualquier subconjunto de \(h\) variables con \(h < p\) tendra distribución normal \(h\)-variada. Es una extensión del la propiedad anterior y se demuestra análogamente.
Sea \({\bf c}\) el vector de constantes, entonces \({\bf c}^t{\bf x} =c_1{\bf x}_1+c_2{\bf x}_2+\cdots c_p{\bf x}_p \sim N_p\left({\bf c}^t\mu,{\bf c}^t \Sigma {\bf c} \right)\).
Sea \(C\) cualquier matriz de constantes, entonces el vector de combinaciones lineales: \(C\cdot {\bf x} \sim N_p\left(C\mu,C\Sigma C^t\right)\).
Al cortar con hiperplanos paralelos al definido por las \(p\) variables que forman la variable vectorial, \({\bf x}\), se obtienen las curvas de nivel, cuya ecuación es: \[ ({\bf x}-\mu)^t \Sigma^{-1} ({\bf x}-\mu) = cte \] Las curvas de nivel son, por tanto, elipsoides, y definen una medida de la distancia de un punto al centro de la distribución. Esta medida es conocida en la literatura como distancia de Mahalanobis.
El vector \(\overline{\bf x}\) tiene una distribución aproximada normal multivariada. \(\overline{\bf x} \sim N_p \left( \mu,\frac{1}{n}\Sigma\right)\), cuando \(n\) es suficientemente grande. Este resultado se conoce como Teorema central del límite.
Importancia de la normal multivariada
La mayoría de las técnicas del análisis multivariado suponen que las observaciones proceden de una población normal multivariada y de aquí precede la gran importancia, puesto que de no cumplirse dicha normalidad, los resultados obtenidos por cualquiera de esas técnicas no serían confiables.
Sin embargo, si la muestra es suficientemente grande, y las técnicas empleadas solamente depende del comportamiento de \(\overline{\bf x}\) o de distancias relacionadas con \(\overline{\bf x}\) de la forma \(\left(\overline{\bf x}-\mu \right)^t \Sigma^{-1}\left(\overline{\bf x}-\mu \right)\), el supuesto de normalidad es menos crucial, debido a los resultados de las propiedades vistas. Sin embargo, la calidad de la inferencia obtenida por estos métodos depende de qué tan cercana esté la verdadera población de la distribución normal multivariada.
Por tanto es importante conocer procedimientos que permitan evaluar la normalidad de un conjunto de datos multivariado. Basados en las propiedades de la distribución normal multivariada, sabemos que todas las combinaciones lineales de las variables son normales y que los contornos de la distribución normal multivariada son elipsoides. Por tanto, en la verificación de la normalidad debería responder a:
¿Las variables marginales son normales?
¿Algunas combinaciones lineales de las variables parecen ser normales?
¿Los diagramas de dispersión de los pares de variables presentan una apariencia elíptica?
Evaluación de la normal univariada
Como hemos visto en los cursos básicos de Estadística, el supuesto de normalidad se puede corroborar tanto gráficamente como analíticamente. A continuación trabajaremos las pruebas existentes tanto gráficamente como analíticamente.
Pruebas gráficas
Trabajaremos con una base de datos simulada, donde las dos primeras variables serán normales y las otras dos no. Luego, aplicaremos las diferentes prueba para analizar que resultados obtenemos.
# simulamos dos variables normales y dos variables betas
n <- 50
mu <- c(0,2)
Sigma <- matrix(c(1,0.9,0.9,1), nrow=2 , byrow = TRUE)
set.seed(1212)
x1 <- rmvnorm(n,mu,Sigma)
x3 <- rbeta(n,5,1)
x4 <- rbeta(n,1,5)
data <- data.frame(cbind(x1,x3,x4))
colnames(data) <- c("x1","x2","x3","x4")
# Estandarizamos las variables
datas <- data.frame(scale(data))
# Visualizaciones para analizar la normalidad
# 1. Boxplot
boxplot(datas)
# 2. Histogramas
par(mfrow = c(2, 2))
x<- seq(-4, 4, length = 100)
hist(datas[,1], prob = TRUE, xlab=" ", main="Histograma x1")
lines(x,dnorm(x),col="red",lwd=2)
hist(datas[,2], prob = TRUE, xlab=" ", main="Histograma x2")
lines(x,dnorm(x),col="red",lwd=2)
hist(datas[,3], prob = TRUE, xlab=" ", main="Histograma x3")
lines(x,dnorm(x),col="red",lwd=2)
hist(datas[,4], prob = TRUE, xlab=" ", main="Histograma x4")
lines(x,dnorm(x),col="red",lwd=2)
# 3. Q-Q plot
par(mfrow = c(2, 2))
qqnorm(datas[,1], pch = 1, main="Q-Q plot x1")
qqline(datas[,1], col = "blue", lwd = 2)
qqnorm(datas[,2], pch = 1, main="Q-Q plot x2")
qqline(datas[,2], col = "blue", lwd = 2)
qqnorm(datas[,3], pch = 1, main="Q-Q plot (x3")
qqline(datas[,3], col = "blue", lwd = 2)
qqnorm(datas[,4], pch = 1, main="Q-Q plot x4")
qqline(datas[,4], col = "blue", lwd = 2)
# 4 Diagrama resumen y su correlación
library(GGally)
ggpairs(datas) 
Pruebas analíticas
A continuación, trabajamos las diferentes pruebas de hipótesis más empleados para analizar la normalidad de los datos. En todos ellos, se consideran como hipótesis: \[\begin{align*} H_0 &: \text{Los datos son normales.} \\ H_a &: \text{Los datos no son normales.} \end{align*}\]
El p-value de estas pruebas indican la probabilidad de obtener una distribución como la observada si los datos proceden realmente de una población con una distribución normal.
Cuando estas pruebas se emplean con la finalidad de verificar las condiciones de métodos paramétricos, por ejemplo, un PCA o un ANOVA, es importante tener en cuenta que, al tratarse de p-values, cuanto mayor sea el tamaño de la muestra más poder estadístico tienen y más fácil es encontrar evidencias en contra de la \(H_0\) (normalidad). Al mismo tiempo, cuanto mayor sea el tamaño de la muestra, menos sensibles son los métodos paramétricos a la falta de normalidad. Por esta razón, es importante no basar las conclusiones únicamente en el p-value, sino también considerar la representación gráfica y el tamaño de la muestra.
La función mvn
proporciona varias pruebas de normalidad univariante ampliamente
utilizadas, incluido “SW” para Shapiro-Wilk, “CVM” para Cramer-von
Mises, “Lillie” para Lilliefors, “SF” para Shapiro-Francia y “AD” para
Anderson-Darling.
# Pruebas de normalidad univariada
library(MVN)
# 1. Shapiro-Wilks
SW <- mvn(datas, univariateTest = "SW",desc=T)
SW## $multivariateNormality
## Test HZ p value MVN
## 1 Henze-Zirkler 1.131717 0.001731672 NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Shapiro-Wilk x1 0.9885 0.9045 YES
## 2 Shapiro-Wilk x2 0.9713 0.2617 YES
## 3 Shapiro-Wilk x3 0.7638 <0.001 NO
## 4 Shapiro-Wilk x4 0.9188 0.0021 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th
## x1 50 2.070284e-18 1 -0.07607341 -2.148317 2.3641033 -0.8128716
## x2 50 9.172784e-17 1 -0.01818184 -1.654874 2.3031836 -0.7333739
## x3 50 -3.030508e-16 1 0.35961005 -3.397328 0.9807097 -0.2238720
## x4 50 4.228605e-17 1 -0.27345313 -1.230481 2.6702863 -0.8049994
## 75th Skew Kurtosis
## x1 0.6506322 0.1651482 -0.6277562
## x2 0.6116685 0.3024705 -0.4992938
## x3 0.5618706 -1.8352584 2.8350040
## x4 0.7829698 0.7469335 -0.4059265Usa esta prueba cuando el tamaño de la muestra sea pequeña \((<30)\).
# 2. Cramer-Von Mises
CVM <- mvn(datas, univariateTest = "CVM",desc=T)
CVM$univariateNormality## Test Variable Statistic p value Normality
## 1 Cramer-von Mises x1 0.0246 0.9119 YES
## 2 Cramer-von Mises x2 0.0302 0.84 YES
## 3 Cramer-von Mises x3 0.7337 <0.001 NO
## 4 Cramer-von Mises x4 0.2124 0.0035 NOSolo imprimimos la salidas de la prueba univariada, ya que las demás salidas se repiten.
# 3. Lilliefors (correccion de Kolmogorov)
L <- mvn(datas, univariateTest = "Lillie",desc=T)
L$univariateNormality## Test Variable Statistic p value Normality
## 1 Lilliefors (Kolmogorov-Smirnov) x1 0.0579 0.9443 YES
## 2 Lilliefors (Kolmogorov-Smirnov) x2 0.0653 0.8584 YES
## 3 Lilliefors (Kolmogorov-Smirnov) x3 0.2440 <0.001 NO
## 4 Lilliefors (Kolmogorov-Smirnov) x4 0.1298 0.0348 NOLa prueba de Kolmogorov-Smirnov permite estudiar si una muestra procede de una población con una determinada distribución (media y desviación típica), no está limitado únicamente a la distribución normal. A pesar de que ser una excelente prueba para contrastar la normalidad, es poco usada, ya que es necesario conocer la media y varianza poblacional, lo que en la mayoría de los casos no es posible. Para solventar este problema, se desarrolló una modificación del Kolmogorov-Smirnov conocida como la prueba Lilliefors. Esta prueba asume que la media y varianza son desconocidas, estando especialmente desarrollada para contrastar la normalidad. Es la alternativa cuando el número de observaciones es mayor de 30.
# 4. Shapiro Francia
SF <- mvn(datas, univariateTest = "SF",desc=T)
SF$univariateNormality## Test Variable Statistic p value Normality
## 1 Shapiro-Francia x1 0.9908 0.9183 YES
## 2 Shapiro-Francia x2 0.9778 0.3907 YES
## 3 Shapiro-Francia x3 0.7614 <0.001 NO
## 4 Shapiro-Francia x4 0.9247 0.0049 NO# 5. Anderson Darling
AD <- mvn(datas, univariateTest = "AD",desc=T)
AD$univariateNormality## Test Variable Statistic p value Normality
## 1 Anderson-Darling x1 0.1843 0.9039 YES
## 2 Anderson-Darling x2 0.2856 0.612 YES
## 3 Anderson-Darling x3 4.1025 <0.001 NO
## 4 Anderson-Darling x4 1.2905 0.0021 NODado los resultados de las pruebas, podemos comprobar que las dos primeras variables son normales, pero las otras no lo cumplen que debería ser.
Evaluación de la normal multivariada
Es importante verificar el cumplimiento de este supuesto para que los
métodos que trabajaremos durante la segunda parte de este curso tengan
validez. En la literatura estadística, existen diferentes pruebas que
permiten verificar la normalidad multivariada. Sin embargo, aun es tema
de investigación. Los criterios para determinar cuál es la prueba más
adecuada que se debe utilizar bajo ciertas condiciones como: tamaño de
muestra y número de variables. Usaremos las pruebas de: Mardia,
Henze-Zirkler, Royston, Doornik-Hansen y Energy que se encuentran dentro
de la librería MVN
# 1. Mardia
Mardia <- mvn(datas, mvnTest = "mardia")
Mardia$multivariateNormality## Test Statistic p value Result
## 1 Mardia Skewness 45.4917531388445 0.000945978399141201 NO
## 2 Mardia Kurtosis 0.266465493921297 0.78988073086202 YES
## 3 MVN <NA> <NA> NOPara esta prueba se espera que los datos cumplan el supuesto de normalidad si logra conseguir la normalidad tanto en la asimetría como en la curtosis.
# 2. Henze-Zirkler
HZ <- mvn(datas, mvnTest = "hz")
HZ$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 1.131717 0.001731672 NO# 3. Royston
Royston <- mvn(datas, mvnTest = "royston")
Royston$multivariateNormality## Test H p value MVN
## 1 Royston 34.24767 3.71359e-07 NO# 4 Doornik-Hansen
DH <- mvn(datas, mvnTest = "royston")
DH$multivariateNormality## Test H p value MVN
## 1 Royston 34.24767 3.71359e-07 NO# 5 Energy
Energy<- mvn(data, mvnTest = "energy")
Energy$multivariateNormality## Test Statistic p value MVN
## 1 E-statistic 1.482825 0.001 NOVeamos que si eliminamos las dos variables que no son normales, los resultados cambian.
datas1 <-datas[,1:2]
# 1 Mardia
Mardia <- mvn(datas1, mvnTest = "mardia")
Mardia$multivariateNormality## Test Statistic p value Result
## 1 Mardia Skewness 3.32523133752659 0.504944531171233 YES
## 2 Mardia Kurtosis 0.233010966187303 0.815752890340689 YES
## 3 MVN <NA> <NA> YES# 2. Henze-Zirkler
HZ <- mvn(datas1, mvnTest = "hz")
HZ$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 0.579085 0.2981517 YES# 3. Royston
Royston <- mvn(datas1, mvnTest = "royston")
Royston$multivariateNormality## Test H p value MVN
## 1 Royston 1.014103 0.4943198 YES# 4 Doornik-Hansen
DH <- mvn(datas1, mvnTest = "royston")
DH$multivariateNormality## Test H p value MVN
## 1 Royston 1.014103 0.4943198 YES# 5 Energy
Energy<- mvn(datas1, mvnTest = "energy")
Energy$multivariateNormality## Test Statistic p value MVN
## 1 E-statistic 0.706207 0.306 YESEjercicio de clase:
Analice la normalidad de la base de datos Tasas,
la cual cuenta con 4 variables (TN: tasa de natalidad,
TM: tasa de mortalidad, EV: esperanza
de vida, EVM: esperanza de vida en mujeres y
EVH: esperanza de vida en hombres), medidas en 194
países alrededor del mundo.