在一些情況下經常需要用到隨機數，而高斯隨機數又是最常用到的。這一篇講一下如何編程生成符合正態分布的高斯隨機數，甚至任何其他分布的隨機數。
我們知道C語言的標准庫函數可以生成符合均勻分布的偽隨機數。那么如何生成符合高斯分布的隨機數呢？我們知道用逆函數法可以由符合(0,1)均勻分布的隨機數得到符合任意分布的隨機數，因此同樣可以得到符合高斯分布的隨機數。簡單證明如下：
設隨機變量u是符合(0,1)之間的均勻分布，那么有
。設隨機變量X的累積分布函數（CDF）為
，其逆函數為
。令隨機變量
。那么
。
因此只要得到所需要的隨機分布的累計分布函數（CDF）的逆函數，就可以簡單的通過逆函數把(0,1)均勻分布的隨機數變換成所需要的隨機數。高斯分布的概率分布函數（PDF）為
，
其累計分布函數（CDF）是PDF的積分
，為從0到1的單調遞增函數。
但是高斯分布的累積分布函數（CDF）不是初等函數，是沒法用解析式表達的。再者符合高斯分布的隨機變量的范圍是從負無窮到正無窮的，是不可能用計算機生成的。即使用浮點數表示，也不是連續的。比如我們有(0,1)之間均勻分布的隨機數，通過高斯分布CDF的逆函數變換到正負無窮，也只是得到一些離散的點。
因此要解決這些問題，首先用計算機生成的隨機數肯定是離散的。比如C語言的rand()函數返回[0,RAND_MAX]之間的偽隨機整數。所以我們也取一定范圍的整數作為生成的隨機高斯數。然后根據高斯分布的性質，離均值越遠，概率越小，通常3σ以外的地方可以近似的認為概率為0。所以可以截斷高斯分布的范圍，讓生成的高斯隨機數位於[μ-3σ,μ+3σ]。這樣CDF的無窮積分就可以由有限的求和運算代替了。算法描述如下：
- 設定高斯隨機數的范圍是[0,2r]，則均值是r，σ=r/3。
- 由PDF計算截斷的概率分布，然后求和計算累積分布。
- 輸入一個均勻分布的隨機數。從小到大搜索高斯累積分布，第一個大於均勻分布隨機數的的累積分布即為對應的高斯隨機數。重復3，即可產生符合同樣高斯分布的一系列隨機數 。
- 如果所需高斯隨機數的范圍有變化，需從第一步重新開始。
再看一下如何生成均勻分布的隨機數。最方便就是用C語言的rand()函數。在很多系統上RAND_MAX是32767，有時候范圍不太夠用，這樣很不方便。這里我用如下代碼生成。
unsigned long _RandomNumber = time(NULL); #define GET_NEXT_RANDOM (_RandomNumber = (_RandomNumber << 7) + (_RandomNumber >> 7))
只要初始數不為0，就可以連續生成一系列的偽隨機數。經過實驗，我覺得可以滿足一般使用。而優點就是，隨機數的范圍就是你設定的隨機數的類型所能取得的范圍。相對於rand()來說，運算速度更快。
生成高斯隨機數的C代碼如下，GaussianRandom()返回一個在[0,2r]的高斯隨機數。為了避免浮點數的比較，計算PDF和CDF都乘以一個大常數轉為整數。
static unsigned int *pGaussianCD = NULL; void GaussCDF(int radius) { unsigned int Weight; int j, n = 0; n = 2*radius + 1; if ((pGaussianCD = realloc(pGaussianCD, sizeof(unsigned int) * n)) == NULL) { printf("memory failure!\n"); exit(-1); } float sigma = radius / 3.0f; float sigma22 = 2*sigma*sigma; float sigma_sqrt2PI = sqrtf(2*PI)*sigma; Weight = (unsigned int)(expf(-(radius*radius)/sigma22) / sigma_sqrt2PI * 65536.0f); *pGaussianCD = Weight; n = 1; for (j = -radius + 1; j <= radius; j++) { Weight = (unsigned int)(expf(-(j*j)/sigma22) / sigma_sqrt2PI * 65536.0f); pGaussianCD[n] = Weight + pGaussianCD[n-1]; n++; } } /* Return a Gaussian random number between [0, 2*r], mean is r. */ unsigned int GaussianRandom(int radius) { static int r = 0, mn, m; int rand, i; if (r != radius) { r = radius; /* recalculate CDF */ GaussCDF(r); mn = 2*r; m = pGaussianCD[mn] + 1; } rand = GET_NEXT_RANDOM % m; for (i = 0; i <= mn; i++) { if (rand <= pGaussianCD[i]) return i; } /* should not reach here */ assert(0); return r; }
產生一個高斯隨機數主要時間都花在搜索輸入的均勻隨機數在高斯累計分布上的位置。產生高斯隨機數的范圍越大，相對花的時間越長。但是我們知道高斯分布的均值附近是產生隨機數概率最高的地方。如果從均值開始向兩側搜索，則有可能最快搜索到，就大大縮短了花費的時間。把最后一個搜索循環修改如下即可。
if (rand <= pGaussianCD[radius]) { i = radius - 1; while (i >= 0) { if (rand > pGaussianCD[i]) return i + 1; i--; } return 0; } else { i = radius + 1; while (i <= mn) { if (rand <= pGaussianCD[i]) return i; i++; } }
我們來看一下產生的高斯隨機數的分布如何，同時看看用移位加產生偽隨機數的方法是否可行。下面左邊的圖是用移位加的方法生成的[0,200]之間的24萬個偽隨機數和24萬個高斯隨機數的分布。作為對比右邊的圖是用rand()生成的，看起來沒什么大的差別。
移位加生成隨機數分布 rand()生成隨機數分布
移位加生成高斯隨機數分布 rand()生成高斯隨機數分布
這樣，我們就可以生成任意分布的隨機數，即使它的CDF不能用初等函數表達。下圖是生成24萬個泊松隨機數的分布，λ=7。因為離開λ的地方的概率下降很快，橫軸拉大了便於觀察。
生成泊松隨機數分布