如何產生正態分布的隨機數？

01-04

我為了這個問題做了個開源項目 miloyip/normaldist-benchmark · GitHub。

這個項目比較不同正態隨機數生成演算法的性能，它測試的函數原型是：

void normaldistf(float* data, size_t n); void normaldist(double* data, size_t n);

它測試了以下實現的性能，並驗證其正確性（mean, SD, Skewness, Kurtosis）。

boxmuller: Box-Muller Transform。
cpp11normal: C++11 的 std::normal_distribution。（在 clang/LLVM 中的實現是 Box-Muller transform。）
clt：利用中心極限定理，把m個均勻分布隨機數相加。
inverse: 利用正態分布 CDF 的反函數做 Inverse transform sampling。
marsagliapolar: Marsaglia polar method。
ziggurat: Ziggurat algorithm，n=128。
null: 作為比較基準，生成均勻分布的隨機數。

一些演算法包含 SSE2 和 AVX 指令集的版本。

注意 clt 演算法的 abs(kurtosis) &>= 0.01，未能通過正確性驗證。其他除了 null 外，都通過正確性驗證。

先看看 float 版本的結果（iMac Corei5-3330S@2.70GHz clang 6.1 64-bit）：

http://rawgit.com/miloyip/normaldist-benchmark/master/result/Corei5-3330S@2.70GHz_mac64_clang6.1.html

可能很多人以為 Zigguart 演算法是最快的，但查表操作、不定數量的分支都不利於 SIMD 並行化。而 Marsaglia polar method 也要做 rejection sampling，所以也會造成分支。Marsaglia polar method 原意是為了改善 Box-Muller 的性能的，但在此測試中反而較慢。Inverse 雖然簡單，但該反函數也需要分支處理中央和尾的部分，形成分支。

在加入 SSE2 和 AVX 加速後（其 PRNG 為簡單 SIMD 並行的 LCG），Box-Muller 有很明顯優勢，SSE2 已是 ziggurat 的 1.79x，AVX 是 2.99x。AVX 實現中使用了SSE2/SSE4.1的指令分拆來做一些操作，要在 AVX2 中才能完全實現 8-way 計算。

相對於維基上 Box-Muller 的例子代碼：

double generateGaussianNoise(double mu, double sigma) { const double epsilon = std::numeric_limits&::min(); const double two_pi = 2.0*3.14159265358979323846;


	static double z0, z1;

	static bool generate;

	generate = !generate;
	if (!generate)

	   return z1 * sigma + mu;

double u1, u2; do { u1 = rand() * (1.0 / RAND_MAX); u2 = rand() * (1.0 / RAND_MAX); } while ( u1 &<= epsilon ); z0 = sqrt(-2.0 * log(u1)) * cos(two_pi * u2); z1 = sqrt(-2.0 * log(u1)) * sin(two_pi * u2); return z0 * sigma + mu; }

我實現的Box-Muller 做了一些優化：

