Discriminador lineal de Fisher

Supongamos que queremos discriminar los datos entre dos poblaciones normales con diferente vector de medias y la misma matriz de covarianzas. Vamos a ver como hacerlo en el caso bivariante ya que así podemos representar el problema gráficamente.
En vez de realizar el problema con datos simulados o reales, lo vamos a tratar con dos distribuciones teóricas concretas.
Así pues, consideremos la distribución de una población normal bivariante con media mu1=(2.5,4) y matriz de covarianzas Sigma=(2,1,1,2) y otra población con la misma matriz de covarianzas, pero de media mu2=(6,3). El gráfico de arriba representa las densidades bivariantes de estas dos poblaciones. Un segmento de puntos une los dos puntos medios de las poblaciones.

La idea de Fisher consistió en hallar una combinación lineal de las variables originales (X1,X2) de la forma w1X1+w2X2 y tal que discrimine "el máximo posible" las dos poblaciones. Él mismo definió el criterio de máxima discriminación como maximizar la razón entre la suma de cuadrados entre grupos y la suma de cuadrados dentro de los grupos sobre la combinación lineal (discriminant scores). La solución viene dada por el vector w (o es proporcional al vector w)

w = Sigma^{-1} * (mu2 - mu1)

Es decir, entre todas las combinaciones lineales de la forma w1X1+w2X2, la que mejor discrimina es justamente esa. En el ejemplo propuesto, el vector solución w se ha dibujado en el gráfico como una flecha desde el punto mu1. Observemos que w tiene en cuenta la diferencia entre las medias pero también la relación de dependencia entre las variables.

Para visualizar este resultado he dibujado las densidades de las dos poblaciones si consideramos el vector w=(1,0), es decir, que la combinación lineal considerada sea simplemente X1.

Esta no es la solución óptima.

La solución exacta es w=(8/3,-11/6)

En el gráfico inicial, la recta que discrimina las dos poblaciones tiene un vector director ortogonal al vector w. Esa recta es la recta de puntos que equidistan de los dos puntos medios según la distancia de Mahalanobis.

El código para hacer todos los cálculos y los gráficos es el siguiente:


mu1 <- c(2.5,4)
mu2 <- c(6,3)
Sigma <- matrix(c(2,1,1,2),ncol=2)

dnormbv <- function(x1,x2,mu,Sigma) {
sigma1 <- sqrt(Sigma[1,1])
sigma2 <- sqrt(Sigma[2,2])
rho <- Sigma[1,2]/(sigma1*sigma2)
cte <- 1/(2*pi*sigma1*sigma2*sqrt(1-rho^2))
cte * exp((-1/2) * (1/(1-rho^2)) * (
(x1-mu[1])^2/sigma1^2 - 2*rho*(x1-mu[1])*(x2-mu[2])/(sigma1*sigma2) + (x2-mu[2])^2/sigma2^2)
)}

x1<-seq(0,6,0.1)
x2<-seq(1,7,0.1)
z <- matrix(rep(0,length(x1)*length(x2)),ncol=length(x2))
for (i in 1:length(x1)) {
for (j in 1:length(x2)) {
z[i,j] <- dnormbv(x1[i],x2[j],mu1,Sigma)
}
}
plot(outer(x1,x2),xlim=c(0,9),ylim=c(0,9),xlab="",ylab="",type="n")
contour(x1,x2,z,nlev=4,add=T)

x1<-seq(2,9,0.1)
x2<-seq(0,7,0.1)
z <- matrix(rep(0,length(x1)*length(x2)),ncol=length(x2))
for (i in 1:length(x1)) {
for (j in 1:length(x2)) {
z[i,j] <- dnormbv(x1[i],x2[j],mu2,Sigma)
}
}
contour(x1,x2,z,nlev=4,add=T)

segments(mu1[1],mu1[2],mu2[1],mu2[2],lty="dotted")

SigmaInv <- matrix(c(2/3,-1/3,-1/3,2/3),ncol=2)

w <- SigmaInv %*% (mu2-mu1)

arrows(mu1[1],mu1[2],mu1[1]+w[1],mu1[2]+w[2],lwd=2)

pm <- (mu1+mu2)/2
a <- w[1] * pm[1] / w[2] + pm[2]
b <- -w[1]/w[2]
abline(a,b,lty="dashed")


ee <- 3.5
x1 <- seq(mu1[1]-ee,mu1[1]+ee,by=0.1)
x2 <- dnorm(x1,mean=mu1[1],sd=sqrt(Sigma[1,1]))
plot.new()
plot.window(xlim=c(-1,9),ylim=c(0,0.5))
axis(1)
axis(2)
lines(x1,x2)
x1 <- seq(mu2[1]-ee,mu2[1]+ee,by=0.1)
x2 <- dnorm(x1,mean=mu2[1],sd=sqrt(Sigma[1,1]))
lines(x1,x2)



E1 <- t(w) %*% mu1
E2 <- t(w) %*% mu2
var.w <- as.numeric(t(w) %*% Sigma %*% w)

ee <- 7
x1 <- seq(E1-ee,E1+ee,by=0.1)
x2 <- dnorm(x1,mean=E1,sd=sqrt(var.w))
plot.new()
plot.window(xlim=c(-6,16),ylim=c(0,0.5))
axis(1)
axis(2)
lines(x1,x2)
x1 <- seq(E2-ee,E2+ee,by=0.1)
x2 <- dnorm(x1,mean=E2,sd=sqrt(var.w))
lines(x1,x2)



Con datos reales o simulados podemos utilizar la función lda() del paquete MASS.

Bibliografía

Lattin, J. et al, Analyzing Multivariate Data, Ed. Brooks/Cole, Belmont (2003).