2018-03-21

メモ

ごちゃごちゃと考えたときのメモ
- 普通のポアソン分布は非負整数を台としてその上に非負実数を配するルールをパラメタ依存に定める。そういう意味で、パラメトリックな確率分布とは、「ある台」に総和(全体の積分)が１になるような確率質量・密度を定めるルールをパラメタ依存に持っているような「確率質量・密度分布というインスタンスの集合」を定めるもの。パラメタの値を特定すれば、それはあるパラメトリックな確率分布の集合の要素を指定することになる。このようにある台の上に確率質量・密度分布を特定するlaw(決まりごと・関数)をdistributionという
- 二項分布は、サンプル数N(と生起確率(p,1-p)と)を定めることで分布インスタンスを特定する。そういう意味で「N個サンプリングする」というときのNもパラメタのひとつ
- では、Poisson-Dirichlet Distributionsはどうなるのか？
  - 台はexchangeable partitionsである。ある自然数(非負整数)n個を取り出したとして、そのタイプごとの内訳が $\sum_{j=1}^n j m_j =n$ となるような $(m_1,...,m_n)$ となるという事象全体がexchangeable partitions
  - 標本を無限にとってもよいが、有限個をサンプルすることもあって、それは二項分布のNが分布インスタンスを特定するのに必要なのと同様
- 節4まででExchangeable partitionsについて説明を加えてきたが、それはPoisson-Dirichlet distributionsの台を説明するためだったことがわかる
- ではPoisson-Dirichlet Distributionsの名前の由来は何だろうか？(これは(今のところ)想像)
  
  自然数の分割だけれど、それを割合的に捉え、総数がNのところを総数が1にしている点でDirichlet
  
  また、自然数の分割であって、その総数が無限大まで広がっていることから、離散確率分布でもある。離散確率分布と言えばポアソン分布があり、そういう意味でPoissonという単語が使われるし、個々のタイプに観測される個数・回数がポアソン分布で説明することも多いという意味でPoissonという単語を使っているのであろうか
- さて。どういうDistributionsがあるだろうか
  - 同一の分布(たとえば正規分布)でもいくつかの異なる表式がある(平均と分散を用いたパラメタ表現や指数型分属表現や…)し
  - 異なる分布と言えるが台は同じ分布もある(ポアソン分布が負の二項分布の特殊形であるというような関係も含めて)しするので
  - いろいろな、分布集合を表すルール(law)があり、それぞれをPoisson-Dirichlet distributionと呼び、別名をつけたりする。また、それぞれがパラメタで表現されていればパラメトリックな分布と言えるだろう
- ではこの節は、上記のような意味での「いろいろなPoisson-Dirichlet distributionS」をどのように説明しているか、と言うと…
- s-Paintbox
  - 降順に実数列を作りその和が1未満となるようにする。1未満の総和の最後の残りをSingletons用とする
  - これにより一意に定まった実数列が得られる
  - ここに(0,1)の数値を発生(乱数を発生させればExchangeable "random" partitionsが得られる)させ、Singletons用の場所ならシングルトンラベルを、そうでなければ、対応セグメントのラベルをつける。これにより任意の自然数個の要素がラベル付けされる
  - これを基本的なラベル付け手順とする
- s-Paitboxなラベル付け手順のもとではどのように降順実数列(足して1未満)を発生させるかの規則(law)がs-Paintboxの作られ方を決める(それがとりもなおさずdxchangeable partitionsを生む)ので、降順実数列生成lawが(Poisson-Dirichlet) distributionを表す
  - Ewens distribution $P(\Pi_n=\{b_1,b_2,...,b_k\}) = \frac{\theta^k}{\theta(\theta+1)...(\theta+n-1)}\prod_{j=1}^k(b_j-1)!$ は、あるnに対し、 $(b_1,...)$ ごとに、 $\theta$ 依存な確率質量分布が得られた。したがって、取りうるすべての $(b_1,...)$ を台(Exchangeable partitions全体の一部)とした離散確率質量分布を定めるlawとしてのdistributionがEwens distributionであるということになる。このとき $\theta$ が単一のパラメタなので $\theta$ をパラメタとするPoisson-Dirichlet distributionと呼ばれることもある
- ディリクレ・多項分布の次元を無限大にしたときの分割がEwens's sampling formula
5. Two-paramter Ewens-Pitam Distributionとしての考え方

2018-03-21

ぱらぱらめくる『The Ubiquitous Ewens Sampling Formula』

ペイパーはこちら
Ewens sampling formulaは遺伝統計学分野から出た正確確率計算式で、ある理想集団におけるアレル頻度パターンが生じる確率の式であるが、それは離散確率過程において応用範囲が広いものであり、遺伝学・生物学の中では、アレル頻度だけでなく中立仮説の下での生物多様性(種数推定)にも使えるし、確率・統計学一般にはノンパラメトリックベイズ推定、組みあわせ確率過程、帰納的に推定するとはどういうことかといった裾野の広いトピックと関連する。また、より数学的に根源的な自然数の素数分解などにも通じている
このペイパーではEwens's sampling formulaが遺伝学分野から登場してきた経緯から始めて、それの応用展開について概説している
1. Introduction
- Ewens's sampling formula
- n本の染色体を観察したとき、観察しうるアレルは1種類からn種類までである。また、アレル別に観察しうる染色体本数は1からnまでである
- 観測本数がj = 1,2,...,nであるアレル数を $m_j$ とすると、 $n = \sum_{j=1}^n j m_j$ が成り立つ
- このとき、ある一定の条件を満足する集団においてn本の染色体の(n/2人：常染色体の場合)遺伝子座位を観察すると $(m_1,...,m_n)$ の生起確率は、あるパラメタ $\theta$ を用いて以下のようになることが知られている。ただし、 $\theta$ は有効集団サイズと変異率との積に比例した値である
- $P(m_1,...,m_n|\theta) = \frac{n!}{\theta(\theta+1)...(\theta+n-1)}\prod_{j=1}^n \frac{\theta^{m_j}}{j^{m_j}m_j!}$
- この関数をRで作成し、確かに $m_1,...$ の場合わけごとの確率の総和が１になることを確かめてみる

