系统仿真（三）：Box Muller法生成正态分布随机量

Box Muller法生成正态分布随机量

——基于到达教学楼时间拟合函数

一、 实验背景

 结合第一次作业到达教学楼时间统计结果，采用原始843人次到达教学楼时间数据统计结果拟合函数，即$N(22.52，9.565)$。利用第二次作业多组合回归随机数生成器CMRG生成符合$N(22.52，9.565)$的正态分布函数.

二、 实验过程

2.1 采用Box-Muller方法

Box-Muller变换是通过服从均匀分布的随机变量，来构建服从高斯分布的随机变量的一种方法。具体的过程描述为：使用CMRG法生成两组服从[0,1]上均匀分布的随机变量$U_1、U_2$

CMRG法生成两组符合N(0,1)的随机数:

%函数参数
a1 = 1403580;a2 = 810728;a3 = 2.^32-209;
b1 = 527612;b2 = 1370589;b3 = 2.^32-22853;
z(1,1) = round(rand*1000000);z(1,2) = round(rand*1000000);z1(1,3) = round(rand*1000000);
z(2,1) = round(rand*1000000);z(2,2) = round(rand*1000000);z(2,3) = round(rand*1000000);
for i=4:randomNumber+3
    z(1,i)=mod(a1*z(1,i-2)-a2*z(1,i-3),a3);
    z(2,i)=mod(b1*z(2,i-1)-b2*z(2,i-3),b3);
    y(i)=mod(z(1,i)-z(2,i),a3);%i-3
    U1(i-3)=y(i)/a3;
end
z(1,1) = round(rand*1000000);z(1,2) = round(rand*1000000);z1(1,3) = round(rand*1000000);
z(2,1) = round(rand*1000000);z(2,2) = round(rand*1000000);z(2,3) = round(rand*1000000);
for i=4:randomNumber+3
    z(1,i)=mod(a1*z(1,i-2)-a2*z(1,i-3),a3);
    z(2,i)=mod(b1*z(2,i-1)-b2*z(2,i-3),b3);
    y(i)=mod(z(1,i)-z(2,i),a3);%i-3
    U2(i-3)=y(i)/a3;
end

为使$X_1$与$X_2$服从均值为0，方差为1的高斯分布，$X_1$、$X_2$满足
$$\begin{array}{l}{x}_1=\rho \cos \theta =-2\ln U_2\cos \left(2\pi U_1\right)\\ x_2=\rho \sin \theta =-2\ln U_2\sin \left(2\pi U_1\right),\end{array}$$
并进行$x^{\prime }=\mu +\sigma x$ 线性变换生成符合$N(22.52，9.565)$的随机数。

生成符合目标正态分布部分代码

mu = 22.52;
sigma = 9.565;
r = sqrt(-2*log(U2));
theta = 2*pi*U1;
X = [r.*cos(theta); r.*sin(theta)];
X1 = repmat(mu, 1, randomNumber) + sigma*X;

2.2 情况一

随机量选取与第一次实际采集843人次相同量

图1 原始数据柱状图

图2 随机数量为843的柱状图 我们可以观看原始数据图形与随机数生成图形，均出现不均匀现象。

2.3 情况二

随机量选取生成30000个，分别生成直方图和密度分布图

图3 随机数量为30000的柱状图

图4 随机数量为30000密度分布图

2.4 统计量检验

分别使用ks检验和卡方检验，均通过

图5 ks检验与卡方检验

三、 实验结论

通过Box Muller法可以生成的符合正态分布的随机变量，由情况一、二综合可知，当数据量相同为843时，原始数据图形与随机数生成图形均出现不均匀现象。当数据量达到30000时，图像较为均匀，故第一次实验出现波动有可能为统计量较小导致。

附录：

clear;clc;close all
randomNumber = 30000; %843
%函数参数
a1 = 1403580;a2 = 810728;a3 = 2.^32-209;
b1 = 527612;b2 = 1370589;b3 = 2.^32-22853;
z(1,1) = round(rand*1000000);z(1,2) = round(rand*1000000);z1(1,3) = round(rand*1000000);
z(2,1) = round(rand*1000000);z(2,2) = round(rand*1000000);z(2,3) = round(rand*1000000);
for i=4:randomNumber+3
    z(1,i)=mod(a1*z(1,i-2)-a2*z(1,i-3),a3);
    z(2,i)=mod(b1*z(2,i-1)-b2*z(2,i-3),b3);
    y(i)=mod(z(1,i)-z(2,i),a3);%i-3
    U1(i-3)=y(i)/a3;
end
z(1,1) = round(rand*1000000);z(1,2) = round(rand*1000000);z1(1,3) = round(rand*1000000);
z(2,1) = round(rand*1000000);z(2,2) = round(rand*1000000);z(2,3) = round(rand*1000000);
for i=4:randomNumber+3
    z(1,i)=mod(a1*z(1,i-2)-a2*z(1,i-3),a3);
    z(2,i)=mod(b1*z(2,i-1)-b2*z(2,i-3),b3);
    y(i)=mod(z(1,i)-z(2,i),a3);%i-3
    U2(i-3)=y(i)/a3;
end

mu = 22.52;
sigma = 9.565;
r = sqrt(-2*log(U2));
theta = 2*pi*U1;
X = [r.*cos(theta); r.*sin(theta)];
X1 = repmat(mu, 1, randomNumber) + sigma*X;

nbin = 100;
hist3(X1', [nbin nbin]);
set(gcf, 'renderer', 'opengl');
figure
hold on;
histogram(X1,70,'Normalization','pdf')
% 正态分布曲线
x=-10:1:40;
norm=normpdf(x,mu,sigma);
plot(x,norm,'linewidth',1)
%ks检验
pd = makedist('normal','mu',mu,'sigma',sigma);
h=kstest(X1,pd);
if h
    kstest='KS检验未通过'
else
    kstest='KS检验通过'
end