正準判別分析

```{r}

library(dplyr)

n=100

data=data.frame(num=1:n)



for(j in 1:3){
  
data_sub=data.frame(num=1:n,y1=rnorm(n,mean=sample(c(1:10),1),sd=3),y2=rnorm(n,mean=sample(c(1:10),1),sd=3)) 



if(j==1){

S1=var(data_sub[,2:3])

G1=apply(data_sub[,2:3],2,mean)

data_sub=data.frame(num=1:n,y1=rnorm(n,mean=5,sd=5),y2=rnorm(n,mean=-3,sd=5)) 
  
}else{
    
if(j==2){
  
S2=var(data_sub[,2:3])

G2=apply(data_sub[,2:3],2,mean)

data_sub=data.frame(num=1:n,y1=rnorm(n,mean=0,sd=5),y2=rnorm(n,mean=0,sd=5))
    
  
}else{
  
  
 S3=var(data_sub[,2:3])
 
 G3=apply(data_sub[,2:3],2,mean)
 
 data_sub=data.frame(num=1:n,y1=rnorm(n,mean=10,sd=5),y2=rnorm(n,mean=-10,sd=5))

  
  
}  
  
  
}

colnames(data_sub)=c("num",paste0("y1",j),paste0("y2",j))

data=left_join(data,data_sub,by="num")


}

n1=n2=n3=n

n=3*n

S=((n1-1)*S1+(n2-1)*S2+(n3-1)*S3)/(n1+n2+n3-3)

x=c(2,3)

D1=t(x-G1)%*%solve(S)%*%(x-G1)

D2=t(x-G2)%*%solve(S)%*%(x-G2)

D3=t(x-G3)%*%solve(S)%*%(x-G3)

a12=t(G1-G2)%*%solve(S);b12=-(1/2)*t(G1-G2)%*%solve(S)%*%(G1+G2)

a23=t(G2-G3)%*%solve(S);b23=-(1/2)*t(G2-G3)%*%solve(S)%*%(G2+G3)

a13=t(G1-G3)%*%solve(S);b13=-(1/2)*t(G1-G3)%*%solve(S)%*%(G1+G3)

library(rgl)

max_value=max(apply(data[,2:ncol(data)],2,max))

min_value=min(apply(data[,2:ncol(data)],2,min))

seq=seq(min_value,max_value,0.1)

data_sub1=data[,1:3]

colnames(data_sub1)=c("num","x1","x2")

data_sub2=data[,c(1,4,5)]

colnames(data_sub2)=c("num","x1","x2")

data_sub3=data[,c(1,6,7)]

colnames(data_sub3)=c("num","x1","x2")

data=rbind(data_sub1,data_sub2,data_sub3)

data_sub_sub1=rbind(data_sub1,data_sub2)

data_sub_sub2=rbind(data_sub2,data_sub3)

data_sub_sub3=rbind(data_sub3,data_sub1)

data_sub1=data_sub_sub1

data_sub2=data_sub_sub2

data_sub3=data_sub_sub3


data_sub1=data_sub1 %>% mutate(value=a12[1]*x1+a12[2]*x2+b12) %>% mutate(sign=ifelse(value>0,1,0))

data_sub2=data_sub2 %>% mutate(value=a23[1]*x1+a23[2]*x2+b23) %>% mutate(sign=ifelse(value>0,1,0))

data_sub3=data_sub3 %>% mutate(value=a13[1]*x1+a13[2]*x2+b13) %>% mutate(sign=ifelse(value>0,1,0))

```

```{r}

plot(data_sub1$x1[data_sub1$sign==1],data_sub1$x2[data_sub1$sign==1],type="p",col=1,xlim=c(min_value,max_value),ylim=c(min_value,max_value),xlab="x1",ylab="x2")

par(new=T)

plot(data_sub1$x1[data_sub1$sign==0],data_sub1$x2[data_sub1$sign==0],type="p",col=2,xlim=c(min_value,max_value),ylim=c(min_value,max_value),xlab="x1",ylab="x2")


```

