Guía4-i,“Coordenadas Discriminantes”Octavio M. Duarte

library(tidyverse)library(knitr)set.seed(42)

"Af,1.38,1.64Af,1.40,1.70Af,1.24,1.72Af,1.36,1.74Af,1.38,1.82Af,1.48,1.82Af,1.54,1.82Af,1.38,1.90Af,1.56,2.08Apf,1.14,1.78Apf,1.20,1.86Apf,1.18,1.96Apf,1.30,1.96Apf,1.26,2.00Apf,1.28,2.00" %>% read_csv( col_names = c(

"especie","x1","x2"

)) -> md4

md4 <- md4 %>% mutate( y1 = x1, y2 = x1 + x2 )

Enunciado

especie x1 x2 y1 y2Af 1.38 1.64 1.38 3.02Af 1.40 1.70 1.40 3.10Af 1.24 1.72 1.24 2.96Af 1.36 1.74 1.36 3.10Af 1.38 1.82 1.38 3.20Af 1.48 1.82 1.48 3.30Af 1.54 1.82 1.54 3.36Af 1.38 1.90 1.38 3.28Af 1.56 2.08 1.56 3.64Apf 1.14 1.78 1.14 2.92Apf 1.20 1.86 1.20 3.06Apf 1.18 1.96 1.18 3.14Apf 1.30 1.96 1.30 3.26Apf 1.26 2.00 1.26 3.26Apf 1.28 2.00 1.28 3.28

1

Prueba para la Igualdad de Matrices de Covarianza

Estamos en el caso H1 de la prueba sobre Σ para dos poblaciones.poblacion <- "especie"variables <- c('y1','y2')

obtenerQi <- function(marcoDeDatos,listVariables,poblacion,colPoblacion) {indices <- which( marcoDeDatos[,quo_name(colPoblacion)] == poblacion )cov(

marcoDeDatos[indices,variables]) * length(indices)

}

obtenerQ <- function(marcoDeDatos,listVariables) {cov(

marcoDeDatos[,variables]) * nrow(marcoDeDatos)

}

todosLosQ <- function(marcoDeDatos,listVariables,colPoblacion) {listaPoblaciones <- marcoDeDatos[,colPoblacion] %>% unique() %>% unlist()listaQ <- map( listaPoblaciones, ~obtenerQi(marcoDeDatos, listVariables,.x, colPoblacion ) )names(listaQ) <- listaPoblacionesreturn(

listaQ)

}

ques <- todosLosQ(md4,variables,"especie")ques

## $Af## y1 y2## y1 0.08820 0.16095## y2 0.16095 0.38560#### $Apf## y1 y2## y1 0.02368 0.04976## y2 0.04976 0.12256mediaUnaPoblacion <- function(marcoDeDatos,listVariables,colPoblacion,poblacion) {

indices <- which( md4[,colPoblacion] == poblacion )return(

colMeans( md4[indices,variables] ))

}

mediaUnaPoblacion(md4,variables,"especie","Af")

## y1 y2## 1.413333 3.217778

2

mediaGeneral <- function(marcoDeDatos,listVariables) {return( colMeans( md4[,variables] ) )

}

obtenerNi <- function(marcoDeDatos,colPoblacion) {listaPoblaciones <- marcoDeDatos[,colPoblacion] %>% unique() %>% unlist()listaN <- map_dbl( listaPoblaciones,

~length( which( marcoDeDatos[,colPoblacion] == .x ) ))

names(listaN) <- listaPoblacionesreturn( listaN )

}

nies <- obtenerNi(md4,"especie")nies

## Af Apf## 9 6obtenerH <- function(marcoDeDatos,listVariables,colPoblacion) {

