重箱の隅過ぎて伝わらないZiggurat法: tailの部分の計算方法がなんでそうなるのか調べてみた。

By | 2018年11月3日 , Last update: 2022年8月7日
0 Comment

はじめに

久々に計算用紙的な記事を書きます。

ちょっと前の記事で、標準正規分布の確率密度関数を使って累積分布関数を求める方法について書きました。

上記の記事を書いた後で、標準正規分布に従う乱数の生成方法を調べてみたところ、Ziggurat法というものにたどり着きました。

Ziggurat法自体の説明や実装については上記のWikipediaのリンクを始めとして、いろいろなところで行われています(MATLAB等でも実装されているようです。)あることと、本サイトの管理人たるpandaが付け足せることはあまりないので、詳細についてはそちら(「参考文献」の節参照。)をご覧いただければと思います。

Ziggurat法では乱数生成の最初の手順において、面積が相等しい長方形を使って単調減少な確率密度関数$y = f(x) (x \ge 0)$と$x$軸及び$y$軸が囲む領域を覆いつくすような領域を考えます(よって長方形及び後述の「layer 0」の面積の総和はもとの領域の面積よりも大きくなります)。

しかし、最も$x$軸寄り、すなわち一番下の領域(以下、「layer 0」と呼びます。)だけは(面積は他の長方形に等しくなるように$y$軸方向の幅を決めるものの、)この近似を行わず、$x$軸及び$x$軸に平行な直線$y = y_1$、$y$軸ならびに$y = f(x) (x \ge x_1, y_1 = f(x_1))$で囲まれた領域(下図の黄色の領域)を考えます。

上記の「layer 0」及び長方形群の中からランダムに1個の領域を選ぶのがZiggurat法による乱数生成の最初の手順となります。

しかしながら、単調減少な確率密度関数が標準正規分布の確率密度関数(の右半分)

\begin{align}
\phi(x) &= \frac{1}{\sqrt{2\pi}}\,e^{-\frac{x^2}{2}}\qquad (x \ge 0)
\end{align}

であり、さらに「layer 0」が選ばれてしまった場合でかつ、その後の手順で$x \gt x_1$となること(すなわちtailの部分になること。上図の(b)の部分に相当します。なお、$x \le x_1$となった場合、すなわち上図の(a)の部分が選択されることに相当しますので、$x$がそのまま生成された乱数となります。)が確定してしまった場合の計算方法がきわめてシンプルな割にどうしてそうなるのかが(少なくともGoogle先生にお伺いを立ててみた限りでは)あまり詳細に記述されている資料が見当たらないようですので、試しに計算してみることにしました。

スポンサーリンク

その筋の論文等に割とよく載っているtailの計算方法

参考文献[2]を始めとした論文や実装を説明した文献などに割とよく載っているtailの計算方法は以下の通りです。

  1. $U_1, U_2 \in [0,1]$を一様分布の乱数として生成する。
  2. $Y = -\log U_1, X = -\displaystyle\frac{\log U_2}{x_1}$を計算し、
    \begin{align}
    X^2 &\lt 2Y \label{eq:tailtest}
    \end{align}

    であれば$X+x_1$を$x$として返す。(\ref{eq:tailtest})が成り立たない場合には、$U_1, U_2$を新たに生成するところからやり直す。

でも、これってなんでそうなるのかって上記の説明だけでは一撃ではわからないですよね?

少なくとも本サイトの管理人たるpandaは最初は理解できませんでした。(´・ω・｀)

見た目が極めて簡単な式ですので、「今日中にコードを書け!!」とか言われたらなんでそうなるのかなんてことは考えずにコードを書いてしまうところです。

そこで、なんでそうなるのか調べてみました。

まず、ぐぐります。


スポンサーリンク

まず手始めに、Google先生に何回かお伺いを立ててみたところ、参考文献[3]及び[4]が出てきました([3]はいろいろなところで引用されているTechnometricsの論文の下書きっぽいテクニカルレポートです)。参考文献[4]については、とりあえず必要なところが記述されているPDFファイルだけをダウンロードすることにします。

なお、本全体に興味のある方は↓をご参照ください。

[3]は3ページで、[4]は2通りの方法が3ページ足らずでそれぞれさらっと書いてありますが、さらっと書かれすぎていてわかりにくいので、[4]に書いてある2通りの方法に沿って試しに計算してみることにします。

