読者です 読者をやめる 読者になる 読者になる

tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

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

thinkbayes 統計 R

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

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

問題

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

バクテリア母集団には全部で$K$種類いるとして、それぞれの存在割合を$f_{i} (1 \leq i \leq K)$とおく。また、$\mathbf{F}_{K} = (f_{i})$とする。

一方データベースには$k$種類のバクテリアが登録され、それぞれ$u_{i}$ 回ずつ、合計$Q$回観測されている。さらに$v_{i} = u_{i} + 1$とおく。$\mathbf{U}_{k} = (u_{i})$、$\mathbf{V_{k}} = (v_{i})$とする。

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

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

尤度は $$ \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}) $$

となり、ディリクレ分布。(ただし$Z=\mathcal{B}(\mathbf{V}) = \dfrac{\prod_{i}\Gamma(v_{i})}{\Gamma(\sum v_{i})}$)

問題を一度整理してみる

$K$、$\mathbf{F}_{K}$、$\mathbf{U}_{k}$がランダムに選択される過程は

$K$

$\mathbf{F}_{K} \sim \mathcal{D}(\mathbf{F}_{k}; \mathbf{1})$

$\mathbf{U}_{k} \sim \mathcal{D}(\mathbf{U}_{k}; \mathbf{F}_{K})$

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

1. $K$の分布について

$P(K| \mathbf{U}_{k})$を知りたいので逆の確率を計算していく。 $$ 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} $$ ここで$\Omega_{K}$は$\mathbf{F}_{K}$の動きうる全範囲とする。これが解析的に求まるならシミュレーションを用いなくても直接$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} $$ ここで尤度比$R(K_{1}|K_{2})$を考えると、$$ \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} $$

よって$\dfrac{\Gamma(Q+K)}{\Gamma(2Q+K)}$について調べればいいと分かる。

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 $$

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

2.1.

以下$K$を固定して考える。

一度にたくさんの種類について考えるのは難しいのでまずは$f_{1}$を固定して考える。$f_{2}+ \cdots +f_{K} = 1-f_{1}$であり、これを満たして$f_{2} \cdots f_{K}$が動く空間を$\Omega_{-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} $$ デルタ関数を用いることで積分区間を簡単にすることができた。

$il = -\kappa$とおくと、 $$ \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} $$

ラプラス変換、逆変換を用いることで簡単に積分できた。他の種類についても同様な考えが成り立つので、$\dfrac{\Gamma (\sum v_{j})}{\Gamma (v_{i})\Gamma(\sum v_{j}-v_{i})} f_{i}^{v_{i}-1} (1-f_{i})^{\sum v_{j}-v_{i}-1}$。

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

実装

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

$\dfrac{\Gamma(Q+K)}{\Gamma(2Q+K)}$ <- (K, 0, 0)

$f_{i}^{u_{i}}$ <- (0, 0, F) ** (0, U, 0)
$(1-f_{i})^{\sum{u_{j}} - u_{i} + K - 2}$ <- (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()

f:id:kuyata:20140707195815j:plain

続きはその2へ。