Distribucion Normal Multivariante

Tema de investigación

Se realiza un estudio de analisis de las variables que influyen en los sueldos de los trabajadores en una empresa.

Las variables a estudiar son la siguientes :

X1: Sueldo de trabajador en soles
X2: Tiempo de Servicio del trabajador en años
X3: Edad del trabajador en años
X4: Años de estudio

Asumimos que las variables provienen de una Distribución Normal Multivariante.

Actividad previa: iniciamos el entorno R y referenciando a librerias:

# Para limpiar el workspace, por si hubiera algun dataset 
# o informacion cargada
rm(list = ls())

#Referenciando el paquete de matrices
library(MASS)

#Referenciando el paquete de multivariado
#install.packages("MVN",dependencies = TRUE )
library(MVN)

#Refenciando el paquete para verificar si matriz es definida positiva
library(matrixcalc)

Distribucion Normal Bivariante

Definimos una Población Normal Bivariada, verificando que la matriz de varianzas y covarianzas sea simétrica y definida positiva.

#----CREANDO LA POBLACION BIVARIADA
#Numero de variables y cantidad de observaciones
N = 500000
p = 2
n = 50 


#Creando el vector de medias poblacional
U = matrix(c(2191.149915 ,4.662904),nrow = p,  ncol = 1)
U

##             [,1]
## [1,] 2191.149915
## [2,]    4.662904

#Creando la matriz de varianzas y covarianzas
V = matrix(c(2711321.073,3083.095,3083.095,19.265652),nrow = p,  ncol = p)
V

##             [,1]       [,2]
## [1,] 2711321.073 3083.09500
## [2,]    3083.095   19.26565

#Validando matriz V sea simetrica y definida postiva
if (identical(V,t(V))==TRUE) {
print("La Matriz es Simétrica")
} else {
print("La Matriz es No es Simétrica")
}

## [1] "La Matriz es Simétrica"

if (is.positive.definite(V) == TRUE ){
  print("La Matriz es Definida Positiva")
} else {
  print("La Matriz es No es Definida Positiva")
}

## [1] "La Matriz es Definida Positiva"

Generamos la Población Norma Bivariada de tamaño N = 500000 y mostramos el histrograma de cada una de variables

# Configuramos la semilla de generación de aleatorios
set.seed(1)

#Generamos la población de tamaño N
T = mvrnorm(N,U, V)

head(T)

##           [,1]      [,2]
## [1,] 3222.6695 10.122491
## [2,] 1888.7639  1.899650
## [3,] 3567.0967 12.188312
## [4,] -435.6548  2.597481
## [5,] 1648.5805  3.154567
## [6,] 3542.1378  8.609249

hist(T[,1],main= 'Histograma  X1', xlab = 'X1', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(T[,1]),col="red")

hist(T[,2],main= 'Histograma  X2', xlab = 'X2', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(T[,2]),col="blue")

Se observa que el histograma de las variables X1 y X2 de la población corresponde a una curva normal.

1. Genere una muestra aleatoria de tamaño 50

Generamos la muestra de tamaño 50 y mostramos la 6 primeras fila

#----TOMANDO UNA MUESTRA LA POBLACION BIVARIADA

set.seed(1)

M = T[sample(1:N,n,replace=F),]
head(M)

##          [,1]      [,2]
## [1,] 3693.467  2.894912
## [2,] 3556.567  2.646271
## [3,] 1815.180  4.514337
## [4,] 5401.018  8.109380
## [5,] 2751.476 13.715344
## [6,] 4789.793  8.244346

hist(M[,1],main= 'Histograma  X1', xlab = 'X1', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(M[,1]),col="red")

hist(M[,2],main= 'Histograma  X2', xlab = 'X2', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(M[,2]),col="blue")

Nótese que los histogramas de las variables X1 y X2 generan curvas con sesgo casi central y platicurticas.

2. Vector de medias, matriz de varianza y convarianza, matriz de correlación y varianza generalizada

#Vector de medias 
X = colMeans(M)
X

## [1] 2260.339743    3.580762

#Matriz de covarianza
S = cov(M)
S

##             [,1]       [,2]
## [1,] 2694302.361 3962.82720
## [2,]    3962.827   27.94453

#Matriz de correlación
R = cor(M)
R

