2014-04-30

次に何が出る？　ディリクレ実践編

本当にできるかやってみよう
まずは、kカテゴリ観察されているときに、元集団のカテゴリ数がkである、と限定した場合
- この場合は、kカテゴリでのディリクレ分布そのものであり、ディリクレ分布に従うとき、n回のうち $a_i-1$ 回観察されたカテゴリの期待値は $\frac{a_i}{\sum_{i=1}^k a_i}$ であると知られているから、これで終わり。
- たとえば、 $n=4$ で $k=3$ で、観測回数が $(2,1,1)$ であるとすれば $a_1=3,a_2=2,a_3=2$ であるから、第1,2,3カテゴリの観測の期待値は $\frac{3}{3+2+2}=\frac{3}{7},\frac{2}{7},\frac{2}{7}$ となる
では、kとK=k+1の２カテゴリ数だけを考慮する場合はどうなるだろうか
- K=k+1という仮説とkという仮説との比は $R(K|k) = \frac{\Gamma(n+K)}{\Gamma(2n+K)}\frac{\Gamma(2n+k)}{\Gamma(n+k)}$ である
- kという仮説の下で観測総数が $n=\sum_{i=1}^k (a_i-1)$ であり、第i番カテゴリの観測回数が $a_i-1$ のとき、第i番カテゴリの観測確率の期待値は $\frac{a_i}{\sum_{i=1}^k a_i}$ 。もちろん、この仮説の下での、k+1番目のカテゴリの観測確率は0。
- 同様にKという仮説の下でのそれは $\frac{b_i}{\sum_{i=1}^{K=k+1} b_i}$ であって、 $i \le k$ ならば $\frac{a_i}{(\sum_{i=1}^k a_i)+1}$ 。 $i=k+1$ ならば $\frac{1}{(\sum_{i=1}^k a_i) +1}$
- たとえば、でで、観測回数がである例でやってみる
  - $a_1=3,a_2=2,a_3=2$ である
  - k=3,K=k+1=4カテゴリでの $R(K|k)$ は、昨日の関数を使って

my.beta <- function(a){
	exp(sum(lgamma(a))-lgamma(sum(a)))
}
#[tex:\frac{(\sum a_i)!}{\prod a_i!}]
my.comb <- function(a){
	exp(lgamma(sum(a)+1)-sum(lgamma(a+1)))
}
#[tex:C(\mathbf{a})\frac{B(2\mathbf{a}-1)}{B(\mathbf{a}}]
my.CBB <- function(a){
	my.like.dir(a)*exp(lgamma(sum(a)-length(a)+1)-sum(lgamma(a)))
}
#組合せ項は比を取る時は不要なので取り除いたもの
my.like.dir <- function(a){
	my.beta(2*a-1)/my.beta(a)
}
# 観測度数bを使って比を上記の比を計算
my.like.dir.obs <- function(b){
	if(length(b)==1){
		return(1)
	}
	a <- b+1
	my.like.dir(a)
}
# サンプル数n、観察カテゴリ数k、想定カテゴリ数Kの場合のk, K 頂点正単体空間の発生確率総積分
my.ratio.dir.integral <- function(n,k,K){
	exp(lgamma(n+K)-lgamma(2*n+K)-lgamma(n+k)+lgamma(2*n+k))
}

my.ratio.dir.integral(n=4,k=3,K=4)

> my.ratio.dir.integral(n=4,k=3,K=4)
[1] 0.6363636

- - 第1,2,3,4カテゴリのそれぞれの仮説での期待確率は

# サンプル数n、観察カテゴリ数k、想定カテゴリ数Kの場合のk, K 頂点正単体空間の発生確率総積分
my.ratio.dir.integral <- function(n,k,K){
	exp(lgamma(n+K)-lgamma(2*n+K)-lgamma(n+k)+lgamma(2*n+k))
}
# k種類観測されたときにK種類を想定した場合の期待値
# aは観測度数に1を加えたベクトル
my.expected.dirichlet <- function(a,K){
	b <- c(a,rep(1,K-length(a)))
	b/sum(b)
}
# vは観測回数のベクトル
my.expected.dirichlet.obs <- function(v,K){
	a <- v+1
	my.expected.dirichlet(a,K)
}
# k=3仮説での第1,2,3,4カテゴリ観測確率期待値
c(my.expected.dirichlet.obs(c(2,1,1),3),0)
# K=4仮説での第1,2,3,4カテゴリ観測確率期待値
my.expected.dirichlet.obs(c(2,1,1),4)

> c(my.expected.dirichlet.obs(c(2,1,1),3),0)
[1] 0.4285714 0.2857143 0.2857143 0.0000000
> my.expected.dirichlet.obs(c(2,1,1),4)
[1] 0.375 0.250 0.250 0.125

- - さて、仮説の事後確率比と併せて、「総合的な第i番観測確率期待値」を算出しよう
    - kとKの比より、kとKの相対値のベクトルが返る方がよいので関数を書いて…

# この比を出す代わりに、相対値をベクトルで…
my.relative.dir.integral <- function(n,K){
	exp(lgamma(n+K)-lgamma(2*n+K))
}
my.relative.dir.integral(4,c(3,4))

