利用反函數及變換抽樣法生成隨機數

2011-09-12 01:00:26王丙參魏艷華

重慶高教研究 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］.