生成符合高斯分布或者其余任意分布的随机数

在一些状况下常常须要用到随机数,而高斯随机数又是最经常使用到的。这一篇讲一下如何编程生成符合正态分布的高斯随机数,甚至任何其余分布的随机数。web

咱们知道C语言的标准库函数能够生成符合均匀分布的伪随机数。那么如何生成符合高斯分布的随机数呢?咱们知道用逆函数法能够由符合(0,1)均匀分布的随机数获得符合任意分布的随机数,所以一样能够获得符合高斯分布的随机数。简单证实以下:算法

设随机变量u是符合(0,1)之间的均匀分布,那么有。设随机变量X的累积分布函数(CDF)为,其逆函数为。令随机变量。那么编程

所以只要获得所须要的随机分布的累计分布函数(CDF)的逆函数,就能够简单的经过逆函数把(0,1)均匀分布的随机数变换成所须要的随机数。高斯分布的几率分布函数(PDF)为dom

函数

其累计分布函数(CDF)是PDF的积分spa

,为从0到1的单调递增函数。3d

 

可是高斯分布的累积分布函数(CDF)不是初等函数,是无法用解析式表达的。再者符合高斯分布的随机变量的范围是从负无穷到正无穷的,是不可能用计算机生成的。即便用浮点数表示,也不是连续的。好比咱们有(0,1)之间均匀分布的随机数,经过高斯分布CDF的逆函数变换到正负无穷,也只是获得一些离散的点。code

所以要解决这些问题,首先用计算机生成的随机数确定是离散的。好比C语言的rand()函数返回[0,RAND_MAX]之间的伪随机整数。因此咱们也取必定范围的整数做为生成的随机高斯数。而后根据高斯分布的性质,离均值越远,几率越小,一般3σ之外的地方能够近似的认为几率为0。因此能够截断高斯分布的范围,让生成的高斯随机数位于[μ-3σ,μ+3σ]。这样CDF的无穷积分就能够由有限的求和运算代替了。算法描述以下:orm

  1. 设定高斯随机数的范围是[0,2r],则均值是r,σ=r/3。
  2. 由PDF计算截断的几率分布,而后求和计算累积分布。
  3. 输入一个均匀分布的随机数。从小到大搜索高斯累积分布,第一个大于均匀分布随机数的的累积分布即为对应的高斯随机数。重复3,便可产生符合一样高斯分布的一系列随机数 。
  4. 若是所需高斯随机数的范围有变化,需从第一步从新开始。

再看一下如何生成均匀分布的随机数。最方便就是用C语言的rand()函数。在不少系统上RAND_MAX是32767,有时候范围不太够用,这样很不方便。这里我用以下代码生成。blog

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。由于离开λ的地方的几率降低很快,横轴拉大了便于观察。

      生成泊松随机数分布

相关文章
相关标签/搜索