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

はじめに

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

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

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

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

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

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

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

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

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

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

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

1. を一様分布の乱数として生成する。
2. を計算し、
   
   (2)
   
   であればをとして返す。(2)が成り立たない場合には、を新たに生成するところからやり直す。

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

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

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

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

まず、ぐぐります。

スポンサーリンク

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

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

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

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

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

その1: のようなタイプの式を使って乱数を生成する方法

ところで、では(4)式が成り立ちます。

(4)式の両辺にをかけると、

となります。

ここで、(5)式の右辺を区間でについて積分すると…

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

そこで、

とおいてで微分してその絶対値をとると…

となって、(5)式の右辺と一致しますので、(7)式をについて解いたは、確率変数の変数変換を行うために使えそうです。

を区間で連続一様分布の乱数とすると、であり、この値における確率密度は(5)式の右辺に一致します。

さらに、区間で連続一様分布の乱数でとは独立の乱数を生成させて、(9)の条件を満たすまでを繰り返し生成させます。

(9)式の両辺をで割ると…

と変形できます。

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

その2: 評価式の不等式の一方の辺が のようなタイプの式に変形できる式を使って乱数を生成する方法

次に、であればであることから、(11)式が成り立つことを使って、前節と同様の議論ができないか考えてみることにします。

ここで、(11)式の右辺を区間でについて積分すると…

と計算できます。前節と違って計算結果は1にはなりませんでしたが、正規化した関数(13)式

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

そこで、

とおくと、は(13)式の右辺と等しくなります。したがって、を区間で連続一様分布の確率変数としたときの確率変数の変数変換の変換式として利用できます。(14)式をについて解くと、

になります。

を生成したら、今度は区間で連続一様分布の乱数でとは独立の乱数を生成させて、(16)の条件を満たすまでを繰り返し生成させます。

(16)式の両辺をで割ると…

両辺の対数をとると…

(18)式に(15)式を代入し、さらに両辺に-2をかけると不等号の向きが変わって、

これは、とおくと、(2)式と一致します(の順番が入れ替わっていますが、単なるnotationの違いですので、気にしない方向でお願いします)。

何とか(2)式にたどりづきました!!

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

ところで、(4)式の右辺にをかけて、さらに不定積分を考えると、

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

そこで、(20)式の右辺をと置くと、なのでとなります。つまり、F(x)は以下のように表せます。

(21)式はでのみ定義されているレイリー分布(切断レイリー分布)の累積分布関数です。

まとめ

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

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

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

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

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

この記事は以上です。

References / 参考文献

• [1] Ziggurat method
• [2] George Marsaglia; Wai Wan Tsang (2000). “The Ziggurat Method for Generating Random Variables”. Journal of Statistical Software. 5 (8).
• [3] George Marsaglia; generating a variable from the tail of the normal distribution
• [4] Luc Devroye; Non-Uniform Random Variate Generation pp.380-382