Box-Muller 與 ziggurat

1. Ziggurat 算法與 Box-muller 算法的效率比較

2. Box-Muller

a. 一般形式 因函數調用較多，速度慢，當u接近0時存在數值穩定性問題

先假設。

用Box-Muller方法，隨機抽出兩個從均勻分布的數字和。然后


那和都是正態分布的。
證明可用極坐標，請參考教科書中的Box-Muller方法。

另外使用反函數，先隨機抽出一個從均勻分布的數字，然后

那是正態分布的。

b.極坐標形式，速度快且有更好的數值魯棒性

double pf_ran_gaussian(double sigma)
{
  double x1, x2, w, r;

  do
  {
    do { r = drand48(); } while (r==0.0);
    x1 = 2.0 * r - 1.0;
    do { r = drand48(); } while (r==0.0);
    x2 = 2.0 * r - 1.0;
    w = x1*x1 + x2*x2;
  } while(w > 1.0 || w == 0.0);

  return(sigma * x2 * sqrt(-2.0*log(w)/w));
}

 3. zigurat 方法

Box–Muller 算法雖然非常快，但是由於用到了三角函數和對數函數，相對來說還是比較耗時的，如果想要更快一點有沒有辦法呢？

當然有，這就是 Ziggurat 算法，不僅可以用於快速生成正態分布，還可以生成指數分布等等。Ziggurat 算法的基本思想是利用拒絕采樣，什么是拒絕采樣呢？

拒絕采樣（Rejection Sampling），有的時候也稱接收 - 拒絕采樣，使用場景是有些函數$\mathrm{p(x)}$

太復雜在程序中沒法直接采樣，那么可以設定一個程序可抽樣的分布$\mathrm{q(x)比如正態分布等等，然后按照一定的方法拒絕某些樣本，達到接近p(x)}$

分布的目的：

具體操作如下，設定一個方便抽樣的函數 $q(x)$，以及一個常量 $k$，使得 $p(x)$ 總在 $kq(x)$ 的下方。（參考上圖）

• $x$$\mathrm{軸方向：從q(x)分布抽樣得到a}$

• $y$$\mathrm{軸方向：從均勻分布(0,kq(a))中抽樣得到u}$

• 如果剛好落到灰色區域：$\mathrm{u>p(a)，}$拒絕；否則接受這次抽樣

• 重復以上過程

證明過程就不細說了，知道怎么用就行了，感興趣的可以看看這個文檔

• Acceptance-Rejection Method

不過在高維的情況下，拒絕采樣會出現兩個問題，第一是合適的 $q$ 分布比較難以找到，第二是很難確定一個合理的 $k$ 值。這兩個問題會造成圖中灰色區域的面積變大，從而導致拒絕率很高，無用計算增加。

采用拒絕采樣來生成正態分布，最簡單直觀的方法莫過於用均勻分布作為 $q(x)$，但是這樣的話，矩形與正態分布曲線間的距離很大，就會出現剛才提到的問題，高效也就無從談起了。

而 Ziggurat 算法高效的秘密在於構造了一個非常精妙的$\mathrm{q(x)}$，看下面這張圖

我們用多個堆疊在一起的矩形，這樣保證陰影部分（被拒絕部分）的始終較小，這樣就非常高效了

簡單對圖作一個解釋：

• 我們用R[i]來表示一個矩形，R[i]的右邊界為x[i]，上邊界為y[i]。S[i]表示一個分割，當i=0時，S[0]=R[0]+tail，其他情況S[i]==R[i]
• 除了R[0]以外，其他每個矩形面積相等，設為定值A。R[0]的面積 = A-tail的面積。這樣保證從任何一個分割中抽樣(x,y)的概率相等
• 當任意選定一個R[i]在其中抽樣(x,y)，若x<x[i+1]，y必然在曲線下方，滿足條件，接受x；若x[i+1]<x<x[i]，則還需要進一步判斷。同樣這里R[0]和tail中的樣本需要進行特殊處理
• 這里為了方便解釋，只用了幾個矩形，在程序實現的時候，一般使用128或256個矩形

可以看出，為了提高效率，Ziggurat 算法中使用了許多技巧性的東西，這在其C代碼實現中更加明顯，使用了與運算和字節等各種小技巧，代碼就不在這里貼了，感興趣可以看看下面幾個版本，C版本的追求的是極致的速度，每個矩形的邊界已經提前計算好了。C#版本中的注釋非常詳細，Java版的基本與C#一致，但是效率一般。

• C
• C#
• Java

4. 總結

Box-muller 算法應對一般的需求足夠了，但是要生成大量服從正態分布的隨機數時，Ziggurat 算法效率會更高一點。

 

參考： https://www.taygeta.com/random/gaussian.html    // Box-Muller的介紹

            https://cosx.org/2015/06/generating-normal-distr-variates   // 對比介紹

            https://www.zhihu.com/question/29971598