使用C语言产生符合高斯分布的随机数基础理论、示例教程

拟采用Box-Muller模型产生符合高斯分布的随机数，Box-Muller模型可以产生均匀分布的随机数，然后将均匀分布的随机数通过方法变为符合其他分布的随机数。

一、Box-Muller使用方法如下：

比如：使用两个等式中的任何一个计算出符合正态分布的随机数Z：[BaiduBaike]

若（0,1 ] 内存在两个独立的随机数U1、U2：

或

其中：

可知，正态值Z有值为0的均值与值为1的标准偏差。

所以将Z映射到X中得：

二、正文

原文链接见http://www.design.caltech.edu/erik/Misc/Gaussian.html

此文主要讲述产生符合高斯分布的伪随机数，文章底部的参考文献也是相关话题的文章，大家也可以深入学习一下。

解决这个问题的方法很多，但是这篇文章仅仅使用一个非常重要的模型。如果我们可以推理出高斯分布的方程式，我们就可以借助我们的方程式通过公式间的基本变换可以变换出我们所需要的分布规律，通过这种分布规律，我们输入随机数，就可以得到符合我们需要的随机数。当然最著名的公式准则就是大家都比较熟悉的Box-Muller方法。Box-Muller方法可以将符合一种分布的随机数变换成一种新的随机数（当然就可以产生高斯分布、均匀分布）。

大多数的基本变换形式如下：

我们使用x1、x2（0到1之间的均匀分布）通过转换得到高斯分布的均值与标准差。当使用这种方式产出成千上万的随机数的时候，这种方式就有很多的弊端，比如说：

（1）、因为需要使用很多计算库，所以计算速度很慢

（2）、当X1特别接近0的时候，得到的数据误差比较大

但是Box-muller方法无论是计算速度还是数据误差都稳定于之前的基本变换方法。

Box-muller的c语言变换形式如下：

         float x1, x2, w, y1, y2;
 
         do {
                 x1 = 2.0 * ranf() - 1.0;
                 x2 = 2.0 * ranf() - 1.0;
                 w = x1 * x1 + x2 * x2;
         } while ( w >= 1.0 );

         w = sqrt( (-2.0 * log( w ) ) / w );
         y1 = x1 * w;
         y2 = x2 * w;

注意：rand()函数可以产生值在[0,1]之间的均匀分布的随机数。由于没有涉及到sin、cos等三角函数，所以这种方法的计算速度是非常快的。

完整的C语言产生高斯分布的程序如下：

/* boxmuller.c           Implements the Polar form of the Box-Muller
                         Transformation

                      (c) Copyright 1994, Everett F. Carter Jr.
                          Permission is granted by the author to use
			  this software for any application provided this
			  copyright notice is preserved.

*/

#include 


extern float ranf();         /* ranf() is uniform in 0..1 */


float box_muller(float m, float s)	/* normal random variate generator */
{				        /* mean m, standard deviation s */
	float x1, x2, w, y1;
	static float y2;
	static int use_last = 0;

	if (use_last)		        /* use value from previous call */
	{
		y1 = y2;
		use_last = 0;
	}
	else
	{
		do {
			x1 = 2.0 * ranf() - 1.0;
			x2 = 2.0 * ranf() - 1.0;
			w = x1 * x1 + x2 * x2;
		} while ( w >= 1.0 );

		w = sqrt( (-2.0 * log( w ) ) / w );
		y1 = x1 * w;
		y2 = x2 * w;
		use_last = 1;
	}

	return( m + y1 * s );
}

博主翻译有点渣，此文将继续完善。