ほくそ笑む

R言語と統計解析について

非心ベータ分布を Stan で推定する (1)

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")

f:id:hoxo_m:20180823223727p:plain

赤線がベータ分布、青線が非心ベータ分布であり、ヒストグラムへの当てはまりは青線の方が良い。

非心ベータ分布はベータ分布の中心をずらした分布になる。 右にずらす場合は 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 はベータ分布の密度関数である。

ここから尤度関数を求めたいが、 j についての無限和の存在が問題になる。 ここで、ベータ分布の密度関数を  j に対するオーダーで書き直すと

\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 はベータ関数であり、そのオーダーとして  \alpha \beta が整数の場合を考えた。

一方、ポアソン分布の密度関数の  j に対するオーダーは次のようになる。

\begin{align} \mathrm{Pois}(j \mid \lambda) &= \frac{\left(\lambda/2\right)^j}{j!} e^{-\lambda} \\ &= \frac{O(\lambda^j)}{O(j!)} \end{align}

すなわち、ベータ分布は  x = 1 のときでも高々  O(j^\beta) でしか増加しないのに対し、ポアソン分布は  O(j!) の速度で減少する。

したがって、非心ベータ分布の密度関数の  j に関して次のことが言える。

  •  j が十分大きいとき項は 0 とみなせる。
  • ポアソン分布の重みが支配的である。

すなわち、0 とみなせる  j についてはポアソン分布の重みのみを考えれば良い。

ポアソン分布の相対的な重みを比較するには次のようにする。

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 

これより、 \lambda \leq 7 に限れば、 j は 0 から 10 までを考えれば十分である。  \lambda の上限を上げると  j の範囲も広がり、そのぶん計算にかかる時間は増加することに注意。

以上によりベータ分布の対数尤度関数は次のように表される。

\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

真の値  \alpha = 10,  \beta = 5,  \lambda = 3 がうまく推定できたようである。

まとめ

非心ベータ分布を Stan で推定する方法を述べた。 上記のコードだと私の PC では 1 chain 80秒程度かかる。

次回はこれを高速化して 1 chain 10秒程度にする方法を書こうと思う。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)


  1. npc は Non-Centrality Parameter の略