c语言学习进阶-C语言程序产生正态分布随机数

C语言程序产生正态分布随机数

（1）给出正态分布随机数的产生办法（算法）
在连续型的随机变量中，最简单而最基本的分布是单位均匀分布，由单位均匀分布抽取的简单例子称为随机数序列，其中每个随机序列个体称为随机数。利用物理方法在计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5。数值方法求正态分布常用中心极限定理，Hasiting有理逼近法以及反函数法。

中心极限定理（大数定理）

（一般$n\ge 10$），产生服从$N\left(\mu ,{\sigma }^2\right)$的算法：

• 产生n个均匀分布随机数：$r_1,r_2,\cdot \cdot \cdot ,r_n$，$E\left(r_i\right)=\frac{1}{2}$， $D\left(r_i\right)=\frac{1}{12}$；
• 计算$x=\frac{\left({\sum }_{i=1}^n{r}_i-\frac{n}{2}\right)}{\sqrt{\frac{n}{12}}}$；
• 计算$y=\sigma x+\mu$,y是服从$N\left(\mu ,{\sigma }^2\right)$分布的随机数。

注意：采用这种方法，随机数的个数n应该尽量大。
给出实现代码LLNgauss.c 如下：

#include <stdio.h>
#include <stdlib.h>
#include <math.h >
#define NSUM 25
double gaussrand()
{
	double x = 0,y=0;
	int i;
	for(i = 0; i < NSUM; i++)//NUSM表示均匀分布随机数个数再次设置为25个
	{
		x += (double)rand() / RAND_MAX;//生成25个均匀分布随机数并将所有随机数求和
	}
	x -= NSUM / 2.0;//求和后的随机数减
	x /= sqrt(NSUM / 12.0);//求得x
	y=sqrt((2.5))*x+1.2;//映射到y服从均值为1.2，方差为2.5的正态随机数
	return y;
 }
 int main()
{
	printf("%lf", gaussrand());
	return 0;
}
//验证
//gcc LLNgauss.c -o LLNgauss.exe
//LLNgauss.exe

中心极限定理生成正态随机数运行结果如图：

Hasiting有理逼近法：

• 构造分布函数的反函数$x=y-a_0+a_y+a_2{y}^2$，$x=y-1+b_{1}y+b_2{y}^2+b_3{y}^3$，这里$y={\left(-\text{2}\ln r\right)}^{\frac{1}{2}}$，$rU\left(\text{0,}1\right)$,$a_0=2.515517$，$a_1=0.802853$，$a_2=0.010328$，$b_1=1.432788$，$b_2=0.189269$，$b_3=0.001308$

反函数产生给定分布的随机数法：

（随机变量Y的概率密度为 $f\left(x\right)$ ）：

• 产生n个均匀分布随机数：$r_1,r_2,\cdot \cdot \cdot ,r_n$ ；
• 从等式$r_i={\int }_{-\infty }^{y_i}f\left(y\right)dy$中所得$y_i,i=\text{1,2,}\cdot \cdot \cdot ,n$即所求 。

Box-Muller法得到服从正态分布的随机数：

• 使用梅森旋转算法产生两个独立的服从[0,1)均匀分布随机数$r_1,r_2$ ；
• 将服从均匀分布的随机数转变为服从正态分布。
  令 $Y_1=\sqrt{-\text{2}\ln r_1}\cos \left(2\pi r_2\right)$ , $Y_2=\sqrt{-\text{2}\ln r_2}\cos \left(2\pi r_1\right)$ ， ${\left(Y_1,Y_2\right)}^T$ 服从二元正态分布，即利用两个独立的遵从均匀分布的随机数得到两个独立的正态分布的随机数。
  $$f\left(y\right)=\frac{1}{\sqrt{2\pi }}\exp \left\{-\frac{y^2}{2}\right\}$$
  正态值$Y_1$、$Y_2$的平均值为0，标准差为1，需要通过 $X=\mu +\left(Y\ast \sigma \right)$ 将 $Y_1$ 、 $Y_2$ 映射到均值为$\mu$标准差为$\sigma$的统计量X。
  $$f\left(x\right)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\frac{{\left(x-u\right)}^2}{2{\sigma }^2}\right\}$$
  给出实现代码BMgauss.c如下：

#include <stdlib.h> 
#include <stdio.h>
#include <math.h > 
double gaussrand( )
{
	static double U, V, y;//定义静态局部变量
	static int phase = 0;
	double z;
	double PI=3.141592654; 
	if(phase == 0)产生两个随机数U和V
	{
		U = rand() / (RAND_MAX + 1.0);
		V = rand() / (RAND_MAX + 1.0);
		z = sqrt(-2.0 * log(U))* sin(2.0 * PI * V);
	}
	else
	{
		z= sqrt(-2.0 * log(U)) * cos(2.0 * PI * V);
	}
	phase = 1 - phase;
	y=sqrt((2.5))*z+1.2;//映射到y服从均值为1.2，方差为2.5的正态随机数
	return y;
 }
