Función de Test de Mardia - Muestras Pequeñas
mardia<-function(datos){
n<-nrow(datos)
k<-ncol(datos)
x<-as.matrix(datos)
media<-colMeans(datos)
matriz.media<-matrix(media,nrow=n,ncol=k,byrow=TRUE)
S<-cov(datos)*(n-1)/n
SI<-solve(S) #Inversa de la matriz de covarianzas
distm<-((x-matriz.media)%*%SI)%*%t((x-matriz.media))
gl<-k*(k+1)*(k+2)/6
c<-(n+1)*(n+3)*(k+1)/(n*((n+1)*(k+1)-6))
AM<-sum(distm^3)/n^2
#Muestras Pequeñas
TAM<-n*c/6*AM
valor_p1<-pchisq(TAM,gl,lower.tail = FALSE)
KM<-sum((diag(distm))^2)/n
TKM<-(KM-k*(k+2))*sqrt(n/(8*k*(k+2)))
TKM2<-TKM^2
valor_p2<-pchisq(TKM2,1,lower.tail = FALSE)
cat("Test de Normalidad Multivariada", "\n" ,
"Ho: Los datos siguen una Distribución Normal Multivariada","\n")
cat("----------------------------------------------------------------","\n")
cat("Prueba de Mardia - Asimetría","\n",
"Asimetría Multivariada =", AM, "\n",
"Valor de TAM =",TAM, "\n",
"Valor p=", valor_p1, "\n")
cat("----------------------------------------------------------------","\n")
cat("Prueba de Mardia - Curtosis","\n",
"Curtosis Multivariado =", KM, "\n",
"Valor de TKM =",TKM2, "\n",
"Valor p=", valor_p2, "\n")
}
cereales<-read.csv("cereales.csv",sep=";")
head(cereales)
## x1 x2 x3 x4 x5
## 1 H 38 1.0 100 2
## 2 H 34 1.0 120 2
## 3 H 21 1.5 110 12
## 4 H 23 1.5 110 9
## 5 H 23 1.5 110 8
## 6 H 51 2.0 180 17
library("MVN")
## Warning: package 'MVN' was built under R version 3.5.2
## sROC 0.1-2 loaded
NM<-mvn(cereales[,2:5])
NM$multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 38.8450836931911 0.00696934780762377 NO
## 2 Mardia Kurtosis 0.216049907122914 0.828948845824358 YES
## 3 MVN <NA> <NA> NO
Ejemplo Muestras Grandes
mardia(cereales[,2:5])
## Test de Normalidad Multivariada
## Ho: Los datos siguen una Distribución Normal Multivariada
## ----------------------------------------------------------------
## Prueba de Mardia - Asimetría
## Asimetría Multivariada = 12.61323
## Valor de TAM = 38.84508
## Valor p= 0.006969348
## ----------------------------------------------------------------
## Prueba de Mardia - Curtosis
## Curtosis Multivariado = 24.80009
## Valor de TKM = 0.04667756
## Valor p= 0.8289488
mardia.df<-function(datos){
n<-nrow(datos)
k<-ncol(datos)
x<-as.matrix(datos)
media<-colMeans(datos)
matriz.media<-matrix(media,nrow=n,ncol=k,byrow=TRUE)
S<-cov(datos)*(n-1)/n
SI<-solve(S) #Inversa de la matriz de covarianzas
distm<-((x-matriz.media)%*%SI)%*%t((x-matriz.media))
gl<-k*(k+1)*(k+2)/6
c<-(n+1)*(n+3)*(k+1)/(n*((n+1)*(k+1)-6))
AM<-sum(distm^3)/n^2
#Muestras Pequeñas
TAM<-n*c/6*AM
valor_p1<-pchisq(TAM,gl,lower.tail = FALSE)
if (valor_p1<0.05){
test1 <- "No Normalidad"
} else {
test1 <- "Normalidad"
}
KM<-sum((diag(distm))^2)/n
TKM<-(KM-k*(k+2))*sqrt(n/(8*k*(k+2)))
TKM2<-TKM^2
valor_p2<-pchisq(TKM2,1,lower.tail = FALSE)
if (valor_p2<0.05){
test2 <- "No Normalidad"
} else {
test2 <- "Normalidad"
}
df<-data.frame(Prueba=c("Prueba de Asimetría","Prueba de Curtosis"),Medida=c(AM,KM),Estadístico=c(TAM,TKM2),valor_p=c(valor_p1,valor_p2),Conclusión=c(test1,test2))
df
}
mardia.df(cereales[,2:5])
## Prueba Medida Estadístico valor_p Conclusión
## 1 Prueba de Asimetría 12.61323 38.84508369 0.006969348 No Normalidad
## 2 Prueba de Curtosis 24.80009 0.04667756 0.828948846 Normalidad
Función de Test de Mardia - Muestras pequeñas y grandes
mardia2<-function(datos){
n<-nrow(datos)
k<-ncol(datos)
x<-as.matrix(datos)
media<-colMeans(datos)
matriz.media<-matrix(media,nrow=n,ncol=k,byrow=TRUE)
S<-cov(datos)*(n-1)/n
SI<-solve(S) #Inversa de la matriz de covarianzas
distm<-((x-matriz.media)%*%SI)%*%t((x-matriz.media))
gl<-k*(k+1)*(k+2)/6
# Medida de Asimetría Multivariada
AM<-sum(distm^3)/n^2
if (n<30){
# Test Multivariado de Asimetría de Mardia - Muestras Pequeñas
c<-(n+1)*(n+3)*(k+1)/(n*((n+1)*(k+1)-6))
TAM<-n*c/6*AM
valor_p1<-pchisq(TAM,gl,lower.tail = FALSE)
} else {
# Test Multivariado de Asimetríıa de Mardia - Muestras grandes
TAM<-n/6*AM
valor_p1<-pchisq(TAM,gl,lower.tail = FALSE)
}
#Medida de Kurtosis Multivariado
KM<-sum((diag(distm))^2)/n
#Test Multivariado de Kurtosis de Mardia
TKM<-(KM-k*(k+2))*sqrt(n/(8*k*(k+2)))
TKM2<-TKM^2
valor_p2<-pchisq(TKM2,1,lower.tail = FALSE)
cat("Test de Normalidad Multivariada", "\n" ,
"Ho: Los datos siguen una Distribución Normal Multivariada","\n")
cat("----------------------------------------------------------------","\n")
cat("Prueba de Mardia - Asimetría","\n",
"Asimetría Multivariada =", AM, "\n",
"Valor de TAM =",TAM, "\n",
"Valor p=", valor_p1, "\n")
cat("----------------------------------------------------------------","\n")
cat("Prueba de Mardia - Curtosis","\n",
"Curtosis Multivariado =", KM, "\n",
"Valor de TKM =",TKM2, "\n",
"Valor p=", valor_p2, "\n")
}
Ejemplo Muestras Grandes
head(mtcars)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
nrow(mtcars)
## [1] 32
NM2<-mvn(mtcars[,1:2])
NM2$multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 5.77902926070842 0.216269437201514 YES
## 2 Mardia Kurtosis -0.929901505898089 0.352422083240033 YES
## 3 MVN <NA> <NA> YES
mardia2(mtcars[,1:2])
## Test de Normalidad Multivariada
## Ho: Los datos siguen una Distribución Normal Multivariada
## ----------------------------------------------------------------
## Prueba de Mardia - Asimetría
## Asimetría Multivariada = 1.083568
## Valor de TAM = 5.779029
## Valor p= 0.2162694
## ----------------------------------------------------------------
## Prueba de Mardia - Curtosis
## Curtosis Multivariado = 6.684921
## Valor de TKM = 0.8647168
## Valor p= 0.3524221
mardia2.df<-function(datos){
n<-nrow(datos)
k<-ncol(datos)
x<-as.matrix(datos)
media<-colMeans(datos)
matriz.media<-matrix(media,nrow=n,ncol=k,byrow=TRUE)
S<-cov(datos)*(n-1)/n
SI<-solve(S) #Inversa de la matriz de covarianzas
distm<-((x-matriz.media)%*%SI)%*%t((x-matriz.media))
gl<-k*(k+1)*(k+2)/6
# Medida de Asimetría Multivariada
AM<-sum(distm^3)/n^2
if (n<30){
# Test Multivariado de Asimetría de Mardia - Muestras Pequeñas
c<-(n+1)*(n+3)*(k+1)/(n*((n+1)*(k+1)-6))
TAM<-n*c/6*AM
valor_p1<-pchisq(TAM,gl,lower.tail = FALSE)
} else {
# Test Multivariado de Asimetríıa de Mardia - Muestras grandes
TAM<-n/6*AM
valor_p1<-pchisq(TAM,gl,lower.tail = FALSE)
}
#Medida de Kurtosis Multivariado
KM<-sum((diag(distm))^2)/n
#Test Multivariado de Kurtosis de Mardia
TKM<-(KM-k*(k+2))*sqrt(n/(8*k*(k+2)))
TKM2<-TKM^2
valor_p2<-pchisq(TKM2,1,lower.tail = FALSE)
if (valor_p1<0.05){
test1 <- "No Normalidad"
} else {
test1 <- "Normalidad"
}
if (valor_p2<0.05){
test2 <- "No Normalidad"
} else {
test2 <- "Normalidad"
}
df<-data.frame(Prueba=c("Prueba de Asimetría","Prueba de Curtosis"),Medida=c(AM,KM),Estadístico=c(TAM,TKM2),valor_p=c(valor_p1,valor_p2),Conclusión=c(test1,test2))
df
}
mardia2.df(mtcars[,1:2])
## Prueba Medida Estadístico valor_p Conclusión
## 1 Prueba de Asimetría 1.083568 5.7790293 0.2162694 Normalidad
## 2 Prueba de Curtosis 6.684921 0.8647168 0.3524221 Normalidad