主成分分析による手書き数字の復元 - INFINITY -数学とかプログラミングとか-

手書き数字について、主成分分析を行うことで、特徴抽出、画像の文字の復元することができる。

https://www.kaggle.com/c/digit-recognizer/data
のtrain.csvデータにheaderにlabel, pixel0 から pixel783まで定義されている。

train.csvを用いて、手書き数字の復元する。

どの数字でも良いのですが、今は0にする。0である手書き数字の標本数は4131である。

Xを $p \times n$ の0である標本ベクトル(縦のp次元ベクトル)を横に標本数分並べる。

つまり、
$X=(\b{x_1},\cdots,\b{x_n})$
とする。

このとき、 $E\{(X-E(X))(X'-E(X'))\}$ の固有ベクトルは主成分ベクトルになる。(詳細は略)

固有値が大きい順に並べて、対応する固有ベクトルが第1主成分,第2主成分,...に対応する。

第1列X[,1]の画像を復元しようとおもったら、

$E:=(\b{e_1},\cdots,\b{e_n})'$ とするとき
$\b{y}=E X[,1]$ より、
$E\b{y}=EE'\b{x}=\b{x}$
が成立するので,
受信側が固有ベクトル $\b{e_i}$ を持っていれば、
送信側は、 $y_i$ を送信すれば、受信側は、復元できることになる。

しかも、 $y_i,\b{e_i}$ はそれほど多く取得しておく必要もない。

なぜなら、以下のように、画像の復元が少ない主成分(今は16,32,64,128)で復元できているからである。

今、n=4131で、m=16,32,64,128について、

$\b{\hat{x}}=E\b{y}=\sum_{i=1}^{m}y_i\b{e_i},(m < n)$
が成立しているからである。

f:id:infinity_th4:20130209035939p:plain
f:id:infinity_th4:20130209045250p:plain
f:id:infinity_th4:20130209045654p:plain
f:id:infinity_th4:20130209050100p:plain

最後に、0の固有ベクトルを25個まで表示する。
f:id:infinity_th4:20130209174108p:plain

所詮、PCAは2次の情報(相関)しか持っていないので、精度は悪い。
今度は、ICAを試すことにする。

お腹が痛くなってきたのでこの辺で。。

githubにコードを上げようと思ったけどやめた。

setwd("./Dropbox/R/handwritten")

getwd()

## Read in data
train <- read.csv("./data/train.csv", header=TRUE)
train<-as.matrix(train)
dim(train)
#42000   785

##Color ramp def.
colors<-c('white','black')
cus_col<-colorRampPalette(colors=colors) #16?i???ŕԂ?

cus_col(256)

## Plot the average image of each digit
par(mfrow=c(1,5),pty='s',mar=c(1,1,1,1),xaxt='n',yaxt='n')
all_img<-array(dim=c(10,28*28)) #NA??10*784?̍s???쐬
#all_img

di0 = train[train[,1]==0,-1]
di0mat=di0[-1,]
dim(di0mat) #4131  784

di0.cor=cor(di0mat)

di0.cor


X=t(di0mat)
dim(X) #p x n

Xvar=function(X){
    return((1/n)*rowSums(X))#各行の和
}

xvar=Xvar(X) 
length(xvar)

xvar#標本平均ベクトル

n=ncol(X)
n
p=nrow(X)
p

#標本分散共分散行列の初期化
S=matrix(0,nrow=p,ncol=p)
#S
#(X[,i]-xvar)%*%t(X[,i]-xvar)
for(i in 1:n){
    S=S+(X[,i]-xvar)%*%t(X[,i]-xvar)
}

S=(1/n)*S

S
dim(S)

eigen(S)$values
eigen(S)$vectors[,1]
length(eigen(S)$vectors[,1]) #p


di=0
m=128
yvec=matrix(0,nrow=m,ncol=1)
for(i in 1:m){
    yvec[i]=t(eigen(S)$vectors[,i])%*%X[,1] #1番目の標本の第i主成分
}


