0
$\begingroup$

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$):

rm(list=ls())
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)
  return(c(val0,val1,p0,p1))
}
classv(-3,m0,var0,m1,var1)
classv(-1,m0,var0,m1,var1)
classv(1,m0,var0,m1,var1)
classv(7.2,m0,var0,m1,var1)
###LDA
fit1 <- lda(y~x1,data = DF,prior=c(0.5,0.5))
pred <- predict(fit1,DF)
###QDA
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:

pred$posterior
          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:

pred2$posterior
          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.

$\endgroup$

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.