my.Ewens.prob <- function(ms,theta,log=FALSE){
	n <- sum(1:length(ms)*ms)
	if(length(ms) < n){
		ms <- c(ms,rep(0,n-length(ms)))
	}
	tmp <- lfactorial(n) -sum(log((0:(n-1))+theta)) + sum(ms * log(theta))-sum(ms*log(1:n))-sum(lfactorial(ms))
	if(log){
		return(tmp)
	}else{
		return(exp(tmp))
	}
}

my.Ewens.prob(c(0,1),1)
my.Ewens.prob(c(2,0),1)

library(partitions)
n <- 6
theta <- 0.1
prts <- parts(n)
prbs <- rep(0,length(prts[1,]))
for(i in 1:length(prts[1,])){
	prt <- prts[,i]
	print(prt)
	tab <- tabulate(prt)
	prbs[i] <- my.Ewens.prob(tab,theta)
}

sum(prbs)

2. 中立モデル(遺伝)
- Wright-Fisher model
  - 人口一定、ランダムメイティング、無限サイト、復元抽出
  - 遺伝子のモデルだが、確率モデルとしてはかなりがっちりとシンプルに徹しているモデル
- 遺伝学的使い道
  - 観察アレル種類数の確率分布も出せる(観察アレル数というのは、n=4のときに2本ずつ2種類、3本と1本との2種類というように、アレル種類数が同じならプールすること)
  - それがわかると、サイズNの母集団が有するアレル種類数と、サイズnの標本集団が有するアレル種類数についても分布がわかり、標本に観察されない未観測アレル種類数などについても予想がつけられる
  - Ewens's formulaからベタに計算することもできるし、第１種スターリング数を使った定義式 $P(K=k) = |S^k_n| \frac{\theta^k}{\theta(\theta+1)...(\theta+n-1)}$ でも計算できる
  - また、このように確率モデルのパラメタ $\theta$ が観察条件 n と観察項目 K(アレルの種類数)とで表されているということは、Kを観察することが $\theta$ を推定するための十分統計量となることを意味しており、それは、種類数の確率分布が指数型分布族形式で表され、分布を定めるパラメタ $\theta$ と観察状態とが絡む項に、観察項目数のみが登場することと関係している(本ペイパーの3.7で示されている)
  - ちなみに。ここで第１種スターリング数が出てくるのは、スターリング数がある場合の数を表しているからであり、その「ある場合の数」とは、『アレル種類数がk』であるような数字の組の場合の数であることが理由である。
  - Rでその両方を書いておく

my.Prob.numallele <- function(n,theta){
	prts <- parts(n)
	prbs <- rep(0,length(prts[1,]))
	for(i in 1:length(prts[1,])){
		prt <- prts[,i]
		print(prt)
		tab <- tabulate(prt)
		prbs[i] <- my.Ewens.prob(tab,theta)
	}
	ks <- apply(prts,2,function(x){length(which(x!=0))})
	K <- sum(prbs*ks)
	ks2 <- 1:n
	prbs2 <- rep(0,n)
	for(i in 1:n){
		prbs2[i] <- sum(prbs[which(ks==i)])
	}
	return(list(K=K,ks=ks,prbs=prbs,ks2=ks2,prbs2=prbs2,prts=prts))
}
library(gmp) # Stirling1()
my.Prob.numallele2 <- function(n,theta){
	ks <- 1:n
	prbs <- rep(0,n)
	for(i in 1:n){
		tmp <- ks[i] * log(theta) - sum(log((0:(n-1)) + theta))
		prbs[i] <- abs(as.integer(Stirling1(n,ks[i]))) * exp(tmp)
	}
	K <- sum(ks*prbs)
	return(list(K=K,ks=ks,prbs=prbs))
}
my.Prob.numallele(4,2.1)
my.Prob.numallele2(4,2.1)

- Ewens's formulaと中華料理店過程
  - このペイパーでは中華料理店過程は少し先の方で出てくるのだが、先取りすると、中華料理店過程では、n+1番目の客は $\frac{\theta}{\theta+n}$ の確率で誰も座っていないテーブルにつき、そのほかの場合には、既に客が座っているテーブルのどれかに着席するような過程である。Ewens's formula がこの過程と同じ確率過程になっていることを以下のようにRの計算的に確認する
  - n染色体でアレル種類数が1,2,...,nになっている確率は上述の方法で算出される
  - n+1染色体での1,2,...,n+1も算出できる。n染色体からn+1染色体への移行は確率事象であり、n+1染色体でn+1種類と言う場合は、n染色体でn種類の場合から、新しい染色体が新しいアレルタイプだった場合のみである。これを使うと、n染色体でn種類だった状況から、n+1染色体でn種類のままである場合が差分として計算できる
  - 次にn+1染色体でn種類である場合は、n染色体でn種類の状態からの推移とn染色体でn-1種類だったところから種類が1つ増える場合とで構成されるから、n染色体でn-1種類だったときに新種が生じる確率が差分で計算できる
  - これを繰り返すことで、n染色体でk=1,2,...,n種類だった状況から、新種が加わる場合と旧種に落ちる場合との確率がすべて列挙できる
  - その新種付加とそれ以外の振り分け割合が $\theta : n$ となる(k=1,2,...,nのすべての段階で)ことを以下のように確かめる

