ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2018-05-15 最小二乗法

[][]固有値分解する

  • (X ¥beta - y)^T (X ¥beta -y)の最小化
  • ((X,y) ¥begin{pmatrix}¥beta¥¥ -1¥end{pmatrix})^T ((X,y) ¥begin{pmatrix}¥beta¥¥ -1¥end{pmatrix})の最小化
  • (¥beta -1)^T (X y)^T (X y) ¥begin{pmatrix} ¥beta ¥¥ -1 ¥end{pmatrix}
  • (¥beta -1)^T M ¥begin{pmatrix} ¥beta ¥¥ -1 ¥end{pmatrix}
  • Mの最小固有値に対応する固有ベクトルが一番拡大率が小さいから、その固有ベクトル方向の(¥beta,-1)が求める解
n <- 100
X <- matrix(rnorm(n*2),ncol=2)
Y <- X %*% c(1,2) + rnorm(n,0,0.01)
library(rgl)

plot3d(X[,1],X[,2],Y)

lm(Y ~ X-1)

XY <- cbind(X,Y)
M <- t(XY) %*% XY

eigenout <- eigen(M)

eigenout[[2]]

eigenout[[2]][,3]/(-eigenout[[2]][3,3])

2018-03-21 ぱらぱらめくる『The Ubiquitous Ewens Sampling Formula』

