高斯随机数产生原理及代码

第一种是近似法生成，由中心极限定理可知独立同分布的多个随机变量和的分布趋近于正态分布。

设期望(expection)为Ex,方差(variance)为Sigma^2，则下面的这个得到的value符合标准正态分布:

sum=X1+X2+...+Xn;

value=(sum-n*Ex)/(sqrt(n)*sigma);

取k个在[0,1]均匀分布的随机变量x1,x2,...,xk，则它们的和近似服从正态分布。

实践中，一般取k=12，（因为D(xi)=1/12），则新的随机变量

y＝x1＋x2＋...＋x12－6可以满足一般精度下的N(0,1)分布的要求(因E(y)=0，D(y)=12*1/12＝1）。

代码如下：

double gauss() /* generate the random number */

{

int i;

double x=0.0;

for(i=0;i<12;i++)

{

x=x+rand()/RAND_MAX;

}

x=x-6.0;

return x;

}

下面给出符合正态分布N(expection,sigma^2)的随机数产生的思想和程序代码.

设X,Y,Z为3个随机变量,且X~N(0,1),Y~N(expection,sigma^2).考虑这个变换: Z=sigma*X+expection;

其期望为expection,方差为sigma^2.因此,我们可以通过上面的处理,由X~N(0,1)来

得到特定期望和方差的正态分布.

具体程序如下(按自己理解写的):

double Gauss(double expection,double sigma)

{

int i;

double x=0.0;

for(i=0;i<12;i++)

{

x=x+rand()/RAND_MAX;

}

x=x-6.0;

x=x*sigma+expection;

return x;

}

//下面的程序有点疑问,该程序摘自网上.

double gauss(doulbe mean,double sigam,long int s)

{

int i;

double x,y;

for(x=0.0;i=0;i<50;i++)

{

x=x+uniform(0.0,1.0,s)

x=x-n/12;

y=mean+sigma*x;

}

}

方法二：Box 和Muller 在1958 年给出了由均匀分布的随机变量生成正态分布的随机变量的算法。

设U1, U2 是区间(0, 1) 上均匀分布的随机变量，且相互独立。令

X1 = sqrt(-2*log(U1)) * cos(2*PI*U2);

X2 = sqrt(-2*log(U1)) * sin(2*PI*U2);

那么X1, X2 服从N(0,1) 分布，且相互独立。等于说我们用两个独立的U(0,1) 随机数

得到了两个独立的N(0,1)随机数。

double carmen_gaussian_random(double mean, double std)

{

const double norm = 1.0 / (RAND_MAX + 1.0);

double u = 1.0 - rand() * norm; /**//* can't let u == 0 */

double v = rand() * norm;

double z = sqrt(-2.0 * log(u)) * cos(2.0 * M_PI * v);

return mean + std * z;

}

方法三：Marsaglia 和Bray 在1964 年提出了一种改进算法，避免使用三角函数。

以下的实现代码用的就是这种改进算法。

// Gaussian Random Number Generator class

// ref. “Numerical Recipes in C++ 2/e”, p.293 ~ p.294

public class GaussianRNG

{

int iset;

double gset;

Random r1, r2;

public GaussianRNG()

{

r1 = new Random(unchecked((int)DateTime.Now.Ticks));

r2 = new Random(~unchecked((int)DateTime.Now.Ticks));

iset = 0;

}

public double Next()

{

double fac, rsq, v1, v2;

if (iset == 0)

{

do{

v1 = 2.0 * r1.NextDouble() - 1.0;

v2 = 2.0 * r2.NextDouble() - 1.0;

rsq = v1*v1 + v2*v2;

}while (rsq >= 1.0 || rsq == 0.0);

fac = Math.Sqrt(-2.0*Math.Log(rsq)/rsq);

gset = v1*fac;

iset = 1;

return v2*fac;

}

else {

iset = 0;

return gset;

}

}

}