LDA

lda

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,])

#散布図を描く
y=sb$class
x=sb[,-7]

pairs(x,main="銀行紙幣の散布図",pch = 21, bg=y+2)

#線形判別分析
train_x=sb_train[,-7]
train_y=sb_train$class
sb_lda=lda(train_x,train_y,prior=c(0.5,0.5))
sb_pred=predict(sb_lda)

#学習用データの結果
table(sb_pred$class,train_y)

#テスト用データの結果
test_x=sb_test[,-7]
test_y=sb_test$class
test_lda=predict(sb_lda,test_x)
table(test_lda$class,test_y)


#2変数に絞った時の散布図
x1=as.matrix(x$bottom)
x2=as.matrix(x$diagonal)
x=cbind(x1,x2)

#散布図を描く
plot(x,xlim=c(7,13),ylim=c(137,143),xlab="bottom",ylab="diagonal",main="銀行紙幣の散布図",col=as.numeric(y)+2)

#等分散の仮定を用いて線形判別分析を行う
sb.lda=lda(x,y,prior=c(0.5,0.5))
pred=predict(sb.lda)
pred

p1=predict(sb.lda)$class
table(y,p1)

mu=apply(sb.lda$means,2,mean)            #列ごとの平均を求め、真偽に関係なく平均を求める
a0=as.vector(sb.lda$scaling)%*%mu
keisu=as.vector(sb.lda$scaling)

c(a0,keisu)    #判別関数の係数を得る
abline(a=a0/keisu[2], b=-keisu[1]/keisu[2])  #散布図に線形判別曲線を描く



#等高線の図示
gd=100
xp=seq(7,13,length=gd)
yp=seq(137,143,length=gd)

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)

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

#線形判別分析では各群団で共通の共分散行列を使用する
tx.v=(hx.v+nx.v)/2

f <- function(xp,yp,hm.v) { dmvnorm(matrix(c(xp,yp), ncol=2), mean=hx.m, sigma=tx.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グラフィクス描画

f <- function(xp,yp,nm.v) { dmvnorm(matrix(c(xp,yp), ncol=2), mean=nx.m, sigma=tx.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)