int main()
{
	printf("%lf", gaussrand());
	return 0;
}
//验证
//gcc BMgauss.c -o BMgauss.exe
//BMgauss.exe

Box-Muller 法生成正态随机数运行结果如下图：

（2）采用命令行参数的办法输入参数，最后将产生结果写到一个结果文本文件当中，每行一个数。
命令行参数格式如：命令均值方差产生随机数的个数输出文档。例如：n.exe 1.2 2.5 8000result.txt 则产生均值为1.2, 方差为2.5 的8000 个数据，并写到指定的文本文件result.txt 中。
设计需求: 采用命令行参数定义均值，方差，输入行数，及输出文件；自定义生成正态分布随机数函数；调用生成正态分布随机数函数生成数据写入文件。
实现思路: 根据LLNgauss.c 定义高斯随机数生成函数（相应的函数参数从命令行参数获取）；主函数定义5 个命令行参数命令、均值、方差、输入行数、输出文件，并且创建文件写入流，循环写入生成的随机数，直到随机数的个数足够为止。
给出实现代码n.c 如下：

#include <stdio.h>
#include <stdlib.h>
#include <math.h >
#define NSUM 25
//自定义正态随机数生成的函数
double gaussrand(double mu,double sigma2)
{
	double x = 0,y=0,Y=0;
	int i;
	for(i = 0; i < NSUM; i++)//NUSM表示均匀分布随机数个数再次设置为25个
	{
		x += (double)rand() / RAND_MAX;//生成25个均匀分布随机数并将所有随机数求和
	}
	x=x- (NSUM / 2.0);//求和后的随机数减
	Y =x/ sqrt(NSUM / 12.0);//求得x
	y=sqrt(sigma2)*Y+mu;//映射到y服从均值为1.2，方差为2.5的正态随机数
	return y;
 }
 int main(int argc,char* argv[])
{
	double mu;
	sscanf(argv[1],"%lf",&mu);argv[3]是字符串格式，从一个字符串中读进双精度类型数据
	double sigma2=atof(argv[2]);//从命令行参数获取方差
	double data=0;//定义数据
	int line=atoi(argv[3]);//定义数据生成行数
	FILE* fp1=fopen(argv[4],"w");
	//检查命令行参数是否正确
	if(argc!=5)
	{
		printf("参数输入不正确！");	
		printf("提示：命令 均值 方差 数据行数 输出文件名\n");
		return -1;
	}
	//产生指定行数数据
	for(int i=0;i<line;i++)
	{
		data=gaussrand(mu,sigma2);//调用自定义函数生成正态随机数数据
		fprintf(fp1, "%lf\n", data);//写入文件
	}
	printf("输出成功\n");
	fclose(fp1);//关闭数据输出流
	if (fclose(fp1)!=0);//关闭文件成功关闭则返回值为0
	{
		printf("文件 %s 关闭错误\n",argv[4]);
	}
}
//验证
//gcc n.c -o n.exe
//n.exe 1.2 2.5 8000 result.txt

命令行参数生成8000 个正态随机数运行结果如下图：

将8000 个正态随机数写入文件result.tct 如下图：

matlab 验证生成随机数文档数据正确性

验证c 语言产生的正态分布的随机数是按参数要求产生的，从而验证程序的正确性。在matlab 中利用统计直方图画出c 语言产生的这些数据的统计值方图，观察是否符合正态分布的概率密度函数的形状。

将result.txt 的计算结果数据，设法读到matlab 中，利用求平均和求方差的函数计算这些结果数据的均值和方差是否达到了要求。
设计需求：将result.txt 的数据读入matlab；计算数据的均值与方差；绘制数据直方图。
matlab 导入数据可在“主页”-菜单栏选择“导入数据”-选择矩阵-导入所选内容，另一种是使用代码data=load(’result.txt’) 或data=textread(’result.txt’)，textread 语句将文本数据赋值给矩阵或data = textscan(fileID,formatSpec,N) 将已打开的文本文件中的数据读取到data 数组读取N 次缺省时读一次。

实现思路：利用matlab 中的textread() 函数将文件数据写入一个矩阵data 中；利用函数mean() 求均值、函数std() 求方差、函数max() 求最大值、函数min() 求最小值、函数rand()、median() 求中位数、函数mode() 求众数（由于随机数小数位数较多在取众数前对数据进行小数位数第二位四舍五入操作）；利用函数histogram() 画出数据的直方图。
matlab 验证代码如下：