##           [,1]      [,2]
## [1,] 1.0000000 0.4567027
## [2,] 0.4567027 1.0000000

#Varianza Generalizada
det(S)

## [1] 59587001

Se puede observar que el vector de medias muestral y la varianza muestral se aproxima a la poblacional.

3. Distancia de Mahalanobis

#Distancia de Mahalanobis
d2M =  mahalanobis(M,X,S) 

d2M

##  [1] 1.11520280 0.98856857 0.18762054 3.66137952 4.09522481 2.41491810
##  [7] 2.70910131 0.45073511 0.63122462 1.52050375 0.08057830 5.37136960
## [13] 4.22788572 2.94342195 0.29869836 3.43480427 0.84886851 0.22911557
## [19] 4.51562444 1.75165558 2.98391768 0.56158553 2.08816149 5.48021547
## [25] 1.45510552 2.56131254 4.17112404 1.19968343 1.39104025 8.11952927
## [31] 0.24971625 2.43437395 0.49110722 1.99030073 0.00900876 2.06269587
## [37] 3.35918597 1.97441073 0.14127664 1.32462682 0.49980885 3.83529424
## [43] 2.42966178 0.79789307 1.48547771 1.75234713 0.15845717 0.52605287
## [49] 0.04886660 0.94126100

4. Test de Multinomalidad: Método Gráfico

#Multinormalidad Test gráfico Q-Q Plot
plot(qchisq(((1:nrow(M)) - 1/2)/nrow(M),df=p),sort( d2M ) )
abline(a=0, b=1,col="red")

Graficamente se visualiza que el cuadrado de la distancia de Mahalanobis se ajusta una distribución Chi2, en consecuencia es muy probable que la muestra provenga de una distribución normal multivariada.

5. Test de Multinomalidad: Método Sesgo de Mardia

#Multinormalidad Test Sesgo de Mardia

mvn(M,subset = NULL,mvn = "mardia", covariance = FALSE,showOutliers = FALSE)

## $multivariateNormality
##              Test         Statistic           p value Result
## 1 Mardia Skewness  1.54300753341728 0.818996887825243    YES
## 2 Mardia Kurtosis -1.08908015020925 0.276118542449808    YES
## 3             MVN              <NA>              <NA>    YES
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk  Column1     0.9900    0.9467    YES   
## 2 Shapiro-Wilk  Column2     0.9769    0.4304    YES   
## 
## $Descriptives
##    n        Mean     Std.Dev     Median          Min        Max
## 1 50 2260.339743 1641.433021 2142.05963 -1223.733866 5841.97542
## 2 50    3.580762    5.286258    3.17326    -7.945419   16.18597
##          25th       75th       Skew   Kurtosis
## 1 1105.026297 3532.06353 0.03492158 -0.6403271
## 2    1.249477    7.35956 0.08734077 -0.1038849

Es Test indica que la muestra proviene de una distribución normal. Tanto el Test de Sesgo indica como el Test de Kurtosis indican que provienen de una distribución normal.

Distribucion Normal Multivariante

#----CREANDO LA POBLACION BIVARIADA
#Numero de variables y cantidad de observaciones
N = 500000
p = 4
n = 250
f = 0.3


#Creando el vector de medias poblacional
U = matrix(c(2191.149915,4.662904,33.434251 ,18.896082 ),nrow = p,  ncol = 1)
U

##             [,1]
## [1,] 2191.149915
## [2,]    4.662904
## [3,]   33.434251
## [4,]   18.896082

#Creando la matriz de varianzas y covarianzas
V = matrix( ,nrow = p,  ncol = p)

V[1,1] = 2711321.073
V[1,2] = 3083.095
V[1,3] = 7904.764
V[1,4] = 1503.618
V[2,1] = 3083.095
V[2,2] = 19.265652
V[2,3] = 20.088258
V[2,4] = 1.881349 
V[3,1] = 7904.764
V[3,2] = 20.088258
V[3,3] = 60.720102
V[3,4] = 4.917222
V[4,1] = 1503.618
V[4,2] = 1.881349
V[4,3] = 4.917222
V[4,4] = 4.185428

V

