遠近関係をサポートベクターマシンに持ち込む

• カーネル・トリックについて昨日書いた
• 簡単にまとめると、動かす前の測度から簡単に計算できるような測度を動かした後の測度として定義すると線形サポートベクターマシンに持ち込める、ということだった
• Rのe1071パッケージのsvm()関数では、動かす前の座標から動かした後の測度を計算してくれるカーネル関数をオプションで指定している
• したがって、オプションにあるカーネル関数は使えるが、それ以外の変換は受け付けてくれない
• 「それ以外の変換」がしたくなったらどうすればいいか、という話
• 変換前の座標から、変換後のペアワイズ距離を出す
• その距離行列からmdsで座標を決めて、それに対してサポートベクターマシンする
• 変換後のペアワイズ距離がどれくらい簡単に計算できるか、鍵

############
############
# のたくり雲作成 ここから
############
############
library(igraph)
library(rgl)
# サンプル数
Ns <- 500
# マーカー数
Nm <- 50
# サンプルのパターン数(群数)
Ns.pt <- Ns
# マーカーのパターン数(群数)
Nm.pt <- 10
# サンプル・マーカーの存在位置を多次元空間酔歩の道として作る
trail <- matrix(rnorm(Ns.pt*Nm.pt),Ns.pt,Nm.pt)
trail <- apply(trail,2,cumsum)
# ３次元分だけ見てみよう
library(rgl)
plot3d(trail[,1:3])
matplot(trail,type="l")
# パターン数(群数)ごとにいくつのサンプル、いくつのマーカーを帰属させるかをランダムに決める
library(MCMCpack)
ps <- rdirichlet(1,rep(1,Ns.pt))
pm <- rdirichlet(1,rep(1,Nm.pt))
ss <- sample(1:Ns.pt,Ns,replace=TRUE,prob=ps)
sm <- sample(1:Nm.pt,Nm,replace=TRUE,prob=pm)
# 行数=サンプル数、列数=マーカー数の行列
M <- trail[ss,sm]
# 少し乱す
M <- jitter(M,100)

############
############
# のたくり雲作成 ここまで
############
############

# 雲の重心からの距離に応じて、0,1のフェノタイプを与える
center <- apply(M,2,mean)
d.from.center <- apply((t(M)-center)^2,2,sum)

phenotype <- rep(0,length(M[,1]))

phenotype[which(d.from.center>mean(d.from.center))] <- 1

# フェノタイプで色分けして３次元分だけプロット
plot3d(M[,1:3],col = phenotype +1)

library(vegan)

# MSTオブジェクトを経由してグラフオブジェクトを作る
make.mst.graph <- function(X){
	mst.x <- spantree(dist(X))
	e.x <- cbind(2:(length(mst.x$kid)+1),mst.x$kid)
	#graph.edgelist(e.x)
	return(list(g = graph.edgelist(e.x),mst=mst.x))
}

# g は最小全域木のグラフオブジェクト
g <- make.mst.graph(M)
# このグラフのペアワイズ点間距離行列を作る
sp.g <- shortest.paths(g$g)

#plot(g$g)

# 最小全域木上のパスで定められた距離行列からMDSする
mds <- cmdscale(sp.g,length(M[1,]))



#mds.data <- cbind(mds,factor(phenotype))
# MDS座標を使って、フェノタイプの判別をサポートベクターマシンする
library(e1071)
svm.out <- svm(phenotype ~ mds)

print(svm.out)
summary(svm.out)

pred <- predict(svm.out,mds)
table(pred,phenotype)

pred <- predict(svm.out,mds,decision.values = TRUE)
attr(pred,"decision.values")

plot(cmdscale(dist(mds)),col=as.integer(phenotype+1))