2013-05-18

なんちゃって平滑化

度数分布平滑化視覚

こちらで『なんちゃってPCA』というのをやった
どういうことかというと、乳幼児(のようにPCAとか行列とか算数とかがわからない生物)は、そんなことを知らなくても、視覚処理をするし、眺めるべき方向を選択することはできる、という話だった
じゃあ、そんな「なんちゃって」な乳幼児が１次元空間の点の標本分布を眺めるとき、どんな風に処理するのだろう、処理して「かいつまむ」としたら、どんな風にするだろう、という話
「なんちゃって乳幼児」は、１次元空間上の密度を知りたい、などとは考えない
視覚刺激をなんらかの系統だった、しかも単純な処理をすることで、「いい感じ」な視神経-脳神経電気刺激パターンを得られるように視覚処理系を訓練していくはず
単純な処理としては、網膜の光刺激で励起する細胞の第１層から始まってそれを多段階的に処理する多層があるだろう。そして、その多層の細胞は「和」と「差」を取ることくらいしかしない
n個の細胞が１列に並び、そこが刺激される。第2層以降は、全部でn層あって、n-1,n-2,...,1個の細胞があるようなパスカルの三角形的な層構造にしてみよう(ちょっと厚すぎるが)
そのうえで、第i層は第i-1層の隣同士の２細胞の刺激の和を計算することにする
ここまでが「和」の処理
次は「差」の処理
「和」の処理によって、１列の長さは短くなりつつ、だんだん平滑化していくが、問題はどこまで短くするか、どこまで滑らかにするか、である
平滑〜滑らか、というのは、「がたぼこしていない」ということだが、「がたぼこ」というのは、どういう状態だろう
隣接する細胞の値の差を符号付きで全部足せば、それは結局、最初の細胞の値と最後の細胞の値の差になるが、細胞間の差の絶対値を足し合わせると、単調のときは、符号付きの和と同じになるが、単調でなくがたぼこすると、大きくなっていくので、この「単調のときよりも大きくなる程度」を定量してやると、なにかしらがたぼこの程度と関係する値になりそうだ。しかも、この処理は、和と正負の反転しか使わない単純処理
さて、こうしてやると、平滑化を進め、その平滑化の程度の評価もできる
残るは、平滑化をやりすぎていないかどうかの判断をするためのルールを入れればよい
層を上げ過ぎると、細胞数が減っているわけなので、がたぼこ度の値も要素数が少ない分だけ小さくなりがち。したがって、単位要素数あたりのがたぼこ度として測ることにしてはどうだろうか。この処理は掛け算・割り算が入るので、ちょっと面倒くさいけれど、神経たちならそれなりにやってくれるレベルと思われる
さて、実際
- ヒストグラム

- density()関数によるカーネル推定プロット

- なんちゃって乳幼児プロット

# 1次元分布を作ろう
x <- c(rnorm(100,10),rnorm(50,80,20),runif(150,30,50))
hist(x)
plot(sort(x))

#網膜細胞を並べよう
X <- seq(from=min(x)-50,to=max(x)+50,length=100)
# パスカルの三角形的に細胞を配置して、細胞ごとに値を持たせる
M <- list()
M[[1]] <- rep(0,length(X)-1)
for(i in 1:length(M[[1]])){
	M[[1]][i] <- length(which(x>=X[i] & x<=X[i+1]))
}
# 層の値の総和で標準化
M[[1]] <- M[[1]]/sum(M[[1]])

# 層を上に上がって行こう
cnt <- 1
loop <- TRUE
while(loop){
	tmp <- length(M[[cnt]])
	M[[cnt+1]] <- M[[cnt]][1:(tmp-1)] + M[[cnt]][2:tmp]
	M[[cnt+1]] <- M[[cnt+1]]/sum(M[[cnt+1]])
	if(length(M[[cnt+1]])==1){
		loop <- FALSE
	}
	cnt <- cnt+1
}
# どんな風に平滑化したかを見てみよう
for(i in 1:length(M)){
	plot(M[[i]],type="l")
}

# 「差」を格納する
Ds <- list()
#DsPM <- Ds
for(i in 1:length(M)){
	Ds[[i]] <- list()
	#DsPM[[i]] <- list()
	Ds[[i]][[1]] <- M[[i]]
	#DsPM[[i]][[1]] <- M[[i]]
	if(length(M[[i]])>1){
		for(j in 2:length(M[[i]])){
			Ds[[i]][[j]] <- diff(Ds[[i]][[j-1]])
			#DsPM[[i]][[j]] <- sign(diff(sign(Ds[[i]][[j-1]])))
		}
	}
}
# 三角錐状のデータ格納になっている
# 以降の処理がしやすいように、値の納め方を変える
Ds.2 <- Ds
for(i in 1:length(Ds)){
	for(j in 1:length(Ds[[i]])){
		#DsPM.2[[j]][[i]] <- DsPM[[i]][[j]]
		Ds.2[[j]][[i]] <- Ds[[i]][[j]]
		#for(k in 1:length(DsPM[[i]][[j]])){
			#DsPM.2[[j]][[i]][k] <- DsPM[[i]][[j]][k]
			#Ds.2[[j]][[i]][k] <- Ds[[i]][[j]][k]
		#}
	}
}
# 平滑度の指標を計算する関数
my.smooth <- function(x){
	(sum(abs(x))-sum(x))
}

#abs.DsPM <- lapply(DsPM,lapply,abs)
#sum.abs.DsPM <- lapply(abs.DsPM,lapply,sum)
#sm.DsPM <- lapply(DsPM,lapply,my.smooth)
#abs.DsPM.2 <- lapply(DsPM.2,lapply,abs)
#sum.abs.DsPM.2 <- lapply(abs.DsPM.2,lapply,sum)
sm.Ds.2 <- lapply(Ds.2,lapply,my.smooth)

#par(ask=TRUE)

#for(i in 1:length(M)){
#	#plot(unlist(sum.abs.DsPM.2[[i]]))
#	plot(unlist(sm.Ds.2[[i]]))
#}
# ぎざぎざ量を観測点数で割る
# 単位要素当たりのぎざぎざがたぼこ度を細胞層ごとに計算してプロット

mean.giza <- unlist(sm.Ds.2[[2]])*(1/(length(unlist(sm.Ds.2[[2]])):1))
plot(mean.giza)

# がたぼこ度はだんだん下がって、あるところで反転する、そこを見つける
diff.mean.giza <- diff(mean.giza)
giza.pm <- sign(diff.mean.giza)
# 最初にがたぼこ度が増加に転じた層
first.posi <- which(giza.pm>0)[1]
# その一つ手前がボトム
kizami <- 1/length(M[[first.posi-1]])
xpt <- ppoints(length(M[[first.posi-1]]),a=0.5)
ylim <- range(M[[first.posi-1]])
plot(xpt,M[[first.posi-1]],main=first.posi-1,type="b",ylim=ylim)
first.posi-1

dev.new()
par(ask=TRUE)
for(i in 1:length(M)){
	kizami <- 1/length(M[[i]])

	xpt <- ppoints(length(M[[i]]),a=0.5)

	plot(xpt,M[[i]],main=i,type="b")
}
plot(density(x))


#plot(unlist(sm.Ds.2[[2]])*(1:length(unlist(sm.Ds.2[[2]]))))