隨機(jī)數(shù)產(chǎn)生原理及實(shí)現(xiàn)

1、電子信息與通信工程學(xué)院實(shí)驗(yàn)報告實(shí)驗(yàn)名稱隨機(jī)數(shù)的產(chǎn)生課程名稱隨機(jī)信號分析姓名顧康學(xué)號U201413323日期6月6日地點(diǎn)南一樓東204成績教師董燕 以上為6種分布的實(shí)驗(yàn)結(jié)果1. 均勻分布 隨機(jī)變量XU(0，1)的一組樣本值的模擬值一般采用某種數(shù)值計(jì)算方法產(chǎn)生隨機(jī)數(shù)序列，在計(jì)算機(jī)上運(yùn)算來得到，通常是利用遞推公式: Xn=f(Xn-1,.,Xn-k)1.1 同余法 Xn+1 = Xn(mod M) Rn=Xn/M R1 R2.Rn 即為（0,1）上均勻分布的隨機(jī)數(shù)列。而上述方法是偽隨機(jī)的，Rn本質(zhì)上是遞推公式給定的周期序列，周期T可看做log（ M）。 解決方法是 ：選擇模擬參數(shù) 并 對序列進(jìn)行統(tǒng)計(jì)

2、檢驗(yàn)。 1.2選擇模擬參數(shù) 1）周期長度取決于Xo,， M的選擇 2）通過選取適當(dāng)?shù)膮?shù)可以改善隨機(jī)數(shù)的性質(zhì) 幾組參考的取值 Xo =1 , =7 , M=1010 Xo =1 , =513 , M=2 *1010 Xo =1 , =517 , M=10121.3對數(shù)列進(jìn)行統(tǒng)計(jì)檢驗(yàn) 對應(yīng)序列能否看作X的獨(dú)立同分布樣本，須檢驗(yàn)其 獨(dú)立性 和 均勻性 for i=2:1:size %同余法 均勻分布 x(i)= mod ( v*x(i-1), M); y(i)=x(i)/M;endsubplot(2,3,1);hist(y,100)ahat,bhat,ACI,BCI=unifit(y)% 以0.9

3、5的置信度估計(jì)樣本的參數(shù) 首先我們的標(biāo)準(zhǔn)是U （0，1），而實(shí)驗(yàn)值，ACI表示ahat的范圍-0.0030,0,BCI表示bhat的范圍1.0000,1.0030。同時樣本的均值和方差分別為0.4932和0.0830，結(jié)論與理論值很接近。該樣本以0.95的可信度服從（0,1）均勻分布。2. 伯努利分布2.1算法原理 若隨機(jī)變量R服從（0,1）， P(X=Xi)=Pi P(0)=0, P(n)=Pi PP(n-1)<R<=P(n)=P(n)-P(n-1)=Pn 令P(n-1)<X<=P(n)=X=Xn 有P(X=Xn)=Pn從理論上講，已經(jīng)解決了產(chǎn)生具有任何離散型隨機(jī)分布

4、的問題。具體執(zhí)行仍有困難，如X取值為無窮時。2.2算法 對于伯努利分布只需用到上述算法最簡單的情形，即取n為2。%伯努利分布k1=0,k2=1p1=0.2,p2=0.8;r=zeros(1,size);for j=1:1:size if y(j)>p1 r(j)=k2; else r(j)=k1; endend subplot(2,3,2);hist(r) title('伯努利分布'); PHAT,PCI=binofit(r,1000)%以0.95的置信度檢驗(yàn)樣本參數(shù) PHAT=0.198，而PCI= 0.195,0.212 ，而我設(shè)置的P=0.2，與實(shí)驗(yàn)結(jié)果十分接近，可

