2011-11-13

駆け足で読む『Lectures on Algebraic Statistics』再び２ 1. Markov Bases

目次はこちら
こちらやこちらも参考
独立な分割表の確率行列のランクは１

Nr<-sample(2:10,1)
Nc<-sample(2:10,1)
library(MCMCpack)

Pr<-rdirichlet(1,rep(1,Nr))
Pc<-rdirichlet(1,rep(1,Nc))

Pmat<-t(Pr) %*% Pc

Pmat

qr(Pmat)$rank

m次元分割表の正確確率検定について考える
- m次元ある
- 各次元に $r_i;i=1,2,...,m$ カテゴリある
- 複合カテゴリは $R=\pi_{i=1}^m r_i$ ある
- この $R$ 種類あるカテゴリの生起確率を $p=(p_i);i=1,2,...,R=\pi_{j=1}^m r_j$ とすれば、 $\sum_{i=1}^{R}p_i=1$ であって、これは、 $\Delta_{R-1}$ 正単体(頂点数が $R$ で $R-1$ 次元空間にある正単体)の内部として表される
- 今、m次元分割表があって、その周辺度数が与えられたとき、 $\Delta_{R-1}$ 上のどこの点なら、OKかには制約があるはず
- その制約を数学的に書き記すことを考える
  - 簡単な例から入る方がわかりやすいので、 $m=2,r_1=3,r_2=2$ の例で考える
  - $\begin{pmatrix}p_{1,1}&p_{1,2}&p_{1,3}\\p_{2,1}&p_{2,2}&p_{2,3}\end{pmatrix}$
  - 今、周辺が決まっていれば、 $p_{i,j}=s_i \times t_j$ なる $s,t$ がある
  - 行列 $p$ の代わりに、その対数をとったものを成分とする行列 $\log(p)=\begin{pmatrix}\log(p_{1,1})&\log(p_{1,2})&\log(p_{1,3})\\\log(p_{2,1})&\log(p_{2,2})&\log(p_{2,3})\end{pmatrix}$ を考える
  - $\log(s_i\times t_j)=\log(s_i)+\log(t_j)$ に気を付ければ
  - $\log(p)=\begin{pmatrix}\log(s_1)+\log(t_1)&\log(s_1)+\log(t_2)&\log(s_1)+\log(t_3)\\\log(s_2)+\log(t_1)&\log(s_2)+\log(t_2)&\log(s_2)+\log(t_3)\end{pmatrix}$
  - さらに $\log(p)=$
  - $\log(s_1)\begin{pmatrix}1&1&1\\0&0&0\end{pmatrix}+\log(s_2)\begin{pmatrix}0&0&0\\1&1&1\end{pmatrix}+\log(t_1)\begin{pmatrix}1&0&0\\1&0&0\end{pmatrix}+\log(t_2)\begin{pmatrix}0&1&0\\0&1&0\end{pmatrix}+\log(t_3)\begin{pmatrix}0&0&1\\0&0&1\end{pmatrix}$
  - と変形できる
  - ここで、0,1でできた行列は成分数が $R=\pi_{i=1}^m r_i$ である。これを、その長さのベクトルに見立てて考えれば
  - $A=\begin{pmatrix}1&1&1&0&0&0\\0&0&0&1&1&1\\1&0&0&1&0&0\\0&1&0&0&1&0\\0&0&1&0&0&1\end{pmatrix}$ と
  - $u=\begin{pmatrix}\log(s_1)&\log(s_2)&\log(t_1)&\log(t_2)&\log(t_3)\end{pmatrix}$ を使って
  - $\log(p)=A^T u$ と表すことができる( $\log(p)$ も長さ $R$ のベクトルにしてある)
  - 周辺度数を等しくするm次元分割表が２つあるとき、その表成分を長さ $R$ のベクトルにしてものは $A^T v = A u$ となる。このような表の集合は $u$ のfiberと呼ぶ
  - 対数が出てきたのは、確率の積を和の計算に変えて、行列計算に持ち込めるから(でもある)
- 途中からの説明は、 $m=2$ に限定していないので、 $m$ がは任意の説明として受け取れるだろう
- さて。実際の観測表が $u=(u_{i_1,i_2,....,i_R);n=\sum_{i=1}^R u_i$ であるとする。これと周辺度数を等しくする表は $A^T \alpha$ を等しくするものとする
- そのような観測をする生起確率は
- - これは、行列とベクトルの計算について、少し気を付けて式変形をすると得られる
- 正確な生起確率が計算できた
- これを、fiberの成分のすべての表の生起確率の和に対する比率にすると、 $\alpha$ の出てこない式が得られる
- 言い換えると、観察テーブルのセルの値と、fiberに含まれるテーブルのセルの値とだけで計算できることになる
- まず、fiberの成分のすべての表の生起確率の和を出そう
  - $P(AU=Au)=\sum_{v \in F(u)} \frac{n!}{\pi_{i=1}^R u_i!}e^{\alpha^T (A u)}$
  - $=n! e^{\alpha^T (A u)} \sum_{v \in F(u)} \frac{1}{\pi_{i=1}^R v_i!}$
  - ただし $F(u)$ はuのfiber
- であるから
- $P\frac{(U=u)}{P(AU=Au)}=\frac{1/(\pi_{i=1}^R u_i!)}{\sum_{v\in F(u)} 1/(\pi_{i=1}^R v_i!)}$
- ここで $A$ として0,1からなる特殊な行列を使ったが、すべてのセルが分割表的に平等な扱いになるような条件を満足すればよく、その条件は、行列Aのすべての列の和が等しいという条件(らしい)
ちょっとRメモ

n<-2
rs<-sample(2:4,n,replace=TRUE)
R<-prod(rs)
# 適当にAを作る
A<-matrix(0,n,R)
N<-100
for(i in 1:R){
	s<-sample(1:n,N,replace=TRUE)
	A[,i]<-c(tabulate(s,nbins=n))
}
A
# 周辺度数条件でAを作る
Arr<-array(0,rs)
address<-which(Arr==0,arr.ind=TRUE)
A<-matrix(0,sum(rs),R)
counter<-1
for(i in 1:n){
	for(j in 1:rs[i]){
		these<-which(address[,i]==j)
		
		A[counter,these]<-1
		counter<-counter+1
	}
	
}

A
g<-rnorm(length(A[,1]))
logp<-t(A)%*%g
p<-exp(logp)
p<-p/sum(p)
p

正確確率は、場合の数が大きすぎてできないが、モンテカルロで近似はできる
- その時に使うのが、マルコフ基底、メトロポリス-ヘイスティングサンプリング
- 周辺度数を守るようにしか動けないので、そのような動きを取りまとめたのがマルコフ基底
- ＯＫな動きのうち、どの動きをするかをランダムに決め、そのうえで、動くか動かないかに関して、位置の相対的な生起確率に基づいて確率的に決める動きが、メトロポリス-ヘイスティング。これをすることで、 $\Delta_{R-1}$ をくまなく、生起確率に応じてサンプリングすることができる
こちらに続く