creathist.m

实现正态分布的伪随机数发生器 实现正态分布的伪随机数发生器

function creathist(mean1,sigma1,mean2,sigma2,scalar)

%实现正态分布的伪随机数发生器
% CREATHIST(mean1,sigma1,mean2,sigma2,scalar)
% 函数作用：生成两个相互独立且满足Gauss分布的随机变量的样本
% 入口参数：mean1－－样本1的均值
% sigma1－－样本1的标准差
% mean2－－样本2的均值
% sigma2－－样本2的标准差
% scalar－－直方图的数目
% 调用格式：1）creathist()--生成两个分别满足N(2,-2)、N(-2,2)且相互独立
% 的随机变量，hist等分成100份。
% 2) creathist(mean1,sigma1,mean2,sigma2)--设置默认hist等分成
% 100份。
% 3) creathist(mean1,sigma1,mean2,sigma2,scalar)

% 出错处理
if nargin==0 % 输入参数为0时的设置
    mean1=2;
    sigma1=-2;
    mean2=-2;
    sigma2=2;
    scalar=100; 
elseif nargin==4 % 输入参数为4时的设置
    scalar=100;
elseif nargin~=5 % 输入参数不是0、4、5时的出错处理
        error('请用help creathist查看函数允许使用的格式');
end

% 下载两个满足U(0,1)的均匀分布的样本
% 每个样本容量为100000，文件大小为1.35M
% 样本数据的生成是采用C语言编程实现
% 算法采用线性同余法
load('data1.txt','-ascii');
load('data2.txt','-ascii');

% 画出直方图，采用hist默认等分成10份
% 也就是说，在样本最小值和最大值范围内等分成10份
% 落在每一份中的样本点总数作为纵坐标绘制直方图
% 这也是为什么要采用hist而不用plot的原因所在
subplot(3,1,1);hist(data1);legend('U1(0，1)样本');
subplot(3,1,2);hist(data2);legend('U2(0，1)样本');

% 利用均匀分布生成两个相互独立且满足N(0,1)的高斯分布的随机变量
% 采用的原理如下：任意分布－－>均匀分布－－>任意分布
% 设X满足均匀分布，Y满足瑞利分布，则Fy(Y)=X
% 两边取反函数运算求出Y的表达式Y=finverse(Fy(X))
% Z=Ycos(2*pi*X)即满足正态分布
x1 = sqrt(-2*log(1-data1)) .* cos(2*pi*data2);
x2 = sqrt(-2*log(1-data1)) .* sin(2*pi*data2);

% 生成满足N(mean,sigma^2)的随机变量
y1=sigma1*x1+mean1;
y2=sigma2*x2+mean2;
% 绘制直方图
subplot(3,1,3);hist(y1,scalar);legend('图形验证');
hold on;
hist(y2,scalar);
hold off;

% 根据样本数据进行验证
fprintf('根据样本数据进行验证\n');

% 方法1：E(xy)=E(x)E(y)
% 说明：对结果采取了适当精确度的判定，本例所选精度为0.01
fprintf('方法1：E(xy)=E(x)E(y)\n');
a=mean(mean(y1.*y2)); % 特别注意mean的使用方法
                           % 明确此处可以执行的理论推导
b=mean(mean(y1))*mean(mean(y2));
fprintf('E(xy)=%f\n',a);
fprintf('E(x)E(y)=%f\n',b);
fprintf('e(n)=|E(xy)-E(x)E(y)|=%f\n',abs(a-b));
if abs(a-b) < 0.01
    fprintf('E(xy)=E(x)E(y)，两个随机变量相互独立\n');
else
    fprintf('E(xy)!=E(x)E(y)，两个随机变量不相互独立\n');
end

% 方法2：cov(x,y)=0
% 说明：此法实质上与1相同，主要想从相关系数的角度来分析；
% 而且使用cov更为简单一些。本例所选精度为0.01，应该与
% 方法1测试结果相同。
fprintf('\n方法2：cov(x,y)=0\n');
c=cov(y1,y2);
fprintf('cov(x,y)=%f\n',c(2));
fprintf('e(n)=|p(x,y)|=%f\n',abs(c(2)));
if abs(c(2)) < 0.01
    fprintf('cov(x,y)=0，两个随机变量相互独立\n');
else
    fprintf('cov(x,y)!=0，两个随机变量不相互独立\n');
end