2013-04-21

全列挙

分割表

１番じゃなきゃシリーズでは、モンテカルロな評価と、正確確率的な評価の両方が面白い
正確確率的な評価には起こりうるすべての表をすべて挙げたい
こちらでそれをうまくやるために図形数とかも考えてみたけれど、図形数として扱いやすいのはユークリッド座標系になっている場合なので、こちらでは、k次元立方格子上に総標本数0,1,...,Nの分割表を対応付け、この立方格子上の変化として扱うことにする
この立方格子の図形数として、総数Nの分割表の総数なども現れる

# k次元、総標本数Nのための立方格子を作る
N <- 4
k <- 8
X <- array(0,rep(N+1,k))
# このアレイXはR的にはベクトルなので、アレイのk次元座標をアレイXのベクトル座標としての対応座標を取りたい
# その逆もハンドリングしたい
num.cds <- function(v,N,k){
	sum(v*(N+1)^(0:(k-1)))+1
}
num.cds.inv <- function(n,N,k){
	tmp.n <- n-1
	ret <- rep(0,k)
	for(i in 1:k){
		ret[i] <- tmp.n %% (N+1)
		tmp.n <- (tmp.n - ret[i])/(N+1)
	}
	ret
}
cds <- num.cds(v,N,k)
v <- sample(0:N,k,replace=TRUE)
cds <- num.cds(v,N,k)
v2 <- num.cds.inv(cds,N,k)
v-v2

# セル数k個のときに総標本数が1増えるとき、増え方はk通り
# 長さkの確率ベクトルが成長を定義するのでそれをディリクレ分布で適当にとる
N <- 10
k <- 8
library(MCMCpack)
r <- c(rdirichlet(1,rep(1,k)))
# 初期状態は、総数0の場合に確定しているとする
# そのときのk次元ベクトルは１つ
current.v <- matrix(rep(0,k),nrow=1)
# そのアレイXをベクトルとみなしたときの座標は１つ
current <- num.cds(current.v,N,k)
X <- array(0,rep(N+1,k))
# その番地の確率が1
X[current] <- 1
# N=0,1,..の状態を格納
Y.v <- list()
Y.p <- list()
# 表のk個の値に対応するベクトルをためる
Y.v[[1]] <- current.v
# 対応する確率をためる
Y.p[[1]] <- X[current]
# ステップごとに推移させる
for(i in 1:N){
	unique.next <- c()
	unique.next.v <- matrix(NA,nrow=0,k)
	for(j in 1:k){
		ippo <- rep(0,k)
		ippo[j] <- 1
		next.v <- t(t(current.v) + ippo)
		#print(next.v)
		next.num <- rep(NA,length(current))
		for(j2 in 1:length(next.num)){
			next.num[j2] <- num.cds(next.v[j2,],N,k)
		}
		tmp.r <- r[j]
		#print(current)
		#print(next.num)
		#print(tmp.r)
		X[next.num] <- X[next.num] + X[current] * tmp.r
		unique.next <- unique(c(unique.next,next.num))
	}
	current <- unique.next
	#print(unique.next)
	print(sum(X[current]))
	current.v <- matrix(0,length(current),k)
	for(j in 1:length(current)){
		#print(current[j])
		current.v[j,] <- num.cds.inv(current[j],N,k)
	}
	Y.v[[i+1]] <- current.v
	Y.p[[i+1]] <- X[current]
}

結局これは、デカルト座標を利用して、継承関係を記録すればよいので(そのうえで継承関係に応じた生起確率を作ればよいので)、継承関係のみをアレイXをベクトルとみたときの要素番号をノードIDとした継承関係を取り出すことにする

