2011-12-11

分割表の位置・自由度・座標変換のまとめ　分割表の正単体・複体座標表現によるカイ二乗検定

分割表外積代数複体代数統計 R

こちらの整理記事
$k$ 次元の表 $t$ がある
$i-$ th( $i=1,2,...,k$ )次元のカテゴリ数を $\mathbf{r}=\{r_i\}$ とする
$t$ のセルの総数は $R=\prod_{i=1}^k r_i$ である
には、複体で表される周辺度数制約がある
- ここで $S_j$ は(正)単体であり、 $S_j$ の頂点は $\{1,2,...,k\}$ の部分集合である
- $\mathbf{S}$ による周辺度数制約とは、 $S_j$ が作る、 $T$ の部分表に相当する周辺度数が与えられていることを意味する
長さ $R$ のベクトル $a$ は表 $t$ のセルの値とする
今、の構成であるが与えられているとき、を同長のベクトルに１対１対応で移す行列()を次のような条件を満足するようにとることができる
- ２つの表 $t_1,t_2$ があり、それらは $S_x$ で表される部分表に相当する周辺度数が等しいとすると、この２表に対応する $a_1,a_2,b_1,b_2$ は $bB_1=Xa_1,b_2=Xa_2$ であって、長さ $R$ のベクトル $b_1,b_2$ の値は、部分表のセルの数 $R(S_x)$ だけ一致する
ここで、 $k,\mathbf{r},\mathbf{S}$ が与えられ、 $\mathbf{S}$ に対応して周辺度数表 $\mathbf{\tau(\mathbf{S})}$ が定まるとき、その条件に合致する表の集合 $\mathbf{T}=\{t|k,\mathbf{r},\mathbf{S},\mathbf{\tau(\mathbf{S})}\}$ が得られる
$t_i,t_j \in \mathbf{T}$ を考え、その差 $d_{i,j}=t_i-t_j$ を取る
対応するベクトル $a_i,a_j,\alpha_{i,j}=a_i-a_j$ , $b_i,b_j,\beta_{i,j}$ は $b_i=Xa_i,b_j=Xa_j,\beta_{i,j}=X\alpha_{i,j}$ の関係にあり、 $b_i,b_j$ は制約 $\mathbf{S}$ が定める要素数が等しいから、 $\beta_{i,j}$ は制約 $\mathbf{S}$ が定める要素数だけ、値が0であるようなベクトルである
このように $t_i,t_j\in \mathbf{T}$ であるとき、 $\beta_{i,j}$ の要素のうち必ず0となるような要素の数を $g$ とすれば、 $f=R-g$ は、 $\mathbf{T}$ の自由度である
$\beta_{i,j}=X\alpha_{i,j},\alpha_{i,j}=X^{-1}\beta_{i,j}$ であるときに
長さ $R$ の正数からなるベクトルと $e\in \mathbf{T}$ を定め、その成分の平方根の逆数を対角成分とする $R\times R$ 対角行列 $E$ を定め、スカラー量 $K^2=(E\alpha_{i,e})^{T}(E\alpha_{i,e})$ を定義する
$\alpha_{i,e}=X^{-1}\beta_{i,e}$ であるから、 $K^2=(EY\beta_{i,e})^{T}(EY\beta_{i,e})$ である。ただし $Y=X^{-1}$
ここで $\forall i$ について $\beta_{i,e}$ の $g$ 個の成分は $0$ であるから、その成分を除去した長さ $f$ のベクトルを $\beta_{i,e}'$ とし、 $Y$ からそれに対応する列ベクトルのみを抽出した $R\times f$ 行列を $Y'$ とすると、 $K^2=(EY'\beta_{i,e}')^{T}(EY'\beta_{i,e}')$
これを書き換えて $K^2=\beta_{i,e}'^{T} (Y'^{T}E^TEY')\beta_{i,e}'$ となるが、 $E^T=E$ であることに注意すれば、 $W=Y'^{T}E^2Y'$ なる $f \times f$ 行列を用いて $K^2=\beta_{i,e}'^{T} W)\beta_{i,e}'$ と表せることがわかる
ここで $W=V\Sigma \Sigma V^{-1}$ と固有値分解することで $K^2=(\Sigma V^{-1}\beta_{i,e}')^T(\Sigma V^{-1}\beta_{i,e}')$ となり、これは長さ $f$ のベクトル $Q=\Sigma V^{-1}\beta_{i,e}'$ のノルムの二乗である
$K^2$ の定義に戻ると、 $e$ が期待値表に対応していると見れば、これは期待値からの逸脱に関するカイ二乗統計量であり、周辺度数制約 $\mathbf{S}$ を満足する表 $t_i$ のセルの値のベクトル $a_i$ の線形代数変換によって得られた自由度の長さのベクトル $Q$ のノルムの二乗として表せることが示された
言い換えると、周辺度数制約 $\mathbf{S}$ を満足する表 $t \in \mathbf{T}$ のうち、 $K^2$ を等しくするものを、自由度 $f$ 次元空間の球表面へと写す座標変換ができたことになる
実行