listaPoblaciones <- marcoDeDatos[,colPoblacion] %>% unique() %>% unlist()listaMedias <- map( listaPoblaciones, ~mediaUnaPoblacion(marcoDeDatos,listVariables,colPoblacion,.x) )names(listaMedias) <- listaPoblacionesmediaGeneral <- mediaGeneral(marcoDeDatos,variables)listaN <- obtenerNi(marcoDeDatos,colPoblacion)listaASumar <- map2( listaMedias, listaN, ~(( .x - mediaGeneral ) %*% t( .x - mediaGeneral )) * .y )matrizH <- listaASumar %>% reduce(`+`)return(matrizH)

}

matrizHP <- obtenerH( md4, variables, "especie" )matrizHP

## y1 y2## [1,] 0.12544000 0.04330667## [2,] 0.04330667 0.01495111gammaEstrella1 <- function(marcoDeDatos,listVariables,colPoblacion) {

listaPoblaciones <- marcoDeDatos[,colPoblacion] %>% unique() %>% unlist()listaQ <- todosLosQ(marcoDeDatos,listVariables,colPoblacion)listaN <- obtenerNi(marcoDeDatos,colPoblacion)

numerador <- map2_dbl(listaQ,listaN,~(det(.x/.y))^(.y/2)

) %>% reduce(`*`)

nTotal <- nrow(marcoDeDatos)denominador <- ( det( reduce( listaQ, `+` ) / nTotal ) )^(nTotal/2)return( numerador / denominador )

}

estadistico <- gammaEstrella1( md4,variables,"especie" )

3

indicesNu <- function(listVariables,colPoblacion) {n <- length(colPoblacion)p <- length(listVariables)k <- colPoblacion %>% unique() %>% length()nu <- p*k+(p*(p+1)*k)/2nu1 <- p*k+(p*(p+1))/2nu2 <- p+(p*(p+1))/2return(

list("nu" = nu,"nu1" = nu1,"nu2" = nu2

))

}

indices <- indicesNu(variables,md4$especie)indices

## $nu## [1] 10#### $nu1## [1] 7#### $nu2## [1] 5pruebaSigmaEnPoblaciones <- function(marcoDeDatos,listVariables,colPoblacion,alfa) {

estadistico <- gammaEstrella1(marcoDeDatos,listVariables,colPoblacion)indices <- indicesNu(listVariables,marcoDeDatos[[colPoblacion]])percentil <- qchisq( p = 1 - alfa, df = indices[["nu"]] - indices[["nu1"]] )gradLib1 <- indices[["nu"]] - indices[["nu1"]]return(

tibble(rechazo = -2*log(estadistico) > percentil,pValor = pchisq( df = gradLib1, q = -2 * log(estadistico) ),estadistico = estadistico,percentil = qchisq( df = gradLib1, p = 1- alfa )

))

}

Evalue si las poblaciones tienen igual matriz de covarianza. Tome α = 0, 01. Enbase al resultado decida si es razonable hacer un gráfico de la primera coordenadadiscriminante.Gráfico preliminar y pueba

md4 %>% ggplot() +aes(x=y1, y=y2,color=especie) +geom_point()

4

2.9

3.1

3.3

3.5

1.2 1.3 1.4 1.5

y1

y2

especie

Af

Apf

Prueba de Hipótesis

pruebaEj4a <- pruebaSigmaEnPoblaciones(md4,variables,"especie",0.01)pruebaEj4a %>% kable()

rechazo pValor estadistico percentilFALSE 0.7304125 0.1404498 11.34487

La prueba donde H0 : Σ1 = Σ2 y H1 : Σ1 6= Σ2 nos da un pvalor de 0.7304125 y por lo tanto indica que esrazonable asumir que las matrices de covarianza son iguales para ambas poblaciones. Esta es la hipótesis quenos habilita a buscar coordenadas discriminantes.

• Dado que no rechazamos la hipótesis, es razonable hacerlo.

Obtención de las coordenadas disciminantes

estimarSigmaDentro <- function(marcoDeDatos,listaVariables,colPoblacion) {listaQ <- todosLosQ(marcoDeDatos,listVariables,colPoblacion)n <- nrow(marcoDeDatos)sDentro <- reduce(listaQ,`+`) /nreturn(sDentro)

}

