はじめに
こんにちわ、研究開発部の島田です。今回は統計的学習で基本的なサンプリングアルゴリズムを一つ紹介します。
統計的学習におけるサンプリング手法はいくつかありますが、大別するとマルコフ連鎖モンテカルロ法(MCMC)を使わないサンプリング手法とMCMCを使ったサンプリング手法に分けられます。
MCMCを使わないサンプリング手法の最も基本的なアルゴリズムとして、棄却サンプリングがあげられます。このアルゴリズムは難しい数式を必要としないので、直感的にもわかりやすいです。
今回はこの棄却サンプリングについて簡単な説明と実装を行います。
棄却サンプリング
ベイズ統計でやりたいことは、ある複雑な事後分布からモンテカルロ法を使ってサンプリングをしたいということになります。
そして棄却サンプリングでは、直接事後分布からサンプリングすることは難しいのでもっと単純な分布(これを提案分布と言います)を利用しようという手法です。
ここで、予測分布を、提案分布
とすると、
上の条件を満たすおよび
を求める必要があります。この条件は予測分布が提案分布を
倍した分布に覆われていることを意味しています。
棄却サンプリングの手順は非常にシンプルで、下記のようになります。
- 提案分布
に従う乱数
を発生
]に従う乱数
を発生
の場合は受理。そうでない場合は棄却。
- 1から3を
回繰り返す。
ベータ分布
今回は予測分布をベータ分布として棄却サンプリングを試してみます。 ベータ分布とは、確率密度関数が下記の式で表される確率分布のことです。
ただし、は正の実数のパラメータです。
また、ベータ分布は区間[0, 1]上の連続型の確率分布であることが特徴です。
では、ベータ分布の形を確認するためにScipyというPythonライブラリを使って確認してみます。
まずは必要なライブラリのインポートです。
import numpy as np import matplotlib.pyplot as plt from scipy import stats from scipy import optimize as opt from scipy.stats import beta, uniform
次にベータ分布を描画してみます。必要なコードはたったこれだけです。
np.random.seed() a, b = 1.8, 2.9 x = np.linspace(beta.ppf(0.001, a, b), beta.ppf(0.999, a, b), 100) plt.plot(x, beta.pdf(x, a, b))
棄却サンプリングの実行
では、いよいよ棄却サンプリングを試してみます。
提案分布はどんな分布であっても構わないのですが、今回はサンプリングが簡単で理解もしやすい一様分布とします。
まずは定数を求める必要があるのですが、下の図からわかるように
はベータ分布の最大値
と求めることが出来ます。
そして、横軸方向について一様分布からランダムな値
を決めます。具体的には0から1の間で一つ値を決定することになります。
次に縦軸方向については]の一様分布に従ってランダムな値
を決めます。
この処理によって横軸と縦軸のスケールの異なる一様分布から生成された二つの乱数によって座標が一点決まることになります。
そしてこの乱数によって求まった値が予測分布よりも下の範囲であれば採用します。とてもシンプルですね!
当然ながらサンプリング数が多ければ多いほど、予測分布に近しい結果になります。
では、実装して結果を見ていきます。
まずはを求めますが、今回はScipyのOptimize関数を使います。ScipyのOptimize関数は目的関数を最適化(最小値もしくは最大値を求める)する便利な関数です。
f = beta(a=a, b=b).pdf res = opt.fmin(lambda x: -f(x), 0.1) M = f(res) print("M : ", M)
にマイナスをつけているのは、ベータ分布の最大値を求めるために最小値問題に帰着させるためです。
結果を見てみると上手くベータ分布の最大値を求められていることがわかります。
次に乱数を生成して提案分布が予測分布に近しい結果になっていくことを確認します。
まずは乱数を1000個生成した場合を見てみます。
N_MC = 1000 y = uniform.rvs(size=N_MC) mq = M * uniform.rvs(size=N_MC) accept = y[mq <= f(y)] plt.hist(accept, normed=True, bins=35, rwidth=0.8, label="rejection sampling") x = np.linspace(beta.ppf(0.001, a, b), beta.ppf(0.999, a, b), 100) plt.plot(x, beta.pdf(x, a, b), label="target distribution") plt.legend()
1000個の乱数だと、まだ少し予測分布(ベータ分布)に対して少しバラツキがありますね。
では、乱数を50000個に増やしてみましょう。
すると、乱数が1000個に比べてベータ分布に対してほぼ同じような曲線を描いていることがわかります。
最後に
今回は基本的なサンプリングアルゴリズムである棄却サンプリングを試してみました。異なる分布から求めたい分布のサンプルを非常にシンプルな考えで得ることが出来て面白いですね。サンプリングアルゴリズムは他にもたくさんあるので、是非試してみたいと思います。