my.singleton.pp <- function(n,theta){
	tmp <- my.Prob.numallele2(n,theta)
	tmp2 <- my.Prob.numallele2(n+1,theta)
	prbs <- tmp$prbs
	prbs2 <- tmp2$prbs

	a <- b <- rep(0,n)
	for(i in n:1){
		if(i == n){
			a[i] <- prbs2[i+1]
			b[i] <- prbs[i]-a[i]
		}else{
			a[i] <- prbs2[i+1]-b[i+1]
			b[i] <- prbs[i]-a[i]
		}
	}
	return(list(a=a,b=b))
	
}
n <- 5
theta <- runif(1)
out <- my.singleton.pp(n,theta)
out[[1]]/(out[[1]]+out[[2]])
theta/(theta+n)

> out[[1]]/(out[[1]]+out[[2]])
[1] 0.1287112 0.1287112 0.1287112 0.1287112 0.1287112
> theta/(theta+n)
[1] 0.1287112

- Coalescent theoryはEwens's formulaと同じWright-Fisherモデルを無限大集団に想定し、無限にある染色体が同祖で統合される過程(時間逆行の過程)に関するものである。中立的統合事象が集団メンバー(染色体)が作るペア数に比例して生じるので統合イベントは時間に関して指数関数表現を取る。そのようにしてできるCoalescent treeについて、ランダムに起きる変異を考慮するとアレル種類数についての確率過程を表したものとなる
3. Ewens's sampling formulaの特徴
- アレル別本数情報(全部でn本のうち、あるアレルがj本あるようなアレルが $m_j$ 種類ある) $\sum_{j=1}^n j m_j= n$ には、次のような推移関係がある
- $m_j' = m_j + 1; m_{j-1}' m_{j-1} - 1; m_i' = m_i (i \ne j,j-1)$
- これはEwens's sampling formulaが使う整数ベクトル
  - $P(m_1,...,m_n|\theta) = \frac{n!}{\theta(\theta+1)...(\theta+n-1)}\prod_{j=1}^n \frac{\theta^{m_j}}{j^{m_j}m_j!}$ がそのformula
- これを、アレル種類数ごとにプールすると第１種スターリング数が出てきて $P(K=k) = |S^k_n| \frac{\theta^k}{\theta(\theta+1)...(\theta+n-1)}$ となる
- さらにこれらとは別にと、k種類のアレルに分ける、という分け方もある。であり、,である。
  - Ewens distributionと呼ばれる以下の式が成り立つ
  - $P(\Pi_n=\{b_1,b_2,...,b_k\}) = \frac{\theta^k}{\theta(\theta+1)...(\theta+n-1)}\prod_{j=1}^k(b_j-1)!$ 、kは種類数
- k=1,2,...,n種類に分ける集合としての分け方を全列挙し、そのそれぞれの生起確率を計算すれば、総和が１になる
- それをRで確認する

my.Ewens.dist <- function(bs,theta,log=FALSE){
	k <- length(bs)
	n <- sum(bs)
	tmp <- k * log(theta) -sum(log((0:(n-1))+theta)) + sum(lfactorial(bs-1))
	if(log){
		return(tmp)
	}else{
		return(exp(tmp))
	}
}
# Ewens formulaと関連付けておく
my.Ewens.dist2 <- function(bs,theta,log=FALSE){
	n <- sum(bs)
	
	tab <- tabulate(bs)
	tmp <- my.Ewens.prob(tab,theta,log=TRUE)
	
	tmp2 <- tmp -lfactorial(n) + sum(tab*lfactorial(1:length(tab))) + sum(lfactorial(tab))
	if(log){
		return(tmp2)
	}else{
		return(exp(tmp2))
	}

}
# 集合としての分割全列挙とその確率
my.list.parts <- function(n,theta){
	lst <- listParts(n)
	ret <- rep(0,length(lst))
	for(i in 1:length(ret)){
		ret[i] <- my.Ewens.dist2(sapply(lst[[i]],length),theta)
	}
	return(list(lst=lst,prob=ret))
}

n <- 4
theta <- 0.3
tmp.out <- my.list.parts(n,theta)
tmp.out$lst
sum(tmp.out$prob)

> tmp.out$lst
[[1]]
[1] (1,2,3,4)

[[2]]
[1] (1,2,4)(3)

[[3]]
[1] (1,2,3)(4)

[[4]]
[1] (1,3,4)(2)

[[5]]
[1] (2,3,4)(1)

[[6]]
[1] (1,4)(2,3)

[[7]]
[1] (1,2)(3,4)

[[8]]
[1] (1,3)(2,4)

[[9]]
[1] (1,4)(2)(3)

[[10]]
[1] (1,2)(3)(4)

[[11]]
[1] (1,3)(2)(4)

[[12]]
[1] (2,4)(1)(3)

[[13]]
[1] (2,3)(1)(4)

[[14]]
[1] (3,4)(1)(2)

[[15]]
[1] (1)(2)(3)(4)

> sum(tmp.out$prob)
[1] 1