###############xxx###########
NFacets<-10
Niter<-10
S1<-S2<-S3<-rep(0,NFacets*Niter)
cnt<-1
for(ii in 1:NFacets){


	# 表の次元・次数の設定
	Nv<-5
	Rs<-sample(2:3,Nv,replace=TRUE)
	Vs<-1:Nv
	ns<-Rs
	# Facesをランダムに作る
	Nf<-3

	Faces<-list()
	maxN<-3
	Ns<-sample(1:maxN,Nf,replace=TRUE)
	for(i in 1:Nf){
		tmpVs<-sample(Vs,Ns[i])
		Faces[[i]]<-as.set(tmpVs)
	}
	# Facetsにする
	Facets<-MakeFacets(Faces)
	# グラフにしてみる

	g<-GraphFacets(Nv,Facets,plot=FALSE)
	##############
	#Facets<-set(1:Nv)
	# 観測表をランダムに作る
	Obs<-array(runif(prod(ns)),ns)
	Obs<-Obs/sum(Obs)

	# 期待値表を作る
m.e.out<-MakeExpected4(Obs,ns,Facets)
Etable<-m.e.out$Etable
X<-m.e.out$X
Z<-m.e.out$Z
#Etable-Etable2
# ObsとFacetsとから、期待値表を得よう

v2<-c(Obs-Etable)
#out.makeX2<-MakeX2(ns,Facets)
#X<-out.makeX2$matrixX
#Z<-out.makeX2$zeros
	P<-X%*%v2
	print(Z)
	Xinv<-solve(X)
	XinvPartial<-Xinv[,which(Z==1)]
	XinvPartialt<-t(Xinv)[which(Z==1),]
	W<-XinvPartialt%*%diag(1/c(Etable))%*%XinvPartial
	eigen.out<-eigen(W)
	V<-eigen.out$vectors
	S<-diag(eigen.out$values)
	
	E<-diag(1/sqrt(c(Etable)))
	Einv<-diag(sqrt(c(Etable)))
	Q<-E%*%Xinv%*%P
	Qz<-diag(sqrt(eigen.out$values))%*%solve(V)%*%P[which(Z==1)]
	S1[cnt]<-sum(Q^2)
	S2[cnt]<-sum(Qz^2)
	S3[cnt]<-sum((v2^2)/c(Etable))
	cnt<-cnt+1
	for(jj in 2:Niter){
		QzRandom<-rnorm(length(Qz))
		QzRandom<-QzRandom/sqrt(sum(QzRandom^2))
		QzRandom<-QzRandom*sqrt(sum(Qz^2))
		PzRandom<-V%*%diag(1/sqrt(eigen.out$values))%*%QzRandom
		PRandom<-rep(0,length(P))
		PRandom[which(Z==1)]<-PzRandom
		QRandom2<-E%*%Xinv%*%PRandom
		DRandom<-Xinv%*%PRandom
		S1[cnt]<-sum(QRandom2^2)
		S2[cnt]<-sum(QzRandom^2)
		S3[cnt]<-sum(DRandom^2/c(Etable))
		cnt<-cnt+1
	}
}

plot(as.data.frame(cbind(S1,S2,S3)))