投针法生成正态分布随机数

用投针法,将“均匀分布随机数”转化为“正态分布随机数”,附 C 语言源代码。

恩,主流程序设计语言都提供生成(伪)随机数的函数,比如 C 的 rand():

The rand() function returns a pseudo-random integer in the range 0 to RAND_MAX inclusive (i.e., the mathematical range [0, RAND_MAX]).

江湖传言,该函数给出的随机数服从均匀分布。果真如此吗……?恩,试一下不就知道了嘛。来十万个。

好极了。不过区间 [0, RAND_MAX] 并不实用。我们用以下函数把它拓展为任意的 [min, max]:

double randUniform(double min, double max){
	return (rand()/(RAND_MAX + 1.0)*((max)-(min))+(min));
}

下面我们将利用这个函数生成服从正态分布的随机数。我们打算用投针法。

原理是,先画出所要正态分布的概率密度曲线,然后拿一个矩形框起来,然后在这个矩形内随机取点(投针)。舍弃那些落在曲线上方的点。剩下的点的 x 坐标是我们要的。

从原理可以看出,本方法可以拓展到任意已知概率密度函数的分布。缺陷是,只能在有限区间内生成随机数,且并不能保证在有限步骤内结束

效果是极好的。比如我们生成 100000 个,保存到文件中,用 Mathematica 画出来,与标准正态曲线放在一起比较:

以下是 C 语言的源代码。

#include 
#include 
#include 
#include 

#define PI 4.0*atan(1.0)
#define E exp(1.0)

//uniform distribution in range [min, max]
double randUniform(double min, double max){
	return (rand()/(RAND_MAX + 1.0)*((max)-(min))+(min));
}
//probability density function of normal distribution
double pdfNormal(double x, double mu, double sigma){
	return 1/(pow(E,pow(-mu+x,2)/(2.0*pow(sigma,2)))*sqrt(2*PI)*sigma);
}
//normal distribution (mu, sigma), in range [-k*sigma, k*sigma]
double randNormal(double mu, double sigma, double k){
	double x, y, ymax=pdfNormal(mu,mu,sigma), xmax=k*sigma;
	while(233){
		x = randUniform(-xmax, xmax);
		y = randUniform(0,ymax);
		if(y
	

发表评论

电子邮件地址不会被公开。 必填项已用*标注