種類数の予測(その１) - 「thinkbayes」第15章

この記事には数式の間違いがあるので注意

thinkbayesの第15章では「belly buttom bacteria」の数を予測する問題を取り扱っている。テキストではシミュレーションによって予測していたが、解析的に解けることが分かったのでメモ。

問題

サンプルから収集したサンプルからシーケンスを調べ、そのデータを分類する。そのデータを元に、全バクテリアはどれくらいの種類があるのか、またそれぞれどれくらいの割合でいるのかを考える。

バクテリア母集団には全部で種類いるとして、それぞれの存在割合をとおく。また、とする。

一方データベースには種類のバクテリアが登録され、それぞれ 回ずつ、合計回観測されている。さらにとおく。、とする。

単純に考えるとバクテリアはデータで得られた分ずつの割合でいると考えることになる。しかし、この考え方ではサンプルできていない未知の種類について考慮していないので発見された種について過剰に、発見されていない種について過小に評価している。また、バクテリアが、同じ1:3:2の割合で現れたと言っても、10, 30, 20回ずつ見つかった場合と10000, 30000, 20000回ずつ見つかった場合では事実の重さが違う。以下では、より正確なバクテリアの種類数と、それぞれの存在する割合を計算する。そのために、それぞれの種類は独立な確率で発見されるという事を仮定し、多項分布にしたがったサンプルが得られると考えることにする。また、種類数もランダムに選択されているものとする。

ディリクレ分布、多項分布について

尤度は $$ \begin{align} P(\mathbf{U}_{k}| \mathbf{F}_{K}) &= &\binom{Q}{u_{1} \cdots u_{k}} \prod_{i} f_{i}^{u_{i}}\\ &= &\dfrac{Q!}{\prod_{i}u_{i}!}\prod_{i} f_{i}^{u_{i}}\\ &= &\mathcal{M}(\mathbf{U}_{k}; \mathbf{F}_{K}) \end{align} $$ より、多項分布。

事後分布は $$ \begin{align} P(\mathbf{F}_{K}| \mathbf{U}_{k}) &\propto &\mathcal{M}(\mathbf{U}_{k}; \mathbf{F}_{K})P(\mathbf{F}_{K})\\ &= &\mathcal{M}(\mathbf{U}_{k}; \mathbf{F}_{K})\mathcal{D}(\mathbf{U}_{k}; \mathbf{1})\\ &\propto &\prod f_{i}^{u_{i}}\prod f_{i}^{\mathbf{1}_{i}-1}\\ &= &\prod f_{i}^{u_{i}+\mathbf{1}_{i}-1} \end{align} $$

より、 $$ P(\mathbf{F}_{K}| \mathbf{U}_{k}) = \dfrac{1}{Z}\prod f_{i}^{v_{i}-1} = \mathcal{D}(\mathbf{F}_{K}; \mathbf{U}_{k}) $$

となり、ディリクレ分布。(ただし)

問題を一度整理してみる

、、がランダムに選択される過程は

↓

↓

の順である。分かっているのはだけなので、を予測しようとすると段階を二つ戻ることになって難しい。「thinkbayes」ではこの問題を、両方についてランダムにサンプリングすることでシミュレーションで解決していた。ここでは解析的に考えることとして、をまず予測し、その分布に従ってを予測する。

1. の分布について

を知りたいので逆の確率を計算していく。 $$ P(\mathbf{U}_{k}|K) = \int_{\mathbf{F}_{K} \in \Omega_{K}} \mathcal{M}(\mathbf{U}_{k}; F_{K}) \mathcal{D}(\mathbf{V}; \mathbf{1}) d\mathbf{F}_{K} $$ ここではの動きうる全範囲とする。これが解析的に求まるならシミュレーションを用いなくても直接を予測できるのでやってみる。

展開して $$ \begin{align} P(U_{k}|K) &= &\int_{\Omega_{K}} \dbinom{Q}{u_{1} \cdots u_{k}} \prod f_{i}^{u_{i}} \dfrac{1}{Z} \prod f_{i}^{u_{i}} d\mathbf{F}_{K}\\ &= &\dfrac{1}{Z} \dbinom{Q}{U_{k}} \int_{\Omega_{K}} \prod f_{i}^{2u_{i}} d\mathbf{F}_{K}\\ &= &\dfrac{1}{Z} \dbinom{Q}{U_{k}} \mathcal{B}(2\mathbf{U}_{k}+\mathbf{1})\int_{\Omega_{K}} \dfrac{1}{\mathcal{B}(2\mathbf{U}_{k}+\mathbf{1})} \prod f_{i}^{2u_{i}+1-1} d\mathbf{F}_{K}\\ &= &\dfrac{1}{\mathcal{B}(\mathbf{V}_{k})} \dbinom{Q}{U_{k}} \mathcal{B}(2\mathbf{U}_{k}+\mathbf{1}) \end{align} $$ ここで尤度比を考えると、$$ \begin{align} R(K_{1}|K_{2}) &= &\dfrac{\mathcal{B}(\mathbf{V}_{K2}) \mathcal{B}(2\mathbf{U}_{K1}+\mathbf{1})} {\mathcal{B}(\mathbf{V}_{K1}) \mathcal{B}(2\mathbf{U}_{K2}+\mathbf{1})}\\ &= &\dfrac{\Gamma(Q+K_{1})\Gamma(2Q+K_{2})} {\Gamma(Q+K_{2})\Gamma(2Q+K_{1})} \end{align} $$