標準正規分布の確率密度関数のtailの部分の乱数は以下のような手順で生成します。

  1. $x \ge x_1$において$e^{-\frac{x^2}{2}} \le g(x)$でそれ以外の$x$では$g(x) = 0$となり、かつ$\displaystyle\lim_{x \to \infty} g(x) = 0$となるような関数$g(x)$を考えます。
  2. \begin{align}
    C\int^{\infty}_{x_1} g(x) dx &= 1
    \end{align}

    となるような定数$C$を求めます。
  3. 確率密度関数$g^{*}(x) = Cg(x)$を考えてそれに従う乱数$x^{*}$を生成します。
  4. $x^{*}$が標準正規分布に従う乱数として採用可能ならその数を採用し、そうでなければ再度$x^{*}$を生成するという手順を繰り返します(参考文献[3]では”rejection technique”と書かれている方法です)。

その1: $xe^{-x/2}$ のようなタイプの式を使って乱数を生成する方法

ところで、$x \ge x_1$では(\ref{eq:exa})式が成り立ちます。

\begin{align}
e^{-\frac{x^2}{2}} &\le \frac{x}{x_1}e^{-\frac{x^2}{2}} \label{eq:exa}
\end{align}

(\ref{eq:exa})式の両辺に$x_1e^{\frac{x_1^2}{2}}$をかけると、
\begin{align}
x_1e^{\frac{x_1^2-x^2}{2}} &\le xe^{\frac{x_1^2-x^2}{2}} \label{eq:dtr}
\end{align}

となります。

ここで、(\ref{eq:dtr})式の右辺を区間$[x_1,\infty)$で$x$について積分すると…

\begin{align}
\int_{x_1}^{\infty}xe^{\frac{x_1^2-x^2}{2}}dx &= \left[ -e^{\frac{x_1^2-x^2}{2}} \right]_{x_1}^{\infty} \nonumber \\
&= 1
\end{align}

と計算できるので、確率密度関数として使えそうです。

そこで、

\begin{align}
U_1 &= e^{\frac{x_1^2-x^2}{2}} \label{eq:firsttransform}
\end{align}

とおいて$x$で微分してその絶対値をとると…
\begin{align}
\left| \frac{dU_1}{dx} \right| &= xe^{\frac{x_1^2-x^2}{2}}
\end{align}

となって、(\ref{eq:dtr})式の右辺と一致しますので、(\ref{eq:firsttransform})式を$x$について解いた$x = \sqrt{x_1^2-\log U_1}$は、確率変数の変数変換を行うために使えそうです。


スポンサーリンク

$U_1$を区間$[0,1]$で連続一様分布の乱数とすると、$x \gt x_1$であり、この値における確率密度は(\ref{eq:dtr})式の右辺に一致します。

さらに、区間$[0,1]$で連続一様分布の乱数で$U_1$とは独立の乱数$U_2$を生成させて、(\ref{eq:firsttest})の条件を満たすまで$U_2$を繰り返し生成させます。

\begin{align}
U_2xe^{\frac{x_1^2-x^2}{2}} &\le x_1e^{\frac{x_1^2-x^2}{2}} \label{eq:firsttest}
\end{align}

(\ref{eq:firsttest})式の両辺を$e^{\frac{x_1^2-x^2}{2}}$で割ると…
\begin{align}
U_2x &\le x_1 \label{eq:firstfinal}
\end{align}

と変形できます。

ものすごくシンプルな式になりましたね。✨

その2: 評価式の不等式の一方の辺が$e^{-(x-a)^2/2}$ のようなタイプの式に変形できる式を使って乱数を生成する方法

次に、$x \in \mathbb{R}$であれば$(x-x_1)^2\ge 0$であることから、(\ref{eq:exb})式が成り立つことを使って、前節と同様の議論ができないか考えてみることにします。

\begin{align}
e^{-\frac{x^2}{2}} &\le e^{\frac{x_1^2}{2}-x_1x} \label{eq:exb}
\end{align}

ここで、(\ref{eq:exb})式の右辺を区間$[x_1,\infty)$で$x$について積分すると…

\begin{align}
\int_{x_1}^{\infty} e^{\frac{x_1^2}{2}-x_1x} &= \left[ -\frac{1}{x_1}e^{\frac{x_1^2}{2}-x_1x} \right]_{x_1}^{\infty} \nonumber \\
&= \frac{1}{x_1}e^{-\frac{x_1^2}{2}}
\end{align}

と計算できます。前節と違って計算結果は1にはなりませんでしたが、正規化した関数(\ref{eq:sdf})式
\begin{align}
f(x) &= x_1e^{x_1^2-x_1x} \label{eq:sdf}
\end{align}

は、確率密度関数として使えそうです。

そこで、

\begin{align}
U_1 &= e^{x_1^2-x_1x} \label{eq:seconduone}
\end{align}

