こんにちは。mutopsyです。この記事では，同一のカウントデータを二項分布でモデリングした場合と幾何分布でモデリングしたときに，推定される確率パラメータの事後分布が一致することを説明します。間違いがあったらご指摘ください。
次の2つの例を考えてみましょう。
- N箱のチョコボールを購入したら銀のエンゼルが1枚当たった。銀のエンゼルが当たる確率θ1を推定せよ。
- 銀のエンゼルが当たるまでチョコボールを購入したらN箱目で当たった。銀のエンゼルが当たる確率θ2を推定せよ。
θ1とθ2はどちらも試行にかかわらず一定の値をとり，かつ各試行は独立であると仮定します。ここではチョコボールを例に挙げていますが，2つの事象 (e.g., 当たるか外れるか) のカウントデータであれば何でも構いません。どちらの文も，「銀のエンゼルがN箱中1箱当たった」という同じ結果を述べていますが，データの収集方法に違いがあります。1番目のケースではNがあらかじめ決まっていた状況で銀のエンゼルが1回当たったという結果が得られていますが，2番目のケースでは銀のエンゼルが1回当たることがあらかじめ決まっていた状況でN箱を要したという結果が得られているという違いです。この違いは確率パラメータの推定結果に影響を与えるのでしょうか。
この例のようにNの決め方 (決まり方) が異なるケースでは，異なる統計モデルを想定するのが自然です。1番目の例のように，Nが所与のときに銀のエンゼルが当たった枚数 (K = 1) のデータが得られた場合には，二項分布を用いて
と表現することができます。一方，2番目の例のように，銀のエンゼルが当たった枚数 (K = 1) が所与のときにNというデータが得られた場合には，幾何分布を用いて
と表現できます。幾何分布は，同一の確率パラメータを持つベルヌーイ分布に従う試行を独立に繰り返し行ったときに初めて成功する (今回の例では銀のエンゼルが当たる) のに要する回数の分布です。二項分布では銀のエンゼルが当たった回数 (K = 1) が説明されるのに対し，幾何分布では銀のエンゼルが当たるまでに要した試行回数 (N) が説明されるのがポイントです。
では，2つのモデルでθ1とθ2の推定結果がどうなるのかを確認してみましょう。今回は，N = {10, 20, 50, 100} の場合のθ1とθ2をStanでベイズ推定してみます。4通りのN と2つのモデル (Is_Censored = 0なら二項分布モデル，Is_Censored = 1なら幾何分布モデル) の組み合わせを用いて8個のθをまとめて推定します。Stanコードは次の通り。rstan2.18.2，Rtools3.5，R3.5.2，Windows10で実行しました。
Stanコード01: 二項分布モデルと幾何分布モデルの比較
data {
int N_ID; //標本の数 = 8
int N[N_ID]; //試行回数 = {10,20,50,100,10,20,50,100}
int K[N_ID]; //成功回数 (今回は1で固定；2以上の場合は負の二項分布になる)
int Is_Censored[N_ID]; //打ち切られたデータか否か = {0,0,0,0,1,1,1,1}
}
parameters{
real<lower=0,upper=1> theta[N_ID]; //当たる確率
}
model{
for (i in 1:N_ID){
if (Is_Censored[i] == 0){
K[i] ~ binomial(N[i],theta[i]); //二項分布
} else if (Is_Censored[i] == 1){
N[i]-K[i] ~ neg_binomial(K[i],theta[i]/(1-theta[i])); //負の二項分布 (K = 1のとき幾何分布と一致)
}
}
}
幾何分布は負の二項分布の特殊なケースなのでモデル上は負の二項分布を使って表現しています (現時点ではStanには幾何分布用の関数はありません)。neg_binomialの第2引数にはθそのものではなくθのオッズを与える必要があります。また，左辺にはNではなくN-1を与えます (Stanで使われている負の二項分布の定義に合わせるため)。このモデルを走らせるためのRコードは以下の通り。
Rコード01: MCMCの実行
library(rstan)
N_ID <- 8
N <- c(10,20,50,100,10,20,50,100)
K <- rep(1,N_ID)
Is_Censored <- c(0,0,0,0,1,1,1,1)
#コンパイル
stanmodel <- stan_model(file="binomgeom.stan")
#データの準備
data <- list(N_ID = N_ID, N = N,
K = K, Is_Censored = Is_Censored)
#サンプリング
fit <- sampling(stanmodel, data=data, seed=123,
chains=4, iter=250500, warmup=500, thin=1)
各々の場合のθの事後分布とMAP推定値は以下の通り。
図1. 幾何分布モデルと二項分布モデルで推定したθの事後分布とそのMAP推定値。
二項分布モデルで推定したθ1を半透明の赤色，幾何分布モデルで推定したθ2を半透明の青色で描画していますが，ほとんど重なっていて紫色にしか見えません。MAPも，多少の誤差はありますがほぼ一致しています。どうやら，どちらのモデルを使っても推定結果は一致ししそうです。
両者が一致することは，事後分布を解析的に求めれば分かります。まず，二項分布の確率質量関数は，Nを所与として
と表されます。K = 1のとき，
となります。データが1つしかないときのθ1の尤度関数は
と表せます。最初のNはθ1とは無関係な定数なので消去しています。最後に，θ1について無情報事前分布を仮定して尤度関数を正規化することで，次の事後分布が得られます。
次に幾何分布ですが，確率質量関数は
と表されます。二項分布モデルの場合と同じように，これをθ2の関数 (尤度関数) とみなして正規化すると，
という事後分布が得られます。先ほどの二項分布モデルの事後分布と一致していることが確認できると思います。実はこの事後分布はどちらも (α, β) = (2, N) のベータ分布と一致します。二項分布と幾何分布の確率質量関数の違いは組み合わせ (C) を掛けるか否かの違いのみで，この部分は確率パラメータに依存しないため，事後分布を計算する過程で消去されるという訳です。ちなみに，データが複数ある場合やKが2以上の場合 (つまり幾何分布ではなく負の二項分布を使う場合) にも事後分布は同じパラメータを持つベータ分布 (α = 1 + 成功回数，β = 1 + 失敗回数) となり一致します (計算は省略)。
という訳で，二項分布モデルを使っても幾何分布モデルを使っても確率パラメータの推定結果は一致することが分かりました。とはいえ，データの持つ情報 (データの収集方法の違い) を正確に表現にするためにもモデル上は両者を区別したほうが良さそうです。