よってについて調べればいいと分かる。

2. 種類の割合について

2.0. ラプラス変換、逆変換の復習

• ラプラス変換
  $$ F(s) = \mathcal{L}(f(t)) \equiv \int_{0}^{\infty} f(t)e^{-st}dt $$

$$ \mathcal{L}(t^{\alpha}) = \dfrac{\Gamma(\alpha + 1)}{s^{\alpha + 1}} $$

• ラプラス逆変換 $$ f(t) = \mathcal{L}^{-1}(F(s)) \equiv \dfrac{1}{2\pi i} \int_{-i\infty}^{i\infty} F(s)e^{st}ds $$

より、 $$ \mathcal{L}(\dfrac{1}{s^{\alpha}}) = \dfrac{t^{\alpha - 1}}{\Gamma(\alpha)} $$

2.1.

以下を固定して考える。

一度にたくさんの種類について考えるのは難しいのでまずはを固定して考える。であり、これを満たしてが動く空間をとおくことにする。 $$ \begin{align} P(f_{1}) &= &\int_{\Omega_{-1}} P(\mathbf{F}_{K}| \mathbf{U}_{k}) d\mathbf{F}_{K}\\ &= &\int_{\Omega_{-1}} \mathcal{D}(\mathbf{F}_{K}; \mathbf{U}_{k}) d\mathbf{F}_{K}\\ &= &\int \cdots \int_{\Omega_{-1}} \dfrac{1}{Z} \prod f_{i}^{u_{i}} df_{2} \cdots df_{K}\\ &= &\int_{0}^{\infty} \cdots \int_{0}^{\infty} \dfrac{1}{Z} \prod f_{i}^{u_{i}} \delta(1-f_{1}-\cdots f(K)) df_{2} \cdots df_{K}\\ &= &\int_{0}^{\infty} \cdots \int_{0}^{\infty} \dfrac{1}{Z} \prod f_{i}^{u_{i}} \Bigl( \int_{-\infty}^{\infty} \dfrac{1}{2\pi} e^{-il(1-\sum f_{i})} dl \Bigr) df_{2} \cdots df_{K}\\ &= &\dfrac{1}{Z} \dfrac{1}{2\pi} \int_{-\infty}^{\infty} e^{-il(1-\sum f_{1})} f_{1}^{u_{1}} \Bigl( \prod \int_{0}^{\infty} f_{i}^{u_{i}}e^{il f_{i}} df_{i} \Bigr) dl \end{align} $$ デルタ関数を用いることで積分区間を簡単にすることができた。

