棄却法によるサンプリングで確率変数を生成する方法をわかりやすく解説

2024年9月17日2024年9月20日

確率分布に従った乱数を生成する方法は多く存在しますが、複雑な分布では直接サンプリングが難しいことがあります。そんなときに役立つのが 棄却法（Rejection Sampling） です。この記事では、棄却法の原理と使い方について、わかりやすく解説します。

この記事では、確率変数は全て連続型のものを想定しています。

標準一様分布に従う確率変数\(U\)と、確率密度関数が\(g\)である確率変数\(Y\)と、特定のルールLから、
確率密度関数が\(f\)である確率変数\(X\)を生成する方法を考えてみます。

補足ですが、標準一様分布に従う確率変数というのは、確率密度関数が
\begin{align*} f = \begin{cases} 1 &(x \in [0, 1]) \\ 0 & (x \notin [0,1]) \end{cases}\end{align*}
である確率変数です。別の言い方をすると、値域が0から1の、一様分布に従う確率変数です。

補足

例えば、四面サイコロ(正四面体の形をしていて、それぞれの面に1から4まで刻まれているサイコロ)が手元になくて困っているとします。
しかし、手元には六面サイコロはあるとします。そこで、
サイコロを振って5, 6が出たらノーカウントにする。
というルールを定めると、実質的にそのサイコロは四面サイコロだと思えることが直感的にわかると思います。

少し状況はことなりますが、このように、確率変数\(Y\)に然るべきルールを適用することで、確率変数\(X\)を生成します。

\(X\)の分布関数を\(F\)で表すことにします。
\(F(x)\)を、標準一様分布の確率密度関数と、何か適当な、確率密度関数が\(g\)だとわかっている\(Y\)で強引に表現してみましょう。
良い表現を得るために、\(f(x) \leq g(x) \quad (\forall x )\)であるような\(g\)を用いることにします。
\begin{align} F(x) &= \int_{-\infty}^x f(\xi) d\xi \\&= \int_{-\infty}^x \frac{f(\xi)}{g(\xi)}g(\xi) d\xi \\&= \int_{-\infty}^{x} \left(\int_0^{\frac{f(\xi)}{g(\xi)}} d \eta \right) g(\xi) d\xi \\&= P(Y \leq x, U \leq \frac{f(Y)}{g(Y)})\end{align}
と、無理やり表現してやります。
(補足ですが、\(f(x) \leq g(x)\)でないと、\(P(U \leq \frac{f(x)}{g(x)}) = \int_0^\frac{f(x)}{g(x)} d \eta \)となりません)

天下り的ですが、それに合わせて、
\begin{align*}1 &= \int_{-\infty}^\infty f(\xi) d\xi \\&= \int_{-\infty}^\infty \frac{f(\xi)}{g(\xi)}g(\xi) d\xi \\&= \int_{-\infty}^{\infty} \left(\int_0^{\frac{f(\xi)}{g(\xi)}} d \eta \right) g(\xi) d\xi \\&= P(U\leq \frac{f(Y)}{g(Y)}) \end{align*}
と表現してやります。

これらを合わせると、
\begin{align*} F(x) &= \int_{-\infty}^{x} \left(\int_0^{\frac{f(\xi)}{g(\xi)}} d \eta \right) g(\xi) d\xi \\&= \frac{\int_{-\infty}^{x} \left(\int_0^{\frac{f(\xi)}{g(\xi)}} d \eta \right) g(\xi) d\xi}{1} \\&= \frac{P(Y \leq x, U \leq \frac{f(Y)}{g(Y)})}{P(U\leq \frac{f(Y)}{g(Y)})} \\&= P(Y \leq x \mid U \leq \frac{f(Y)}{g(Y)})\end{align*}
となります。

棄却法もどき