close all;
clear all;
clc;
data=textread('result.txt');
mu=mean(data);
sigma2=(std(data))*(std(data));
MAX=max(data);
MIN=min(data);
RAND=sort(data,1,'descend');
MEDIAN=median(RAND);
ROUND=roundn(RAND,-1);%向小数点后1位四舍五入
MODE=mode(ROUND);
histogram(data,(-10:10));
fprintf('mu=%f\n',mu);
fprintf('sigma2=%f\n',sigma2);
fprintf('MAX=%f\n',MAX);
fprintf('MIN=%f\n',MIN);
fprintf('MEDIAN=%f\n',MEDIAN);
fprintf('MODE=%f\n',MODE);

matlab 计算数据均值mu、方差σ2，最大值MAX，最小值MIN，中位数MEDIAN、众数MODE 如图:

数据分布的直方图如图:

c 语言验证生成随机数文档数据正确性

用c 语言实现上一步matlab 中求平均和求方差的计算过程，将result.txt 的数据读到c 语言的程序中，然后进行计算，通过命令行打印出来计算结果的均值和方差。
例如js.exe result.txt 将显示类似如下面的信息: 平均值=1.245 方差=2.5 数据个数=8000数据中的最大值=4.663 数据中的最小值=-1.008
设计需求：将存数据的文件名设置为命令行参数；采用c 语言实现对文本数据操作；编程实现求数据的平均值，方差，标准差，数据个数，数据中的最大值，数据中的最小值，数据的中值与众数。
实现思路：设置一个命令行参数定义文件名，并给出参数输入报错的提示；创建文件读入流读取文件数据并赋值给变量；采用简单四则运算计算数据统计值；采用冒泡排序实现数的排序，找中位数众数。
c 语言验证代码如下：

#include<stdio.h>
#include<math.h>
#include<string.h>
#include<stdlib.h>
int main(int argc,char*argv[])
{
	//读文件操作
	FILE*fp1;
	if(argc!=2)
	{
		printf("参数输入不正确！");
		printf("提示：命令 输入文件");
	}
	fp1=fopen(argv[1],"r");
	if(fp1==NULL)
	{
		printf("文件打开失败！");
		return 0;
	}
	fclose(fp1);
	//计行操作
	fp1=fopen(argv[1],"r");
	int c,lines=0;
	while((c=fgetc(fp1))!=EOF)
	{
		if (c=='\n')
		{
			lines++;
		}
	}
	lines=lines+1;
	fclose(fp1);
	printf("数据个数lines=%d\n",lines);
	//赋值给数组
	int i;
	fp1=fopen(argv[1],"r");
	double data[lines];
	for( i=0;i<lines;i++)
	{
		fscanf(fp1,"%lf",&data[i]);
		//printf("data[%d]=%lf\t",i,data[i]);
	}
	fclose(fp1);
	//计算均值
	double mu;
	double sum1=0;
	for( i=0;i<lines;i++)
	{
		sum1=sum1+data[i];
	}
	mu=sum1/lines;
	printf("均值mu=%lf\n",mu);
	//计算方差
	double sigma2;
	double sum2=0;
	for( i=0;i<lines;i++)
	{
		sum2=sum2+pow((data[i]-mu),2.0);
	}
	sigma2=sum2/lines;
		printf("方差sigma2=%lf\n",sigma2);
	//求标准差
	double sigma;
	sigma=sqrt(sigma2);
	printf("标准差sigma=%lf\n",sigma);
	//数组排序输出最大值与最小值
	double a[lines];
	int j;
	double temp;//定义冒泡空盒
	for( i=0;i<lines;i++)
	{
	a[i]=data[i];	
	}
	for(i=0;i<lines-1;++i)//比较lines-1轮
	{
		for(j=0;j<lines-1-i;++j)//每轮比较lines-1-i次
		{
			if(a[j]<=a[j+1])//降序
			{
				temp=a[j];
				a[j]=a[j+1];
				a[j+1]=temp;
			}	
		}
	}
	//for( i=0;i<lines;i++)
	//{
	//	printf("a[%d]=%lf\n",i,a[i]);
	//}
	printf("最大值MAX=%lf\n",a[0]);
	printf("最小值MIN=%lf\n",a[lines-1]);
	//求中位数
	double MEDIAN;//定义中位数
	if (lines%2==0)
	{
		MEDIAN=	(a[(lines/2)-1]+a[(lines/2)+1])/2;
	}
	else
	{
		MEDIAN=	a[(lines+1)/2];
	}
	printf("中位数=%lf\n",MEDIAN);
	
}
	
//验证
//gcc js.c -o js.exe
//js.exe result.txt

c 语言计算数据均值mu、方差σ2，最大值MAX，最小值MIN，中位数MEDIAN 如下图:

c语言根据给出区间宽度对随机数画直方图