sDentroP <- estimarSigmaDentro( md4, variables, "especie" )sDentroP

5

## y1 y2## y1 0.007458667 0.01404733## y2 0.014047333 0.03387733estimarSigmaEntre <- function(marcoDeDatos,listaVariables,colPoblacion) {

n <- nrow(marcoDeDatos)matrizH <- obtenerH(marcoDeDatos,listaVariables,colPoblacion) / nreturn( matrizH )

}

coordenadasDiscriminantes <- function(marcoDeDatos, listaVariables, colPoblacion) {sEntre <- estimarSigmaEntre(marcoDeDatos,listaVariables,colPoblacion)sDentro<- estimarSigmaDentro(marcoDeDatos,listaVariables,colPoblacion)matrizC <- chol(sDentro)inversaC <- solve(matrizC)matrizB <- t(inversaC) %*% sEntre %*% inversaCdescompB <- eigen(matrizB)## indicesNoNulas <- which( descompB$values != 0 )## listaVec <- map( indicesNoNulas, ~ inversaC %*% (descompB$vectors)[,.x] )listaVectores <- map( seq(1,length(descompB$values)), ~descompB$vectors[,.x] )listaDirecciones <- map( listaVectores, ~ inversaC %*% .x )columnasProyecciones <- map(listaDirecciones, ~ (marcoDeDatos[,listaVariables] %>% as.matrix()) %*% .x )return(

list(proyecciones = columnasProyecciones,direcciones = listaDirecciones)

)}

coordenadasDiscriP <- coordenadasDiscriminantes(md4,variables,"especie")coordenadasDiscriP

## $proyecciones## $proyecciones[[1]]## [,1]## [1,] -4.5943182## [2,] -4.3080252## [3,] -1.7286619## [4,] -3.3223048## [5,] -2.8412235## [6,] -4.3315828## [7,] -5.2257984## [8,] -2.0620703## [9,] -2.9916224## [10,] 0.3460623## [11,] 0.2309999## [12,] 1.5030133## [13,] -0.2854179## [14,] 0.7003024## [15,] 0.4022305#### $proyecciones[[2]]## [,1]## [1,] -16.06545## [2,] -16.52713

6

## [3,] -15.99176## [4,] -16.61435## [5,] -17.20235## [6,] -17.61590## [7,] -17.86404## [8,] -17.70764## [9,] -19.58893## [10,] -15.95718## [11,] -16.71060## [12,] -17.25950## [13,] -17.75576## [14,] -17.84298## [15,] -17.92569###### $direcciones## $direcciones[[1]]## [,1]## y1 -24.643008## y2 9.739415#### $direcciones[[2]]## [,1]## y1 2.180558## y2 -6.316100proyeccionesP <- coordenadasDiscriP$proyeccionesdireccionesP <- coordenadasDiscriP$direcciones

sEntreP <- estimarSigmaEntre(md4,variables,"especie")sEntreP

## y1 y2## [1,] 0.008362667 0.0028871111## [2,] 0.002887111 0.0009967407matrizC <- chol(sDentroP)matrizC

## y1 y2## y1 0.08636357 0.16265345## y2 0.00000000 0.08614631inversaC <- solve(matrizC)

matrizB <- t(inversaC) %*% sEntreP %*% inversaCmatrizB

## y1 y2## y1 1.121201 -1.728890## y2 -1.728890 2.665946descompB <- eigen(matrizB)dir1 <- ( inversaC %*% (descompB$vectors)[,1] )dir1

## [,1]

7

## y1 -24.643008## y2 9.739415dir2 <- ( inversaC %*% (descompB$vectors)[,2] )dir2

## [,1]## y1 2.180558## y2 -6.316100normaVector <- function(vector) {

norm( x = matrix( vector, ncol=1 ),type="2" )}

colProy <- (md4[,variables] %>% as.matrix()) %*% dir1

