ぱらぱらめくる『The Ubiquitous Ewens Sampling Formula』 - ryamadaの遺伝学・遺伝統計学メモ

ペイパーはこちら
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. 組みあわせ、代数、数論での位置