2012-01-23

周辺度数を計算する

周辺度数は「特定の軸」についての足し算である
Rでは、arrayという形式のオブジェクトにapply()関数を使って、特定の軸の組に足し算(sum())を作用させることで得られる

k<-3
r<-sample(2:3,k,replace=TRUE)
print(r)
R<-prod(r)
print(R)
A<-array(1:R,r)
print(A)

apply(A,c(1),sum)
apply(A,c(2),sum)
apply(A,c(3),sum)
apply(A,c(1,2),sum)
apply(A,c(2,3),sum)
apply(A,c(3,1),sum)
apply(A,c(1,2,3),sum)

手計算との一致を確認するのがよい

> k<-3
> r<-sample(2:3,k,replace=TRUE)
> print(r)
[1] 2 3 2
> R<-prod(r)
> print(R)
[1] 12
> A<-array(1:R,r)
> print(A)
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

> 
> apply(A,c(1),sum)
[1] 36 42
> apply(A,c(2),sum)
[1] 18 26 34
> apply(A,c(3),sum)
[1] 21 57
> apply(A,c(1,2),sum)
     [,1] [,2] [,3]
[1,]    8   12   16
[2,]   10   14   18
> apply(A,c(2,3),sum)
     [,1] [,2]
[1,]    3   15
[2,]    7   19
[3,]   11   23
> apply(A,c(3,1),sum)
     [,1] [,2]
[1,]    9   12
[2,]   27   30
> apply(A,c(1,2,3),sum)
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

軸数kについて、周辺度数の取り方は $2^k-1$ ある。すべてのセルの和も周辺度数の一つとみなせば、取り方は $2^k$ ある
これは、集合 $\{1,2,...,k\}$ のべき集合の要素数であって、べき集合の要素が、周辺度数の取り方を規定しているとも言い換えられる
そのことを使って、任意の表のすべての周辺度数を計算する処理を以下のように作る

# 集合を扱うパッケージをインストールする
# install.packages("sets")
# 集合を扱うパッケージを読み込む
library(sets)
# 集合sを作る
s<-as.set(1:k)
# 見てみる
print(s)
# sのべき集合を作る
powset<-set_power(s)
# 見てみる
print(powset)

# 周辺度数を格納する容器を作る
marginals<-list()
# べき集合の要素数のループを回す
# ループの回数をカウントするための変数cnt
cnt<-1
for(i in powset){
	# べき集合の要素は集合なので、使い勝手が悪いから、整数のベクトルにする
	elems<-as.integer(i)
	# 見てみる
	print(elems)
	# 空集合の場合とそうでない場合とを分ける
	if(length(elems)==0){# 空集合の場合
		# 要素のリスト(空)と周辺度数=総和とを格納する
		marginals[[cnt]]<-list(elems=c(),marginals=sum(A))
	}else{# 空集合でない場合
		# 要素のリストと周辺度数とを格納する
		marginals[[cnt]]<-list(elems=elems,marginals=apply(A,elems,sum))
	}
	# ループ数のカウンタを一つ上げる
	cnt<-cnt+1
}
# できたものを見てみる
marginals

> library(sets)
> # 集合sを作る
> s<-as.set(1:k)
> # 見てみる
> print(s)
{1L, 2L, 3L}
> # sのべき集合を作る
> powset<-set_power(s)
> # 見てみる
> print(powset)
{{}, {1L}, {2L}, {3L}, {1L, 2L}, {1L, 3L}, {2L, 3L}, {1L, 2L, 3L}}
> 
> # 周辺度数を格納する容器を作る
> marginals<-list()
> # べき集合の要素数のループを回す
> # ループの回数をカウントするための変数cnt
> cnt<-1
> for(i in powset){
+ # べき集合の要素は集合なので、使い勝手が悪いから、整数のベクトルにする
+ elems<-as.integer(i)
+ # 見てみる
+ print(elems)
+ # 空集合の場合とそうでない場合とを分ける
+ if(length(elems)==0){# 空集合の場合
+ # 要素のリスト(空)と周辺度数=総和とを格納する
+ marginals[[cnt]]<-list(elems=c(),marginals=sum(A))
+ }else{# 空集合でない場合
+ # 要素のリストと周辺度数とを格納する
+ marginals[[cnt]]<-list(elems=elems,marginals=apply(A,elems,sum))
+ }
+ # ループ数のカウンタを一つ上げる
+ cnt<-cnt+1
+ }
integer(0)
[1] 1
[1] 2
[1] 3
[1] 1 2
[1] 1 3
[1] 2 3
[1] 1 2 3
> # できたものを見てみる
> marginals
[[1]]
[[1]]$elems
NULL

[[1]]$marginals
[1] 78


[[2]]
[[2]]$elems
[1] 1

[[2]]$marginals
[1] 36 42


[[3]]
[[3]]$elems
[1] 2

[[3]]$marginals
[1] 18 26 34


[[4]]
[[4]]$elems
[1] 3

[[4]]$marginals
[1] 21 57


[[5]]
[[5]]$elems
[1] 1 2

[[5]]$marginals
     [,1] [,2] [,3]
[1,]    8   12   16
[2,]   10   14   18


[[6]]
[[6]]$elems
[1] 1 3

[[6]]$marginals
     [,1] [,2]
[1,]    9   27
[2,]   12   30


[[7]]
[[7]]$elems
[1] 2 3

[[7]]$marginals
     [,1] [,2]
[1,]    3   15
[2,]    7   19
[3,]   11   23


[[8]]
[[8]]$elems
[1] 1 2 3

[[8]]$marginals
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

この処理を関数化する

# 上記の処理を関数にする
CalcAllMarginals<-function(Table,k){
	s<-as.set(1:k)
	powset<-set_power(s)
	marginals<-list()
	cnt<-1
	for(i in powset){
		elems<-as.integer(i)
		if(length(elems)==0){
			marginals[[cnt]]<-list(elems=c(),marginals=sum(Table))
		}else{
			#print(apply(Table,elems,sum))
			marginals[[cnt]]<-list(elems=elems,marginals=apply(Table,elems,sum))
		}
		cnt<-cnt+1
	}
	marginals
}

# 使ってみる
CalcAllMarginals(A,k)