これは、\(Y= y, U = u\)と実際に値を出してみた時に、
\begin{align*} u > \frac{f(y)}{g(y)} \end{align*}
であった場合にはノーカウントとしたときの、確率変数\(Y\)の分布関数です。
これが\(X\)の分布関数\(F\)と一致しているため、
\(Y= y, U = u\)と実際に値を出してみた時に、
\begin{align*} u > \frac{f(y)}{g(y)} \end{align*}
であった場合にはノーカウントとするというルールを設けて、それ以外の場合はノーカウントとしなければ、その出力は確率変数\(X\)として振る舞うということです。

ここで、この手続きを洗練させます。
\(Y= y, U = u\)と実際に値を出してみた時に、\begin{align*} u > \frac{f(y)}{g(y)} \end{align*}
だとノーカウントとしましたが、ノーカウントが多すぎると時間の無駄なので、できるだけカウントされるようにしましょう。

\(g\)がクソでかだと、\(\frac{f}{g}\)がほとんど\(0\)になってしまい、ほぼほぼノーカウントが発生してしまうので、\(g\)を適当に定数\(c\)倍して、\(f \leq c g\)が成り立つ範囲でできるだけ小さくします。
つまり、
\begin{align*} f(x) \leq c g(x) \quad ( \forall x) \end{align*}
が成り立つ最小の\(c\)を求め、ルールを改変します。

\(g\)を\(cg\)に置き換えると、（式変形は全く同じなので読み飛ばして大丈夫です）

\begin{align} F(x) &= \int_{-\infty}^x f(\xi) d\xi \\&= \int_{-\infty}^x \frac{f(\xi)}{cg(\xi)}cg(\xi) d\xi \\&= \int_{-\infty}^{x} \left(\int_0^{\frac{f(\xi)}{cg(\xi)}} d \eta \right) cg(\xi) d\xi \\&= P(Y \leq x, U \leq \frac{f(Y)}{cg(Y)})\end{align}

\begin{align*}1 &= \int_{-\infty}^\infty f(\xi) d\xi \\&= \int_{-\infty}^\infty \frac{f(\xi)}{cg(\xi)}cg(\xi) d\xi \\&= \int_{-\infty}^{\infty} \left(\int_0^{\frac{f(\xi)}{cg(\xi)}} d \eta \right) cg(\xi) d\xi \\&= P(U\leq \frac{f(Y)}{cg(Y)}) \end{align*}

\begin{align*} F(x) &= \int_{-\infty}^{x} \left(\int_0^{\frac{f(\xi)}{cg(\xi)}} d \eta \right) cg(\xi) d\xi \\&= \frac{\int_{-\infty}^{x} \left(\int_0^{\frac{f(\xi)}{cg(\xi)}} d \eta \right) cg(\xi) d\xi}{1} \\&= \frac{P(Y \leq x, U \leq \frac{f(Y)}{cg(Y)})}{P(U\leq \frac{f(Y)}{cg(Y)})} \\&= P(Y \leq x \mid U \leq \frac{f(Y)}{cg(Y)})\end{align*}

となるので、

棄却法

\(Y= y, U = u\)と実際に値を出してみた時に、
\begin{align*} u > \frac{f(y)}{cg(y)} \end{align*}
であった場合にはノーカウントとするというルールを設けて、それ以外の場合はノーカウントとしなければ、その出力は確率変数\(X\)として振る舞い、かつ、ノーカウントとなる確率も最小化されています。

このルールにより\(X\)を再現する方法を棄却法といいます。

また、補足ですが、\(c\)倍をするという操作により、最初に選ぶ確率変数\(Y\)は\(f \leq g\)を満たしていなくても良いことがわかります。

つまり、まとめると、

STEP

確率密度関数が\(f\)である確率変数\(X\)の生成をしたい。

STEP

標準一様分布に従う確率変数\(U\)と、確率密度関数が\(g\)である確率変数\(Y\)を用意する。

最初に選ぶ確率変数\(Y\)の確率密度関数\(g\)は\(f \leq g\)を満たしていなくても良い

STEP

\(f(x) \leq c g(x) \quad (\forall x)\)を満たす最小の\(c\)を求める。

STEP

\(U = u, Y = y\)と、値を実際に出力する。