> my.relative.dir.integral(4,c(3,4))
[1] 0.0001984127 0.0001262626

- - これを相対割合(足して1)にして計算すれば

v <- c(2,1,1)
n <- sum(v)
k <- length(v)
K <- k+1
pr <- my.relative.dir.integral(n,c(k,K))
pr
pk <- c(my.expected.dirichlet.obs(v,k),rep(0,K-k))
pk
pK <- my.expected.dirichlet.obs(v,K)
pK
pr. <- pr/sum(pr)
pr.
pkK <- rbind(pk,pK)
pkK
tmp <- pkK * pr.
tmp
apply(tmp,2,sum)
sum(apply(tmp,2,sum))

> v <- c(2,1,1)
> n <- sum(v)
> k <- length(v)
> K <- k+1
> pr <- my.relative.dir.integral(n,c(k,K))
> pr
[1] 0.0001984127 0.0001262626
> pk <- c(my.expected.dirichlet.obs(v,k),rep(0,K-k))
> pk
[1] 0.4285714 0.2857143 0.2857143 0.0000000
> pK <- my.expected.dirichlet.obs(v,K)
> pK
[1] 0.375 0.250 0.250 0.125
> pr. <- pr/sum(pr)
> pr.
[1] 0.6111111 0.3888889
> pkK <- rbind(pk,pK)
> pkK
        [,1]      [,2]      [,3]  [,4]
pk 0.4285714 0.2857143 0.2857143 0.000
pK 0.3750000 0.2500000 0.2500000 0.125
> tmp <- pkK * pr.
> tmp
        [,1]       [,2]       [,3]       [,4]
pk 0.2619048 0.17460317 0.17460317 0.00000000
pK 0.1458333 0.09722222 0.09722222 0.04861111
> apply(tmp,2,sum)
[1] 0.40773810 0.27182540 0.27182540 0.04861111
> sum(apply(tmp,2,sum))
[1] 1

- じゃあ、k,K=k+1に限定せず、Kをどんどん増やそう。大きすぎるKは実質的に寄与割合が0に収束するので、適当な大きさのKまでを計算しよう

# v を観測したときにk=length(v)から、k,k+1,....とやる関数
# 最大Kと事前確率ベクトルはデフォルト値を与える
my.exp.dirichlet.multi <- function(v,maxK=length(v)+1000,pre=rep(1/(maxK-length(v)+1),maxK-length(v)+1)){
	# vは観測回数ベクトル。0より大の値
	k <- length(v)
	K <- k:maxK
	# 総観測数
	n <- sum(v)
	# 仮説の事後確率
	pr <- my.relative.dir.integral(n,K)
	pr. <- pr/sum(pr)
	pr.. <- pre * pr.
	pr.. <- pr../(sum(pr..))
	# 仮説別期待値
	pkK <- matrix(0,length(K),maxK)
	for(i in 1:length(K)){
		pkK[i,1:K[i]] <- my.expected.dirichlet.obs(v,K[i])
	}
	tmp <- pkK * pr..
	expected.prob <- apply(tmp,2,sum)
	return(list(v=v,n=n,k=k,maxK=maxK,K=K,pre=pre,post=pr..,pkK=pkK,pr.matrix=tmp,expected.prob=expected.prob))
}

v <- c(2,1,1)
tmp.out <- my.exp.dirichlet.multi(v)

- これを使って、k=3、観測回数(2,1,1)の場合を見てみよう
- まず、k=3,4,5,...の事後確率は？
  - すごく大きいKまでと、比較的小さいKの範囲でプロット。観察されたカテゴリだけ、というのは30％と見積もられており、Kの大きさに連れて確率は減衰しているのがわかる

par(mfcol=c(1,2))
plot(tmp.out$K,tmp.out$post)
plot(tmp.out$K[1:20],tmp.out$post[1:20])
par(mfcol=c(1,1))

- どのカテゴリが起きるかの期待値は？
  - 上と同じくとても大きなKまでと、小さ目のKでのプロット。2回観察されたカテゴリが最も確率が高く、3割強だが、これは2/(2+1+1)=0.5よりはずいぶん小さい。1回観察された２つのカテゴリは等確率でそれに続く。それ以外は「未観測カテゴリたち」である。

par(mfcol=c(1,2))
plot(1:tmp.out$maxK,tmp.out$expected.prob,type="h")
plot(1:20,tmp.out$expected.prob[1:20],type="h")
par(mfcol=c(1,1))

- 未観測カテゴリはたくさんの種類を想定しているので、それが「次」に出るときには「どれを見ても「初出」としか感じないので、併せて１カテゴリ〜初出カテゴリ、とすると、確かに第１カテゴリは高確率だが、第１、２，３の既出カテゴリにそれほど差はなく、初出カテゴリの確率は、第1カテゴリのそれより低いが、第２，３よりは高い、という結果になる

# 未観測カテゴリを合算する
plot(1:(tmp.out$k+1),tmp.out$expected.prob.obs.unobs,type="b",ylim=c(0,1))