2014-01-26

Isotonic Regressionと半順序

Isotonic Regression 半順序グラフ R quadprog igraph

n次元空間に点が散在しているとする
この点にはスカラー値が観察されている
このスカラー値は、n次元空間に面をなすスカラー値の分布の観察値である
観測点における、スカラー値の大小関係がわかっているものがあるときに、その大小関係制約を入れて、推定値を求めたい
その推定値は、ある評価関数を最小にするものとする
これはIsotonic regressionのアルゴリズムに乗せることができる
n次元上の点集合と、スカラー値不等式制約があるとする
すべての点ペアについて、あり得るすべてのスカラー値不等式制約を有向グラフで表すことができる
さらに、その不等式制約には冗長性があるので、それを排除することもできる
Rでやってみる

library(quadprog)
library(geometry)
library(rgl)
library(gtools)
library(igraph)
library(vegan)

# 次元
d <- 2
# 点の数
n.pt <- 100
#X <- matrix(runif(n.pt*d)-0.5,n.pt,d)
# d次元座標
X <- matrix(rnorm(n.pt*d),ncol=d)
#Y <- apply(X^0.3,1,sum)
#Y <- X[,1]^2+X[,2]^0.2
# 値のベクトル
Y <- 1/(1+1/exp(X[,1]*2)) * 1/(1+1/exp(X[,2]*0.8))
#Y <- 1/(1+1/exp((X[,1]+X[,2])*5))
# 乱雑項を入れる
Y.jit <- Y + rnorm(n.pt,0,0.1)
# ベルヌーイ事象にする
Y.jit <- rep(0,n.pt)
for(i in 1:n.pt){
	if(runif(1) < Y[i]){
		Y.jit[i] <- 1
	}
}

plot3d(cbind(X,Y.jit))

# 値の面を描くためにドロネー分割する
d.X <- delaunayn(X)

# ドロネー分割は、三角形(四面体…)の頂点組を返すが
# そこから重複の無い辺リストを作る関数
ed.list.delaunay <- function(d.X){
	n <- length(d.X[1,])
	comb <- combinations(n,2)
	tmp <- matrix(0,0,2)
	for(i in 1:length(comb[,1])){
		tmp <- rbind(tmp,cbind(d.X[,comb[i,1]],d.X[,comb[i,2]]))
	}
	tmp <- as.data.frame(tmp)
	tmp <- unique(tmp)
	tmp
}
# ドロネー分割の重複の無い辺のリスト
ed.list <- ed.list.delaunay(d.X)

plot(X)
segments(X[ed.list[,1],1],X[ed.list[,1],2],X[ed.list[,2],1],X[ed.list[,2],2])

plot3d(cbind(X,Y))
segments3d(X[c(t(ed.list)),1],X[c(t(ed.list)),2],Y[c(t(ed.list))])

points3d(cbind(X,Y.jit))

# ドロネーの辺について、その「値」の大小で向きを決める
ed.list.2 <- as.matrix(ed.list)
for(i in 1:length(ed.list.2[,1])){
	if(Y[ed.list.2[i,1]] >= Y[ed.list.2[i,2]]){
		ed.list.2[i,1] <- ed.list[i,2]
		ed.list.2[i,2] <- ed.list[i,1]
	}
}

g <- graph.empty(n.pt)
g <- graph.edgelist(ed.list.2)
plot(g,layout=X,vertex.size=0.2)

# すべての点のペアについて、isotonic regressionの不等式関係を確認する
# すべての座標軸の値で大小関係がそろっているときには、不等式関係を与える
D <- matrix(0,n.pt,n.pt)
for(i in 1:n.pt){
	tmp <- t(t(X)-X[i,])
	tmp.2 <- apply(tmp >= 0,1,prod)
	D[i,] <- tmp.2
}
diag(D) <- 0
# 不等式関係の２要素を列とする行列
uneq.ed <- which(D==1,arr.ind=TRUE)

g2 <- graph.empty(n.pt)
g2 <- add.edges(g2,t(uneq.ed))
#g2 <- graph.edgelist(uneq.ed)
plot(g2,layout=X,vertex.size=0.2)

