Análisis Multivariado
Pilar
2023-01-20
X=data.frame(V1=c (1,1,2,3 ),V2=c (- 1,0,3,0), V3=c(3,3,0,1))
X## V1 V2 V3
## 1 1 -1 3
## 2 1 0 3
## 3 2 3 0
## 4 3 0 1sapply(X, mean)## V1 V2 V3
## 1.75 0.50 1.75#Varianza muestral
sapply(X, var )## V1 V2 V3
## 0.9166667 3.0000000 2.2500000 #Desviación t í p i c a muestral
sapply (X, sd )## V1 V2 V3
## 0.9574271 1.7320508 1.5000000 #Varianza
var_ver=function(x){var(x)*(length(x)-1)/length(x)}
sapply(X,var_ver)## V1 V2 V3
## 0.6875 2.2500 1.6875#Desviación t í p i c a
sd_ver=function(x){sqrt(var_ver(x))}
sapply(X,sd_ver )## V1 V2 V3
## 0.8291562 1.5000000 1.2990381Código tipificación tabla de datos
scale(X)## V1 V2 V3
## [1,] -0.7833495 -0.8660254 0.8333333
## [2,] -0.7833495 -0.2886751 0.8333333
## [3,] 0.2611165 1.4433757 -1.1666667
## [4,] 1.3055824 -0.2886751 -0.5000000
## attr(,"scaled:center")
## V1 V2 V3
## 1.75 0.50 1.75
## attr(,"scaled:scale")
## V1 V2 V3
## 0.9574271 1.7320508 1.5000000apply(scale(X),2,mean)## V1 V2 V3
## -4.163336e-17 0.000000e+00 0.000000e+00apply(scale(X),2,var)## V1 V2 V3
## 1 1 1scale(X, center=TRUE, scale = FALSE)## V1 V2 V3
## [1,] -0.75 -1.5 1.25
## [2,] -0.75 -0.5 1.25
## [3,] 0.25 2.5 -1.75
## [4,] 1.25 -0.5 -0.75
## attr(,"scaled:center")
## V1 V2 V3
## 1.75 0.50 1.75X=data.frame(V1=c (1,1,2,3 ),V2=c (- 1,0,3,0), V3=c(3,3,0,1))
X## V1 V2 V3
## 1 1 -1 3
## 2 1 0 3
## 3 2 3 0
## 4 3 0 1#Covarianza entre l a s dos primeras columnas 9 #Muestral
cov(X$V1,X$V2)## [1] 0.5cov(X)## V1 V2 V3
## V1 0.9166667 0.500000 -1.083333
## V2 0.5000000 3.000000 -2.166667
## V3 -1.0833333 -2.166667 2.250000n=dim(X)[1]
A<-((n-1/n))*cov(X)
A## V1 V2 V3
## V1 3.4375 1.875 -4.0625
## V2 1.8750 11.250 -8.1250
## V3 -4.0625 -8.125 8.4375#Cálculo c o r r e l a c i o n e s
#Matriz c o r r e l a c i o n e s
cor (X)## V1 V2 V3
## V1 1.0000000 0.3015113 -0.7543365
## V2 0.3015113 1.0000000 -0.8339504
## V3 -0.7543365 -0.8339504 1.0000000#Código dos variables en función plot
iris.pet=iris[,c("Petal.Length","Petal.Width")]
plot(iris.pet, pch=20,xlab="Largo", ylab="Ancho")
Diagramas de dispersión con plot
plot(iris,1:4)
plot(iris,1:4, col=iris$Species)
###############################
Distribucion Normal Multivariante
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.
library(MASS)
library(MVN)## Warning: package 'MVN' was built under R version 4.2.2library(matrixcalc)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 = 500 000 y mostramos el histograma 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.609249hist(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.
Genere una muestra aleatoria de tamaño 50
Generamos la muestra de tamaño 50 y mostramos la 6 primeras fila
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.244346hist(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.
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] 59587001Distancia 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.94126100Test de Multinomalidad: Método Gráfico
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.
Test de Multinomalidad: Método 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 Anderson-Darling Column1 0.1555 0.9523 YES
## 2 Anderson-Darling Column2 0.5089 0.1895 YES
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th
## 1 50 2260.339743 1641.433021 2142.05963 -1223.733866 5841.97542 1105.026297
## 2 50 3.580762 5.286258 3.17326 -7.945419 16.18597 1.249477
## 75th Skew Kurtosis
## 1 3532.06353 0.03492158 -0.6403271
## 2 7.35956 0.08734077 -0.1038849Es 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 norma ##########################################
La distribución multinomial en R
Hay tres candidatos para alcalde; el candidato A recibe el 10% de los votos, el candidato B recibe el 40% y el candidato C recibe el 50%.
¿Cuál es la probabilidad de que 3 personas elijan al candidato A, 4 elijan al candidato B y 5 elijan al candidato C de una muestra aleatoria de 10 votantes?
dmultinom(x=c(3,4,5), prob = c(0.1,0.4,0.5))## [1] 0.022176La probabilidad de que precisamente 3 personas votaran por A, 4 por B y 5 por C es 0.022176.
Ejercicio 2
Suponga que una urna contiene cinco canicas amarillas, tres canicas rojas y dos canicas rosas.
¿Cuál es la probabilidad de que las cuatro bolas sean amarillas si elegimos al azar cuatro bolas de la urna con reemplazo?
dmultinom(x=c(4,0,0), prob = c(0.5,0.3,0.2))## [1] 0.0625Es 0,0625 veces más probable que las cuatro bolas sean amarillas.
Digamos que dos alumnos compiten en un juego de ajedrez. Las probabilidades de que el estudiante A gane un juego en particular son 0,6, las probabilidades de que el estudiante B gane un juego en particular son 0,2 y las probabilidades de empate en un juego en particular son 0,2.
¿Cuál es la probabilidad de que el jugador A gane cuatro veces, el jugador B gane cuatro veces y empaten dos veces si juegan diez juegos?
dmultinom(x=c(4,4,2), prob = c(0.6,0.2,0.2))## [1] 0.02612736Alrededor del 0,02612736 por ciento de las veces, el jugador A gana cuatro veces, el jugador B gana cuatro veces y empatan dos veces. ###################################
Dos jugadores de ajedrez tienen la probabilidad de que el jugador A gane es 0,40, el jugador B gane es 0,35, el juego termine en empate es 0,25.
La distribución multinomial se puede usar para responder preguntas como: “Si estos dos jugadores de ajedrez jugaron 12 juegos, ¿cuál es la probabilidad de que el jugador A gane 7 juegos, el jugador B gane 2 juegos y los 3 juegos restantes sean empate?”
dmultinom(x=c(7,2,3), prob = c(0.4,0.35,0.25))## [1] 0.02483712Encuestas de opinión sobre las elecciones
En un pueblo pequeño, el 40% de los votantes elegibles prefieren al candidato A, el 10% prefieren al candidato B, el 50% no tienen preferencia.
Muestra aleatoriamente 10 votantes elegibles. ¿Cuál es la probabilidad de que 4 prefieran al candidato A, 1 prefiera al candidato B y 5 no tengan preferencia?
dmultinom(x=c(4,1,5), prob = c(0.4,0.10,0.50))## [1] 0.1008Tirar un dado N=100 veces
one.dice <- function(){
dice <- sample(1:6, size = 1, replace = TRUE)
return(dice)
}
one.dice()## [1] 4par(mfrow=c(2,2))
for (i in 1:4){
sims <- replicate(100, one.dice())
table(sims)
table(sims)/length(sims)
plot(table(sims), xlab = 'Event', ylab = 'Frequency')
}
Tirar un dado N=10000 veces
par(mfrow=c(2,2))
for (i in 1:4){
sims <- replicate(10000, one.dice())
table(sims)
table(sims)/length(sims)
plot(table(sims), xlab = 'Event', ylab = 'Frequency')
}
# rmultinom(n, size, prob)
my_prob <- c(0.2,0.3,0.1,0.4)
number_of_experiments <- 10
number_of_samples <- 10 #numero de muestras
experiments <- rmultinom(n=number_of_experiments, size=number_of_samples, prob=my_prob)
experiments## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 3 3 1 1 3 3 1 0 2
## [2,] 5 3 1 3 3 2 1 5 3 3
## [3,] 2 2 1 2 1 0 0 1 3 1
## [4,] 2 2 5 4 5 5 6 3 4 4df=data.frame(experiments)/number_of_samples
par(mfrow=c(2,5))
for(i in 1:10) {
barplot(df[,i],ylim=c(0,1))
}
número_de_muestras <- 10000
experiments <- rmultinom(n=number_of_experiments, size=number_of_samples, prob=my_prob)
experiments## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 2 2 0 3 1 1 3 4 1 1
## [2,] 3 3 0 0 0 3 0 2 0 2
## [3,] 3 1 4 1 0 1 4 1 0 1
## [4,] 2 4 6 6 9 5 3 3 9 6df=data.frame(experiments)/number_of_samples
par(mfrow=c(2,5))
for(i in 1:10) {
barplot(df[,i],ylim=c(0,1))
}