生成符合高斯分佈或者其餘任意分佈的隨機數

在一些狀況下常常須要用到隨機數，而高斯隨機數又是最經常使用到的。這一篇講一下如何編程生成符合正態分佈的高斯隨機數，甚至任何其餘分佈的隨機數。

咱們知道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的無窮積分就能夠由有限的求和運算代替了。算法描述以下：

1. 設定高斯隨機數的範圍是[0,2r]，則均值是r，σ=r/3。
2. 由PDF計算截斷的機率分佈，而後求和計算累積分佈。
3. 輸入一個均勻分佈的隨機數。從小到大搜索高斯累積分佈，第一個大於均勻分佈隨機數的的累積分佈即爲對應的高斯隨機數。重複3，便可產生符合一樣高斯分佈的一系列隨機數 。
4. 若是所需高斯隨機數的範圍有變化，需從第一步從新開始。

再看一下如何生成均勻分佈的隨機數。最方便就是用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。由於離開λ的地方的機率降低很快，橫軸拉大了便於觀察。

      生成泊松隨機數分佈