利用反函数及变换抽样法生成随机数

2011-09-12 01:00王丙参魏艳华

重庆高教研究 2011年5期

王丙参，魏艳华，张云

(天水师范学院数学与统计学院，甘肃天水 741001)

用随机模拟方法解决实际问题时，首先要解决的是随机数的产生方法，或称r.v的抽样方法.目前关于随机数的生成的文献很多［1－3］，但基本没有系统介绍反函数及变换抽样法的.Sas软件是最专业的统计软件，可处理各种数理统计问题，进行数据分析.本文结合Sas软件研究U［0，1］随机数的生成，运用反函数及变换抽样法生成各种非均匀随机数，并分析了它们的优缺点.

1 ［0，1］上均匀分布随机数的生成

定义1［4］若 r.v X 的cdf为 F(x)，则 X 的一个样本值称为一个F随机数，若F(x)有pdf f(x)时，也称为f随机数.U［0，1］的n个独立样本称为n个均匀随机数，简称随机数.

非均匀随机数是由 U［0，1］随机数经相应运算而产生的，因此U［0，1］随机数的质量决定了非均匀随机数的质量，生成U［0，1］随机数主要有3 种方法［5］:

(1)手工方法，即采用掷骰子、抽签、抽牌等，许多彩票的发行至今仍采用这种方法.

(2)物理方法，在计算机上安装一台物理随机数发生器，把具有随机性质的物理过程变换为随机数，这样可得到真正的随机数且取之不尽.但缺点是速度慢，无法重复且对统计模拟带来不可验证性，加上物理随机数发生器需经常检查和维修，因而大大降低了这种方法的使用价值.

(3)利用数学方法生成伪随机数.

定理 1［4］设 Xi，i=1，2… ，i.i.d 于 B(1，0.5)，则Y=≜0.X1X2… ～ U(0，1).

此定理给出了产生U［0，1］随机数的方法.目前应用最广泛的随机数发生器是线性同余发生器(LCG)

其中M为模数，a为乘子，c为增量且都为非负整数.当c≠0，上式称为混合同余发生器;当c=0时，称为乘同余发生器，此时当模为素数时，称它为素数模乘同余发生器.伪随机数具有以下特点:

1)近似性［1］.由LCG算法可知，随机数不相互独立，每个数都依赖于前一个数，但如果这种相关关系弱到一定程度就可以认为是相互独立的.我们用离散数据，…，1代替线段［0，1］，设k为表示一个整数值的二进制位数，则m=2k，若m足够大，这种近似是完美的.

2)周期性.理想随机数序列的周期应该是无限长的，但由一定算法产生的伪随机数必然有一定的周期性，从而导致随机数出现一定的相关性和重复性.周期短的直接后果是:按统计模型模拟的结果不可靠，即模拟样本也存在周期性，不符合模型的假定，但周期与计算机的精度直接相关，精度决定了最长周期.

Sas系统中生成U［0，1］随机数的函数有:UNIFORM(seed)其乘子为16 807，模为231的乘同余发生器和一个64位数的搅乱表形成的组合发生器，seed必须是常数，它一般是0，或5位、6位、7位的奇数.RANUNI(seed)其乘子为397 204 094，模为231－1的素数发生器，seed必须是小于模231－1 的任何常数［6］.

2 利用反函数及变换抽样法生成非均匀随机数

如果cdf F(x)严格单调，U ～ U［0，1］，则P(F-1(U)≤x)=P(U≤F(x))=F(x)，即F-1(u)是一个F随机数，其中 u ～ U［0，1］.但很多cdf并非严格单调如离散型r.v，不存在逆函数，故可定义广义的逆函数

F-1(u)=inf{x:F(x)≥u}，

并有:

定理2［5］如果F(x)是cdf，u ～ U［0，1］，则F-1(u)=inf{x:F(x)≥u}是一个F随机数.

显然，若X～G(x)，则Y=F-1(G(X))～F(x)，即由已知G(x)随机数可以得到∀F(x)随机数.为叙述方便，文中的u及ui都是均匀随机数.

设离散r.v cdf为F(x)，pdf为f(xi)=P(X=xi)=pi，将［0，1］分为一些互不相交的子区间，使第i个子区间Ji的长度为pi.任取一u～U［0，1］，若 u∈Ji，令 z=xi，则 z ～ F ，Sas命令为 rantbl(seed，p1，…，pn).若 F(x)为 {1，2，…，n}上离散均匀r.v，由于

则 z= ［nu］+1为{1，2，…，n}上离散均匀随机数;

若X ～Ge(p)，即

f(k)=P(X=k)=qk-1p，k≥1，q=1-p，

由于

因此

若 X ～ P(λ)，Uki.i.d于 U(0，1)，则 Tk=-ln Uk是i.i.d于Exp(1)，那么

{Nt=sup{k:τk=T1+ … +Tk≤ t}，t≥0}为强度为1的齐次p.p，Nλ=m⇔τm≤λ ＜ τm+1，即则满足的m是泊松随机数.

若X ～ P(λ)，即有

Sas命令为ranpoi(seed，lambda)，于是有产生随机数X的算法:

(1)置 k:=0，pk=e-λ和 p(k)=0;

(2)产生 u ～ U［0，1］;

(3)令pk+1=和p(k+1)=p(k)+pk+1;

(4)若p(k)＜ u≤p(k+1)，令x=k+1;否则令k:=k+1，返回(3).

如果0=t0＜t1＜… ＜tn=1，ti-ti-1=pi，i≤n，F(x)= ∑ni=1piFi(x)，其中Fi(x)为cdf，u为随机数，Z≜.若ti-1＜U≤ ti，则

若 F(x)=1-e-λx，则 z=- λ-1ln(1-u)是Exp(λ)随机数.n个独立的均值为θ的指数随机数的和为 Γ(n，θ)随机数.注意 RANEXP(seed)产生 Exp(1)随机数，y=RANEXP(seed)/λ，则产生 Exp(λ)随机数;若 y=α－β*LOG(RANEXP(seed))，则y为具有位置参数α和尺度参数为β的极值分布随机数;若

y=FLOOR(－RANEXP(seed)/LOG(p))，则y～Ge(p);若x1～Γ(n，1)，x2～Γ(m-1)且相互独立，则是Beta(n，m)随机数;若，则～ F(m，n).

注意若 x=rangam(seed，α)，则 y=x/β 为 Γ(α，β)随机数，即Sas系统中默认β=1.

若 X ～ tr(a，b，m)，pdf

当a=0，b=1，它称为标准三角形分布.若X1～ tr(0，1，c)，0 ≤ c≤1 ，则有 a+(b-a)X1～tr(a，b，m)，其中 m=a+(b-a)c.故只需考虑怎样生成tr(0，1，c)的随机数.Sas系统中三角分布随机数函数为Rantri(seed，h)，其中0＜h＜ 1，即 tr(0，1，h)，tr(a，b，c)随机数 y 可由tr(0，1，h)随机数的线性变换得到，即:y=(ba)*Rantri(seed，h)+a，h=(c-a)/(b-a)，c∈［a，b］.