##             [,1]        [,2]        [,3]        [,4]
## [1,] 2711321.073 3083.095000 7904.764000 1503.618000
## [2,]    3083.095   19.265652   20.088258    1.881349
## [3,]    7904.764   20.088258   60.720102    4.917222
## [4,]    1503.618    1.881349    4.917222    4.185428

#Validando matriz V sea simetrica y definida postiva
if (identical(V,t(V))==TRUE) {
print("La Matriz es Simétrica")
} else {
print("La Matriz es No es Simétrica")
}

## [1] "La Matriz es Simétrica"

if (is.positive.definite(V) == TRUE ){
  print("La Matriz es Definida Positiva")
} else {
  print("La Matriz es No es Definida Positiva")
}

## [1] "La Matriz es Definida Positiva"

Se verifica que la matriz de varianza y covarianza s simetrica y definida postiva.

Generamos una población normal multivariado de tamano N = 500000 y mostramos el histograma de cada una de las variables.

#Generamos la población de tamaño N
set.seed(1)

T = mvrnorm(N,U, V)

hist(T[,1],main= 'Histograma  X1', xlab = 'X1', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(T[,1]),col="red")

hist(T[,2],main= 'Histograma  X2', xlab = 'X2', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(T[,2]),col="blue")

hist(T[,3],main= 'Histograma  X3', xlab = 'X3', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(T[,3]),col="green")

hist(T[,4],main= 'Histograma  X4', xlab = 'X4', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(T[,4]),col="yellow")

Observamos que las variables X1, X2, X3 y X4 se ajustan a una curva normal.

6. Genere una muestra aleatoria de tamaño 250

#----TOMANDO UNA MUESTRA LA POBLACION BIVARIADA
set.seed(1)

M = T[sample(1:N,n,replace=F),]
head(M)

##            [,1]      [,2]     [,3]     [,4]
## [1,]   688.8541 -2.226566 25.05052 17.70339
## [2,]   825.7563  2.533415 23.38260 16.32519
## [3,]  2567.1176  5.396072 34.89823 19.09926
## [4,] -1018.7171  3.638471 22.60512 17.87042
## [5,]  1630.7713  9.617298 44.40410 19.04934
## [6,]  -407.4966  7.333977 24.63084 16.69613

hist(M[,1],main= 'Histograma  X1', xlab = 'X1', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(M[,1]),col="red")

hist(M[,2],main= 'Histograma  X2', xlab = 'X2', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(M[,2]),col="blue")

hist(M[,3],main= 'Histograma  X3', xlab = 'X3', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(M[,3]),col="red")

hist(M[,4],main= 'Histograma  X4', xlab = 'X4', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(M[,4]),col="blue")

Se observa que las variables X1, X2, X3 y X4 se aproximan a una curva normal sin llegar a serlo.

7. Vector de medias, matriz de varianza y convarianza, matriz de correlación y varianza generalizada

#Vector de medias 
X = colMeans(M)
X

## [1] 2238.18857    4.93132   33.55510   19.01999

#Matriz de covarianza
S = cov(M)
S

##             [,1]        [,2]        [,3]        [,4]
## [1,] 2295195.555 2004.830954 6159.164084 1173.723213
## [2,]    2004.831   16.028448   16.120738    1.333562
## [3,]    6159.164   16.120738   56.276155    3.513265
## [4,]    1173.723    1.333562    3.513265    4.088677

#Matriz de correlación
R = cor(M)
R

##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.3305386 0.5419381 0.3831461
## [2,] 0.3305386 1.0000000 0.5367561 0.1647313
## [3,] 0.5419381 0.5367561 1.0000000 0.2316100
## [4,] 0.3831461 0.1647313 0.2316100 1.0000000

G = det(S)
G

## [1] 3612572274

El vector de medias y varianza muestral se aproxima a la poblacional.

8. Distancia de Mahalanobis

#Distancia de Mahalanobis
d2M =  mahalanobis(M,X,S)
head(d2M)

## [1] 3.4515009 2.9854564 0.0552654 5.2213606 4.3882178 5.8336460

9. Test de Multinomalidad: Método Gráfico

#Multinormalidad Test gráfico Q-Q Plot
plot(qchisq(((1:nrow(M)) - 1/2)/nrow(M),df=p),sort( d2M ) )
abline(a=0, b=1,col="red")

