概率特性仿真實驗與程序-Matlab仿真-隨機數生成-負指數分布-k階愛爾蘭分布-超指數分布

使用Java中的SecureRandom.nextDouble()生成一個0~1之間的隨機浮點數，然后使用反函數法生成一個符合指數分布的隨機變量（反函數求得為x=?ln(1?R)λ）。指數分布的參數λ為getExpRandomValue函數中的參數lambda。生成一個指數分布的隨機變量的代碼如下，后面都將基于該函數生成一組負指數分布、K階愛爾蘭分布、2階超指數分布隨機變量，然后將生成的隨機數通過matlab程序進行仿真，對隨機數的分布特性進行驗證。

public static double getExpRandomValue(double 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 double[size];while(--size>=0) { array[size] = 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的區間為統計單位，統計各區間內隨機數出現的頻數，求出在各區間的概率，繪制圖形，與參數為λ 的真實負指數分布曲線進行對比。以下為matlab代碼。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%測試以λ=0.2為參數的負指數分布 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%randomValues = load ('d:/exp.txt');%從文件導入生成的隨機數 X = 1:1:80;%以長度為1的區間為統計單位，統計1~80內的隨機數頻數m = mean(randomValues);%計算平均值，如果生成的隨機數正確，均值應=1/λ=1/0.2=5 d = var(randomValues);%計算方差，方差應=1/λ^2=1/(0.2^2)=25all_count = length(randomValues);%隨機數個數，方面后面將頻數轉成概率 [f,xout] = hist(randomValues, X);%按區間統計頻數 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.96761，約等于25，與實際情況相符。此外，通過matlab統計的概率密度函數曲線與真實曲線基本重合（其中在0-1之間沒有重合的原因是，實際情況是在0-1之間有無數個點，而matlab統計時以1為一個區間進行統計，只生成了一個統計項，而這無數個點的概率全部加到1點處，因此兩條線沒有重合，而且1點處的值遠大于實際值，如果統計單位劃分越細，0-1之間的擬合度更高），表明生成的隨機數符合負指數分布。

生成一組參數為lambda（λ）的k階愛爾蘭分布的隨機變量

通過下面的函數生成一組λ參數為lambda的k階愛爾蘭分布隨機變量，其中size表示隨機變量的個數，k表示階數。由于k階愛爾蘭分布是k個相同lambda的負指數分布的串聯，因此可以將連續k個負指數分布的隨機變量相加成為一個愛爾蘭分布的隨機變量，從而生成愛爾蘭分布的隨機變量，如下面程序所示。通過該函數生成之后，可以將這些隨機值保存在文件中，以備分析和驗證，比如保存在erlang_k.txt文件中，供下面介紹的matlab程序分析。

public static double[] genErlang(int size, double lambda, int k){double[] array = new double[size];while(--size>=0) { for(int i = 0; i<k; i++)array[size] += getExpRandomValue(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，方差應等于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;%以長度為1的區間為統計單位，統計1~80內的隨機數頻數m_2 = mean(randomValues_2);%計算平均值，如果生成的隨機數正確，均值應=K/λ=2/0.2=10 d_2 = var(randomValues_2);%計算方差，方差應=K/λ^2=2/(0.2^2)=50m_4 = mean(randomValues_4);%計算平均值，如果生成的隨機數正確，均值應=K/λ=4/0.2=20 d_4 = var(randomValues_4);%計算方差，方差應=K/λ^2=4/(0.2^2)=100m_8 = mean(randomValues_8);%計算平均值，如果生成的隨機數正確，均值應=K/λ=8/0.2=40 d_8 = var(randomValues_8);%計算方差，方差應=1/λ^2=8/(0.2^2)=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, 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('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，以上結果都與實際情況符合。

生成一組2階超指數分布的隨機變量

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

public static double[] genHyperExp(int size, double lambda1, double lambda2){double a1 = 0.3;//a1:進入分支1的概率 因此a2=1-a1=0.7double[] array = new double[size];while(--size>=0) { if(SecureRandom.nextDouble()>a1)array[size] = getExpRandomValue(lambda2);elsearray[size] = getExpRandomValue(lambda1);}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/λ1^2+α2/λ2^2)-(α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的區間為統計單位，統計1~80內的隨機數頻數m = mean(randomValues);%計算平均值，如果生成的隨機數正確，均值應=α1/λ1+α2/λ2=0.3/0.2+0.7/0.5=2.9 d = var(randomValues);%計算方差，方差應=2*(α1/λ1^2+α2/λ2^2)-(α1/λ1+α2/λ2)^2=12.19all_count = length(randomValues);%隨機數個數，方面后面將頻數轉成概率 [f,xout] = hist(randomValues, X);%按區間統計頻數 for i=1:length(X)f(i) = f(i)/all_count;%頻數轉概率end;plot(X,f); grid on; % 顯示格線 title_str = sprintf('均值:%d 方差:%d', m, d); title(title_str);

如下圖所示，均值為2.896629，約等于2.9，方差為12.17702，約等于12.19，以上結果與實際情況符合。

總結

以上是生活随笔為你收集整理的概率特性仿真实验与程序-Matlab仿真-随机数生成-负指数分布-k阶爱尔兰分布-超指数分布的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。