ベータ分布

\(\alpha\) と \(\beta\) のスライダーを動かしてベータ分布を変化させよう。

定義

ベータ分布 (beta distribution) は以下の式で表される連続確率分布。確率変数は \(0 \lt x \lt 1\) の範囲を取る。\[ P(x; \alpha, \beta) = \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)}, \,\,\, \alpha \gt 0, \beta\gt 0 \] ここで \(B(\alpha,\beta)\) はベータ関数である。\[ B(\alpha,\beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} \]

\(x\) の区間に 0 と 1 を含むか (つまり閉区間か開区間か) については論議があり、ベータ分布は \(\alpha, \beta \lt 0\) において発散することから区間に含めない立場を取ることもある。実数演算においてはオーバーフローを防ぐ目的で \(0 \lt x \lt 1\) とする。

アルゴリズム

確率関数

ベータ分布の確率関数は指数とガンマ関数を含むため \(\alpha\), \(\beta\) の値が極端に大きいときに数値演算上のオーバーフローを起こしやすため一度対数化してから指数演算を行う。\[ \begin{align*} P(x) & = {\rm exp}\left[\log\left\{\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)}\right\}\right] \\ & = {\rm exp}\left\{(\alpha-1)\log x + (\beta-1)\log(1-x) \right. \\ & ~~~~~ \left. + \log\Gamma(\alpha) + \log\Gamma(\beta) - \log\Gamma(\alpha+\beta)\right\} \end{align*} \] ただし \(x = 0, 1\) においてそれぞれ \(\log x, \log (1-x)\) が発散することからコーナーケースを考慮しなければならない。\(x\) を実数全域で考えたとき以下の条件に分けられる。\[ P(x) = \left\{ \begin{array}{ll} {\rm exp}\left\{\log P(x)\right\} & (0 \lt x \lt 1) \\ 0 & ({\rm otherwise}) \end{array} \right. \]

対数関数 log(Double.MaxPositiveMin) の結果はオーバーフローしないことから機械イプシロンは考慮しなくても良いため実数対実数の同値比較 == で判断しても問題はない。

ここで gammaln() はガンマ関数で実装した対数ガンマ関数である。

def beta(alpha:Double, beta:Double)(x:Double):Double = {
  if(x <= 0.0 || x >= 1.0) 0.0 else math.exp(
    (alpha - 1.0) * math.log(x) + (beta - 1.0) * math.log(1.0 - x)
      + gammaln(alpha) + gammaln(beta) - gammaln(alpha + beta)
  )
}

以下に REPL での実行結果を記す。

scala> for(i <- 0 to 10) println(f"$i%.1f: ${beta(1,1)(i * 0.1)}%f")
0.0: 0.000000
1.0: 1.000000
2.0: 1.000000
3.0: 1.000000
4.0: 1.000000
5.0: 1.000000
6.0: 1.000000
7.0: 1.000000
8.0: 1.000000
9.0: 1.000000
10.0: 0.000000

scala> for(i <- 0 to 10) println(f"$i%.1f: ${beta(2,2)(i * 0.1)}%f")
0.0: 0.000000
1.0: 0.015000
2.0: 0.026667
3.0: 0.035000
4.0: 0.040000
5.0: 0.041667
6.0: 0.040000
7.0: 0.035000
8.0: 0.026667
9.0: 0.015000
10.0: 0.000000

scala> for(i <- 0 to 10) println(f"$i%.1f: ${beta(0.5,0.5)(i * 0.1)}%f")
0.0: 0
1.0: 10.471976
2.0: 7.853982
3.0: 6.855517
4.0: 6.412749
5.0: 6.283185
6.0: 6.412749
7.0: 6.855517
8.0: 7.853982
9.0: 10.471976
10.0: 0