[][][][][]ぱらぱらめくる『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,...の場合わけごとの確率の総和が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からベタに計算することもできるし、第1種スターリング数を使った定義式P(K=k) = |S^k_n| ¥frac{¥theta^k}{¥theta(¥theta+1)...(¥theta+n-1)}でも計算できる
      • また、このように確率モデルのパラメタ¥thetaが観察条件 n と観察項目 K(アレルの種類数)とで表されているということは、Kを観察することが¥thetaを推定するための十分統計量となることを意味しており、それは、種類数の確率分布が指数型分布族形式で表され、分布を定めるパラメタ¥thetaと観察状態とが絡む項に、観察項目数のみが登場することと関係している(本ペイパーの3.7で示されている)
      • ちなみに。ここで第1種スターリング数が出てくるのは、スターリング数がある場合の数を表しているからであり、その「ある場合の数」とは、『アレル種類数が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
    • これを、アレル種類数ごとにプールすると第1種スターリング数が出てきてP(K=k) = |S^k_n| ¥frac{¥theta^k}{¥theta(¥theta+1)...(¥theta+n-1)}となる
    • さらにこれらとは別に¥{b_1,b_2,...,b_k}と、k種類のアレルに分ける、という分け方もある。¥sum_{i=1}^k b_i = nであり、¥sum_{j=1}^n m_j = k,¥sum_{j=1}^n j m_j = nである。
      • 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種類に分ける集合としての分け方を全列挙し、そのそれぞれの生起確率を計算すれば、総和が1になる
    • それを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)!)と書けて、パラメタ¥thetaNumber(¥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(パラメタ数2、もしくは1)
  • Exchangeable partitionsを台とする確率分布を考える
    • n個の内訳が¥sum_{j=1}^n j m_j =nとなるような(m_1,...,m_n)を確率的に定めるのがExchangeable partitionsを台とする確率分布であり、そのような確率分布に従って、ランダムに内訳を作るとき、それをExchangeable random partitionsと言う
  • どんな分布があるかというと、Poisson-Dirichlet Distributionと言う分布(確率分布をパラメトリックに決める規則(law))がある
    • 2パラメタの分布だが、1パラメタで表現できるものが含まれ、それらにいろいろ説明がついている(そのためにかえって混乱してしまうかもしれない)
      • 1パラメタのPoisson-Dirichlet Distributionの方には、2パラメタを1パラメタにする方法によって異なる規則が含まれているが、そのうちのひとつに対応する確率質量分布がEwens's sampling formula(¥alpha=0)であり。2パラメタのPoisson-Dirichlet Distributionの方は、その確率質量分布としてEwens-Pitman formula(distribution)が知られている
    • もちろん、そのほかにも、何かしらの確率過程を考案すれば、それがExchangeable random partitionsを生成することはあるし、それを規則としているということで分布名をつけてもよいが、Ewens-Pitman Distributionは前述のSelf-similarityとかNon-interferenceとか指数型分布族とか言うよい性質を持っていることから取り立てられている(のだと思う)
    • 1変数をやってから2変数をやるのもよいが(ペイパーはそうなっている)、2変数を先に与えて、1変数を眺めなおすほうが(自分には)わかりやすいので、そのようにしてやる
    • Ewens-Pitman sampling formula
      • ある2つのはパラメタ¥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の¥alpha=0の場合
      • 式を比較してみる
      • 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であった。それに対応する分割分布((b_1,...,b_{K});¥sum_{i=1}^{K} b_i = n,ただし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)

    • 確率過程的な意味づけ
      • 式としては上記の通りだが、いくつかの確率過程としての説明がある
      • ¥alpha=0,¥theta > 0
        • この場合は、いわゆるEwens sampling formulaに対応する場合
        • 遺伝学でのいわゆる集団サイズ一定、ランダムメイティング、非抽出サンプリング、無限サイトモデルのときにn染色体のアレルタイプ別観測数の分布を表す
        • それは、中華料理店過程で、新規テーブルに着席する確率を¥frac{¥theta}{¥theta+n}とする場合にも相当する(壷のモデルも同じこと)
        • そして、(-alpha,-alpha,....)というパラメタのディリクレ分布から発生した頻度ベクトルに対して、n標本をサンプリングしたときの多項分布について、多項を第1、第2、…のように項の違いを意識するのではなく、項が揃っていれば第何番目の項かは無視するとしたときの、多項分布を考えるときに、この項数を無限大にしたときが、やはりこれに対応することが示せる。ペイパーでは、そのあたりの式変形を示しているが、ディリクレ分布に観測を上乗せした確率の積分の式で、その中にディリクレ分布が隠れていることとそれを全体に積分すると1になることと、項を区別しないために、設定としての項数kと観察された項数pについてk(k-1)...(k-p+1)の重複分の考慮などに注意すれは式変形を追うことはできる
    • 他方、Ewens-Pitman sampling formulaは2パラメタである
      • これに対応する中華料理店過程を考えると、「具体的な過程」のイメージは湧く
      • すでに客が居るテーブル数がKだとしたとき、新しいテーブルに着く確率を¥frac{¥theta+K¥alpha}{¥theta+n}とし、¥frac{n_i-¥alpha}{¥theta+n}の確率ですでに客が居る席に着くという過程(Pitman-Yor過程)を考える
      • この¥theta,¥alphaに対応する確率質量分布が2パラメタの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. 組みあわせ、代数、数論での位置

メモ

  • ごちゃごちゃと考えたときのメモ
    • 普通のポアソン分布は非負整数を台としてその上に非負実数を配するルールをパラメタ依存に定める。そういう意味で、パラメトリックな確率分布とは、「ある台」に総和(全体の積分)が1になるような確率質量・密度を定めるルールをパラメタ依存に持っているような「確率質量・密度分布というインスタンスの集合」を定めるもの。パラメタの値を特定すれば、それはあるパラメトリックな確率分布の集合の要素を指定することになる。このようにある台の上に確率質量・密度分布を特定する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-20 ポアソン点過程・分割・ノンパラメトリックベイズ

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

  • 動機
    • 色々動機はあるかもしれないが
    • 多数のもの・無限個あるかもしれないもののタイプ分けが興味の対象
    • クラスタ数不定な状況でのクラスタリング
    • そのための確率モデル
    • その確率モデルの下での生起確率・事前確率・尤度・事後確率
  • モデル
    • 具体的な説明から始めよう
      • 中華料理店過程
        • 中華料理店に客が1人ずつ入ってきて着席する
        • 着席には確率的ルールがあると仮定する
        • 客は着席風景を眺め、誰も座っていないテーブルをある確率で選ぶ。その確率はすでに座っている客の数が多いほど小さくなる(¥frac{¥alpha}{¥alpha + n}、ただしn+1-番目の客の場合)
        • 誰も座っていないテーブルに座ることはせず、誰かが座っている席に座るとしたら、座っているテーブルの人数に比例して座る(同じことだが、既に座っている人の中から等確率で1人を選んでその人と同じテーブルに座る)。その確率は、¥frac{n_i}{¥alpha + n}である。ただしn_iはテーブル-iの人数。
        • もちろん¥frac{¥alpha}{¥alpha + n}+¥sum_i ¥frac{n_i}{¥alpha + n} = 1
    • ルール¥frac{¥alpha}{¥alpha + n},¥frac{n_i}{¥alpha + n}を見直す
      • 壺モデル
        • ¥alphaは正の実数として考えているが、これを1とすると、Polya urnモデル。黒いボールが1個からスタートし、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分布を作るときに多数を定めてその総和が1になるようにすることと同じなので、ポアソン・ディリクレ過程と言う(らしい)
    • もう一つ。Size-biased的にでもpoisson-dirichlet的にでもよいが、何かしらの分割をしたとしよう。無限の分割点・無限の線分長ができてしまう。それが目的だから仕方ないが、やはり無限個の要素は扱いにくい。有限個の要素であたかも無限個の要素を扱えるのが、Kingmanのpaintbox。これは単位線分を有限線分に分割し、そのうちの1つを除いたすべての有限長線分は、タイプ別の割合を表し、1つの線分は無限要素を持つ全体集合の中でSingletonになっているものたちの集合であるとする考え方。これをすることで、有限個の、有限割合を有するタイプ(そのタイプには無限個の要素が帰属する)と、無限個のタイプ(ただしすべてSingleton)とに分割することになる。そしてこの一見特殊な無限タイプへの分割が任意のExchangeable random partitionsに1対1対応することがしられている(Kingman's correspondence)ので、これも便利。1つの便利さは、有限分割して、それに応じて無限にタイプ別標本を発生させられること

2018-03-19 Ewens sampling formula

[]Ewens sampling formula

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

2018-03-17 Random Exchangeable Partitions

[][][][]Random Exchangeable Partitions

  • 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

[][][][][]Kingman's theorem、Random Exchangeable Partitions

  • 無限大(N ¥to ¥infty)のRandom Exchangeable Partitions
      • 限大にするとちょっと厄介
  • こんな方法(KingmanのPaintboxの方法)というのがある
    • 1,2,...,Nという数列を長さ1の線分に見立てて、それを分割する
    • ただし、Nは無限大なのでこの線分上には無限個の自然数が並んでいる
    • これを有限個に区切って、そこには線分の長さに比例した数の自然数が帰属するとする。ただし、全体が無限個の自然数に対応するから、有限の長さの線分はどんなに短くても無限個の自然数が対応する
    • このようにすると、無限個の自然数を有限個にしか分割できない
  • 分割数を無限にする
    • 無限個にある分割が有限の長さに対応するとそれを合算すると長さ1の線分に納まらない
    • なので、無限個の分割には無限小の長さを対応付ける
    • これを実現するために、長さ1の線分を有限個に分割し、そのうちの1つは、自然数が1個帰属する小さな分割の集まりとする。それ以外の分割線分には長さに応じた無限個の自然数が帰属するとする
    • こうすることで、無限個の¥frac{1}{¥infty}の長さの線分の集合があることを表現できる
    • Exchangeableであることを考慮すれば、長さ1の線分を有限個に適当に分割し、無限分割用の線分を取り除ける。それ以外の線分は長さの順にならべても一般性を失わない
    • そんな風にして無限個への分割を作成し
    • さらに、無限個の自然数を割り当てることにする。ただし、無限個の自然数の割り当ては無限に続けなくてはならないので、適当な数までランダムに割り付けることにする
# Kingman's paint box

k <- 10
library(MCMCpack)
r <- rdirichlet(1,rep(1,k))
r. <- cumsum(r)

n <- 10^5
u <- runif(n)
my.label <- function(x){
	tmp <- which(r. > x)
	tmp[1]	
}
lb <- unlist(sapply(u,my.label))

plot(r,c(table(lb)))
abline(0,n,col=2)
  • このKingmanのPaintboxの方法は単位線分のタイリングとも称されるが、exchangeable random partitionとKingmanのpaintbox タイリングとは1対1対応することが知られている

[][][][]確率質量分布で表すRandom Exchangeable Partitions

  • Kingmanのpaintbox〜単位線分タイリングがrandom exchangeable partitionsの表現であることがわかったが、実際、どんなrandom exchangeale partitionsが現れるのかは確率事象なので
  • 確率事象として生成されるときにどんな確率事象なのか、そのときにどのような分布になるのか、ということが気にかかる
  • 前の記事では有限個線分の個数を指定し、一様ディリクレ乱数によって線分の長さを定めて生成した。この方法だと、有限個線分の指定個数がかなり大きな制約条件とはなっているが、これも、一つの生成確率事象ルールではある
  • 以下では、いくつかの方法を記述する
  • どんな条件があるかと言えば、¥sum_{i=0}^¥infty p_i=1;p_i ¥ge 0。確率が0より大とすれば確実に無限個に分割。確率が0のものが無限個あるとすれば、有限個への分割もこの定義でカバーされる
  • 自然数・非負整数を大とする確率質量分布は使える
    • 有限個の項に対応する多項分布の確率質量に比例した単位線分分割は、Kingmanのタイリングの無限小分割を持たないバージョンとなる
k <- 10
sort(rdirichlet(1,rep(1,k)))
    • ポアソン分布の確率質量は0,1,...,無限大のそれぞれに正の確率質量が与えられその和は1
plot(dpois(0:100,3),type="h")
plot(sort(dpois(0:100,3)),type="h")
    • ポアソン分布よりも分散が大きい確率質量分布に負の二項分布がある
plot(dnbinom(0:100,3,0.3),type="h")
plot(sort(dnbinom(0:100,3,0.3)),type="h")
    • ディリクレ分布、ポアソン分布、負の二項分布のいずれも、ガンマ分布を使って次のように表現できる
      • k項のディリクレ乱数は、尺度パラメターが1のガンマ分布に従うk個の乱数¥Gamma_1,...,¥Gamma_kを発生し¥frac{¥Gamma_i}{¥sum_{i=1}^k ¥Gamma_i}とすることで生成できる。このとき¥sum_{i=1}^k ¥Gamma_iは有限の値を取ることが知られている。このときのガンマ分布の形状パラメタベクトルがディリクレ分布の形状パラメタになっている(らしい)
      • ポアソン分布の場合は1個のガンマ分布乱数をその期待値として指定することでガンマ分布と紐づく
      • 負の二項分布の場合は、ガンマ乱数を発生させ、その乱数に対応するポアソン分布から乱数を一つ取り出すことを繰り返すとその結果が負の二項分布となると知られています。負の二項分布のsizeとprobとのパラメタとガンマ分布のshape,rateパラメタとの関係はshape = size, rate = prob/(1-prob)です)こちら)
