2013-05-19

なんちゃって平滑化その２

度数分布平滑化視覚

こちらで『なんちゃって度数分布平滑化』というのをやった
多次元に拡張しよう
やり方は同じ。乳幼児の学習過程レベルの処理に限定する
多次元の点分布を「感知」し、それを隣接細胞の刺激を順次足し合わせていく
また、順次、差分を取る。ただし多次元に上がったので一次の差分は１地点あたり次元数の方向の偏差分になる
また、１次元度数分布では、２次の差分も取ったが、多(n)次元に上げると $n^2$ の偏差分が必要になる
偏差分の成分が $n^2$ と大きくなることからわかるように、各点の勾配情報はn方向に関する、単調増・単調減・極大・極小の４通りについてn方向組み合わせになってくる
まずは、和をとって平滑化することと、２次の偏差分までとってみよう
次元方向の２次の偏差分の正負入れ替わりで最適化してみる
２次元程度なら視細胞数がそこそこだが、次元が上がると、素子数が多くなってコンピュータ上では問題が生じるが…

# データの座標によらず、視野内に収めるための補正
vision.scale <- function(x,scl){
	(x-min(x))/(max(x)-min(x)) * scl + (1-scl)/2
}
# 観察対象をきちんと視野内に収めるために、周辺部にはシグナルがないように調整する
vision.standard <- function(x,scl = 0.9){
	rg <- apply(x,2,range)
	tmp <- apply(x,2,vision.scale, scl)
	return(list(st.x = tmp,rg=rg))
}
# 1辺n.cell数の多次元視覚器にシグナルが入ることは、多次元度数分布を作ることだが
# floor()関数を使って、その作業をする
# これで視覚第１層のシグナルパターンが得られる
multi.dim.vision <- function(x,n.cell=50,scl=0.5){
	d <- length(x[1,])
	tmp.out <- vision.standard(x,scl=scl)
	rg <- tmp.out$rg
	X <- matrix(floor(tmp.out$st.x * n.cell),ncol=d)
	X.val <- apply(t(X) * n.cell^(0:(d-1)),2,sum)+1
	vs <- array(0,rep(n.cell,d))
	for(i in 1:length(X[,1])){
		vs[X.val[i]] <- vs[X.val[i]] + 1
	}
	vs
}
# 1層ずつ隣接シグナルを加算する
one.step.array <- function(a){
	dm <- dim(a)
	d <- length(dm)
	n <- dm[1]
	shifts <- as.matrix(expand.grid(rep(list(0:1),d)))
	base.address <- as.matrix(expand.grid(rep(list(0:(n-2)),d)))
	ret <- array(0,rep(n-1,d))
	for(i in 1:length(shifts[,1])){
		tmp.shift <- t(t(base.address) + shifts[i,])
		tmp.shift.2 <- apply(t(tmp.shift) * n^(0:(d-1)),2,sum)+1
		ret <- ret + a[tmp.shift.2]
	}
	ret
}
# 一階差分
diff.partial <- function(a){
	dm <- dim(a)
	d <- length(dm)
	n <- dm[1]
	diff.list <- list()
	for(i in 1:d){
		diff.list[[i]] <- array(apply(a,i,diff),rep(n-1,d))
	}
	diff.list
}
# 二階差分
diff2.partial <-function(diff.a){
	dm <- dim(diff.a[[1]])
	d <- length(dm)
	n <- dm[1]
	ret <- list()
	for(i in 1:d){
		ret[[i]] <- diff.partial(diff.a[[i]])
	}
	ret
}
# 二階差分
diff2.partial.2 <-function(a){
	dm <- dim(a)
	d <- length(dm)
	n <- dm[1]
	diff.list <- list()
	for(i in 1:d){
		diff.list[[i]] <- array(apply(a,i,diff),rep(n-1,d))
		diff.list[[i]] <- array(apply(diff.list[[i]],i,diff),rep(n-1,d))
	}
	diff.list
}

# n次元分布を作ろう
n <- 3
x1 <- c(rnorm(100,10),rnorm(50,80,20),runif(150,30,50))
x2 <- c(rnorm(50,5),rnorm(150,20,10),runif(100,4,50))
x3 <- c(rnorm(150,40),rnorm(50,10,100),runif(100,40,50))
x <- cbind(x1,x2,x3)
x <- cbind(x1,x2)
n <-2
library(rgl)
#plot3d(x)
#x <- matrix(rnorm(300),100,3)
m.out <- multi.dim.vision(x,n.cell=100)
plot(x)
dev.new()
#image(m.out)
# コントラストがきついので対数スケールでプロットしてみる
# ３方向からのプロット
#image(apply(log(m.out+1),c(1,2),sum))
#image(apply(log(m.out+1),c(2,3),sum))
#image(apply(log(m.out+1),c(2,3),sum))

# １層ずつ平滑化
tmp.m.out <- m.out
#tmp.m.out.2 <- one.step.array(tmp.m.out)

par(ask=TRUE)
n.step <- 20
X.series <- Diff1.series <- Diff2.series <- Diff2.series.2 <- list()
v.series <- rep(0,n.step)
for(i in 1:n.step){
	X.series[[i]] <- one.step.array(tmp.m.out)
	X.series[[i]] <- X.series[[i]]-min(X.series[[i]])
	X.series[[i]] <- X.series[[i]]/max(X.series[[i]])
	tmp.m.out <- X.series[[i]]
	# １階差分
	Diff1.series[[i]] <- diff.partial(X.series[[i]])
	# ２階差分
	Diff2.series[[i]] <- diff2.partial(Diff1.series[[i]])
	# ２階差分の２
	Diff2.series.2[[i]] <- diff2.partial.2(X.series[[i]])
	# ２階差分のがたぼこ評価
	#v.series[i] <- (sum(abs(unlist(Diff2.series[[i]]))) - sum(unlist(Diff2.series[[i]])))/length(unlist(Diff2.series[[i]]))
	v.series[i] <- (sum(abs(unlist(Diff2.series.2[[i]]))) - sum(unlist(Diff2.series.2[[i]])))/length(unlist(Diff2.series.2[[i]]))

	#v.series[i] <- (sum(abs(unlist(Diff1.series[[i]]))) - sum(unlist(Diff1.series[[i]])))/length(unlist(Diff1.series[[i]]))
	
	#image(apply(log(X.series[[i]]+1),1:2,sum))
	image(log(X.series[[i]]+1))
}
# １階、２階の偏差分
#diff.x <- diff.partial(x)
#diff2.x <- diff2.partial(diff.x)

choice <- which(v.series==min(v.series))[1] -1

image(log(X.series[[choice]]+1))
dev.new()
par(ask=TRUE)
for(i in 1:n.step){
	image(log(X.series[[i]]+1),main=i)
}