QDA

qda

#MASS読み込み
library(MASS)
library(mvtnorm)

sb=read.csv("sbnote.csv")

#学習用データの区分
sb_train=rbind(sb[1:50,],sb[101:150,])
sb_test=rbind(sb[51:100,],sb[151:200,])
test_x=sb_test[,-7]
test_y=sb_test$class

#二次判別用
sb_qda=qda(class~.,sb_train)
table(predict(sb_qda)$class,sb_train$class)

test_qda=predict(sb_qda,test_x)
table(test_qda$class,test_y)

#等分散の仮定を置かずに２変数に絞って判別分析を行う
gd=100
xp=seq(7,13,length=gd)
yp=seq(137,143,length=gd)

x2=cbind(sb[,4],sb[,6])
names(x2)=c("bottom","diagonal")
y=sb$class

sb.qda=qda(x2,y,prior=c(0.5,0.5))
pred=predict(sb.qda)

con=expand.grid(bottom=xp,diagonal=yp)

qpred=predict(sb.qda,con)
qp=qpred$posterior[,1]-qpred$posterior[,2]

plot(x2,pch=21,col=y+2,xlim=c(7,13),ylim=c(137,143),xlab="bottom",ylab="diagonal",main="２次判別の例")
contour(xp,yp,matrix(qp,gd),add=T,levels=0,col=4)

#誤識別したデータを図示
pairs(test_x,main="銀行紙幣の散布図",pch = 21, bg=4+as.numeric(test_qda$class)-test_y)



#等高線の図示
plot(x2,pch=21,col=y+2,xlim=c(7,13),ylim=c(137,143),xlab="bottom",ylab="diagonal",main="２次判別の例")
contour(xp,yp,matrix(qp,gd),add=T,levels=0,col=4)

honmono=subset(sb,sb$class==0)
nisemono=subset(sb,sb$class==1)

hx1=as.matrix(honmono$bottom)
hx2=as.matrix(honmono$diagonal)
hx=cbind(hx1,hx2)
hx.m=apply(hx,2,mean)
hx.v=var(hx)

f <- function(xp,yp,hm.v) { dmvnorm(matrix(c(xp,yp), ncol=2), mean=hx.m, sigma=hx.v) }
  # 上の sigma.positive を分散共分散行列とする密度関数
z <- outer(xp, yp, f,hm.v)
contour(xp, yp, z, add=T,levels=c(0.2,0.1,0.05,0.02,0.005,0.001,0.0001),col =6 )
  #共分散が正である場合の3Dグラフィクス描画

nx1=as.matrix(nisemono$bottom)
nx2=as.matrix(nisemono$diagonal)
nx=cbind(nx1,nx2)
nx.m=apply(nx,2,mean)
nx.v=var(nx)

f <- function(xp,yp,nm.v) { dmvnorm(matrix(c(xp,yp), ncol=2), mean=nx.m, sigma=nx.v) }
  # 上の sigma.positive を分散共分散行列とする密度関数
z <- outer(xp, yp, f,nm.v)
contour(xp, yp, z, add=T,levels=c(0.2,0.1,0.05,0.02,0.005,0.001,0.0001),col = 5)
  #共分散が正である場合の3Dグラフィクス描画

par(new=T)
plot(hx.m[1],hx.m[2],xlim=c(7,13),ylim=c(137,143),xlab="",ylab="",pch = "重心",font=2,col = 6,ps = 40)
par(new=T)
plot(nx.m[1],nx.m[2],xlim=c(7,13),ylim=c(137,143),xlab="",ylab="",pch = "重心",font=2,col = 5,ps = 40)