- Exchangeabilityは、ラベルパーミュテーションによってブロック帰属相互関係が変わらないような時に同一視えきるような分割の性質
- Consistency under subsamplingは、帰属要素を減らした分割が入れ子関係になることを言う
- Self-similarityはフラクタルのようにどこを取っても分割として同じ性質を持つこと
- Noninterferenceは分割要素の関係が対等・平等
- 指数型分布族である。 $P(\Pi_n=\{b_1,b_2,...,b_k\}) = \frac{\theta^k}{\theta(\theta+1)...(\theta+n-1)}\prod_{j=1}^k(b_j-1)!$ は $P(\Pi_n=\pi) = exp(Number(\pi) \log{\theta} - \sum_{j=0}^{n-1}\log{\theta+j}) \times \prod_{b \in \pi} (Number(b)-1)!)$ と書けて、パラメタ $\theta$ は $Number(\pi)$ とだけ絡んでいるから、 $Number(\pi)$ が十分統計量
- ポアソン過程とつないで考えることもできて、それをLogarithmic combinatorial structuresと言う。
  - Ewens's sampling formulaが $P(m_1,...,m_n|\theta) = \frac{n!}{\theta(\theta+1)...(\theta+n-1)}\prod_{j=1}^n \frac{\theta^{m_j}}{j^{m_j}m_j!}$ であったことを思い出そう
  - 今、期待値が $\frac{\theta}{j};j=1,2,...$ であるようなポアソン分布に従う乱数 $Y1,Y2,...$ があるとする。ポアソン分布は $\frac{\lambda^k}{k!}e^{-\lambda}$ なので、この $\lambda$ が $\frac{\theta}{j}$ とすると、 $P((Y_1,...,Y_n)=(m_1,...,m_n)|\sum_{j=1}^n j Y_j=n)$ という確率は、n個のポアソン分布の確率質量関数の積なので $\prod_{j=1}^n \frac{\theta^{m_j}}{j^{m_j}m_j!}e^{-\frac{\theta}{j}}$ に比例する。 $\prod e^{-\frac{\theta}{j}}$ の部分はmによらないから無視するとして、残りの部分がEwens's sampling formulaと一致することから、Ewens's sampling formulaはn個の乱数かポアソン乱数であって、その期待値が特定のセットになっている特別な場合であることがわかる
  - これを一般化すると、多数の0,1,2,...の値をとる乱数を考え、ただしその乱数の確率質量分布は自由にすることにした分割の生起確率を考えることができる。これが、乱数の積としてのランダム分割になる
  - ちなみに、期待値 $\theta/j$ のポアソン乱数 $Y_j$ がとる値は、n本の染色体の観察でj本の同じタイプの染色体が観察されるタイプ数に対応していることでEwens's sampling formulaとつながっている
4. Poisson-Dirichlet Distribution(パラメタ数２、もしくは１)
Exchangeable partitionsを台とする確率分布を考える
- n個の内訳が $\sum_{j=1}^n j m_j =n$ となるような $(m_1,...,m_n)$ を確率的に定めるのがExchangeable partitionsを台とする確率分布であり、そのような確率分布に従って、ランダムに内訳を作るとき、それをExchangeable random partitionsと言う
どんな分布があるかというと、Poisson-Dirichlet Distributionと言う分布(確率分布をパラメトリックに決める規則(law))がある
- ２パラメタの分布だが、１パラメタで表現できるものが含まれ、それらにいろいろ説明がついている(そのためにかえって混乱してしまうかもしれない)
  - １パラメタのPoisson-Dirichlet Distributionの方には、２パラメタを１パラメタにする方法によって異なる規則が含まれているが、そのうちのひとつに対応する確率質量分布がEwens's sampling formula( $\alpha=0$ )であり。２パラメタのPoisson-Dirichlet Distributionの方は、その確率質量分布としてEwens-Pitman formula(distribution)が知られている
- もちろん、そのほかにも、何かしらの確率過程を考案すれば、それがExchangeable random partitionsを生成することはあるし、それを規則としているということで分布名をつけてもよいが、Ewens-Pitman Distributionは前述のSelf-similarityとかNon-interferenceとか指数型分布族とか言うよい性質を持っていることから取り立てられている(のだと思う)
- １変数をやってから２変数をやるのもよいが(ペイパーはそうなっている)、２変数を先に与えて、１変数を眺めなおすほうが(自分には)わかりやすいので、そのようにしてやる
- Ewens-Pitman sampling formula
  - ある２つのはパラメタ $\theta,\alpha$ で定義された確率過程が分割を定める規則になっているとする。n標本を抽出したところ、同種の標本が $j=1,2,...,n$ であるような種が $m_j$ 種類あったとする。 $\sum_{j=1}^n j m_j=n$ であり、 $k=\sum_{j=1}^n m_j$ は観察された種類の数
  - このとき、 $(m_1,...,m_n)$ となる確率が以下のようになるとされ、その確率計算式をEwens-Pitman sampling formulaと言う
  - $P((m_1,...,m_n)|\alpha,\theta) = \frac{n!}{f(\theta,n)}\prod_{l=1}^{k}(\theta+(l-1)\alpha)\prod_{j=1}^n \frac{f(1-\alpha,j-1)^{m_j}}{(j!)^{m_j} m_j!}$ ,ただし $f(\theta,n) = \theta (\theta+1) ... (\theta + n-1)$
- Ewens sampling formulaはEwens-Pitman sampling formulaのの場合
  - 式を比較してみる
  - $P^{Ewens}((m_1,...,m_n)|\theta) = \frac{n!}{f(\theta,n)}\prod_{j=1}^n \frac{\theta^{m_j}}{j^{m_j} m_j!}$
  - $\alpha=0$ のとき
  - Ewens-Pitman formulaの $\prod_{l=1}^{k}(\theta+(l-1)\alpha) = \prod_{l=1}^k \theta = \theta^k = \theta^{\sum_{j=1}^n m_j}=\prod_{j=1}^n \theta^{m_j}$ であり
  - $\frac{f(1-\alpha,j-1)^{m_j}}{(j!)^{m_j}} = \frac{f(1,j-1)^{m_j}}{(j!)^{m_j}} = \frac{1}{j^{m_j}}$ であることから、一致することがわかる
