# TRAIN A LINEAR DISCRIMINANT CLASSIFIER. Takes the training inputs (as an # N by p matrix) and the training classes (vector of length N, containing # integers from 1 to K) as arguments. Returns the constant terms (a) and the # coefficients (b) for the K discriminant functions. lda.train <- function (train.x, train.y) { N <- nrow(train.x) p <- ncol(train.x) K <- max(train.y) # Find the estimates of mu and pi. mu <- matrix(NA,K,p) pi <- numeric(K) for (k in 1:K) { mu[k,] <- apply(train.x[train.y==k,,drop=F],2,mean) pi[k] <- mean(train.y==k) } # Find the estimate of the common Sigma from deviations from the estimated # means. Sigma <- matrix(0,p,p) for (k in 1:K) { d <- train.x[train.y==k,,drop=FALSE] for (j in 1:ncol(d)) { d[,j] <- d[,j] - mu[k,j] } Sigma <- Sigma + t(d) %*% d } Sigma <- Sigma / (N-K) # Compute the discriminant functions. Sigma.inv <- solve(Sigma) a <- numeric(K) b <- matrix(0,K,nrow(Sigma)) for (k in 1:K) { a[k] <- log(pi[k]) - 0.5 * t(mu[k,]) %*% Sigma.inv %*% mu[k,] b[k,] <- Sigma.inv %*% mu[k,] } # Return the coefficients of the discriminant functions. list (a=a, b=b) } # USE LINEAR DISCRIMINANT TO PREDICT CLASSES. Takes as arguments the estimates # from lda.train and the matrix of test inputs (M by p). Returns a vector # of predicted classes (of length M). lda.predict <- function (ld, test.x) { a <- ld$a b <- ld$b K <- length(a) pred <- numeric(nrow(test.x)) for (i in 1:nrow(test.x)) { x <- test.x[i,] max.disc <- -Inf for (k in 1:K) { d <- a[k] + t(x) %*% b[k,] if (d>max.disc) { max.disc <- d pred[i] <- k } } } pred }