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

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

2018-07-07 ぱらぱらめくる『Mathematical Theory of Bayesian Statistics』

[][][]ぱらぱらめくる『Mathematical Theory of Bayesian Statistics』

  • Prefaceより
    • (1) 3つ組(真の分布、統計モデル、事前分布)の間の関係を任意の場合について取り扱えるような枠組みが欲しい
    • (2) それにより、真の分布が不明な時にも、統計モデルと事前分布とのペアの適切さについて評価できるようにしたい
    • (3) さらにその結果、統計モデルと事前分布ペアとのペアの最適なものを選んでベイズ推定に使えるようにしたい
    • 従来は、事後分布に正規分布を仮定することによって、どんな統計モデルにどんな事前分布を取るかを決めていたが、もっと現実に即したものにしたい
  • 目次
  • 細目次
    • 1 ベイズ統計の基礎
      • 1.1 ベイズ統計
      • 1.2 確率分布
      • 1.3 真の分布
      • 1.4 モデル、事前分布、事後分布
      • 1.5 事後分布の例
      • 1.6 推定と一般化("Generalization loss")
      • 1.7 Marginal LikelihoodとPartition Function
      • 1.8 条件付きで独立な場合
    • 2 統計モデル
      • 2.1 正規分布
      • 2.2 多項分布
      • 2.3 線形回帰
      • 2.4 ニューラルネットワーク
      • 2.5 有限な混合正規分布
      • 2.6 ノンパラメトリックな(無限な?)混合正規分布
    • 3 観測量のベイズにおける基礎的な式
      • 3.1 真とモデルとの関係の定式化
      • 3.2 Normalized Observables
      • 3.3 キュムラント母関数
      • 3.4 基本的なベイズ理論
    • 4 真の分布と事後分布とが合っている場合
      • 4.1 Partitio FunctionのDivision
      • 4.2 漸近的な自由エネルギー
      • 4.3 漸近的な損失
      • 4.4 Asymptotic Expansionsの証明
      • 4.5 点推定
    • 5 事後分布が正規分布の場合
      • 5.1 標準的な式
      • 5.2 State Density Function
      • 5.3 Asymptotic Free Energy
      • 5.4 Renormalized Posterior Distribution
      • 5.5 条件付きで独立な場合
    • 6 一般化した場合
      • 6.1 ベイズ分解(Bayesian Decomposition)
      • 6.2 Resolution o Singularities
      • 6.3 General Asymptotic Theory
      • 6.4 Maximum A Posteriori Method
    • 7 MCMC
      • 7.1 メトロポリス法
      • 7.2 Gibbsサンプラー
      • 7.3 Numerical Approximation of Observables
    • 8 情報量基準
      • 8.1 モデル選択
      • 8.2 ハイパーパラメタの最適化
    • 9 各論
      • 9.1 最適とは
      • 9.2 ベイズによる仮説検定
      • 9.3 ベイズによるモデル比較
      • 9.4 Phase Transition
      • 9.5 Discovery Process
      • 9.6 階層ベイズ
    • 10 確率論の基礎事項
      • 10.1 デルタ関数
      • 10.2 Kullback-Leibler Divergence
      • 10.3 Probability State
      • 10.4 Empirical Process
      • 10.5 Convergence of Expected Values
      • 10.6 ディリクレ過程の混合

