y_hat_qda <- predict(train_qda, mnist_27$true_p, type = "prob")[,2]

grid.arrange(plot_cond_prob(),
             plot_cond_prob(y_hat_qda, title = "QDA estimate"),
             ncol = 2)

Linear Discriminant Analysis (LDA)

The LDA is similar to the QDA, however we assume that the correlation structure is the same for all classes, therefore we utilize the mean valor for the standard deviations and correlations between the outcomes.

params <- mnist_27$train %>%
  group_by(y) %>%
  summarize(avg_1 = mean(x_1), avg_2 = mean(x_2),
            sd_1 = sd(x_1), sd_2 = sd(x_2),
            r = cor(x_1, x_2))
params <- params %>% mutate(sd_1 = mean(sd_1), sd_2 = mean(sd_2), r = mean(r))
params

A tibble: 2 × 6

y	avg_1	avg_2	sd_1	sd_2	r
<fct>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
2	0.1287798	0.2831542	0.07102743	0.08129729	0.4279554
7	0.2341133	0.2881862	0.07102743	0.08129729	0.4279554

train_lda <- train(y ~., method = "lda", data = mnist_27$train)
train_lda$finalModel

Call:
lda(x, grouping = y)

Prior probabilities of groups:
      2       7 
0.47375 0.52625 

Group means:
        x_1       x_2
2 0.1287798 0.2831542
7 0.2341133 0.2881862

Coefficients of linear discriminants:
          LD1
x_1 15.509380
x_2 -4.998317

y_hat_lda <- predict(train_lda, mnist_27$test)
confusionMatrix(data = y_hat_lda, reference = mnist_27$test$y)

Confusion Matrix and Statistics

          Reference
Prediction  2  7
         2 82 26
         7 24 68
                                         
               Accuracy : 0.75           
                 95% CI : (0.684, 0.8084)
    No Information Rate : 0.53           
    P-Value [Acc > NIR] : 1.266e-10      
                                         
                  Kappa : 0.4976         
                                         
 Mcnemar's Test P-Value : 0.8875         
                                         
            Sensitivity : 0.7736         
            Specificity : 0.7234         
         Pos Pred Value : 0.7593         
         Neg Pred Value : 0.7391         
             Prevalence : 0.5300         
         Detection Rate : 0.4100         
   Detection Prevalence : 0.5400         
      Balanced Accuracy : 0.7485         
                                         
       'Positive' Class : 2              
                                         

confusionMatrix(data = y_hat_lda, reference = mnist_27$test$y)$overall["Accuracy"]
cat("F1 Score:", F_meas(data = y_hat_lda, reference = mnist_27$test$y))

Accuracy: 0.75

F1 Score: 0.7663551

y_hat_lda <- predict(train_lda, mnist_27$true_p, type = "prob")[,2]

grid.arrange(plot_cond_prob(),
             plot_cond_prob(y_hat_lda, title = "LDA estimate"),
             ncol = 2)

The lack of flexibility of the LDA does not permit the method to capture the non-linearity of the problem. By applying the same standard deviations and correlations, the boundary is a line, as with logistic regression.

In the case of 784 predictors, it becomes unpratical to use QDA, due to overfitting and the number of correlations that must be calculated.
LDA would have to estimate over 600,000 parameters and QDA would have to estimate over 6,000,000 parameters.

QDA and LDA for 3 classes: 1, 2, 7

We will look up on the QDA and LDA models with 3 outcomes. The purpose here is to see how these techniques relies on multivariate normal distributions and how the accuracy gets lower when there are more outcomes and the predictors do not clearly follow a normal distribution.
We sample from these 3 outcomes, create training and test subsets and do some preprocessing:

ind <- sample(which(mnist$train$labels %in% c(1,2,7)), 2000)
x <- mnist$train$images[ind,]
y <- mnist$train$labels[ind]
index_train <- createDataPartition(y, p=0.8, list = FALSE)

# Binarize the data
x <- x > 200

# Proportion of dark pixels in the two quadrants
x <- cbind(rowSums(x[ ,upper_left_ind])/rowSums(x),
           rowSums(x[ ,lower_right_ind])/rowSums(x))

# Train and test sets
train_set <- data.frame(y = factor(y[index_train]),
                        x_1 = x[index_train,1],
                        x_2 = x[index_train,2])

test_set <- data.frame(y = factor(y[-index_train]),
                       x_1 = x[-index_train,1],
                       x_2 = x[-index_train,2])