### criterio del discriminante lineal
### (equivale a menor distancia de Majalanobis si hay homocedasticidad)
discriminante.lineal <- function (X, tipo, nuevo)
  ## clasifica a "nuevo" en alguno de los niveles de "tipo"
  {
    p    <- ncol(X)                 # número de variables de clasificación
    Ni   <- table (tipo)            # tamaños de muestra
    library (Matrix)                # para diagonal por bloques con "bdiag"
    ## matriz que calcula las medias por grupos:
    D    <- as.matrix (bdiag (lapply (Ni, function (ni) matrix (1/ni, ni, ni))))
    n    <- sum (Ni)                # tamaño total de la muestra
    Jn   <- matrix(1,n,n)/n         # matriz que calcula media global de las variables
    I    <- diag(rep(1,n))          # matriz identidad
    H    <- I - Jn                  # matriz de centrado
    T    <- t(X) %*% H %*% X        # matriz de variabilidad total
    B    <- t(X) %*% (D - Jn) %*% X # matriz de variabilidad entre grupos
    W    <- t(X) %*% (I - D)  %*% X # matriz de variabilidad intra grupos
    ## lista de matrices para luego poner en diagonal por bloques:
    bloques <- lapply (Ni, function (ni) matrix(1/ni,1,ni))
    ## matriz que calcula las medias por grupos:
    Dbis    <- as.matrix (bdiag (bloques))
    dimnames (Dbis) [[1]] <- names (bloques)
    Med  <- Dbis %*% X              # matriz de medias por grupos
    S    <- W / n                   # matriz de varianzas y covarianzas
    A    <- 2 * solve(S) %*% t(Med)
    C    <- diag (Med %*% solve(S) %*% t(Med))
    V    <- nuevo %*% A - C
    asigna.indice <- which.max(V)
    asigna.nombre <- dimnames (V) [[2]] [asigna.indice]
    asigna.nombre
  }

### factor discriminante de Fisher
discriminante.factor <- function (X, tipo, nuevo)
  {
    p      <- ncol(X)                 # número de variables de clasificación
    Ni     <- table (tipo)            # tamaños de muestra
    library (Matrix)                # para diagonal por bloques con "bdiag"
    ## matriz que calcula las medias por grupos:
    D      <- as.matrix (bdiag (lapply (Ni, function (ni) matrix (1/ni, ni, ni))))
    n      <- sum (Ni)                # tamaño total de la muestra
    Jn     <- matrix(1,n,n)/n         # matriz que calcula media global de las variables
    I      <- diag(rep(1,n))          # matriz identidad
    H      <- I - Jn                  # matriz de centrado
    T      <- t(X) %*% H %*% X        # matriz de variabilidad total
    B      <- t(X) %*% (D - Jn) %*% X # matriz de variabilidad entre grupos
    W      <- t(X) %*% (I - D)  %*% X # matriz de variabilidad intra grupos
    autoW  <- eigen (W)
    invW12 <- autoW$vectors %*% diag (1 / sqrt (autoW$values)) %*% t(autoW$vectors)
    auto   <- eigen (invW12 %*% B %*% invW12)
    V      <- auto$vectors
    landa  <- auto$values
    varexp <- landa / sum(landa)
    U      <- invW12 %*% V              # vectores solución
    F      <- X %*% U                   # puntuaciones factoriales
    bloques<- lapply (Ni, function (ni) matrix(1/ni,1,ni))
    Dbis   <- as.matrix (bdiag (bloques))
    dimnames (Dbis) [[1]] <- names (bloques)
    MedF   <- Dbis %*% F                # medias en los factores
    ## calcular distancias a las medias:
    Fnuevo <- nuevo %*% U
    Dist <- sapply (levels (tipo), function (t) sum ((Fnuevo-MedF[t,])^2))
    asigna.indice.F <- which.min (Dist)
    asigna.nombre.F <- names (asigna.indice.F)
    asigna.nombre.F
  }

### datos
library (datasets)              # para hacer que "iris" sea accesible desde Rcmdr
## equivale lo anterior a Rcmdr > Datos > ...en paquetes > Leer... > datasets > iris
X    <- as.matrix(iris[,-5])    # variables de clasificación
tipo <- iris[,5]                # variable de grupo

### validación cruzada "sacando uno"

mean (sapply (1:nrow(X),
              function (i)
              {
                Xi   <- X[-i,]
                tipoi <- tipo[-i]
                tipo[i] == discriminante.lineal (Xi, tipoi, X[i,])
              }))

mean (sapply (1:nrow(X),
              function (i)
              {
                Xi   <- X[-i,]
                tipoi <- tipo[-i]
                tipo[i] == discriminante.factor (Xi, tipoi, X[i,])
              }))