2013-06-03

ヒストグラムと情報量推定

確率密度分布情報量 R

昨日の記事で分割表の生起確率と情報量との関係、２つの変数が生起するときの確率密度分布の情報量と２つの変数の非独立性の話を書いた
じゃあ、多変量の確率密度分布があると、その情報量は、構成変数間(変数のセットの間でもよい)の非独立性の指標になることになるわけで
たいていの場合、標本データからその情報量を推定してやりたくなる
Wiki記事
このWiki記事にもあるように、大雑把にヒストグラムを取ってそれに基づいて情報量推定する方法と、密度分布関数を推定してそれに基づいて情報量推定する話などがある
ヒストグラムを使う方法の場合、ビンの刻みを変えると推定情報量も変わる
すべての点がばらばらのビンに入れば $N \times \frac{1}{N}\log{\frac{1}{N}}=-\log{N}$ である
すべての点が１つのビンに入れば $1 \times 1 \log{1}=0$ である
今、一様分布だとしてビンをk個に分けるとすると、 $k \times \frac{1}{k} \log{\frac{1}{k}}=-\log{k}$ となる
一様分布からのずれを問題にするなら、分けるビン数に関して標準化するとよいだろう
逆に、いくつに分けるかに迷ったら、標準化後の情報量が大きい状態がある意味でよいヒストグラムとも言えそうだ
そんなことをRでやってみる(多次元)

# データの座標によらず、視野内に収めるための補正
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))
}
# 視野全体のピクセル数を指定して、該当する超立方体の番地にする
round.coords <- function(x,n.cell=50){
	d <- length(x[1,])
	matrix(floor(x * n.cell),ncol=d)
}

mult.dim.hist <- function(x,n.cells=50,scl=0.7){
	d <- length(x[1,])
	x.st <- vision.standard(x,scl=scl)
	x <- x.st$st.x

	#Ls <- 1/(c(16,8,4,2))
	d <- length(x[1,])
	round.x <- count.x <- membership.x <- cluster.x <- represent.x <- infos <- st.infos <- list()
	dist.mat.x <- list()
	trees.out <- list()
	dist.mat.x[[1]] <- NULL
	round.x[[1]] <- x
	N <- length(x[,1])
	membership.x[[1]] <- 1:N
	count.x[[1]] <- rep(1,N)
	cluster.x[[1]] <- 1:N
	represent.x[[1]] <- 1:N
	infos[[1]] <- -log(N)
	st.infos[[1]] <- infos[[1]]/(-d*log(min(diff(sort(x)))))

	for(i in 1:length(n.cells)){
		if(length(matrix(round.x[[i]],ncol=d)[,1])>1){
			if(i==1){
				tmp <- round.coords(x,n.cells[i])
			}else{
				tmp <- round.coords(round.x[[i]]/n.cells[i-1],n.cells[i])
			}
			
			n <- length(tmp[,1])
			membership.x[[i+1]] <- rep(0,n)
			count.x[[i+1]] <- vector()
			represent.x[[i+1]] <- vector()
			tobechecked <- rep(1,n)
			while(sum(tobechecked)>0){
				current.checker <- which(tobechecked==1)[1]
				current.tobechecked <- which(tobechecked==1)[-1]
				#if(length(current.tobechecked)>0){
					diffs <- matrix(t(tmp[current.tobechecked,])-tmp[current.checker,],nrow=d)
					diff.all <- apply(diffs,2,sum)
					zeros <- which(diff.all==0)
					membership.x[[i+1]][c(current.checker,current.tobechecked[zeros])] <- max(membership.x[[i+1]])+1
					
					tobechecked[c(current.checker,current.tobechecked[zeros])] <- 0
					represent.x[[i+1]] <- c(represent.x[[i+1]],current.checker)
					count.x[[i+1]] <- c(count.x[[i+1]],sum(count.x[[i]][c(current.checker,current.tobechecked[zeros])]))
			}
			round.x[[i+1]] <- matrix(tmp[represent.x[[i+1]],],ncol=d)
			dist.mat.x[[i+1]] <- as.matrix(dist(round.x[[i+1]],method="manhattan"))
			#trees.out[[i+1]] <- one.direction.2(dist.mat.x[[i+1]],round.x[[i+1]])
			tmp.prob <- count.x[[i+1]]/sum(count.x[[i+1]])
			infos[[i+1]] <- sum(tmp.prob*log(tmp.prob))
			st.infos[[i+1]] <- infos[[i+1]]/(-log(n.cells[i])*d)
		}else{
			break
		}
	}
	return(list(round.x = round.x,count.x = count.x, membership.x = membership.x, represent.x = represent.x, infos = infos, st.infos = st.infos ))
}

x <- runif(100)
x <- matrix(x,ncol=1)
n.cells <- 50:10

d <- 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.st <- vision.standard(x,scl=0.7)
x <- x.st$st.x
n <- length(x[,1])

m.out <- mult.dim.hist(x,n.cells)

plot(unlist(m.out$st.infos))