#图3.1
#左图
rm(list=ls())
x=seq(from=-5.6,to=14,by=0.01)
y1=dnorm(x,0,2)   #计算密度
y2=dnorm(x,8,2)   #计算密度
plot(x,y1,type="l",xlab="", ylab="",axes=F,xlim=c(-5.8,14), ylim=c(-0.03,0.4))
lines(x,y2,type="l")
rect(-5.8, 0, 14, 0)
rect(0, 0, 0, 0.198)
rect(4, 0, 4, 0.025)
rect(8, 0, 8, 0.198)
text(4,-0.015,expression(bar(italic(μ))),cex=0.7)
text(0,-0.015,expression(italic(μ)[1]),cex=0.7)
text(8,-0.015,expression(italic(μ)[2]),cex=0.7)
text(0,0.225,expression(italic(X)[1]),cex=0.7)
text(8,0.225,expression(italic(X)[2]),cex=0.7)
x=seq(from=4,to=7,by=0.01)
y1=dnorm(x,0,2)
x[1]=4;y1[1]=0
polygon(x, y1-0.0001,col="#dddddd")
#右图
rm(list=ls())
x=seq(from=-5.6,to=14,by=0.01)
y1=dnorm(x,0,2)   #计算密度
y2=dnorm(x,8,2)   #计算密度
plot(x,y2,type="l",xlab="", ylab="",axes=F,xlim=c(-5.8,14), ylim=c(-0.03,0.4))
lines(x,y1,type="l")
rect(-5.8, 0, 14, 0)
rect(0, 0, 0, 0.198)
rect(5, 0, 5, 0.065)
rect(8, 0, 8, 0.198)
text(5,-0.015,expression(italic(ξ)),cex=0.7)
text(0,-0.015,expression(italic(μ)[1]),cex=0.7)
text(8,-0.015,expression(italic(μ)[2]),cex=0.7)
text(0,0.225,expression(italic(X)[1]),cex=0.7)
text(8,0.225,expression(italic(X)[2]),cex=0.7)
x=seq(from=5,to=0,by=-0.01)
y2=dnorm(x,8,2)
x[1]=5;y2[1]=0
polygon(x, y2-0.0001,col="#dddddd")
x=seq(from=-5.6,to=14,by=0.01)
y1=dnorm(x,0,2)   #计算密度
lines(x,y1,type="l")


#案例3.1精神病距离判别
DA1<-function
　   (TrnX1, TrnX2, TstX = NULL, var.equal = FALSE){if (is.null(TstX) == TRUE) TstX<-rbind(TrnX1,TrnX2)
　   if (is.vector(TstX) == TRUE)  TstX<-t(as.matrix(TstX)) else if (is.matrix(TstX) != TRUE)
　      TstX<-as.matrix(TstX)
　   if (is.matrix(TrnX1) != TRUE) TrnX1<-as.matrix(TrnX1)
　   if (is.matrix(TrnX2) != TRUE) TrnX2<-as.matrix(TrnX2); nx<-nrow(TstX)
　   blong<-matrix(rep(0, nx), nrow=1, byrow=TRUE, dimnames=list("blong", 1:nx))
　   mu1<-colMeans(TrnX1); mu2<-colMeans(TrnX2) 
　   if (var.equal == TRUE  || var.equal == T){S<-var(rbind(TrnX1,TrnX2))
　      w<-mahalanobis(TstX, mu2, S)-mahalanobis(TstX, mu1, S)} else{S1<-var(TrnX1); S2<-var(TrnX2)
　      w<-mahalanobis(TstX, mu2, S2)-mahalanobis(TstX, mu1, S1)}
　   for (i in 1:nx){if (w[i]>0) blong[i]<-1 else blong[i]<-2}; blong}
classX1<-matrix(c(22,20,23,23,17,24,23,18,22,19,20,20,21,13,20,19,20,18,20,23,23,25,23,21,23,
　     6,14,9,1,8,9,13,18,16,18,17,31,9,13,14,18,11,17,7,6,23,9,5,12,7),ncol=2)
classX2<-matrix(c(24,19,11,6,9,10,3,15,14,20,8,20,14,3,10,22,11,6,20,20,15,2,10,13,12,
　     38,36,43,60,32,17,17,56,43,8,46,62,36,12,51,22,30,30,61,43,48,53,43,19,4),ncol=2)
DA1(classX1, classX2, var.equal=TRUE)　#样本协方差相同

# 样本协方差不同
DA1(classX1, classX2)

#输入待测样本
NEW<-matrix(c(28,36),ncol=2)
DA1(classX1, classX2, NEW, var.equal=TRUE)

#案例3.2 Fisher Iris数据
source("DA2.R")
   X<-iris[,1:4];G<-gl(3,50);DA2(X,G)

#案例3.3 海南板块的数据判别
read.csv("hnpb.csv",header=T)->DA
X<-DA[,3:6]
G<-gl(3,6)
DA2(X,G)

#案例3.4用Fisher判别解例3.1
source("DAf.R")
classX1<-matrix(c(22,20,23,23,17,24,23,18,22,19,20,20,21,13,20,19,20,18,20,23,23,25,23,21,23,
　     6,14,9,1,8,9,13,18,16,18,17,31,9,13,14,18,11,17,7,6,23,9,5,12,7),ncol=2)
classX2<-matrix(c(24,19,11,6,9,10,3,15,14,20,8,20,14,3,10,22,11,6,20,20,15,2,10,13,12,
　     38,36,43,60,32,17,17,56,43,8,46,62,36,12,51,22,30,30,61,43,48,53,43,19,4),ncol=2)
  DAf(classX1, classX2)

#案例3.5 交通设施和教育传媒板块股票的判别分析
#距离判别
source("DA1.R")
read.csv("jt.csv",header=F)->DAjt
read.csv("jy.csv",header=F)->DAjy
#样本协方差相同
DA1(DAjt, DAjy, var.equal=TRUE)　
# 样本协方差不同
DA1(DAjt, DAjy)

#协方差检验
read.csv("jt.csv",header=F)->DAjt
read.csv("jy.csv",header=F)->DAjy
DAjt1=cov(DAjt)
DAjy1=cov(DAjy)
var.test(DAjt1, DAjy1)

#输入待测样本
source("DA1.R")
read.csv("njt.csv",header=F)->DAnjt
DA1(DAjt, DAjy, DAnjt)
read.csv("njy.csv",header=F)->DAnjy
DA1(DAjt, DAjy, DAnjy)


#费歇判别
source("DAf.R")
read.csv("jt.csv",header=F)->DAjt
read.csv("jy.csv",header=F)->DAjy
DAf(DAjt, DAjy)

#输入待测样本
read.csv("njt.csv",header=F)->DAnjt
DAf(DAjt, DAjy, DAnjt)
read.csv("njy.csv",header=F)->DAnjy
DAf(DAjt, DAjy, DAnjy)

####多日数据

#距离判别
source("DA1.R")
read.csv("jtd.csv",header=F)->DAjt
read.csv("jyd.csv",header=F)->DAjy
#样本协方差相同
DA1(DAjt, DAjy, var.equal=TRUE)　
#样本协方差不同
DA1(DAjt, DAjy)　

#费歇判别
source("DAf.R")
read.csv("jtd.csv",header=F)->DAjt
read.csv("jyd.csv",header=F)->DAjy
DAf(DAjt, DAjy)