1、概率特性仿真實驗與程序-Matlab仿真-隨機數生成-負指數分布-k階愛爾蘭分布-超指數分布使用Java中的SecureRandom.nextDouble()生成一個01之間的隨機浮點數，然后使用反函數法生成一個符合指數分布的隨機變量（反函數求得為）。指數分布的參數為getExpRandomValue函數中的參數lambda。生成一個指數分布的隨機變量的代碼如下，后面都將基于該函數生成一組負指數分布、K階愛爾蘭分布、2階超指數分布隨機變量，然后將生成的隨機數通過matlab程序進行仿真，對隨機數的分布特性進行驗證。public static double getExpRandomValue(d

2、ouble lambda)return (-1.0/lambda)*Math.log(1-SecureRandom.nextDouble(); 生成一組參數為lambda（）的負指數分布的隨機變量通過下面的函數生成一組參數為lambda的隨機變量，其中size表示隨機變量的個數。通過該函數生成之后，可以將這些隨機值保存在文件中，以備分析和驗證，比如保存在exp.txt文件中，供下面介紹的matlab程序分析。public static double genExp(int size, double lambda)double array = new doublesize; while(-size

3、>=0) arraysize = getExpRandomValue(lambda); return array;通過genExp(1000000, 0.2)生成1000000個參數為0.2的隨機變量，然后保存到exp.txt中，然后使用下面的matlab程序對這些隨機數的性質進行驗證，如果這些隨機數符合=0.2的負指數分布，則其均值應為1/，即1/0.2=5，其方差應為1/2=1/(0.2*0.2)=25。然后對這些隨機數的概率分布進行統計分析，以長度為1的區間為統計單位，統計各區間內隨機數出現的頻數，求出在各區間的概率，繪制圖形，與參數為的真實負指數分布曲線進行對比。下圖為matla

4、b代碼%測試以=0.2為參數的負指數分布%randomValues = load ('d:/exp.txt');%從文件導入生成的隨機數X = 1:1:80;%以長度為1的區間為統計單位，統計180內的隨機數頻數m = mean(randomValues);%計算平均值，如果生成的隨機數正確，均值應=1/=1/0.2=5d = var(randomValues);%計算方差，方差應=1/2=1/(0.22)=25all_count = length(randomValues);%隨機數個數，方面后面將頻數轉成概率f,xout = hist(randomValues, X);%按

5、區間統計頻數for i=1:length(X) f(i) = f(i)/all_count;%頻數轉概率end;Y = 0.2*exp(-1*0.2*X);%畫出=0.2的負指數分布概密函數曲線plot(X,f,X,Y,'r');%與隨機生成的概密函數曲線對比grid on;%顯示格線legend('統計曲線','實際曲線');%圖形注解title_str = sprintf('參數:0.2 均值:%d 方差:%d', m, d);title(title_str);如下圖所示，均值為4.996423，約等于5，方差為24.9676

6、1，約等于25，與實際情況相符。此外，通過matlab統計的概率密度函數曲線與真實曲線基本重合（其中在0-1之間沒有重合的原因是，實際情況是在0-1之間有無數個點，而matlab統計時以1為一個區間進行統計，只生成了一個統計項，而這無數個點的概率全部加到1點處，因此兩條線沒有重合，而且1點處的值遠大于實際值，如果統計單位劃分越細，0-1之間的擬合度更高），表明生成的隨機數符合負指數分布。生成一組參數為lambda（）的k階愛爾蘭分布的隨機變量通過下面的函數生成一組參數為lambda的k階愛爾蘭分布隨機變量，其中size表示隨機變量的個數，k表示階數。由于k階愛爾蘭分布是k個相同lambda的負

7、指數分布的串聯，因此可以將連續k個負指數分布的隨機變量相加成為一個愛爾蘭分布的隨機變量，從而生成愛爾蘭分布的隨機變量，如下面程序所示。通過該函數生成之后，可以將這些隨機值保存在文件中，以備分析和驗證，比如保存在erlang_k.txt文件中，供下面介紹的matlab程序分析。public static double genErlang(int size, double lambda, int k) double array = new doublesize; while(-size>=0) for(int i = 0; i<k; i+) arraysize += getExpRan

8、domValue(lambda); return array;通過genErlang(1000000, 0.2, 2)、genErlang(1000000, 0.2, 4)、genErlang(1000000, 0.2, 8)分別生成1000000個參數為0.2的2、4、8階愛爾蘭隨機變量，然后分別保存到erlang_2.txt、erlang_4.txt、erlang_8.txt中，然后使用下面的matlab程序對這些隨機數的性質進行驗證，驗證的方法與上面相同，對于k=2，則其均值應為k/，即2/0.2=10，其方差應為k/2=2/(0.2*0.2)=50；同理，對于k=4，均值應等于20，方

9、差應等于100；對于k=8，均值應等于40，方差應等于200。下圖為matlab代碼。%測試以=0.2為參數，K分別為2、4、8的愛爾蘭分布%randomValues_2 = load ('d:/erlang_2.txt');%從文件導入生成的=0.2 K=2的隨機數randomValues_4 = load ('d:/erlang_4.txt');%從文件導入生成的=0.2 K=4的隨機數randomValues_8 = load ('d:/erlang_8.txt');%從文件導入生成的=0.2 K=8的隨機數X = 1:1:80;%以長度為

10、1的區間為統計單位，統計180內的隨機數頻數m_2 = mean(randomValues_2);%計算平均值，如果生成的隨機數正確，均值應=K/=2/0.2=10d_2 = var(randomValues_2);%計算方差，方差應=K/2=2/(0.22)=50m_4 = mean(randomValues_4);%計算平均值，如果生成的隨機數正確，均值應=K/=4/0.2=20d_4 = var(randomValues_4);%計算方差，方差應=K/2=4/(0.22)=100m_8 = mean(randomValues_8);%計算平均值，如果生成的隨機數正確，均值應=K/=8/0

11、.2=40d_8 = var(randomValues_8);%計算方差，方差應=1/2=8/(0.22)=200all_count_2 = length(randomValues_2);%隨機數個數，方面后面將頻數轉成概率f_2,xout_2 = hist(randomValues_2, X);%按區間統計頻數for i=1:length(X) f_2(i) = f_2(i)/all_count_2;%頻數轉概率end;all_count_4 = length(randomValues_4);%隨機數個數，方面后面將頻數轉成概率f_4,xout_4 = hist(randomValues_4

12、, X);%按區間統計頻數for i=1:length(X) f_4(i) = f_4(i)/all_count_4;%頻數轉概率end;all_count_8 = length(randomValues_8);%隨機數個數，方面后面將頻數轉成概率f_8,xout_8 = hist(randomValues_8, X);%按區間統計頻數for i=1:length(X) f_8(i) = f_8(i)/all_count_8;%頻數轉概率end;plot(X,f_2,'r',X,f_4,'g',X,f_8,'b');str1 = sprintf(

13、'k:2 m:%d d:%d', m_2, d_2);str2 = sprintf('k:4 m:%d d:%d', m_4, d_4);str3 = sprintf('k:8 m:%d d:%d', m_8, d_8);legend(str1,str2,str3); %圖形注解 如下圖所示，k=2時，均值為9.992167，約等于10，方差為49.93048，約等于50；k=4時，均值為20.00298，約等于20，方差為100.4140，約等于100；k=8時，均值為40.03118，約等于40，方差為200.4146，約等于200，以上結果

14、都與實際情況符合。生成一組2階超指數分布的隨機變量通過下面的函數生成一組參數分別為lambda1和lambda2的2階超指數分布隨機變量，其中size表示隨機變量的個數，lambda1和lambda2表示兩個負指數分布的參數，這里指定進入分支1的概率為1，進入分支2的概率為2。由于2階超指數分布是2個參數分別為lambda1和lambda2的負指數分布的并聯，且以一定概率進入各分支，因此可以根據概率隨機的從兩個參數不同的負指數分布中抽取一個隨機變量作為一個超指數分布的隨機變量，如下面程序所示。通過該函數生成之后，可以將這些隨機值保存在文件中，以備分析和驗證，比如保存在hyper_exp.txt

15、文件中，供下面介紹的matlab程序分析。public static double genHyperExp(int size, double lambda1, double lambda2)double a1 = 0.3;/a1:進入分支1的概率 因此a2=1-a1=0.7double array = new doublesize; while(-size>=0) if(SecureRandom.nextDouble()>a1) arraysize = getExpRandomValue(lambda2); else arraysize = getExpRandomValue(la

16、mbda1); return array;通過genHyperExp(1000000, 0.2, 0.5)生成1000000個參數分別為0.2和0.5，1=0.3、2=0.7的超指數分布隨機變量，然后保存到hyper_exp.txt中，使用下面的matlab程序對這些隨機數的性質進行驗證，驗證的方法與上面相同，如果生成的隨機數正確，均值應=1/1+2/2=0.3/0.2+0.7/0.5=2.9，方差應=2*(1/12+2/22)-(1/1+2/2)2=12.19。下圖為matlab代碼。%測試1=0.2、1=0.5、1=0.3、2=0.7的2階超指數分布%randomValues = load ('d:/hyper_exp.txt');%從文件導入生成的隨機數X = 1:1:80;%以長度為1的區間為統計單位，統計180內的隨機數頻數m = mean(randomValues);%計算平均值，如果生成的隨機數正確，均值應=1/1+2/2=0.3/0.2+0.7/0.5=2.9d = var(randomValues);%計算方差，方差應=2*(1/12+2/22)-(1/1+2/2)2=12.19all_count = length(rando