md42 <- md4md42$proyecciones1 <- proyeccionesP[[1]] %>% unlist()md42$proyecciones2 <- proyeccionesP[[2]] %>% unlist()

md42 %>% kable()

especie x1 x2 y1 y2 proyecciones1 proyecciones2Af 1.38 1.64 1.38 3.02 -4.5943182 -16.06545Af 1.40 1.70 1.40 3.10 -4.3080252 -16.52713Af 1.24 1.72 1.24 2.96 -1.7286619 -15.99176Af 1.36 1.74 1.36 3.10 -3.3223048 -16.61435Af 1.38 1.82 1.38 3.20 -2.8412235 -17.20235Af 1.48 1.82 1.48 3.30 -4.3315828 -17.61590Af 1.54 1.82 1.54 3.36 -5.2257984 -17.86404Af 1.38 1.90 1.38 3.28 -2.0620703 -17.70764Af 1.56 2.08 1.56 3.64 -2.9916224 -19.58893Apf 1.14 1.78 1.14 2.92 0.3460623 -15.95718Apf 1.20 1.86 1.20 3.06 0.2309999 -16.71060Apf 1.18 1.96 1.18 3.14 1.5030133 -17.25950Apf 1.30 1.96 1.30 3.26 -0.2854179 -17.75576Apf 1.26 2.00 1.26 3.26 0.7003024 -17.84298Apf 1.28 2.00 1.28 3.28 0.4022305 -17.92569

Haga un gráfico de las dos primeras coordenadas discriminantes,¿Qué observaen la segunda coordenada?Sabemos que la segunda coordenada corresponde al autovalor 0 y por lo tanto no aporta información.

Gráficamente esto se observa considerando que en este sistema de coordenadas los centroides de las respectivaspoblaciones están alineados para esta segunda coordenada.varsProy <- c('proyecciones1','proyecciones2')

centroidePAf <- md42[ which( md42$especie == 'Af' ) ,varsProy] %>% colMeans %>% enframe() %>% spread(name,value)centroidePApf <- md42[ which( md42$especie == 'Apf' ) ,varsProy] %>% colMeans %>% enframe() %>% spread(name,value)

centroides <- tibble(proyecciones1 = c(centroidePAf$proyecciones1, centroidePApf$proyecciones1),proyecciones2 = c(centroidePAf$proyecciones2, centroidePApf$proyecciones2),

8

especie = c('Af','Apf'))

md42 %>% ggplot() +aes( x = proyecciones1, y = proyecciones2 , color = especie ) +geom_point() +geom_point(data = centroides, color = "blue", size = 5) +geom_hline( yintercept = centroides[[1,'proyecciones2']], linetype="dashed")

−19

−18

−17

−16

−4 −2 0

proyecciones1

proy

ecci

ones

2

especie

Af

Apf

Esto hace bien evidente la propiedad fundamental de las coordenadas discriminantes, son las que evidencianuna mejor separación entre los datos.

Proyección en Diversas Direcciones

A modo comparativo, realizamos una plasmación unidimensional sobre las dos coordenadas naturales, lacoordenada discriminante y su coordenada ortogonal.md42 %>% gather(y1,y2,proyecciones1,proyecciones2,key='variable',value='magnitud') %>%

ggplot() +aes( x = magnitud, y = 0, color = especie ) +geom_point() +facet_grid(.~variable, scales="free")

## Warning: attributes are not identical across measure variables;## they will be dropped

9

proyecciones1 proyecciones2 y1 y2

−4 −2 0 −19 −18 −17 −16 1.2 1.3 1.4 1.5 2.9 3.1 3.3 3.5−0.050

−0.025

0.000

0.025

0.050

magnitud

0

especie

Af

Apf

Se puede observar:

• Proyectando sobre la coordenada discriminante las especies forman dos cúmulos separados.• Proyectando sobre la coodenada ortogonal a la discriminante, tenemos mínima distinción entre los