- 確率質量分布がEwens-Pitman sampling formulaであった。それに対応する分割分布(,ただしKは観測される種類数、のように分けるわけ方)の確率質量の式は
  - $P(\Pi_n = \pi) = \frac{f(\theta/\alpha,K)}{f(\theta,n)} \prod_{j=1}^K -f(-\alpha,b_j)$
  - この式に対して $\alpha=0$ を考えようとすると困る。その場合には
  - $\alpha= -\frac{\theta}{k}$ と置くことにし、 $k \to \infty$ を考えることとする
  - このとき $P(\Pi_n = \pi) = \frac{f(-k,K)}{f(-k \alpha,n)} \prod_{j=1}^K -f(-\frac{\theta}{k},b_j)$ となる
  - ここで、 $f(-k,K)$ は $k \to \infty$ で $(-k)^K$ とみなせる
  - その $(-k)$ を、 $-f(-\frac{\theta}{k},b_j)$ に掛けると、第一要素は $\theta$ を残し、残りの $b_j-1$ 要素は $1,2,...,b_j-1$ となることから
  - $P(\Pi_n = \pi) = \frac{f(-k,K)}{f(-k \alpha,n)} \prod_{j=1}^K -f(-\frac{\theta}{k},b_j) = \frac{\theta^K}{f(\theta,n)}\prod_{j=1}^K (b_j-1)!$ となる。これはEwens版の式に他ならない
  - なお、 $\alpha <0 , \theta = -k \alpha; k=1,2,...$ または
  - $0 \ge \alpha \ge 1, \theta > \alpha$ がパラメタのとりうる範囲(らしい)
- Rで書いておく

my.ascending <- function(theta,n,log=TRUE){
	ret <- sum(log(abs((1:n) + theta -1)))
	
	if(!log){
		ret <- exp(ret)
	}
	return(ret)
}
my.Ewens.Pitman.dist <- function(b,theta,alpha,log=TRUE){
	n.b <- length(which(b!=0))
	n <-sum(b)
	if(alpha == 0){
		ret <- n.b * log(theta) - my.ascending(theta,n) + sum(lfactorial(b-1))
	}else{
		tmp <- 0
		for(i in 1:length(b)){
			tmp <- tmp + my.ascending(-alpha,b[i])
		}
		ret <- my.ascending(theta/alpha,n.b) - my.ascending(theta,n) + tmp
	}
	if(!log){
		ret <- exp(ret)
	}
	return(ret)
}


my.Ewens.Pitman.dist2 <- function(b,theta,alpha,log=TRUE){
	n.b <- length(which(b!=0))
	n <-sum(b)
	if(alpha == 0){
		ret <- n.b * log(theta) - my.ascending(theta,n) + sum(lfactorial(b-1))
	}else{
		tmp <- 1
		for(i in 1:n.b){
			tmp <- tmp * (theta/alpha + i-1)
		}
		for(i in 1:n){
			tmp <- tmp / (theta + i-1)
		}
		for(i in 1:n.b){
			tmp <- tmp * (-1)
			for(j in 1:b[i]){
				tmp <- tmp * (-alpha + j -1)
			}
		}
		ret <- log(tmp)
	}
	if(!log){
		ret <- exp(ret)
	}
	return(ret)
}


my.Ewens.Pitman.dist(c(3,1,1),1.2,-0.3)
my.Ewens.Pitman.dist2(c(3,1,1),1.2,-0.3,log=TRUE)

bs <- c(50,2,1,1,1)
theta <- 0.2
my.Ewens.dist(bs,theta,log=TRUE)
my.Ewens.Pitman.dist2(bs,theta,0,log=TRUE)

ks <- length(bs):100
alphas <- -theta/ks
out <- rep(0,length(ks))
for(i in 1:length(ks)){
	out[i] <- my.Ewens.Pitman.dist(bs,theta,alphas[i],log=TRUE)
}

plot(out)
out.lim <- my.Ewens.dist(bs,theta,log=TRUE)

rng <- range(c(out),out.lim)

plot(out,ylim = rng + c(-0.5,0.5),type="l")
abline(h=out.lim,col=2)

- 確率過程的な意味づけ
  - 式としては上記の通りだが、いくつかの確率過程としての説明がある
  - - この場合は、いわゆるEwens sampling formulaに対応する場合
    - 遺伝学でのいわゆる集団サイズ一定、ランダムメイティング、非抽出サンプリング、無限サイトモデルのときにn染色体のアレルタイプ別観測数の分布を表す
    - それは、中華料理店過程で、新規テーブルに着席する確率を $\frac{\theta}{\theta+n}$ とする場合にも相当する(壷のモデルも同じこと)
    - そして、 $(-alpha,-alpha,....)$ というパラメタのディリクレ分布から発生した頻度ベクトルに対して、n標本をサンプリングしたときの多項分布について、多項を第１、第２、…のように項の違いを意識するのではなく、項が揃っていれば第何番目の項かは無視するとしたときの、多項分布を考えるときに、この項数を無限大にしたときが、やはりこれに対応することが示せる。ペイパーでは、そのあたりの式変形を示しているが、ディリクレ分布に観測を上乗せした確率の積分の式で、その中にディリクレ分布が隠れていることとそれを全体に積分すると１になることと、項を区別しないために、設定としての項数kと観察された項数pについて $k(k-1)...(k-p+1)$ の重複分の考慮などに注意すれは式変形を追うことはできる
