2012-01-23

周辺度数の制約によって期待値を計算する

log-linear modelにおける期待値は次のように考える
k次元表の各セルは、そのセルが、第1,2,...,k番目の尺度において、それぞれ $c_1,c_2,...,c_k$ 番目のカテゴリであるというように指定することができる
これを $\mathbf{c}=(c_1,...,c_k)$ と表すことにする
他方、周辺度数制約の複体表現において、あるFacetを考えると、このFacetに対応する周辺度数が計算できる
表の任意のセル $\mathbf{c}$ は、任意のFacet,Fの周辺度数の中に、対応する値がある。これを $m^{(F)}(\mathbf{c})$ と書くこととする
このとき、log-linear modelでは各セルの期待値 $e(\mathbf{c})$ が $e(\mathbf{c}) \propto \prod_{F} m^{(F)}(\mathbf{c})$ を満足し、 $\sum_{\mathbf{c}} e(\mathbf{c})$ が表のサンプル総数になるように調整したもの、とする
これを関数にしてみよう

LogLinearExp<-function(Table,r,Facets){
	# すべての要素が1である多次元表を作る
	# 期待値表の素
	E<-array(1,r)
	# 表の番地格納した行列を作る
	address<-which(E==1,arr.ind=TRUE)
	# 集合として入力されるFacetsをベクトルのリストにするための格納庫
	FacetList<-list()
	# 各Facetの周辺度数の格納庫
	margList<-list()
	# 周辺度数の番地の格納庫
	margAddress<-list()
	# すべてのFacetsについてループを回す
	# ループ回数のカウンタ
	cnt<-1
	for(i in Facets){
		# Facetの集合を整数ベクトルにする
		FacetList[[cnt]]<-as.integer(i)
		#print(FacetList[[cnt]])
		# このFacetの周辺度数を計算する
		# Facetの要素の数が1であっても、行列として格納されるようにas.matrix()関数を用いている
		if(length(FacetList[[cnt]])==1){
			margList[[cnt]]<-as.matrix(apply(Table,FacetList[[cnt]],sum))
		}else{
			margList[[cnt]]<-as.array(apply(Table,FacetList[[cnt]],sum))
		}
		
		# 周辺度数行列の番地を格納
		margAddress[[cnt]]<-which(margList[[cnt]]<Inf,arr.ind=TRUE)
		#print(margList[[cnt]])
		#print(margAddress[[cnt]])
		# ループ回数を一つ増やす
		cnt<-cnt+1
	}
	# 期待値表のすべてのセルを順次処理するためのループ
	for(i in 1:length(address[,1])){
		#print(address[i,])
		# すべてのFacetに関する処理のループ
		for(j in 1:length(FacetList)){
			# このFacetが集計している軸について、セルの番地を取り出す
			short.address<-address[i,FacetList[[j]]]
			# 軸数が１のときにも行・列の値が必要なのでそのための処理をする
			if(length(short.address)==1)short.address<-c(short.address,1)
			#print(short.address)
			# Facetの周辺度数のうち、該当するものがどれかを調べる
			diff<-apply((t(margAddress[[j]])-short.address)^2,2,sum)
			#print(diff)
			marg.address<-which(diff==0)
			#print(apply(margAddress[[j]],1,identical,short.address))
			#marg.address<-which(apply(margAddress[[j]],1,identical,short.address))
			#marg.address<-which(identical(margAddress[[j]][,1:length(margAddress[[j]][1,])],short.address))
			#print(marg.address)
			#print(margList[[j]][marg.address])
			# すべてのFacetについて積を取る
			E[i]<-E[i]*margList[[j]][marg.address]
		}
	}
	# 期待値表の総数が、元の表の総数と一致するように補正
	#print(sum(Table))
	E<-E/sum(E)*sum(Table)
	# 返す
	E
}

使ってみよう

# 表の次元を与える
k<-3
# カテゴリ数を適当に取る
r<-sample(2:4,k,replace=TRUE)
R<-prod(r)
# きれいな表を作る
A<-array(1:R,r)
# Facetsを作る
s<-as.set(1:k)
powset<-set_power(s)
# Facesの候補をpowsetから適当に選ぶ
# 大きいFaceが選ばれても、その部分集合のすべてが選ばれるようにはしていないので、ここで作っているFacesはFaceのリストにはなっていないことに注意Faces<-list()
cnt<-1
thres<-0.1
Faces<-list()
for(i in powset){
	if(runif(1)<thres){
		Faces[[cnt]]<-i
		cnt<-cnt+1
	}
}
Facets<-MakeFacets(Faces)
print(Facets)
# log-linear modelでの期待値表を作る
A.e<-LogLinearExp(A,r,Facets)
# AとA.eの全周辺度数を計算する
# Aの周辺度数
m.A<-CalcAllMarginals(A,k)
# Aの期待値表の周辺度数
m.A.e<-CalcAllMarginals(A.e,k)
# すべての周辺度数について、Facetに含まれるかを示しつつ
# Aの周辺度数とA.eの周辺度数をプリントアウトする
for(i in 1:length(m.A)){
	in.facet<-FALSE
	for(j in Facets){
		if(as.set(m.A[[i]]$elems) <= j){
			in.facet<-TRUE
		}
	}
	if(in.facet){
		print("in facet")
	}else{
		print("not in facet")
	}
	print(m.A[[i]])
	print(m.A.e[[i]])
}