Skip to content

Instantly share code, notes, and snippets.

@ac00std ac00std/lda
Last active Jun 2, 2019

Embed
What would you like to do?
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.