- 他方、Ewens-Pitman sampling formulaは２パラメタである
  - これに対応する中華料理店過程を考えると、「具体的な過程」のイメージは湧く
  - すでに客が居るテーブル数がKだとしたとき、新しいテーブルに着く確率を $\frac{\theta+K\alpha}{\theta+n}$ とし、 $\frac{n_i-\alpha}{\theta+n}$ の確率ですでに客が居る席に着くという過程(Pitman-Yor過程)を考える
  - この $\theta,\alpha$ に対応する確率質量分布が２パラメタのPoisson-Dirichlet distribution
  - n標本観察して、何種類がその中に観察されるかの期待値が計算できたりする。nを無限に飛ばしたときにどんな値になるかといえば、有限次元Dirichlet-Multinomに対応する場合は有限の値に収束するし、 $\alpha=0$ は無限種類であるから、nに応じて単調増の値が示される。それらは $\alpha \le 0$ の場合だが、 $\alpha > 0$ の場合には、別の収束関数が認められる
  - Gamma and Stable Subordinatorsの節では、inhomogenous poisson point processに結びつける話が書かれている。ある確率過程に従って何かが生じ、それに基づいた確率過程が起きるとする。最初の確率過程がsubordinator。Wiki記事
  - Size-biased sampling のところでは、だんだん小さく区切る、という方法での分割生成方法が示され、そのひとつとしてStick-breaking法が提示されている。そこで使うベータ分布のパラメタとEwens-Pitman distributionsとの関係が書かれている
6. ノンパラメトリック・ベイズ
- Ewens-Pitman distributionはstick-breakingで作れたり、Subordinatorを用いて作れたりするので、computational Bayesinanに向いている
- 応用例がいろいろある
7. Induvtive inference (帰納的推定)での位置
8. 組みあわせ、代数、数論での位置

2018-03-20

ポアソン点過程・分割・ノンパラメトリックベイズ

ポアソン点過程分割ノンパラメトリックベイズ

動機
- 色々動機はあるかもしれないが
- 多数のもの・無限個あるかもしれないもののタイプ分けが興味の対象
- クラスタ数不定な状況でのクラスタリング
- そのための確率モデル
- その確率モデルの下での生起確率・事前確率・尤度・事後確率
モデル
- 具体的な説明から始めよう
  - 中華料理店過程
    - 中華料理店に客が１人ずつ入ってきて着席する
    - 着席には確率的ルールがあると仮定する
    - 客は着席風景を眺め、誰も座っていないテーブルをある確率で選ぶ。その確率はすでに座っている客の数が多いほど小さくなる( $\frac{\alpha}{\alpha + n}$ 、ただし $n+1$ -番目の客の場合)
    - 誰も座っていないテーブルに座ることはせず、誰かが座っている席に座るとしたら、座っているテーブルの人数に比例して座る(同じことだが、既に座っている人の中から等確率で１人を選んでその人と同じテーブルに座る)。その確率は、 $\frac{n_i}{\alpha + n}$ である。ただし $n_i$ はテーブル-iの人数。
    - もちろん $\frac{\alpha}{\alpha + n}+\sum_i \frac{n_i}{\alpha + n} = 1$
- ルール,を見直す
  - 壺モデル
    - $\alpha$ は正の実数として考えているが、これを1とすると、Polya urnモデル。黒いボールが１個からスタートし、urn(壺)から取り出すにあたり、黒の場合は新しい色(タイプ)のボールと黒ボールを戻し、黒以外のボールの場合は取り出したボールと同色(同タイプ)のボールとの2個を戻すことの繰り返しをするときに壺の中のボールの色別割合が中華料理店過程と同じになる
  - $:alpha$ を自然数とすれば、壺に入れておく黒ボールの個数となる
  - 一般化を進めたPiman-Yor モデル
    - ルール $\frac{\alpha}{\alpha + n}$ , $\frac{n_i}{\alpha + n}$ は、行動1によってタイプ数を増やし、行動1と行動2との割合を試行回数nによって次第に変化させ、タイプ数を増やさない場合に既存タイプ数を増やすにあたり、行動2に割り振られた確率を既存タイプ数別に割り振り、全体の確率を1にする式になっている
    - このルールを守りさえすれば、適当に確率モデルを作ることができる
    - $\frac{\alpha + k \theta}{\alpha + n}$ , $\frac{n_i-\theta}{\alpha + n}$ ,ただしkはその時点でのタイプ数(i=1,2,...,k)という割り振りも条件を満足する。ここでは中華料理店過程に $\theta$ というパラメタを加えて新規タイプの発生確率を変化させている( $\theta>0$ なら新規タイプが中華料理店過程より増えやすい)
