# LDA by hand gives different results than MASS

I'm performing LDA by hand in R by following this formula:

$$\delta_k(x) = x^T\Sigma^-1\mu_k-0.5*\mu_{k}^T\Sigma^-1\mu_k+log(\pi_k)$$

where $$\pi_k$$ is the proportion of the data that's in the group $$k$$, $$\Sigma$$ is the covariance matrix that's assume to be constant through the different groups and $$\mu_k$$ is a vector with means of each predictor for class $$k$$.

I wrote the following code to implement this formula which classify each variable based on the maximum $$\delta_k(x)$$ obtain

library(ISLR)
data("Default")

Y = Default$default Xmat = Default[,-c(1:2)] Xmat = t((t(Xmat)-colMeans(Xmat))/apply(Xmat,2,sd)) Xmat = as.matrix(Xmat) ## 1. Se estima la media de cada variable por cada variable predictora Medias_hat = apply(Xmat,2, function(x) aggregate(x,by = list(Y),FUN = mean)[,2]) rownames(Medias_hat) = levels(Y) ## 2. Se estima la matriz de covarianza Sigma_hat = cov(Xmat) ## 3. Se estima pi pi_hat = table(Y)/length(Y) ## 4. Prediccion temp_class = sapply(1:ncol(Medias_hat), function(x){ Xmat%*%solve(Sigma_hat)%*%Medias_hat[x,]- as.numeric(1/2*t(Medias_hat[x,])%*%solve(Sigma_hat)%*%Medias_hat[x,])+log(pi_hat[x]) }) classificacion = apply(temp_class,1,which.max) classificacion = ifelse(classificacion ==1,"No","Yes")  The problem i'm having is that using the package MASS I get different results and I don't understand why. The code with this package is as follow: LdaMultiPred = lda(Y~as.matrix(Xmat)) predict(LdaMultiPred,as.data.frame(Xmat))$class


So I figure out that the covariance matrix to use is not the one estimated over the hole data set but a covariance matrix that it is weighted for each class. The formula, as in The Elements of Statistical Learning is:

$$\widehat{\Sigma} = \sum_{k=1}^K\sum_{g_i=k}(x_i-\widehat{\mu_k})(x_i-\widehat{\mu_k})^T/(N-K)$$

I used the following R code to estimate the matrix:

temp_cov = NULL
for(i in 1:(1+ncol(Xmat))){
temp = NULL
for(j in 1:length(levels(Y))){
pos_nivel = Y == levels(Y)[j]
if(i>ncol(Xmat)){
temp = c(temp,(t(Xmat[pos_nivel,1]-Medias_hat[j,1])%*%(Xmat[pos_nivel,2]-Medias_hat[j,2]))
/(nrow(Xmat)-length(levels(Y))))
} else {
temp = c(temp,(t(Xmat[pos_nivel,i]-Medias_hat[j,i])%*%(Xmat[pos_nivel,i]-Medias_hat[j,i]))
/(nrow(Xmat)-length(levels(Y))))
}
}
temp_cov = c(temp_cov,sum(temp))
}
Sigma_hat = diag(temp_cov[1:ncol(Xmat)])
Sigma_hat[lower.tri(Sigma_hat)] = temp_cov[-c(1:ncol(Xmat))]
Sigma_hat[upper.tri(Sigma_hat)] = temp_cov[-c(1:ncol(Xmat))]


With this change the code from above yields the same result as the lda function from MASS package