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


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){
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))

1 Answer 1


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]
        temp = c(temp,(t(Xmat[pos_nivel,1]-Medias_hat[j,1])%*%(Xmat[pos_nivel,2]-Medias_hat[j,2]))
      } else {
        temp = c(temp,(t(Xmat[pos_nivel,i]-Medias_hat[j,i])%*%(Xmat[pos_nivel,i]-Medias_hat[j,i]))
    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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.