```{r}


plot(data_sub2$x1[data_sub2$sign==1],data_sub2$x2[data_sub2$sign==1],type="p",col=3,xlim=c(min_value,max_value),ylim=c(min_value,max_value),xlab="x1",ylab="x2")

par(new=T)

plot(data_sub2$x1[data_sub2$sign==0],data_sub2$x2[data_sub2$sign==0],type="p",col=4,xlim=c(min_value,max_value),ylim=c(min_value,max_value),xlab="x1",ylab="x2")

```

```{r}

plot(data_sub3$x1[data_sub3$sign==1],data_sub3$x2[data_sub3$sign==1],type="p",col=5,xlim=c(min_value,max_value),ylim=c(min_value,max_value),xlab="x1",ylab="x2")

par(new=T)

plot(data_sub3$x1[data_sub3$sign==0],data_sub3$x2[data_sub3$sign==0],type="p",col=6,xlim=c(min_value,max_value),ylim=c(min_value,max_value),xlab="x1",ylab="x2")

```

```{r}

B=n1*(G1-apply(data[,2:3],2,mean))%*%t((G1-apply(data[,2:3],2,mean)))+n2*(G2-apply(data[,2:3],2,mean))%*%t((G2-apply(data[,2:3],2,mean)))+n3*(G3-apply(data[,2:3],2,mean))%*%t((G3-apply(data[,2:3],2,mean)))

W=(n1-1)*S1+(n2-1)*S2+(n3-1)*S3

y=solve(W)%*%B

plot=as.matrix(data[,2:3])%*%y

plot=data.frame(plot,sign=0)

plot$sign[1:100]=1

plot$sign[101:200]=2

plot$sign[201:300]=3



```

```{r}

plot(plot$y1[plot$sign==1],plot$y2[plot$sign==1],col=1,type="p",xlim=c(min(plot$y1),max(plot$y1)),ylim=c(min(plot$y2),max(plot$y2)),xlab="y1",ylab="y2")

par(new=T)

plot(plot$y1[plot$sign==2],plot$y2[plot$sign==2],col=2,type="p",xlim=c(min(plot$y1),max(plot$y1)),ylim=c(min(plot$y2),max(plot$y2)),xlab="y1",ylab="y2")

par(new=T)

plot(plot$y1[plot$sign==3],plot$y2[plot$sign==3],col=3,type="p",xlim=c(min(plot$y1),max(plot$y1)),ylim=c(min(plot$y2),max(plot$y2)),xlab="y1",ylab="y2")

```

```{r}

n=100

data=data.frame(x1=rnorm(n,sample(c(0:10),1),sample(c(1:10),1)),x2=rnorm(n,sample(c(0:10),1),sample(c(1:10),1)),x3=rnorm(n,sample(c(0:10),1),sample(c(1:10),1)),x4=rnorm(n,sample(c(0:10),1),sample(c(1:10),1)),x5=rnorm(n,sample(c(0:10),1),sample(c(1:10),1)),x6=rnorm(n,sample(c(0:10),1),sample(c(1:10),1)))

data=data %>% mutate(sign=kmeans(as.matrix(data),2)$cluster)

data1=data %>% filter(sign==1) %>% select(-sign)

data2=data %>% filter(sign==2) %>% select(-sign)

n1=nrow(data1);n2=nrow(data2)

mu=function(dat){
  
return(apply(dat,2,mean))  
  
}

S=function(dat1,dat2){
  
  return(((nrow(dat1)-1)*cov(dat1)+(nrow(dat2)-1)*cov(dat2))/(nrow(dat1)+nrow(dat2)-2))
  
}

f=c()

N=nrow(data);

for(j in 1:ncol(data1)){
  
dat1=data1[,-j];dat2=data2[,-j]

r=ncol(dat1)

D1=(mu(dat1)-mu(dat2))%*%solve(S(dat1,dat2))%*%t(t((mu(dat1)-mu(dat2))))

D2=(mu(data1)-mu(data2))%*%solve(S(data1,data2))%*%t(t((mu(data1)-mu(data2))))

F_value=(N-r)*(D2-D1)/((N*(N+2))/(nrow(dat1)*nrow(dat2))+D1)

f=c(f,F_value)
  
}

qf(1-0.15,1,N-r)

```