I am using linear discriminant analysis (LDA) and quadratic discriminant analysis (QDA) in R. I am working on a small data set of 4 observations and two variables y class variable with outcome 0 and 1 and a continuous variable x1. The definition for LDA states that in calculating the probability for each class it is assumed that variables, in this case x1, come from a gaussian distribution. So we would have for LDA: $$p(y_i=0|x_i)=\frac{p(x_i|y_i=0)p(y_i=0)}{p(x_i)}$$ $$p(y_i=1|x_i)=\frac{p(x_i|y_i=1)p(y_i=1)}{p(x_i)}$$ Where $p(y_i=0)=\pi$ and $p(y_i=0)=1-\pi$ are prior probabilities. We can throw away $p(x_i)$ and we would have: $$p(y_i=0|x_i)\;\alpha\;p(x_i|y_i=0)\pi$$ $$p(y_i=1|x_i)\;\alpha\;p(x_i|y_i=1)(1-\pi)$$ Considering gaussian distribution we would end up with: $$p(y_i=0|x_i)\;\alpha\;\pi \frac{1}{\sqrt{2\pi\sigma_0^2}}e^{\frac{-(x_i-\mu_0)^2}{2\sigma_0^2}}$$ $$p(y_i=1|x_i)\;\alpha\;(1-\pi) \frac{1}{\sqrt{2\pi\sigma_1^2}}e^{\frac{-(x_i-\mu_1)^2}{2\sigma_1^2}}$$ Then this values can be normalized and parameters $\mu$ and $\sigma^2$ can be computed for each class using MLE. Using the previous definitions I replicated on R with the next information (Moreover, I obtained lda and qda for the dataset and $\pi=0.5$):

x1 <- c(-3,-1,1,7.2)
y <- c(0,1,1,0)
DF <- data.frame(x1,y)
#MLE for each class
m0 <- mean(DF$x1[DF$y==0])
m1 <- mean(DF$x1[DF$y==1])
var0 <- var(DF$x1[DF$y==0])
var1 <- var(DF$x1[DF$y==1])
#Compute class for each data point
classv <- function(x,mu,varv,mu2,varv2)
  val0 <- 0.5*dnorm(x,mu,sqrt(varv))
  val1 <- 0.5*dnorm(x,mu2,sqrt(varv2))
  p0 <- val0/(val0+val1)
  p1 <- val1/(val0+val1)
fit1 <- lda(y~x1,data = DF,prior=c(0.5,0.5))
pred <- predict(fit1,DF)
fit2 <- qda(y~x1,data = DF,prior=c(0.5,0.5))
pred2 <- predict(fit2,DF)

My question arises in this point. The estimations computed by hand using classv for each data point are the next (Last two elements of each result are probability for class 0 and 1 respectively):

> classv(-3,m0,var0,m1,var1)
[1] 0.02153879 0.01486629 0.59164248 0.40835752
> classv(-1,m0,var0,m1,var1)
[1] 0.02521621 0.10984782 0.18669820 0.81330180
> classv(1,m0,var0,m1,var1)
[1] 0.02733657 0.10984782 0.19926879 0.80073121
> classv(7.2,m0,var0,m1,var1)
[1] 0.0215387886994 0.0000003318246 0.9999845943289 0.0000154056711

But when I check the fit1 object (LDA) I get the next:

          0         1
1 0.4219232 0.5780768
2 0.4602378 0.5397622
3 0.4990281 0.5009719
4 0.6173124 0.3826876

Which are totally different from the ones computed with the classv function. On the other side, when I analyse the posterior probabilities from fit2 (QDA) I get:

          0             1
1 0.5916425 0.40835752195
2 0.1866982 0.81330180161
3 0.1992688 0.80073120619
4 0.9999846 0.00001540567

That are the same as the ones computed from classv function which considers the definition for LDA. So, I would like help to clarify if I am doing wrong with LDA or QDA because I should get the same results from fit1 using the manual implementation from classv but actually I am getting the results that a QDA produces. Many thanks for your help.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.