5、見該樣本的性質(zhì)較好。該樣本以0.95的可信度服從0.2的伯努利分布。3. 正態(tài)分布3.1算法 設(shè)有兩個在（0,1）上獨(dú)立均勻分布的隨機(jī)數(shù)R1 ,R1;作如下變換 Y1= （-2R1）½ COS(2R2） Y2= （-2R1）½ SIN(2R2）其逆變換為R1=exp(- (Y1²+Y2²)/2)R2=1/2 arctag(Y2/Y1)可導(dǎo)出Y1，Y2的聯(lián)合分布函數(shù)f(Y1,Y2)= 1/2 exp(- (Y1²+Y2²)/2) 故知Y1，Y2相互獨(dú)立且服從N（0,1）再作變換 Xi = Yi+,可得到服從N(，)的X3.2代碼%正態(tài)隨

6、機(jī)分布y=rand(1,size);z=rand(1,size); m=sqrt(-2*log(y).*cos(2*pi.*z);n=sqrt(-2*log(y).*sin(2*pi.*z);t=m,n;subplot(2,3,3);hist(t,100)title('正態(tài)分布');muhat,sigmahat,muci,sigmaci=normfit(t)%以0.95置信度檢驗(yàn)樣本參數(shù)實(shí)驗(yàn)結(jié)果：muhat=0.0060, sigmahat=1.0080; Muci=-0.0382,0.0502 ,sigmaci=0.9777,1.0402;理論上假設(shè)的是標(biāo)準(zhǔn)正態(tài)分布，則該樣本

7、以0.95的可信度服從標(biāo)準(zhǔn)正態(tài)分布。4. 指數(shù)分布4.1算法 已知指數(shù)分布的分布函數(shù)為 F=1-exp(-x),F0,1。而利用服從U(0，1)的變量Y, 帶入方程左邊反解出x： X= -log(1-y)/ 即我們得到了用均勻分布產(chǎn)生指數(shù)分布的方法4.2代碼 %指數(shù)分布y=rand(1,size); lamda=3; x=-log(1-y)/lamda; subplot(2,3,4);hist(x,30)title('指數(shù)分布')muhat,muci=expfit(x)%以0.95的置信度檢驗(yàn)樣本參數(shù)理論上取為3，而實(shí)驗(yàn)值 muhat=0.3204, muci=0.3015,0

8、.3402 ,可得到該樣本以0.95的可信度服從E(3). 5. 泊松分布5.1算法 利用一組在（0,1）獨(dú)立且同均勻分布的變量X，首先理解泊松分布的到達(dá)間隔服從指數(shù)分布： 在一次獨(dú)立實(shí)驗(yàn)中 Xk exp(-)作為判斷成功條件，輸出為k，即X的次數(shù)。5.2代碼%泊松分布 lamda=20;p=exp(-lamda); y=zeros(1,100); for cnt=1:1:100 i=0; q=1; while q>=p q=q*rand(1); i=i+1; end y(cnt)=i-1;endsubplot(2,3,5);hist(y,25)title('泊松分布')

9、 lamdahat,lamdaci=possifit(y)%以0.95的置信度檢驗(yàn)樣本參數(shù)注意到lamda初值被賦為20，實(shí)驗(yàn)值 lamdahat=20.1900, Lamdaci=19.3093,21.0707.所以樣本以0.95的可信度服從P(20).6. 幾何分布6.1算法原理 先由幾何分布的定義出發(fā) P(X=k)= q(k-1) * p;則可利用與產(chǎn)生泊松分布隨機(jī)數(shù)相似的方法，從均勻分布入手。 產(chǎn)生一組（0,1）獨(dú)立地均勻分布隨機(jī)數(shù)，分別與P比較。如果Xi<P,則繼續(xù)X（i+1）；但如果XkP,則輸出K（代表比較的次數(shù)）6.2代碼%幾何分布 y=zeros(1,100);p=0.

10、6;q=0.4;test=0;cnt=1; while cnt<=100 if test=1 x=rand(1,size); end for i=1:1:size if x(i)>p y(cnt)=i; cnt=cnt+1; test=1; break end endendsubplot(2,3,6);hist(y,20)title('幾何分布')phat,pci=mle(y,'distribution','geometric')%極大似然估計(jì)函數(shù)，以0.95置信度預(yù)測p 我給P的初值為0.6，但是值得注意的現(xiàn)象是，多次運(yùn)算得到的phat都在0.3左右。當(dāng)改變P值后，phat仍然會分布在另一個值附近。通過pci的范圍我們可以判斷樣本以0.95的可信度屬于幾何分布，但概率P并不與初值一致。 我認(rèn)為這可能是算法導(dǎo)致的，并無法保證初值P仍