理解根据给出区间宽度对随机数画直方图的原理，例如对学生成绩以每10 分一个分数段统计得分人数进行直方图绘制的原理。对result.txt 中的所有数据（可把这些数据看作学生的得分成绩），编写c 语言命令行参数程序，生成指定区间宽度的直方图绘制数据，然后通过origin 等绘图工具绘制出直方图。与matlab 中绘制的直方图进行对比。
设计需求：将输入数据文件与输出数据文件作为命令行参数fig.exe result.txt result.txt；将存随机数的的文件作为输入文件；采用c 语言实现对文本数据操作；编程实现数据分段；将分段数据写入新文件result1.txt 中。
实现思路：设置三个命令行参数（数据分段步长、输入文件文件名、输出文件名文件名），并给出参数输入报错的提示；创建文件读出与写入流读取文件数据并赋值给数组；找数据区间端点，采用判断语句实现数据区间内命令行指定的步长参数对数据分段；在文件中输出两列数一列是区间，另外一列是随机数的个数。

c 语言实现分段统计随机数代码如下：

#include<stdio.h>
#include<math.h>
#include<string.h>
#include<stdlib.h>
int main(int argc,char*argv[])
{
	//读文件操作
	FILE*fp1,*fp2;
	if(argc!=4)
	{
		printf("参数输入不正确！");
		printf("提示：命令 步长 输入文件");
	}
	fp1=fopen(argv[2],"r");
	if(fp1==NULL)
	{
		printf("文件打开失败！");
		return 0;
	}
	fclose(fp1);
	//计行操作
	fp1=fopen(argv[2],"r");
	int c,lines=0;
	while((c=fgetc(fp1))!=EOF)
	{
		if (c=='\n')
		{
			lines++;
		}
	}
	lines=lines+1;
	fclose(fp1);
	printf("数据个数lines=%d\n",lines);
	//赋值给数组
	int i;
	fp1=fopen(argv[2],"r");
	double data[lines];
	for( i=0;i<lines;i++)
	{
		fscanf(fp1,"%lf",&data[i]);
		//printf("data[%d]=%lf\t",i,data[i]);
	}
	fclose(fp1);
	//数组排序输出最大值与最小值
	double a[lines];
	int j;
	double temp;//定义冒泡空盒
	for( i=0;i<lines;i++)
	{
	a[i]=data[i];	
	}
	for(i=0;i<lines-1;++i)//比较lines-1轮
	{
		for(j=0;j<lines-1-i;++j)//每轮比较lines-1-i次
		{
			if(a[j]<=a[j+1])//降序
			{
				temp=a[j];
				a[j]=a[j+1];
				a[j+1]=temp;
			}	
		}
	}
	printf("最大值MAX=%lf\n",a[0]);
	printf("最小值MIN=%lf\n",a[lines-1]);
	double MAX = ceil(a[0]);//区间最大值
	double MIN = floor(a[lines-1]);//区间最小值
	double area = MAX - MIN;;//区间范围
	double step = atof(argv[1]);//区间要用命令行参数
	double num = area / step;//算出需要几个区间
	printf("区间最大值:%.0f\n", MAX);
	printf("区间最小值:%.0f\n", MIN);
	//int MAX = (int)MAX;
	//int MIN = (int)MIN;
	//int step = (int)step;
	for(double j=MIN;j<MAX;j=j+step)//区间遍历
	{
		int counts=0;//初始化存放区间数目值
		for(int i=0;i<lines-1;i++)//随机数元素遍历
		{
			if(a[i]>=j && a[i]<=j+step)//判断落入区间
			{
				counts++;//计数
			}
		}
		fp2= fopen(argv[3], "a");
		printf("区间%lf-%lf随机数个数：%d\n",j,j+step,counts);
		//写入文件
		//int x = (int)j+step;
		fprintf(fp2, "%.1f %d\n", j,counts);
		fclose(fp2);
	}
}
//验证
//gcc fig.c -o fig.exe
//fig.exe 1 result.txt result1.txt

c 语言实现分段统计随机数代码运行结果如图：

Origin 绘图可导入dat 文件也可导入txt 文件，导入时只需把数据文件拖入Book1 中，也可采用菜单栏的文件导入，导入数据后再绘图菜单栏选择直方图、散点图、折线图。
采用origin 以生成的result1 为数据绘制统计直方图及曲线界面如下图：

origin 以生成的result1 为数据绘制统计直方图及曲线如图：

总结：

在做此任务时还是遇到了一些困难：
在生成随机数时，采用中心极限定理计算，在验证时验证出了问题，发现数据误差较大，返工断点检查，检查出问题是在定义变量类型不够严谨，最后根据数据输出要求修改了数据类型，输出了较准确的结果。
生成随机数程序每次运行时会产生一个随机数，当在命令行运行多次时，应当产生服从正态分布的不同的随机数，但是实际发现每次产生的随机数都是同一个数，而为何将生成随机数嵌入循环多次生成随机数时又和预期一致产生了不同的服从正态分布的随机数。