とおくと、 $$ \begin{align} P(f_{1}) &= &\dfrac{1}{Z} \dfrac{1}{2\pi i} \int_{-i\infty}^{i\infty} e^{\kappa(1-f_{1})} f_{1}^{u_{1}} \Bigl( \prod_{j\neq 1} \int_{0}^{\infty} f_{j}^{u_{j}} e^{-\kappa f_{j}} df_{j} \Bigr) d\kappa\\ &= &\dfrac{1}{Z} \dfrac{1}{2\pi i} \int_{-i\infty}^{i\infty} e^{\kappa(1-f_{1})} f_{1}^{u_{1}} \Bigl( \prod_{j\neq 1} \mathcal{L}(f_{j}^{u^{j}}) \Bigr) d\kappa\\ &= &\dfrac{1}{Z} \dfrac{1}{2\pi i} \int_{-i\infty}^{i\infty} e^{\kappa(1-f_{1})} f_{1}^{u_{1}} \prod_{j\neq 1} \dfrac{\Gamma (u_{j}+1)}{\kappa^{u_{j}+1}} d\kappa\\ &= &\dfrac{1}{Z} f_{1}^{u_{1}} \prod_{j\neq 1} \Gamma(u_{j}+1) \dfrac{1}{2\pi i} \int_{-i\infty}^{i\infty} e^{(1-f_{1})\kappa} \dfrac{1}{\kappa^{\sum u_{k}+K-(1+u_{1})}} d\kappa\\ &= &\dfrac{1}{Z} f_{1}^{u_{1}} \prod_{j\neq 1} \Gamma(v_{j}) \dfrac{1}{2\pi i} \int_{-i\infty}^{i\infty} e^{(1-f_{1})\kappa} \dfrac{1}{\kappa^{\sum v_{k}-v_{1}}} d\kappa\\ &= &\dfrac{1}{Z} f_{1}^{u_{1}} \prod_{j\neq 1} \Gamma(v_{j}) \mathcal{L}^{-1} \Bigl(\dfrac{1}{\kappa^{\sum v_{k}-v_{1}}} \Bigr)\\ &= &\dfrac{1}{Z} f_{1}^{v_{1}-1} \prod_{j\neq 1} \Gamma(v_{j}) \dfrac{(1 - f_{1})^{\sum v_{k}-v_{1}-1}}{\Gamma(\sum v_{k}-v_{1})}\\ &= &\dfrac{\Gamma(\sum v_{i})}{\prod \Gamma(v_{i})} f_{1}^{v_{1}-1} \prod_{j\neq 1} \Gamma(v_{j}) \dfrac{(1 - f_{1})^{\sum v_{k}-v_{1}-1}}{\Gamma(\sum v_{k}-v_{1})}\\ &= &\dfrac{\Gamma (\sum v_{i})}{\Gamma(v_{1}) \Gamma(\sum v_{i}-v_{1})} f_{1}^{v_{1}-1} (1-f_{1})^{\sum v_{i}-v_{1}-1} \end{align} $$

ラプラス変換、逆変換を用いることで簡単に積分できた。他の種類についても同様な考えが成り立つので、。

さらに、について考えているので、 $$ f_{i}^{u_{i}}(1-f_{i})^{\sum{u_{j}} - u_{i} + K - 2} $$ についてのみ調べればいいと分かる。

実装

の三次元のベクトルを扱うことになるがそのまま実装してしまうといたる所ループだらけになる。そこで計算される項の次元について先に考察し、計算量を減らす。またouter関数の恩恵にも預かる。

<- (K, 0, 0)

<- (0, 0, F) ** (0, U, 0)
<- (0, 0, F) ** ((0, U, 0) + (K, 0, 0))

(K, 0, 0) * (0, 0, F) ^ (0, U, 0) * (0, 0, F) ^ ((0, U, 0) + (K, 0, 0))

= (0, 0, F) ^ (0, U, 0) * {(K, 0, 0) * (0, 0, F) ^ ((0, U, 0) + (K, 0, 0))}

第3項の計算が厄介であるがこれはUでループを回すことにする。第1項にはKの次元が含まれておらず、第2項でKについての和を取ってしまってから足し合わせる。

誤差を防ぐための正規化はFの要素について変化がなければ、どのタイミングで行なってもいい。

U <- c(10, 30, 60, 1)
U <- c(U, 0)
Fs <- seq(0.0000001, 1-0.0000001, length.out=100)
Ks <- seq(length(U) - 1, 100)

# 第1項
pmf.1 <- outer(U, Fs, function(x,y) y^x)

lpmf.Ks <- lgamma(sum(U) + Ks) - lgamma(2 * sum(U) + Ks)
# 第2項
pmf.2 <- matrix(nrow=length(U), ncol=length(Fs))
for (i in seq(length(U))){
  #↓KsとFsの要素数が同じになると正しい結果が得られないのに注意
  tmp <- lpmf.Ks + outer((sum(U) - U[i] + Ks - 2), log(1 - Fs))
  # max(tmp)でも良いが精度の維持のため
  tmp <- exp(tmp - apply(tmp, 1, max))
  pmf.2[i,] <- apply(tmp, 2, sum)
}

pmf <- pmf.1 * pmf.2
# Uごとのpmfの和が1になるように調整
pmf <- pmf / apply(pmf, 1, sum)
cmf <- t(apply(pmf, 1, cumsum))

# 概形の確認
par(mfcol=c(2,1))
matplot(t(pmf), type='l', lty=1)
matplot(t(cmf), type='l', lty=1)

# プロット
png('estimaton.png')
for (i in seq(nrow(pmf))){
  plot(Fs, pmf[i,], xlim=c(0, 1), ylim=c(0, 1), ylab="", type='l')
  par(new=T)
}
title(ylab='pmf')
for (i in seq(nrow(pmf))){
  plot(Fs, cmf[i,], xlim=c(0, 1), ylim=c(0, 1), ylab="", type='l')
  par(new=T)
}
title(ylab='cmf')dev.off()

続きはその2へ。