N <- 4
k <- 8
#r <- c(rdirichlet(1,rep(1,k)))
TableHistory <- function(N,k){
	current.v <- matrix(rep(0,k),nrow=1)
	current <- num.cds(current.v,N,k)
	Y.v <- list()
	Y.parent <- Y.parent.id <- list()
	Y.v[[1]] <- current.v
	Y.parent[[1]] <- current
	Y.parent.id[[1]] <- 1
	Z <- Z.id <- list()
	for(i in 1:N){
		unique.next <- c()
		unique.next.v <- matrix(NA,nrow=0,k)
		Z[[i]] <- Z.id[[i]] <- matrix(NA,length(current),k)
		for(j in 1:k){
			ippo <- rep(0,k)
			ippo[j] <- 1
			next.v <- t(t(current.v) + ippo)
			next.num <- rep(NA,length(current))
			for(j2 in 1:length(next.num)){
				next.num[j2] <- num.cds(next.v[j2,],N,k)
			}
			unique.next <- unique(c(unique.next,next.num))
			Z[[i]][,j] <- next.num
		}
		current <- unique.next
		current.v <- matrix(0,length(current),k)
		for(j in 1:length(current)){
			current.v[j,] <- num.cds.inv(current[j],N,k)
		}
		Y.v[[i+1]] <- current.v
		Y.parent[[i+1]] <- current
		Y.parent.id[[i+1]] <- 1:length(current)
	}
	for(i in 1:length(Z.id)){
		for(j in 1:length(Y.parent.id[[i+1]])){
			Z.id[[i]][which(Z[[i]] == Y.parent[[i+1]][j])] <- j
		}
	}
	return(list(N=N,k=k,coords=Y.v,parent=Y.parent,parent.id = Y.parent.id,offspring=Z,offspring.id = Z.id))
}
N <- 4
k <- 8
th.out <- TableHistory(N,k)

# A function for Decision_Dice_beta.2.3() :Decision rule 5
Decision_beta.2 <- function(x){
	if(x[1] < x[4]){
		x <- x[4:1]
	}
	ret <- 0
	common <- -log(x[3]+x[4])+lgamma(x[1]+x[2])+lgamma(x[3]+x[4]+1)-lgamma(x[1])-lgamma(x[2])-lgamma(sum(x)-1)
	for(j in x[3]:(x[3]+x[4]-1)){
		tmp <- lgamma(x[1]+j)+lgamma(sum(x[2:4])-1-j)-lgamma(j+1)-lgamma(x[3]+x[4]-j)
		ret <- ret + exp(tmp+common)
	}
	return(ret)
}


# th.out$Nは最大標本数、th.out$kは次元数、th.out$coordsはテーブルのセルの値のベクトル、th.out$parentは継承元のID番号、th.out$offspringは継承元からk方向別の継承先のID番号

p1 <- 0.5
pA <- c(0.6,0.3)
pB <- c(0.6,0.3)

p.history <- list()
p.history[[1]] <- 1
for(i in 2:(N)){
	p.history[[i]] <- rep(0,length(th.out$parent.id[[i]]))
	tmp.beta.prob.1 <- apply(matrix(th.out$coords[[i-1]][,1:4]+1,ncol=4),1,Decision_beta.2)
	tmp.beta.prob.2 <- apply(matrix(th.out$coords[[i-1]][,5:8]+1,ncol=4),1,Decision_beta.2)
	r1 <- tmp.beta.prob.1 * p1 * pA[1]
	r2 <- tmp.beta.prob.1 * p1 * (1-pA[1])
	r3 <- (1-tmp.beta.prob.1) * p1 * pA[2]
	r4 <- (1-tmp.beta.prob.1) * p1 * (1-pA[2])
	r5 <- tmp.beta.prob.2 * (1-p1) * pB[1]
	r6 <- tmp.beta.prob.2 * (1-p1) * (1-pB[1])
	r7 <- (1-tmp.beta.prob.2) * (1-p1) * pB[2]
	r8 <- (1-tmp.beta.prob.2) * (1-p1) * (1-pB[2])
	rs <- cbind(r1,r2,r3,r4,r5,r6,r7,r8)
	for(j in 1:length(rs[1,])){
		print("--")
		print(p.history[[i]][th.out$offspring.id[[i-1]][,j]])
		print(rs[,j])
		print("##")
		p.history[[i]][th.out$offspring.id[[i-1]][,j]] <- p.history[[i]][th.out$offspring.id[[i-1]][,j]] + rs[,j]*p.history[[i-1]][th.out$parent.id[[i-1]]]
	}
}
lapply(p.history,sum)