とおくと、$\left| \displaystyle\frac{dU_1}{dx} \right|$は(\ref{eq:sdf})式の右辺と等しくなります。したがって、$U_1$を区間$[0,1]$で連続一様分布の確率変数としたときの確率変数の変数変換の変換式として利用できます。(\ref{eq:seconduone})式を$x$について解くと、
\begin{align}
x &= x_1 – \frac{\log U_1}{x_1} \label{eq:secondx}
\end{align}

になります。

$x$を生成したら、今度は区間$[0,1]$で連続一様分布の乱数で$U_1$とは独立の乱数$U_2$を生成させて、(\ref{eq:secondcheck})の条件を満たすまで$U_2$を繰り返し生成させます。

\begin{align}
U_2x_1e^{x_1^2-x_1x} &\lt x_1e^{\frac{x_1^2-x^2}{2}} \label{eq:secondcheck}
\end{align}

(\ref{eq:secondcheck})式の両辺を$x_1e^{x_1^2-x_1x}$で割ると…
\begin{align}
U_2 &\lt e^{-\frac{(x_1-x)^2}{2}} \label{eq:secondstyle}
\end{align}

両辺の対数をとると…
\begin{align}
\log U_2 &= -\frac{(x_1-x)^2}{2} \label{eq:thirdstyle}
\end{align}

(\ref{eq:thirdstyle})式に(\ref{eq:secondx})式を代入し、さらに両辺に-2をかけると不等号の向きが変わって、
\begin{align}
-2\log U_2 &\gt \left(\frac{\log U_1}{x_1}\right) \label{eq:finalstyle}
\end{align}

これは、$Y = -\log U_2, X = -\displaystyle\frac{\log U_1}{x_1}$とおくと、(\ref{eq:tailtest})式と一致します($U_1, U_2$の順番が入れ替わっていますが、単なるnotationの違いですので、気にしない方向でお願いします)。

何とか(\ref{eq:tailtest})式にたどりづきました!! ✨✨

おまけ:切断レイリー分布

ところで、(\ref{eq:exa})式の右辺に$x_1e^{\frac{x_1^2}{2}}$をかけて、さらに不定積分を考えると、

\begin{align}
\int xe^{\frac{x_1^2-x^2}{2}}dx &= -e^{\frac{x_1^2-x^2}{2}} + C \label{eq:idi}
\end{align}

($C$は積分定数です。)となります。

そこで、(\ref{eq:idi})式の右辺を$F(x)$と置くと、$F(x_1)=0$なので$C=1$となります。つまり、F(x)は以下のように表せます。

\begin{align}
F(x) &= 1-e^{\frac{x_1^2-x^2}{2}} \label{eq:trd}
\end{align}

(\ref{eq:trd})式は$x \ge x_1$でのみ定義されているレイリー分布(切断レイリー分布)の累積分布関数です。

まとめ

標準正規分布に従う乱数を生成する際にこの記事で取り扱う計算が行われる確率は実は非常に小さい(標準正規分布, $x$軸及び$y$軸に囲まれた領域を256個に分割する場合で0.03%程度)のですが、Ziggurat法が標準正規分布に厳密に従う乱数を生成する方法であるため、理解しておいても損にはならないと思います。

ただ、確率密度関数の計算を扱う分野ではよくあることなのかどうかわかりませんが、それまでの式を定数倍した式が説明なしにいきなり登場したり、「〜は$e^{-x^2/2}$に比例する。」といった記述が突然登場したりすることがよくあります。特に後者のような書き方をされると、何をかけたらいいのかがすぐにはわからないので、読み解くのにかなり苦労します。

正直なところ、最初はあまりにも理解できなかったので、この記事はボツにしようかとも思ってました。(´・ω・`)

この記事を書くにあたり、国内外のWebサイトの記述を調べてみましたが、見た目は簡単なせいか(\ref{eq:tailtest})式の比較を採用している例が多いようです。(\ref{eq:firstfinal})式の両辺を二乗すると$x$の平方根の計算はこの比較後に行う実装とすることも可能であることと、自然対数の計算は1回だけで済みますので、(\ref{eq:firstfinal})式(+比較式の両辺を二乗する方法)の方がもしかすると処理速度が速いのかもしれません。

また、ここまでに書いた計算の方法または理解が正しいのかどうかについての保証はできませんが、何かの際の参考にしていだけると幸いです。

この記事は以上です。

References / 参考文献

Related Posts / 関連ページ

夏が来れば思い出す、Euler角、四元数(1): ちょっと計算してみた。 整数次の変形Bessel関数の母関数とVon Mises分布の関係。 グラフを描画するプログラムをscratchから書かねばならなくなった時に、軸の目盛り(tick)の間隔をいい感じに設定する方法。 要素ごとに計算して理解する正規方程式。