用计算机产生符合某中概率密度分布的高伪随机数
高伪随机数产生器及程序介绍
在做这个题目之前,我们需要知道伪随机数与真随机数有什么本质上的差别。真随机数是基于系统外部状态或物理上被认为随机的状态做出的随机数,比如对微观粒子的运动随机性进行描述时,得到的数的序列,就称其为真随机数。而伪随机数实际上是一种序列,是计算机根据种子所算出的所有的随机数字,而且在一定时间后(最大为模值)出现重复,就像用计算机模拟掷骰子的过程,该过程是均匀分布,但实际出现的数字是随机的,而用计算机模拟出来的结果在一定时间间隔后会重复,即具有周期性,我们可以用两张图来解释。
大气噪声 模拟大气噪声
图片来源于http://boallen.com/random-numbers.html
可以看出,右图在模拟大气噪声时产生了条纹,即具有周期性。这是伪随机数的局限。
事实上,在观测过程中由于物理偏差的存在,要想得到真随机数也是不可能的,而在应用问题中,比如说超均匀的分布比真随机数的均匀分布更合乎实际需要,这种伪随机数在实际应用中有着极为重要的作用,这也是我们为什么要研究它的原因。
接下来,我们列举一些符合某种特定概率密度分布的伪随机数生成器:
rand 均匀分布的随机数生成器
bertrand 贝塔分布的随机数生成器
binornd 二项分布的随机数生成器
normrnd 正态(高斯)分布的随机数生成器
exprnd 指数分布的随机数生成器
poissrnd 泊松分布的随机数生成器
raylrnd 瑞利分布的随机数生成器
gamrnd 伽玛分布的随机数生成器
geornd 几何分布的随机数生成器
chi2rnd 卡方分布的随机数生成器
frnd f分布的随机数生成器
hygernd 超几何分布的随机数生成器
lognrnd 对数正态分布的随机数生成器
nbinrnd 负二项分布的随机数生成器
ncfrnd 非中心f分布的随机数生成器
nctrnd 非中心t分布的随机数生成器
ncx2rnd 非中心卡方分布的随机数生成器
unidrnd 离散均匀分布的随机数生成器
unifrnd 连续均匀分布的随机数生成器
weibrnd 威布尔分布的随机数生成器
这些生成器都集成在matlab中的统计工具箱里面,接下来我们介绍一下指数分布的随机数生成器(exprnd)以及它的用法。
指数分布的概率密度函数为:
在matlab命令窗口里输入open exprand, 打开相应的m文件:
function r = exprnd(mu,varargin)
%EXPRND Random arrays from exponential distribution.
% R = EXPRND(MU) returns an array of random numbers chosen from the
% exponential distribution with mean parameter MU. The size of R is
% the size of MU.
%
% R = EXPRND(MU,M,N,...) or R = EXPRND(MU,[M,N,...]) returns an
% M-by-N-by-... array.
%
% See also EXPCDF, EXPFIT, EXPINV, EXPLIKE, EXPPDF, EXPSTAT, RANDOM.
% EXPRND uses the inversion method.
% References:
% [1] Devroye, L. (1986) Non-Uniform Random Variate Generation,
% Springer-Verlag.
% Copyright 1993-2009 The MathWorks, Inc.
if nargin < 1
error(message('stats:exprnd:TooFewInputs'));
end
[err, sizeOut] = statsizechk(1,mu,varargin{:});
if err > 0
error(message('stats:exprnd:InputSizeMismatch'));
end
% Return NaN for elements corresponding to illegal parameter values.
mu(mu < 0) = NaN;
% Generate uniform random values, and apply the exponential inverse CDF.
r = -mu .* log(rand(sizeOut)); % == expinv(u, mu)
其中,mu为均值,varargin为序列长度,前面的信息可以忽略,我们看一下r的命令行:
r = -mu .* log(rand(sizeOut)); % == expinv(u, mu)
从该命令行可以看出,生成一个具有指数分布的r矩阵,首先是使用rand命令生成一列规定序列长度的均匀分布的随机数,这些随机数定义为其随机变量的概率,再通过其log函数的值乘以均值最后生成该r矩阵,其产生原理与利用累加指数函数的逆函数x=expinv(u, mu)产生随机数相似,即x=F(p|μ)=-μln(1-p).
实际操作时,我们使用exprnd(mu,1,m)来产生一个一维序列长度为m的伪随机序列,接下来用matlab模拟该过程,在命令窗口输入:
R=exprnd(50,1,100);
plot((1:100),R)
mean=sum(R)/100
title(‘m=100’)
得到的结果是:
mean= 44.3180
可以看出,这里伪随机数的均值误差比较大
再在命令窗口输入以下代码:
R=exprnd(50,1,10000);
plot((1:10000),R)
mean=sum(R)/10000
title(‘m=10000’)
得到的结果如下:
mean= 50.5747
可见随机数的吻合程度已经比较高了。
经过排序后的r如下图所示,更直观地展示了指数分布的伪随机序列:
我们可以用均值的吻合程度来简单判断一下伪随机数的精确度。通过对不同的m取值刻画伪随机数的均值与理论值相比较,得到下图:
由该图可以得出这样一个结论:
1.要想获得高伪随机数,必须保证m(即样本个数足够大)
2.无论m取多大(有限),由于伪随机数的周期性,总是与理论值存在一定的误差。