Graficamente se visualiza que el cuadrado de la distancia de Mahalanobis se aproxima una distribución Chi2, en consecuencia es muy probable que la muestra provenga de una distribución normal multivariada.

10. Test de Multinomalidad: Método Sesgo de Mardia

#Multinormalidad Test gráfico Q-Q Plot
mvn(M,subset = NULL,mvn = "mardia", covariance = FALSE,showOutliers = FALSE)

## $multivariateNormality
##              Test          Statistic           p value Result
## 1 Mardia Skewness   10.6908018656011 0.953850190417407    YES
## 2 Mardia Kurtosis -0.976249029472967 0.328941071829227    YES
## 3             MVN               <NA>              <NA>    YES
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk  Column1     0.9953    0.6483    YES   
## 2 Shapiro-Wilk  Column2     0.9950    0.5961    YES   
## 3 Shapiro-Wilk  Column3     0.9928    0.2621    YES   
## 4 Shapiro-Wilk  Column4     0.9939    0.4029    YES   
## 
## $Descriptives
##     n       Mean     Std.Dev      Median          Min        Max
## 1 250 2238.18857 1514.990282 2275.283940 -1568.809816 6036.80154
## 2 250    4.93132    4.003554    5.282085    -6.573834   16.11905
## 3 250   33.55510    7.501743   33.815368    17.006538   54.88107
## 4 250   19.01999    2.022048   19.062111    14.408722   24.41366
##          25th       75th         Skew   Kurtosis
## 1 1155.774426 3273.49718 -0.007446467 -0.3599870
## 2    2.227946    7.60682 -0.112476801 -0.2425099
## 3   28.249977   38.63460  0.130933546 -0.3212622
## 4   17.598079   20.52005  0.050267033 -0.5221880

Es Test indica que la muestra proviene de una distribución normal.Tanto en el Test de Sesgo y el de Test de Kurtosis indica que si.

11. Contamine el 30% de los datos de la muestra

#Contaminando el f=30% de la muestra
a = matrix(,nrow = n*f, ncol = p)
a[,1] = sample(1:1000,1,replace=FALSE)
a[,2] = sample(1:5,1,replace=FALSE)
a[,3] = sample(1:30,1,replace=FALSE)
a[,4] = sample(1:10,1,replace=FALSE)

head(a)

##      [,1] [,2] [,3] [,4]
## [1,]  472    3   27   10
## [2,]  472    3   27   10
## [3,]  472    3   27   10
## [4,]  472    3   27   10
## [5,]  472    3   27   10
## [6,]  472    3   27   10

#seleccionado las filas a contamina
set.seed(1)
i= sample(1:n,n*f,replace=FALSE)
i

##  [1] 249  68 167 129 162 215  43  14 210 187  51 225  85  21 106 182  74
## [18]   7  73  79 213  37 105 217 110 165  34 236 126  89 172 207  33  84
## [35] 163  70 234  42 166 111 148 156  20  44 121  87 245 206  40 220  25
## [52] 119 198 122  39 179 230 134  24 160 243 130  45 146  22 115 104 204
## [69] 183 184 103  75  13 202 176

#Copiamos la matriz de datos
MC = M

  #Contaminamos
  MC[i,] = M[i,] + a

#Verificamos la correcta contaminación
head(MC - M)

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]    0    0    0    0
## [3,]    0    0    0    0
## [4,]    0    0    0    0
## [5,]    0    0    0    0
## [6,]    0    0    0    0

hist(MC[,1],main= 'Histograma  X1', xlab = 'X1', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(MC[,1]),col="red")

hist(MC[,2],main= 'Histograma  X2', xlab = 'X2', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(MC[,2]),col="blue")

hist(MC[,3],main= 'Histograma  X3', xlab = 'X3', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(MC[,3]),col="red")

hist(MC[,4],main= 'Histograma  X4', xlab = 'X4', ylab = 'Densidad',probability = TRUE,breaks = 100)
lines(density(MC[,4]),col="blue")

Se puede observar que algunas variables contaminadas son bimodales o asimetricas, por lo tanto distan de ser normales.

12. Para la nueva BD responda las preguntas de 7 a 11

#Vector de medias 
XC = colMeans(MC)
XC

## [1] 2379.78857    5.83132   41.65510   22.01999

#Matriz de covarianza
SC = cov(MC)
SC