z<-array(X[,1],dim=c(28,28)) #28*28?̉摜?ɂ??邽?߁B
z<-z[,28:1] ##right side up #??28??????1???ɁA??27??????2???ɁA....
image(1:28,1:28,z,main="original data",col=cus_col(256))
n
restorevec=matrix(0,nrow=p,ncol=1)
for(i in 1:m){
    restorevec=restorevec+yvec[i]*eigen(S)$vectors[,i] #1番目の標本の第i主成分
    if(i==16||i==32||i==64||i==128){
        z<-array(restorevec,dim=c(28,28)) #28*28?̉摜?ɂ??邽?߁B
        z<-z[,28:1] ##right side up #??28??????1???ɁA??27??????2???ɁA....
        image(1:28,1:28,z,main=i,col=cus_col(256))
    }
}

#2番目の標本の復元###################
m=128
z<-array(X[,2],dim=c(28,28)) #28*28?̉摜?ɂ??邽?߁B
z<-z[,28:1] ##right side up #??28??????1???ɁA??27??????2???ɁA....
image(1:28,1:28,z,main="original data",col=cus_col(256))


yvec=matrix(0,nrow=m,ncol=1)
for(i in 1:m){
    yvec[i]=t(eigen(S)$vectors[,i])%*%X[,2] #2番目の標本の第i主成分
}

restorevec=matrix(0,nrow=p,ncol=1)
for(i in 1:m){
    restorevec=restorevec+yvec[i]*eigen(S)$vectors[,i] #1番目の標本の第i主成分
    if(i==16||i==32||i==64||i==128){
        z<-array(restorevec,dim=c(28,28)) #28*28?̉摜?ɂ??邽?߁B
        z<-z[,28:1] ##right side up #??28??????1???ɁA??27??????2???ɁA....
        image(1:28,1:28,z,main=i,col=cus_col(256))
    }
}

#3番目の標本の復元###################
m=128
z<-array(X[,3],dim=c(28,28)) #28*28?̉摜?ɂ??邽?߁B
z<-z[,28:1] ##right side up #??28??????1???ɁA??27??????2???ɁA....
image(1:28,1:28,z,main="original data",col=cus_col(256))


yvec=matrix(0,nrow=m,ncol=1)
for(i in 1:m){
    yvec[i]=t(eigen(S)$vectors[,i])%*%X[,3] #3番目の標本の第i主成分
}

restorevec=matrix(0,nrow=p,ncol=1)
for(i in 1:m){
    restorevec=restorevec+yvec[i]*eigen(S)$vectors[,i] #1番目の標本の第i主成分
    if(i==16||i==32||i==64||i==128){
        z<-array(restorevec,dim=c(28,28)) #28*28?̉摜?ɂ??邽?߁B
        z<-z[,28:1] ##right side up #??28??????1???ɁA??27??????2???ɁA....
        image(1:28,1:28,z,main=i,col=cus_col(256))
    }
}

#4番目の標本の復元###################
m=128
z<-array(X[,4],dim=c(28,28)) #28*28?̉摜?ɂ??邽?߁B
z<-z[,28:1] ##right side up #??28??????1???ɁA??27??????2???ɁA....
image(1:28,1:28,z,main="original data",col=cus_col(256))


yvec=matrix(0,nrow=m,ncol=1)
for(i in 1:m){
    yvec[i]=t(eigen(S)$vectors[,i])%*%X[,4] #4番目の標本の第i主成分
}

restorevec=matrix(0,nrow=p,ncol=1)
for(i in 1:m){
    restorevec=restorevec+yvec[i]*eigen(S)$vectors[,i] #1番目の標本の第i主成分
    if(i==16||i==32||i==64||i==128){
        z<-array(restorevec,dim=c(28,28)) #28*28?̉摜?ɂ??邽?߁B
        z<-z[,28:1] ##right side up #??28??????1???ɁA??27??????2???ɁA....
        image(1:28,1:28,z,main=i,col=cus_col(256))
    }
}

#0の固有ベクトルを25個まで表示する##########
par(mfrow=c(5,5),pty='s',mar=c(1,1,1,1),xaxt='n',yaxt='n')
l=5*5
for(i in 1:l){
    z<-array(eigen(S)$vectors[,i],dim=c(28,28)) #28*28?̉摜?ɂ??邽?߁B
    z<-z[,28:1] ##right side up #??28??????1???ɁA??27??????2???ɁA....
    image(1:28,1:28,z,main=i,col=cus_col(256))
}