template & void boxmuller(T* data, size_t count) { assert(count % 2 == 0); static const T twopi = T(2.0 * 3.14159265358979323846);

LCG& r; for (size_t i = 0; i &< count; i += 2) { T u1 = 1.0f - r(); // [0, 1) -&> (0, 1] T u2 = r(); T radius = std::sqrt(-2 * std::log(u1)); T theta = twopi * u2; data[i ] = radius * std::cos(theta); data[i + 1] = radius * std::sin(theta); } }

直接寫入兩個結果（這個與 API 設計相關）
把 [0, 1) 的隨機數轉換成 (0, 1]，避免了 log(0) 的出現，不需要做 rejection。
在 SSE2/AVX 版本中，利用 sincos() 函數同時計算 sin 和 cos。

如果目標平台有 SIMD 指令集，建議使用 Box-Muller，否則可考慮 Zigguart。

如果有其他改進或演算法，歡迎 pull request。

謝邀。

先假設 $Z sim mathcal{N}(0,1)$ 。

用Box-Muller方法，隨機抽出兩個從 $[0,1]$ 均勻分布的數字 $u$ 和 $v$ 。然後

$z_1 = sqrt{-2 log u} cos 2pi v$

$z_2 = sqrt{-2 log u} sin 2pi v$

那 $z_1$ 和 $z_2$ 都是正態分布的。

證明可用極坐標，請參考教科書中的Box-Muller方法。

另外使用反函數，先隨機抽出一個從 $[0,1]$ 均勻分布的數字 $u$ ，然後

$z=sqrt{2} ext{erf}^{-1}(2u-1)$

那 $z$ 是正態分布的。

Python實現：

import numpy as np from scipy.special import erfinv


def boxmullersampling(mu=0, sigma=1, size=1):

    u = np.random.uniform(size=size)

    v = np.random.uniform(size=size)

    z = np.sqrt(-2*np.log(u))*np.cos(2*np.pi*v)

    return mu+z*sigma

def inverfsampling(mu=0, sigma=1, size=1): z = np.sqrt(2)*erfinv(2*np.random.uniform(size=size)-1) return mu+z*sigma

最簡單的：rejection sampling，思路很簡單，也很容易實現，但效率較差
較複雜的：inverse CDF，直接利用累積分布函數（CDF）的反函數生成隨機數，但計算中牽扯到比較複雜的誤差函數erf（非初等函數）
更好的：Box-Muller演算法，在很長時間內都是生成正態分布隨機數的"標準"演算法。Box-Muller演算法的特點是效率高，並且計算過程比較簡單（只用到了初等函數）。參見：Box-Muller transform
目前最好的（相較於其它實用演算法）：ziggurat演算法，效率很高，很多現代的編程語言都使用了這一演算法。ziggurat並不是人名，其含義是「金字形神塔」，不是埃及那個金字塔，而是古代蘇美爾人建造的類金字塔結構的神壇：神壇由多層平台構成，每層平台都呈矩形、卵形或正方形，且自下而上面積逐漸減小。ziggurat演算法實際上是一種改進的、包含查表操作的rejection sampling。參見：Ziggurat algorithm

--------------------------------------------------------------------------------------------------------------------------------------

樓上知友 @何史提的答案中已經涵蓋了inverse CDF和Box-Muller演算法，這裡只簡要介紹一下ziggurat演算法：

先考慮最簡單的rejection sampling，其方法如下圖所示，首先生成在藍色矩形框中均勻分布的隨機點 $(x,y)$ ，然後捨棄（即拒絕，reject）所有未落在正態分布鐘形曲線以下的點，剩餘點的橫坐標 $x$ 就構成一個服從相應正態分布的隨機數序列。rejection sampling的效率之所以低下，就是由於很多生成的隨機點因為落在鐘形曲線以上而被白白捨棄。

如果有一種方法能夠直接產生在鐘形曲線下方面積內均勻分布的隨機點，那麼就不存在捨棄隨機點的情況，也就能夠避免上面所說的「浪費」問題，這正是ziggurat演算法的基本思路。具體來說，ziggurat演算法的內容如下：

% 以下內容均遵從MATLAB語言的規定，數組及向量的下標索引從1開始編號

正態分布的鐘形曲線左右對稱，這裡只取右半邊考慮，如下圖所示。首先將鐘形曲線以下的面積等分為 $n$ 份，包括 $(n-1)$ 個水平矩形和最下方一直延伸至無限大的尾巴，即 $(n-1)$ 個矩形及最下方尾巴的面積都相等。這些層疊的矩形看起來像古代蘇美爾人建造的金字形神塔，這就是ziggurat這個名稱的由來。圖中取 $n=8$ ，實踐中 $n$ 可能達到 $200$ 以上。設各水平矩形右邊界的橫坐標（即各圓圈標記的橫坐標）為 $z_k (k=1,2,dots,n)$ ，其中 $z_1=0$ 。當 $n$ 確定後， $z_k$ 也就確定了，可以通過數值方法求解超越方程得到 $z_k$ 的值並製表記錄。這裡稱長度方向上與上一層矩形重疊的部分為當前矩形的核心區，下圖中各核心區的右邊界就是虛線邊，顯見最上面的矩形沒有核心區。

圖片出自Numerical Computing with MATLAB一書的第九章Random Numbers。順便提一下，此書作者就是MATLAB的創始者、MathWorks公司的主席和首席數學家Cleve B. Moler。

定義 $sigma_kequiv z_{k-1}/z_k (k=2,3,dots,n)$ ，即相鄰兩層矩形水平方向長度的比值，規定 $sigma_1equiv0$ 。顯然，依照所要求的參數（ $mu,,sigma$ ）生成正態分布時，所有的 $z_k$ 、 $sigma_k$ 都只需要在初始化時計算一次並製表記錄即可，之後不需要再重複計算，除非所要求生成的正態分布的參數（ $mu,,sigma$ ）發生改變。

初始化結束後就可以開始生成正態分布隨機數了：首先生成一個 $[1,n]$ 區間內均勻分布的隨機整數 $j$ 和一個 $(-1,1)$ 區間內均勻分布的隨機浮點數 $u$ 。然後檢查 $|u|<sigma_j$ 是否成立，若成立則 $uz_j$ 就在第 $j$ 段核心區內，顯然也就在鐘形曲線下方，故可將 $uz_j$ 作為正態分布的一個採樣點返回。相應的偽代碼為

j = randint(1, n); % [1,n]區間內均勻分布的隨機整數 u = 2 * rand() - 1; % （-1,1）區間內均勻分布的隨機浮點數 if abs(u) &< sigma[j] return u * z[j]; % u * z[j]在第j段核心區內 end

由此可見，使用ziggurat演算法時只需要生成一個隨機整數和一個隨機浮點數，做一次條件判斷和一次乘法運算，外加兩次查表即可得到一個服從正態分布的隨機數，整個過程中沒有開根號、求對數、算三角函數等複雜運算。

在上面偽代碼中，那個if條件判斷在絕大多數情況下都為true，剩下為false的就是被reject的情況，具體又包含三種可能：(1). $j=1$ ，最上面的矩形沒有核心區；(2). $j=2sim(n-1)$ ， $uz_j$ 恰好落在核心區之外（即右側包含曲線的小矩形中）；(3). $j=n$ ， $uz_j$ 落在底部核心區之外的尾巴中。在這三種情況下就需要依據Box-Muller演算法利用已經生成的均勻分布隨機數進行額外運算以生成正態分布隨機數。容易看出， $n$ 越大，出現這三種需要額外運算的情況的可能性就越低。Numerical Computing with MATLAB一書中給出的數據是當 $n=128$ 時，需要額外運算的可能性就小於3%，所以這部分額外運算對ziggurat演算法的總體效率影響很小。

ziggurat演算法適用於連續生成大量服從同一正態分布的隨機數的情形，在這種情況下ziggurat演算法的效率遠高於Box-Muller演算法。反之，如果只是需要少量正態分布隨機數，那麼考慮到計算 $z_k$ 、 $sigma_k$ 的開銷，ziggurat演算法的總體表現未必優於Box-Muller演算法。

ziggurat演算法的基本思想最早由美國數學家兼計算機科學家George Marsaglia於20世紀60年代提出，在之後的歲月中該演算法不斷得到改進和發展，目前常用的版本大多源自於George Marsaglia和香港數學家曾衛寰（Wai Wan Tsang）在2000年發表的論文The Ziggurat Method for Generating Random Variables。

最後需要說明兩點：

實際使用的ziggurat演算法會比上面介紹的版本更加複雜。如果要在實際項目（特別是與加密有關的項目）中使用正態分布隨機數，強烈建議直接使用標準庫或是經廣泛檢驗的第三方庫所提供的函數，而不是照著上面的介紹自己隨便寫一個湊合著用。
ziggurat演算法實際上並不局限於生成正態分布隨機數，只要一個連續分布的概率密度函數（PDF）圖象是一條下降的曲線或是像正態分布這樣的鐘形曲線，理論上都可以使用ziggurat演算法生成。

剛好我的專欄裡面有一個現成的code：Julia call C - 碼農經濟學 - 知乎專欄

MATLAB大法好

R＝normrnd(MU,SIGMA,m,n)：生成m×n的服從正態分布的隨機數矩陣。

其中：MU為返回均值，SIGMA為標準差

如果要畫直方圖的話可以用bar函數。

如果你的要求不高，比如只是應用在遊戲里，那麼可以簡單地累加若干個均勻分布的隨機數。實際上，累加三個隨機數就可以形成看起來很平滑的分布。我從《遊戲編程精粹7》學到的這個技巧。

取兩個（或多個）隨機數的平均值即可產生非常近似於正態分布的結果。龍與地下城就廣泛使用了這一方法。（就是那個2d6）

我記得統計之都有篇文章講過這個，一個是Box-Muller直接算，一個是 CDF 反演，一個我忘了，最後一個我比較喜歡，是 rejection sample

鏈接在此

Editor: 漫談正態分布的生成

最簡單粗暴的方法：下載r，寫一行r代碼：y&<--dnorm（x，0，1）。這是標準正態分布。把0和1可以替換成想要的mean和standard deviation

如何生成隨機數（上） ? Free Mind

python

random.gauss(aver,sigma)

SAS中使用 RANNOR(seed )或者NORMAL(seed)

在獲得均勻分布隨機數的基礎上, 進行 Box-Muller 變換, 一次計算可獲得兩個正態分布的隨機數.

Box-Muller transform: https://en.wikipedia.org/wiki/Box-Muller_transform

下面是 C# 代碼的實現:

/// & /// 產生正態分布的隨機數 /// & /// &

正態分布的平均值, 即 N(μ, σ^2) 中的 μ & /// &

正態分布的標準差, 即 N(μ, σ^2) 中的 σ & /// & 返回正態分布的隨機數. 理論值域是 μ±∞; 落在 μ±σ, μ±2σ, μ±3σ 的概率依次為 68.26%, 95.44%, 99.74% & public float Normal(float averageValue, float standardDeviation) { return averageValue + standardDeviation * (float) ( Math.Sqrt(-2.0 * Math.Log(1.0 - NextDouble())) * Math.Sin((Math.PI + Math.PI) * NextDouble()) ); }

其中 NextDouble 是產生一個 [0, 1) 區間的均勻分布隨機數. 這個方法丟棄了另一個結果, 並且一次調用導致隨機數種子變化兩次.

老生常談的 Box-Muller 方法不過實現有些不同，不用三角函數。

參考 Generating Gaussian Random Numbers

// return a random double between zero and 1 inline double RandFloat() {return ((rand())/(RAND_MAX+1.0));}


//returns a random number with a normal distribution

inline double RandGaussian(double mean = 0.0, double standard_deviation = 1.0)

{

  double x1, x2, w, y1;

  static double y2;

  static int use_last = 0;
  if (use_last)		        /* use value from previous call */

  {

    y1 = y2;

    use_last = 0;

  }

  else

  {

    do

    {

      x1 = 2.0 * RandFloat() - 1.0;

      x2 = 2.0 * RandFloat() - 1.0;

      w = x1 * x1 + x2 * x2;

    }

    while ( w &>= 1.0 );

w = sqrt( (-2.0 * log( w ) ) / w ); y1 = x1 * w; y2 = x2 * w; use_last = 1; } return( mean + y1 * standard_deviation ); }

1,先產生均勻分布的隨機數(線性取余等)

2,把均勻分布變換到正態分布(經典的有Box-Muller Transform，傅里葉變換等)

用R，輸入：x&<-rnorm(50,60,5) 再輸入x 就可以得到50個均值60，標準差是5的隨機數。

java Random.nextGaussian()

《概率論與數理統計》一書中似乎給了一種辦法

不知道matlab行不行……

順便，那還叫隨機數么？

用均勻分布的隨機數生成器多隨幾個，取平均數就行