1. はじめに
ベータ分布は範囲が固定されている連続データの分布を柔軟に表現できるため便利である。 非心ベータ分布 (Noncentral Beta Distribution) はベータ分布を一般化した確率分布であり柔軟性はさらに上がる。
R ではベータ分布に従う乱数を生成する rbeta()
関数に ncp
引数1を指定することで非心ベータ分布に従う乱数を生成できる。
x <- rbeta(2000, 10, 5, ncp = 3) hist(x, probability = TRUE) curve(dbeta(x, 10, 5), add = TRUE, col = "red") curve(dbeta(x, 10, 5, ncp = 3), add = TRUE, col = "blue")
赤線がベータ分布、青線が非心ベータ分布であり、ヒストグラムへの当てはまりは青線の方が良い。
非心ベータ分布はベータ分布の中心をずらした分布になる。
右にずらす場合は Type I、左にずらす場合は Type II と呼ばれる。
rbeta()
は Type I にしか対応していないが、Type I と Type II は単に左右反転しただけなので 1 - x
を取ることで Type II を表せる。
本記事では Stan を使って非心ベータ分布のパラメータを推定する方法について述べる。
2. 対数尤度関数
Stan には非心ベータ分布を推定するための組み込み関数はない。
したがって、target +=
記法を用いてモデルを記述する。
そのために、非心ベータ分布の対数尤度関数を求める必要がある。
定義より、非心ベータ分布の密度関数は次である。
\begin{align} \mathrm{NCBeta}(x \mid \alpha, \beta, \lambda) = \sum_{j = 0}^\infty \mathrm{Pois}(j \mid \lambda) \mathrm{Beta}(x \mid \alpha + j, \beta) \end{align}
ただし、Pois はポアソン分布の密度関数、Beta はベータ分布の密度関数である。
ここから尤度関数を求めたいが、 についての無限和の存在が問題になる。 ここで、ベータ分布の密度関数を に対するオーダーで書き直すと
\begin{align} \mathrm{Beta}(x \mid \alpha + j, \beta) &= \frac{x^{\alpha+j-1}(1-x)^{\beta-1}}{\mathrm{B}(\alpha + j, \beta)} \\ &= \frac{O(x^j)}{O\left(\frac{(\alpha+j-1)!(\beta - 1)!}{(\alpha + j + \beta - 1)!}\right)} \\ &= \frac{O(x^j)}{O\left(\frac{1}{j^\beta}\right)} \\ &= O(x^j j^\beta) \end{align}
となる。ただし B はベータ関数であり、そのオーダーとして と が整数の場合を考えた。
一方、ポアソン分布の密度関数の に対するオーダーは次のようになる。
\begin{align} \mathrm{Pois}(j \mid \lambda) &= \frac{\left(\lambda/2\right)^j}{j!} e^{-\lambda} \\ &= \frac{O(\lambda^j)}{O(j!)} \end{align}
すなわち、ベータ分布は のときでも高々 でしか増加しないのに対し、ポアソン分布は の速度で減少する。
したがって、非心ベータ分布の密度関数の に関して次のことが言える。
- が十分大きいとき項は 0 とみなせる。
- ポアソン分布の重みが支配的である。
すなわち、0 とみなせる についてはポアソン分布の重みのみを考えれば良い。
ポアソン分布の相対的な重みを比較するには次のようにする。
lambda <- 7 j <- rpois(10000, lambda / 2) round(table(j) / max(table(j)), 2)
j 0 1 2 3 4 5 6 7 8 9 10 11 12 13 0.14 0.49 0.81 1.00 0.82 0.58 0.37 0.18 0.08 0.03 0.01 0.00 0.00 0.00
これより、 に限れば、 は 0 から 10 までを考えれば十分である。 の上限を上げると の範囲も広がり、そのぶん計算にかかる時間は増加することに注意。
以上によりベータ分布の対数尤度関数は次のように表される。
\begin{align} \log\mathrm{NCBeta}(x \mid \alpha, \beta, \lambda) &= \log\left( \sum_{j = 0}^{10} \mathrm{Pois}(j \mid \lambda) \mathrm{Beta}(x \mid \alpha + j, \beta) \right) \\ &= \log \sum_{j=0}^{10} \exp\left( \log \mathrm{Pois}(j \mid \lambda) + \log \mathrm{Beta}(x \mid \alpha + j, \beta) \right) \end{align}
これは log_sum_exp の形になっていることに注意しよう。
3. Stan による実装
非心ベータ分布の対数尤度関数がわかったので、これを Stan で実装すると次のようになる。
data { int<lower=1> N; real<lower=0, upper=1> values[N]; } parameters { real<lower=1> alpha; real<lower=1> beta; real<lower=0, upper=7> lambda; } model { real logliks[11]; for (i in 1:N) { for (j in 0:10) { logliks[j + 1] = poisson_lpmf(j | lambda / 2); logliks[j + 1] += beta_lpdf(values[i] | alpha + j, beta); } target += log_sum_exp(logliks); } }
これを ncbeta01.stan
という名前でファイルに保存する。
推定を実行する R コードは以下である。
set.seed(314) alpha <- 10 beta <- 5 ncp <- 3 N <- 200 x <- rbeta(N, alpha, beta, ncp) standata <- list(N = N, values = x) # Type I # standata <- list(N = N, values = 1 - x) # Type II library(rstan) # options(mc.cores = 3) model <- stan_model("ncbeta01.stan") fit <- sampling(model, data = standata, chains = 3) print(fit)
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat alpha 8.93 0.03 1.15 6.70 8.18 8.91 9.70 11.29 1160 1 beta 4.57 0.01 0.42 3.80 4.28 4.56 4.86 5.43 1491 1 lambda 3.17 0.06 1.97 0.17 1.45 3.04 4.81 6.73 1255 1 lp__ 150.20 0.04 1.25 147.12 149.60 150.52 151.13 151.65 928 1
真の値 , , がうまく推定できたようである。
まとめ
非心ベータ分布を Stan で推定する方法を述べた。 上記のコードだと私の PC では 1 chain 80秒程度かかる。
次回はこれを高速化して 1 chain 10秒程度にする方法を書こうと思う。
StanとRでベイズ統計モデリング (Wonderful R)
- 作者:松浦 健太郎
- 出版社/メーカー: 共立出版
- 発売日: 2016/10/25
- メディア: 単行本
-
npc
は Non-Centrality Parameter の略↩