grupos.• Las otras dos coordenadas muestran grados de separación intermedios.

Realice un gráfico de los puntos en las coordenadas originales y grafique la rectaa ·

{x− x1−x2

2

}= 0.

Hallar los puntos de la curva implícita.

vector1 <- (dir1 / normaVector(dir1) ) %>% unlist() %>% c()vector2 <- c(1,1)proyeccion2en1 <- (vector1 / (normaVector(vector1))^2 ) * c( vector2 %*% vector1 )proyeccion2en1

## [1] 0.5230759 -0.2067302vectorResto <- vector2 - proyeccion2en1vectorParam <- vectorResto / vectorResto[[1]]

listaMediasP <- map( md4$especie %>% unique() %>% unlist() , ~mediaUnaPoblacion(md4,variables,"especie",.x) )mediaMedias <- ( reduce(listaMediasP,`+`) / 2 )

trazarRecta <- function(t) {punto <- c(t) * vectorParam + mediaMedias

10

names(punto) <- c('y1','y2')return( punto )

}

puntosL <- seq(-0.5,0.5,length.out = 15)puntosRecta <- map(puntosL,trazarRecta)colP1 <- puntosRecta %>% map_dbl(~.x['y1'])colP2 <- puntosRecta %>% map_dbl(~.x['y2'])

coordenadasRecta <- tibble(y1 = colP1,y2 = colP2

)

md43 <- md42md43$P1 <- colP1md43$P2 <- colP2

Recta Discriminante

ggplot(data=md43) +aes( x = y1, y = y2, color = especie ) +geom_point() +geom_line( aes( x = P1, y = P2), colour = "purple" )

2.0

2.5

3.0

3.5

4.0

4.5

1.00 1.25 1.50 1.75

y1

y2

especie

Af

Apf

11

Comprobar que la distancia euclídea entre los promedios de cada grupo ex-presados en la primer coordenada coincide con la distancia de Mahalanobisentre los promedios x1 y x2 expresados en las variables originales, es decir√

(x1 − x2)′ S−1 (x1 − x2).

difMedias <- listaMediasP %>% reduce(`-`)listaQP <- todosLosQ(md4,variables,"especie")

matrizS <- (1 / (nrow(md4) - length( unique(md4$especie) ) ) ) * reduce(listaQP, `+`)

(t(difMedias) %*% solve(sDentroP) %*% difMedias) %>% sqrt()

## [,1]## [1,] 3.972377mediaProyAf <- md43[which( md43$especie == "Af" ),"proyecciones1"] %>% colMeans()mediaProyApf <- md43[which( md43$especie == "Apf" ),"proyecciones1"] %>% colMeans()

difProms <- mediaProyApf - mediaProyAf

difProms

## proyecciones1## 3.972377

12

1. Estadístico Usado

γ∗1 =

∏kj=1

∣∣∣Qi

ni

∣∣∣ni2

∣∣Un

∣∣n2

−2 ln (γ∗1 )D→ χ2

ν−ν1

2. Grafique los puntos originales y la recta a′[x− x1+x2

2

]= 0.

R : a′[x− x1 + x2

2

]= 0

a = S−1 · (x1 + x2)

Dado que estamos en el plano, el vector ortogonal a a′ es único y define larecta que queremos en forma explícita.Necesitamos además un punto que pertenezca a al misma, el más simple deobtener es x1+x2

2 .b/b · a = 0,

Sea u ∈ R2/u ∦ ab = a− Proya (u)(0; x1+x2

2

)∈ R

(x1;x2) = b · t+ x1 + x2

2

(x1;x2) = b · t+ o

R (t) =

{x1 (t) = b1t+ o1

x2 (t) = b2t+ o2

Podemos obtener una ecuación explícita:{x1−o1b1

= tx2−o2b2

= t

x2 − o2b2

=x1 − o1b1

x2 − o2 =b2b1x1 −

b2b1o1

R : x2 =b2b1x1 −

b2b1o1 + o2

1