表現
- 壺モデル・中華料理店過程・Pitman-Yorモデルは、壺やテーブル着席と言った、「具体的なイメージ」を伴った確率過程の表現
- 「具体的なイメージ」では第i番目が加わると…という状況であり、かつそれが無限に続けられる状況である。それによって無限自然数分割を扱っている
- それを「全体を無限に分割し続ける」という過程に置き換えたのが、Stick-breaking 過程のように、単位線分に分割点を入れていく表現。この場合は無限に分割を続けられるように、「だんだん小さな線分を作っていく」というやり方を採用している。この「だんだん小さく」という制約と関係するのが"Size-biased"という考え方。
- また、どのように分割を作成するか、と言う点で、だんだん小さくしないやり方もある。たくさんの正の乱数を発生させ、それの和が1になるように乱数の総和で標準化する、という方法。このとき注意するべきは正の乱数の総和が発散しないような乱数発生法を採用すること。そのような発生法として、正の実数直線に正の値を取る関数を定め、その関数の値は、その点における平均生起回数を表し、その平均生起回数に応じてポアソン分布的に値が発生する、というモデルを作ることにする。場所によって平均生起回数パラメタが異なるポアソン点過程なので、inhomogeneous Poisson processと言う。これによって、この関数が「生ぜしめる分割」を定めていることになる。この点生起過程も無限に続けられる。ポアソン点過程を使いつつ、その総和によって「割合化」しているが、この割合化するところは、Dirichlet分布を作るときに多数を定めてその総和が１になるようにすることと同じなので、ポアソン・ディリクレ過程と言う(らしい)
- もう一つ。Size-biased的にでもpoisson-dirichlet的にでもよいが、何かしらの分割をしたとしよう。無限の分割点・無限の線分長ができてしまう。それが目的だから仕方ないが、やはり無限個の要素は扱いにくい。有限個の要素であたかも無限個の要素を扱えるのが、Kingmanのpaintbox。これは単位線分を有限線分に分割し、そのうちの１つを除いたすべての有限長線分は、タイプ別の割合を表し、１つの線分は無限要素を持つ全体集合の中でSingletonになっているものたちの集合であるとする考え方。これをすることで、有限個の、有限割合を有するタイプ(そのタイプには無限個の要素が帰属する)と、無限個のタイプ(ただしすべてSingleton)とに分割することになる。そしてこの一見特殊な無限タイプへの分割が任意のExchangeable random partitionsに1対1対応することがしられている(Kingman's correspondence)ので、これも便利。１つの便利さは、有限分割して、それに応じて無限にタイプ別標本を発生させられること

2018-03-19

Ewens sampling formula

Ewens sampling formulaは、集団のハプロタイプ生成状況にある条件を課したときに、標本の頻度分布の正確確率に関する式
これは集団遺伝学の領域では、最も単純な条件でのCoalescent過程が生成する頻度分布になる
Coalescent過程は頻度分布のみを定めるだけではなく、それをもたらすに至る点変異の集積状態・その時系列状況についてもモデル化したもの
このEwens sampling formulaは遺伝的変異とか時間とはをすっ飛ばして、たんたんと中華料理店過程で生じるテーブル別人数分布(確率的整数分割の一つ)も説明する
変異・アレル頻度が中華料理店過程と頻度分布・単位線分分割・整数分割の極限として同一になるのは、変異・アレル頻度生成が、新規変異とCoalescentという2種類の確率過程の合成確率過程である一方、中華料理店過程も、新しいテーブルを選ぶ・すでにあるテーブルを確率的に選ぶという２つの確率的過程である点での一致だからのようだ(2パラメタPoisson-Dirichlet過程として同じ種類になる)

2018-03-17

Poisson 点過程で作る Random Exchangeable Partitions

Random Exchangeable Partitions Poisson Dirichlet過程 Coalescent Size-biased Pick

Stick breaking process/中華料理店過程でExchangeable random partitionsが作れることを前の記事で書いた
同じことを別の作り方として表現できる
正の値をランダムに発生させ、その総和が1になるように標準化すれば、足して1になる多数の(無限の)正の数の集合ができる。これは今、欲しい値の集合である
ただし、これができるには、総和が有限になるようにランダムな値を発生させることである
この $(0,\infty)$ に値をランダムに発生させつつ、その総和の極限が発散しないような作り方をPoisson 点過程を用いた単位線分の分割作成法と言う
この発散しないような発生方法を、 $(0,\infty)$ に定義された確率密度関数に基づく乱数発生とみなせば、確率密度関数を指定することでRandom Exchangeable Partitionsが定義できる
たとえば、指数分布で発生させたとする。指数乱数はポアソン過程の間隔の分布なので、指数乱数の総和は、ポアソン過程で発生する事象の最後のイベントが起きる時刻に相当する。しかしながら、今、イベント数は無限回にしようとしているので、その総和は無限大になる。このことからわかるのは、指数分布は単調減少する確率密度分布だが、それに従う乱数を使ったのでは、今の用途には不適切だということ
$P(x) = \frac{1}{x}e^{-x}$ のような分布は指数分布よりも減衰が早い分布なので、これなら(多分)使える
さて。前の記事で、Poisson-Dirichlet Processについて書いた。このPoisson-Dirichlet Processは２つのパラメタ $\alpha,\theta$ で定まる過程であったが、これに対応する確率密度分布があって、それにしたがって乱数を発生させ、それの総和を取って除するという処理で単位線分タイルの長さを決めると、Poisson-Dirichlet Processの分割分布になることが知られている
そのような関数として $x^{-\alpha-1}$ とか $\theta x^{-1} e^{-x}$ などがある
さらに、そのような確率密度関数に対応するProbability Generating Functionalのラプラス変換の式は、 $\alpha,\theta$ に関して連続して定義できることも知られている(こちらの2ページから3ページにかけて)
[tex(0,\infty)]の範囲に点が発生する過程であって、位置によって点の発生確率が異なるので、これをinhomogeneous Poisson point processと言う

2018-03-17

Random Exchangeable Partitions

Random Exchangeable Partitions Poisson Dirichlet過程 Coalescent Size-biased Pick

Random Exchangeable Partitions
- Partition(分割)を考える
  - 何を分割するのかが問題になる。ある正の整数Nを分割する。このとき $\{1,2,...,N\}$ という集合を排他的な部分集合に分割する、という考え方もあるが、それだと、「1から正の整数Nまでの整数集合を分割する」と表現する方がより正確であることを考えると、「正の整数Nを分割する」というときには、それを $\{n_1,...,n_k\};\sum n_i = N$ という分割を指すと考えるたい。いわゆる整数分割である
- Exchangeable(交換可能)とは
  - 長さNの数字列をいろいろに分割した集合を考えた時に、分割した後で、1,2,...,Nの数の並びをシャッフルしてできる分割の集合が、元の分割の集合とシャッフル後の分割の集合とが分布として同じになるとき、この分割の作り方とその作り方によって作られる分割の集合が、「交換可能〜交換しても分布として同じ」と言う
- Randomとは
  - 上記で、分割の作り方のことを書いたがその作り方が、ランダムである・確率過程的である、ということ
有限のRandom Exchangeable Partitionsを悉皆列挙するとすれば…

library(partitions)
N <- 3
# 整数分割
prts <- parts(N)
prts
# 整数分割ごとに分割画分の成分を表示
for(i in 1:length(prts[1,])){
	tmp <- prts[,i]
	tmp. <- cumsum(tmp)
	tmp.. <- tmp.[which(tmp!=0)]
	print("---")
	print((1:N)[1:tmp..[1]])
	if(length(tmp..)>1){
		for(j in 2:length(tmp..)){
			print((1:N)[(tmp..[j-1]+1):tmp..[j]])
		}
	}
}
# 順列
prms <- perms(N)
# ある順列でシャッフルした場合の成分列挙
# 分割画分の成分は亜集合として同一

for(i in 1:length(prts[1,])){
	tmp <- prts[,i]
	tmp. <- cumsum(tmp)
	tmp.. <- tmp.[which(tmp!=0)]
	print("---")
	print((prms[,2])[1:tmp..[1]])
	if(length(tmp..)>1){
		for(j in 2:length(tmp..)){
			print((prms[,2])[(tmp..[j-1]+1):tmp..[j]])
		}
	}
}

> prts
          
[1,] 3 2 1
[2,] 0 1 1
[3,] 0 0 1

[1] "---"
[1] 1 2 3
[1] "---"
[1] 1 2
[1] 3
[1] "---"
[1] 1
[1] 2
[1] 3


[1] "---"
[1] 1 3 2
[1] "---"
[1] 1 3
[1] 2
[1] "---"
[1] 1
[1] 3
[1] 2

2018-03-17

確率的に単位線分を分割し続けて作る Random Exchangeable Partitions

Random Exchangeable Partitions Poisson Dirichlet過程 Coalescent Size-biased Pick

長さ１の単位線分をあるルールで確率的に分割していけば、それもRandom Exchangeable Partitionsとなる
Poisson-Dirichlet Processと呼ばれる方法がその一つで、よく研究されている
何度でも分割し続けるルールとして、単位線分から出発して、分割点を１点とり二分する。生じた２つの線分のうち片方(たとえば1に近い方)はそれ以上いじらないことにして、もう片方をまた二分する。これはいつまでも繰り返せる。この二分点の取り方として、ベータ乱数を使う方法は確率過程として扱いやすく、この方法をPoisson-Dirichlet Processと言う
ベータ乱数を発生させる元となるベータ分布としてとするという。ただしベータ乱数の発生に当たり、i回目ごとに用いるベータ分布のパラメタがに応じて変化するように設計されていることに注意。また、この定義の中で,の２つの場合は特殊形として、単純なモデルとして取り上げられることがある
- また $\alpha=0,\theta=1$ の場合は、 $Beta(1,1)$ になり、一様乱数での分割の繰り返しになる
Rでやってみる

theta <- runif(1)
alpha <- runif(1)

n.iter <- 100
beta <- rep(0,n.iter)
for(i in 1:n.iter){
	beta[i] <- rbeta(1,1-alpha,theta + i * alpha)
}
P <- rep(0,n.iter)
P[1] <- beta[1]
for(i in 2:n.iter){
	P[i] <- (1-sum(P[1:(i-1)])) * beta[i]
}

plot(P)

S <- cumsum(P)

plot(S,type="b",pch=20,cex=0.1,ylim=c(0,1))

- 特に $\alhap=0$ の場合は

theta <- runif(1)
alpha <- 0

n.iter <- 10
beta <- rep(0,n.iter)
for(i in 1:n.iter){
	beta[i] <- rbeta(1,1-alpha,theta + i * alpha)
}
P <- rep(0,n.iter)
P[1] <- beta[1]
for(i in 2:n.iter){
	P[i] <- (1-sum(P[1:(i-1)])) * beta[i]
}

plot(P)

S <- cumsum(P)

plot(S,type="b",pch=20,cex=0.1,ylim=c(0,1))

- また、 $\theta=0$ の場合は

theta <- 0
alpha <- runif(1)

n.iter <- 10
beta <- rep(0,n.iter)
for(i in 1:n.iter){
	beta[i] <- rbeta(1,1-alpha,theta + i * alpha)
}
P <- rep(0,n.iter)
P[1] <- beta[1]
for(i in 2:n.iter){
	P[i] <- (1-sum(P[1:(i-1)])) * beta[i]
}
plot(P)

S <- cumsum(P)

plot(S,type="b",pch=20,cex=0.1,ylim=c(0,1))

Poisson-Dirichlet Processの別の見方
- $\alpha=0,\theta$ での分割は、長さ１の線分を確率的ルールで分割していくものと読み取れる。Stick-breaking processと称される理由である
- 次のように中華料理店過程がある意味で(*)同じ分割をもたらす
- 中華料理店過程ではi+1番目の客が、i番目までの客の着席状態に応じて、新しいテーブルに1人で座るか、既に客の居るテーブルのどれかに座るかを決める、というルールで着席状態が変化していく過程である
- の場合には
  - $\frac{\theta}{i+\theta}$ の確率で新しいテーブルに着席し、 $\frac{m_j}{i + \theta}$ の確率でj番目のテーブルに着席するという確率過程のときに現れる整数分割(とそれに対応する単位線分分割)と同じであると言う。ただし $m_j$ はi人まで着席が終わったときのj番目のテーブルの人数であり、 $\sum_{j=1}^k m_j = i$ である。ここでkはi番目までの客で着席者の居るテーブルの数である
- の場合には
  - $\frac{k\alpha}{i}$ の確率で新しいテーブルに着席し、 $\frac{m_j-\alpha}{i}$ の確率でj番目のテーブルに着席する
- この中華料理店過程で得られるランダムな分割はexchangeableなランダム分割になっており、size-biased ordering of the asymptotic block frequeiciesがstick-breakingで得られるものになることが知られている。"size-biased"というのは、「ある誰かしらを特定して、その人が座っているテーブルの人数割合」を調べる場合に使う用語である。そのようにすると、テーブルの人数割合は、割合の大きいテーブルが割合に比例して選ばれるので、そのようにして得られる割合の分布は、割合で重みづいた分布となるのだが、この分割の分野では、この点を考慮することがよくある