N <- 10^4
a <- c(0.1,0.2,0.3)
Rdirichlet <- rdirichlet(N,a)
Rgamma <- cbind(rgamma(N,a[1],1),rgamma(N,a[2],1),rgamma(N,a[3],1))
Rgamma <- Rgamma/apply(Rgamma,1,sum)
plot(sort(Rdirichlet[,1]),sort(Rgamma[,1]))
N <- 10^5
sz <- 4
pr <- 0.3
sp <- sz
rt <- pr/(1-pr)
Rnb <- rnbinom(N,sz,pr)
RgammaPois <- rep(0,N)
for(i in 1:N){
	tmp <- rgamma(1,sp,rt)
	RgammaPois[i] <- rpois(1,tmp)
}

plot(sort(Rnb),sort(RgammaPois))

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

  • 長さ1の単位線分をあるルールで確率的に分割していけば、それもRandom Exchangeable Partitionsとなる
  • Poisson-Dirichlet Processと呼ばれる方法がその一つで、よく研究されている
  • 何度でも分割し続けるルールとして、単位線分から出発して、分割点を1点とり二分する。生じた2つの線分のうち片方(たとえば1に近い方)はそれ以上いじらないことにして、もう片方をまた二分する。これはいつまでも繰り返せる。この二分点の取り方として、ベータ乱数を使う方法は確率過程として扱いやすく、この方法をPoisson-Dirichlet Processと言う
  • ベータ乱数を発生させる元となるベータ分布としてBeta(1-¥alpha,¥theta+i¥alpha);¥theta ¥ge 0, 0 ¥le ¥alpha < 1,i=1,2,...とするという。ただしベータ乱数の発生に当たり、i回目ごとに用いるベータ分布のパラメタが¥alphaに応じて変化するように設計されていることに注意。また、この定義の中で¥theta=0,¥alpha>0,¥theta>0,¥alpha=0の2つの場合は特殊形として、単純なモデルとして取り上げられることがある
    • また¥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での分割は、長さ1の線分を確率的ルールで分割していくものと読み取れる。Stick-breaking processと称される理由である
    • 次のように中華料理店過程がある意味で(*)同じ分割をもたらす
    • 中華料理店過程ではi+1番目の客が、i番目までの客の着席状態に応じて、新しいテーブルに1人で座るか、既に客の居るテーブルのどれかに座るかを決める、というルールで着席状態が変化していく過程である
    • ¥alpha=0,¥thetaの場合には
      • ¥frac{¥theta}{i+¥theta}の確率で新しいテーブルに着席し、¥frac{m_j}{i + ¥theta}の確率でj番目のテーブルに着席するという確率過程のときに現れる整数分割(とそれに対応する単位線分分割)と同じであると言う。ただしm_jはi人まで着席が終わったときのj番目のテーブルの人数であり、¥sum_{j=1}^k m_j = iである。ここでkはi番目までの客で着席者の居るテーブルの数である
    • ¥alha>0,¥theta=0の場合には
      • ¥frac{k¥alpha}{i}の確率で新しいテーブルに着席し、¥frac{m_j-¥alpha}{i}の確率でj番目のテーブルに着席する
    • この中華料理店過程で得られるランダムな分割はexchangeableなランダム分割になっており、size-biased ordering of the asymptotic block frequeiciesがstick-breakingで得られるものになることが知られている。"size-biased"というのは、「ある誰かしらを特定して、その人が座っているテーブルの人数割合」を調べる場合に使う用語である。そのようにすると、テーブルの人数割合は、割合の大きいテーブルが割合に比例して選ばれるので、そのようにして得られる割合の分布は、割合で重みづいた分布となるのだが、この分割の分野では、この点を考慮することがよくある

[][][][]Poisson 点過程で作る Random Exchangeable Partitions

  • 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は2つのパラメタ¥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と言う