# Linear discriminant analysis against quadratic discriminant analysis behavior in R

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.