# 有向グラフとそのエッジリストを与えて
# 半順序関係に注意してエッジ数を減らす関数

my.poset <- function(g,el){
	n.pt <- length(V(g))
	# この有向グラフの２点間の到達可能性を行列にする
	reachable <- shortest.paths(g,mode="out")
	reachable[which(reachable==Inf)] <- 0
	reachable <- sign(reachable)
	# 半順序を表すグラフにしたいので、
	# 「取り去って」も半順序関係が変わらないエッジを除去する
	removed <- c()
	for(i in 1:length(el[,1])){
		tmp.ed.list <- el[-unique(c(i,removed)),]
		tmp.g2 <- graph.empty() + vertices(V(g))
		tmp.g2 <- add.edges(tmp.g2,t(tmp.ed.list))
		#tmp.g2 <- graph.edgelist(tmp.ed.list)
		tmp.reachable <- shortest.paths(tmp.g2,mode="out")
		tmp.reachable[which(tmp.reachable==Inf)] <- 0
		tmp.reachable <- sign(tmp.reachable)
		if(sum(reachable-tmp.reachable)==0){
			removed <- unique(c(removed,i))
		}
	}
	if(length(removed)>0){
		uneq.ed2 <- uneq.ed[-unique(removed),]
	}else{
		uneq.ed2 <- uneq.ed
	}
	
	g2 <- graph.empty() + vertices(V(g))
	g2 <- add.edges(g2,t(uneq.ed2))
	#g2 <- graph.edgelist(uneq.ed)
	#plot(g2,layout=X,vertex.size=0.2)
	return(list(g=g2,ed.list=uneq.ed2))
}

poset.out <- my.poset(g2,uneq.ed)
plot(poset.out[[1]],layout=X,vertex.size=0.2)
uneq.ed2 <- poset.out[[2]]

# 半順序を満足する２つの制約でIsotonic regressionをする
# 半順序を満足する２つの制約とは、半順序に不要なエッジを取り去る前か、後かの２つ
Amat <- matrix(0,length(uneq.ed[,1]),n.pt)
for(i in 1:length(uneq.ed[,1])){
	Amat[i,uneq.ed[i,1]] <- -1
	Amat[i,uneq.ed[i,2]] <- 1
}
Amat2 <- matrix(0,length(uneq.ed2[,1]),n.pt)
for(i in 1:length(uneq.ed2[,1])){
	Amat2[i,uneq.ed2[i,1]] <- -1
	Amat2[i,uneq.ed2[i,2]] <- 1
}

Bmat <- diag(rep(1,n.pt))
ABmat <- rbind(Amat,-Bmat,Bmat)
ABmat2 <- rbind(Amat2,-Bmat,Bmat)
#ABmat <- rbind(Amat,-Bmat)
#bvec1 <- c(rep(0,length(Amat[,1])))
bvec <- c(rep(0,length(Amat[,1])),rep(-1,n.pt),rep(0,n.pt))
bvec2 <- c(rep(0,length(Amat2[,1])),rep(-1,n.pt),rep(0,n.pt))
#bvec <- c(rep(0,length(Amat[,1])),rep(-1,n.pt))
meq <- 0

Dmat <- diag(rep(1,n.pt))
dvec <- Y.jit
#out <- solve.QP(Dmat,dvec,t(Amat),bvec=bvec1,meq=0)
out <- solve.QP(Dmat,dvec,t(ABmat),bvec=bvec,meq=0)
out2 <- solve.QP(Dmat,dvec,t(ABmat2),bvec=bvec2,meq=0)

Y.est <- out[[1]]
Y.est2 <- out2[[1]]
# 半順序の制約の冗長度はあってもなくても結果は同じ
plot(Y.est,Y.est2)

plot3d(cbind(X,Y.est2))
segments3d(X[c(t(ed.list)),1],X[c(t(ed.list)),2],Y.est2[c(t(ed.list))])
points3d(cbind(X,Y),col=2)
points3d(cbind(X,Y.jit),col=3)