[][][]1 ベイズ統計の基礎 ぱらぱらめくる『Mathematical Theory of Bayesian Statistics』

  • 真の分布の事前分布があって、それぞれがある観測をもたらす確率がわかっていれば、観測から真の分布の事後分布が計算できる、というのがベイズの定理
  • 残念ながら、事前分布も正確には解らないし、観測をもたらす確率も解っていない(のでモデルを入れる)から、答えがきっちり決まらない〜 ill-posed problem
  • 仕方がないので、方法を選んで、結果を出してみて、方法について検討する、という手順を踏んで、何かしら正しそうな推定結果に至ることになる
  • 真の分布がわからなければ、これを頑張ってもうまく行かない様に見えるが、そうではない
  • いかなる真分布・統計モデル・事前分布のトリオに対しても、ある数学的な決まりがあって、それを使うと、統計モデルと事前分布のペアの取り方の適切さについて検討できる
  • 本書はその話
  • Training Loss とGeneralization Loss
    • 標本に事前分布とモデルとをあてはめて分布推定をしたら、その推定の良さが知りたい(q(x)は真の分布、p(x|X^n)は推定分布。
      • T_n = -¥frac{1}{n}¥log{p(X_i|X^n): Training Loss
      • G_n = -¥int q(x) ¥log{p(x|X^n) dx: Generalization Loss
      • G_n-S = KL(q(x)||p(x|X^n)), Sは真の分布のエントロピー。Generalization Loss からエントロピーを引いたものが、真と推定とのKL divergence。Generalization Error
      • C_n = -¥frac{1}{n} ¥sum_{i=1}^n ¥log(p(X_i|(X^n -Xi)): Cross Varidation Loss (Leave-one-out type)
        • C_n - S_n: Cross Validation Error
      • Widely Applicabe Information Criterion (WAIC)
        • W_n = T_n +¥frac{1}{n} ¥sum_{i=1}^n V_n¥[¥log(p(X_i|w)¥], V_w¥[¥]はPosterior Variane
        • W_n - S_n: WAIC error
        • W_nは漸近的にCross Validation Lossに相当するという
    • Errors
      • G_n-S,C_n-S_n,W_n-S_nが3つのエラー
      • S,S_nともモデル事前分布によらない
      • Losses であるG_n,C_n,W_nよりも、エラーであるG_n-S,C_n-S_n,W_n-S_nの方が分散が小さいので、エラーを比較することも多い
    • Marginal Likelihood or Partition Function. パラメタで表された事前分布と、パラメタによって決まる個々の変数の条件付き確率分布との積に分離する
    • 条件付き独立な場合に、Cross Validation Loss ではうまく行かないことが、Widely Applicable Information Criterion(WAIC)ではうまく行くことが示せる

[][][]3 Basic Formula of Bayesian Observables 観測量のベイズにおける基礎的な式 ぱらぱらめくる『Mathematical Theory of Bayesian Statistics』

  • (1) 真の分布と統計モデルの関係の定義をする
  • (2) 観測変数のベイズ流定義とのそのnormalized version
  • (3) ベイズ予測のキュムラント母関数の定義
  • (4) ベイズ統計の基礎的な理論事項をキュムラント母関数を用いて証明する(それにより、真の分布とモデルと予測推定の数学的基礎を定める)
  • Model は Truthの近くをカバーするけれど、必ずしもModelにTruthが含まれるわけではない
  • 各種エラー・ロス関数等は、最適解を基準にすると、パラメータ最適解を含まない関数になる。それがnormalized 関数。average and empirical log loss functionを考えることにより、ベイズ流観測変数とnormalized 観測変数とは単純な関係になる
  • キュムラント母関数を作るのだそうだ
  • 級数展開することで、各種関数の値の近似値とその精度に関する情報が解析的に示せる

[][][]5 Standard Posterior Distribution 事後分布が正規分布の場合 ぱらぱらめくる『Mathematical Theory of Bayesian Statistics』

  • モデルが真の分布を含み、推定にあたり、漸近的に正規分布を想定できる場合は簡単。ベイズ推定と最尤推定とがほぼ一致するから
  • 真の分布は正規分布とはみなせないことも多く、その場合は、ベイズ推定のあてはまりがよくなる
  • p(w) ¥propto exp(-n w_1^{2k_1} w_2^{2k_2} ... w_d^{2k_d})と一般化できる(そうだ)
  • この章では、この表現ができるならば、free energy, generalization lossについて有用なことが証明できることを示す
  • これは、ランダム標本から、その母分布を推定するにあたり、ある確率モデルがwでパラメタライズされているときに、その推定にあたり、wの事後分布の形はこのような形(中心があって、2k_i乗的に減衰するようになっていて、しかも、標本数nによって、その事後分布の集中の良さが指定できる、という式である

[][][]6 General Posterior Distribution 一般化した場合 ぱらぱらめくる『Mathematical Theory of Bayesian Statistics』

  • Algebraic geometric transformにより、standard form(5章の形)に持ち込めることを示す
  • p(w) = ¥sum Standard formとなれば、5章で示したStandard formの色々な良い性質がそのまま使える(らしい)

2018-07-06 ぱらぱらめくる『Mathematical Theory of Bayesian Statistics』

[][][]ぱらぱらめくる『Mathematical Theory of Bayesian Statistics』

  • Prefaceより
    • (1) 3つ組(真の分布、統計モデル、事前分布)の間の関係を任意の場合について取り扱えるような枠組みが欲しい
    • (2) それにより、真の分布が不明な時にも、統計モデルと事前分布とのペアの適切さについて評価できるようにしたい
    • (3) さらにその結果、統計モデルと事前分布ペアとのペアの最適なものを選んでベイズ推定に使えるようにしたい
    • 従来は、事後分布に正規分布を仮定することによって、どんな統計モデルにどんな事前分布を取るかを決めていたが、もっと現実に即したものにしたい
  • 目次
    • 1 Definition of Bayesian Statistics
    • 2 Statistical Models
    • 3 Basic Formula of Bayesian Observables
    • 4 Regular Posterior Distribution
    • 5 Standard Posterior Distribution
    • 6 General Posterior Distribution
    • 7 Markov Chain Monte Carlo
    • 8 Information Criteria
    • 9 Topics in Bayesian Statistics
    • 10 Basic Probability Theory

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つの便利さは、有限分割して、それに応じて無限にタイプ別標本を発生させられること