##             [,1]        [,2]       [,3]        [,4]
## [1,] 2326526.738 2283.376613 8469.37845 1985.809675
## [2,]    2283.377   18.303570   35.34664    8.177476
## [3,]    8469.378   35.346641  218.05757   60.941188
## [4,]    1985.810    8.177476   60.94119   24.435641

#Matriz de correlación
RC = cor(MC)
RC

##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.3499094 0.3760209 0.2633733
## [2,] 0.3499094 1.0000000 0.5594930 0.3866691
## [3,] 0.3760209 0.5594930 1.0000000 0.8348596
## [4,] 0.2633733 0.3866691 0.8348596 1.0000000

#varianza Global
GC = det(SC)
GC

## [1] 37813649762

#Distancia de Mahalanobis
d2MC =  mahalanobis(MC,XC,SC) 
head(d2MC)

## [1] 3.8119923 1.9655878 0.4384603 5.6893067 3.1971435 5.8565515

#Multinormalidad Test gráfico Q-Q Plot
plot(qchisq(((1:nrow(MC)) - 1/2)/nrow(MC),df=p),sort( d2MC ) )
abline(a=0, b=1,col="red")

#Multinormalidad Test Sesgo de Mardia
mvn(MC,subset = NULL,mvn = "mardia", covariance = FALSE,showOutliers = FALSE)

## $multivariateNormality
##              Test         Statistic            p value Result
## 1 Mardia Skewness   32.334380087388 0.0398630569740468     NO
## 2 Mardia Kurtosis -1.72584981794673  0.084374434663089    YES
## 3             MVN              <NA>               <NA>     NO
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk  Column1     0.9953  0.6428      YES   
## 2 Shapiro-Wilk  Column2     0.9960  0.7689      YES   
## 3 Shapiro-Wilk  Column3     0.9287  <0.001      NO    
## 4 Shapiro-Wilk  Column4     0.8974  <0.001      NO    
## 
## $Descriptives
##     n       Mean     Std.Dev      Median          Min        Max
## 1 250 2379.78857 1525.295623 2402.578874 -1568.809816 6508.80154
## 2 250    5.83132    4.278267    6.027674    -6.573834   17.15741
## 3 250   41.65510   14.766773   37.528879    17.546708   75.73119
## 4 250   22.01999    4.943242   20.343486    14.505751   34.41366
##          25th        75th        Skew   Kurtosis
## 1 1257.241562 3356.558444  0.05433157 -0.3020213
## 2    3.185285    8.707033 -0.06662191 -0.2269379
## 3   29.565471   53.599932  0.57951416 -0.8315993
## 4   18.218991   26.725507  0.69431288 -0.8196847

Es Test Mardia indica que la muestra contaminada no se ajusta a una Distribucion Nornal Multivariada.

13. Conclusiones y comentarios generales

Para generar una Población Normal Multivariada se debe verificar que la matriz de varianza y covariaza sea simetrica y definida positiva.
El Test de Mardia es más potente cuando el tamaño de la muestra es más grande.
El Test gráfico no es 100% concluyente cuando la d2M se ajusta a la Chi2.
El test de Mardia requiere que tanto es Test de Sesgo como el Test de Kurtosis No Rechacen Ho para concluir que la muestra proviene de una Distribución Normal Multivariado.

Distribucion Normal Multivariante

Alfredo Zea Valenzuela

07 de mayo del 2019

Tema de investigación

Distribucion Normal Bivariante

1. Genere una muestra aleatoria de tamaño 50

2. Vector de medias, matriz de varianza y convarianza, matriz de correlación y varianza generalizada

3. Distancia de Mahalanobis

4. Test de Multinomalidad: Método Gráfico

5. Test de Multinomalidad: Método Sesgo de Mardia

Distribucion Normal Multivariante

6. Genere una muestra aleatoria de tamaño 250

7. Vector de medias, matriz de varianza y convarianza, matriz de correlación y varianza generalizada

8. Distancia de Mahalanobis

9. Test de Multinomalidad: Método Gráfico

10. Test de Multinomalidad: Método Sesgo de Mardia

11. Contamine el 30% de los datos de la muestra

12. Para la nueva BD responda las preguntas de 7 a 11

13. Conclusiones y comentarios generales