set.seed(2^17-1)
library(MASS)
trn_X = rbind(mvrnorm(1000, c( 3,  3), matrix(c(1, 0, 0, 1), nrow=2)),
              mvrnorm(1000, c(-3, -3), matrix(c(1, 0, 0, 1), nrow=2)))
trn_y = factor(c(rep(1, 1000), rep(0, 1000)))
tst_X = matrix(c(1, 1, -1, -1), nrow = 2, byrow = T)

lda.model = lda(trn_X, trn_y)
lda.model
lda.predictions = predict(lda.model, tst_X)
lda.predictions

# implementation
n = nrow(trn_X)
p = ncol(trn_X)
prior = as.vector(table(trn_y)) / n
k = length(prior)
means = tapply(trn_X, list(rep(trn_y, p), col(trn_X)), mean)
centered = trn_X - means[trn_y,]
scaling = diag(1 / sqrt(diag(var(centered))), , p)
X = sqrt(1 / (n - k)) * centered %*% scaling
decomposition = svd(X, nu = 0)
rank = sum(decomposition$d > .0001)
scaling = scaling %*% decomposition$v[, 1:rank] %*% diag( 1 / decomposition$d[1:rank], , rank)
X = sqrt((n * prior) / (k - 1)) * scale(means, center = apply(prior %*% means, 2, sum), scale = F) %*% scaling
decomposition = svd(X, nu = 0)
rank = sum(decomposition$d > .0001 * decomposition$d[1])
scaling = scaling %*% decomposition$v[, 1:rank]
model = list(prior = prior, means = means, scaling = scaling)
model

# project new data [the MASS implementation of LDA uses a projection]
X = scale(tst_X, center = apply(prior * means, 2, sum), scale = F) %*% scaling
projected.means = scale(means, center = apply(prior * means, 2, sum), scale = F) %*% scaling
dist <- matrix(0.5 * apply(projected.means ^ 2, 1, sum) - log(prior), nrow(X), length(prior), byrow = TRUE) - X %*% t(projected.means)
dist <- exp( -(dist - apply(dist, 1L, min, na.rm=TRUE)))
posterior <- dist / drop(dist %*% rep(1, k))
posterior