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

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

public static double getExpRandomValue(double lambda) {return (-1.0/lambda)*Math.log(1-SecureRandom.nextDouble()); }

生成一組參數(shù)為lambda（λ）的負(fù)指數(shù)分布的隨機變量

通過下面的函數(shù)生成一組λ參數(shù)為lambda的隨機變量，其中size表示隨機變量的個數(shù)。通過該函數(shù)生成之后，可以將這些隨機值保存在文件中，以備分析和驗證，比如保存在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個λ參數(shù)為0.2的隨機變量，然后保存到exp.txt中，然后使用下面的matlab程序?qū)@些隨機數(shù)的性質(zhì)進行驗證，如果這些隨機數(shù)符合λ=0.2的負(fù)指數(shù)分布，則其均值應(yīng)為1/λ，即1/0.2=5，其方差應(yīng)為1/λ2=1/(0.2?0.2)=25。然后對這些隨機數(shù)的概率分布進行統(tǒng)計分析，以長度為1的區(qū)間為統(tǒng)計單位，統(tǒng)計各區(qū)間內(nèi)隨機數(shù)出現(xiàn)的頻數(shù)，求出在各區(qū)間的概率，繪制圖形，與參數(shù)為λ 的真實負(fù)指數(shù)分布曲線進行對比。以下為matlab代碼。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%測試以λ=0.2為參數(shù)的負(fù)指數(shù)分布 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%randomValues = load ('d:/exp.txt');%從文件導(dǎo)入生成的隨機數(shù) X = 1:1:80;%以長度為1的區(qū)間為統(tǒng)計單位，統(tǒng)計1~80內(nèi)的隨機數(shù)頻數(shù)m = mean(randomValues);%計算平均值，如果生成的隨機數(shù)正確，均值應(yīng)=1/λ=1/0.2=5 d = var(randomValues);%計算方差，方差應(yīng)=1/λ^2=1/(0.2^2)=25all_count = length(randomValues);%隨機數(shù)個數(shù)，方面后面將頻數(shù)轉(zhuǎn)成概率 [f,xout] = hist(randomValues, X);%按區(qū)間統(tǒng)計頻數(shù) for i=1:length(X)f(i) = f(i)/all_count;%頻數(shù)轉(zhuǎn)概率end;Y = 0.2*exp(-1*0.2*X);%畫出λ=0.2的負(fù)指數(shù)分布概密函數(shù)曲線 plot(X,f,X,Y,'r');%與隨機生成的概密函數(shù)曲線對比 grid on;%顯示格線 legend('統(tǒng)計曲線','實際曲線');%圖形注解title_str = sprintf('參數(shù):0.2 均值:%d 方差:%d', m, d); title(title_str);

如下圖所示，均值為4.996423，約等于5，方差為24.96761，約等于25，與實際情況相符。此外，通過matlab統(tǒng)計的概率密度函數(shù)曲線與真實曲線基本重合（其中在0-1之間沒有重合的原因是，實際情況是在0-1之間有無數(shù)個點，而matlab統(tǒng)計時以1為一個區(qū)間進行統(tǒng)計，只生成了一個統(tǒng)計項，而這無數(shù)個點的概率全部加到1點處，因此兩條線沒有重合，而且1點處的值遠(yuǎn)大于實際值，如果統(tǒng)計單位劃分越細(xì)，0-1之間的擬合度更高），表明生成的隨機數(shù)符合負(fù)指數(shù)分布。

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

通過下面的函數(shù)生成一組λ參數(shù)為lambda的k階愛爾蘭分布隨機變量，其中size表示隨機變量的個數(shù)，k表示階數(shù)。由于k階愛爾蘭分布是k個相同lambda的負(fù)指數(shù)分布的串聯(lián)，因此可以將連續(xù)k個負(fù)指數(shù)分布的隨機變量相加成為一個愛爾蘭分布的隨機變量，從而生成愛爾蘭分布的隨機變量，如下面程序所示。通過該函數(shù)生成之后，可以將這些隨機值保存在文件中，以備分析和驗證，比如保存在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個 λ 參數(shù)為0.2的2、4、8階愛爾蘭隨機變量，然后分別保存到erlang_2.txt、erlang_4.txt、erlang_8.txt中，然后使用下面的matlab程序?qū)@些隨機數(shù)的性質(zhì)進行驗證，驗證的方法與上面相同，對于k=2，則其均值應(yīng)為k/λ，即2/0.2=10，其方差應(yīng)為k/λ2=2/(0.2?0.2)=50；同理，對于k=4，均值應(yīng)等于20，方差應(yīng)等于100；對于k=8，均值應(yīng)等于40，方差應(yīng)等于200。下圖為matlab代碼。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%測試以λ=0.2為參數(shù)，K分別為2、4、8的愛爾蘭分布 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%randomValues_2 = load ('d:/erlang_2.txt');%從文件導(dǎo)入生成的λ=0.2 K=2的隨機數(shù) randomValues_4 = load ('d:/erlang_4.txt');%從文件導(dǎo)入生成的λ=0.2 K=4的隨機數(shù) randomValues_8 = load ('d:/erlang_8.txt');%從文件導(dǎo)入生成的λ=0.2 K=8的隨機數(shù) X = 1:1:80;%以長度為1的區(qū)間為統(tǒng)計單位，統(tǒng)計1~80內(nèi)的隨機數(shù)頻數(shù)m_2 = mean(randomValues_2);%計算平均值，如果生成的隨機數(shù)正確，均值應(yīng)=K/λ=2/0.2=10 d_2 = var(randomValues_2);%計算方差，方差應(yīng)=K/λ^2=2/(0.2^2)=50m_4 = mean(randomValues_4);%計算平均值，如果生成的隨機數(shù)正確，均值應(yīng)=K/λ=4/0.2=20 d_4 = var(randomValues_4);%計算方差，方差應(yīng)=K/λ^2=4/(0.2^2)=100m_8 = mean(randomValues_8);%計算平均值，如果生成的隨機數(shù)正確，均值應(yīng)=K/λ=8/0.2=40 d_8 = var(randomValues_8);%計算方差，方差應(yīng)=1/λ^2=8/(0.2^2)=200all_count_2 = length(randomValues_2);%隨機數(shù)個數(shù)，方面后面將頻數(shù)轉(zhuǎn)成概率 [f_2,xout_2] = hist(randomValues_2, X);%按區(qū)間統(tǒng)計頻數(shù) for i=1:length(X)f_2(i) = f_2(i)/all_count_2;%頻數(shù)轉(zhuǎn)概率end;all_count_4 = length(randomValues_4);%隨機數(shù)個數(shù)，方面后面將頻數(shù)轉(zhuǎn)成概率 [f_4,xout_4] = hist(randomValues_4, X);%按區(qū)間統(tǒng)計頻數(shù) for i=1:length(X)f_4(i) = f_4(i)/all_count_4;%頻數(shù)轉(zhuǎn)概率end;all_count_8 = length(randomValues_8);%隨機數(shù)個數(shù)，方面后面將頻數(shù)轉(zhuǎn)成概率 [f_8,xout_8] = hist(randomValues_8, X);%按區(qū)間統(tǒng)計頻數(shù) for i=1:length(X)f_8(i) = f_8(i)/all_count_8;%頻數(shù)轉(zhuǎn)概率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，以上結(jié)果都與實際情況符合。

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

通過下面的函數(shù)生成一組 λ 參數(shù)分別為lambda1和lambda2的2階超指數(shù)分布隨機變量，其中size表示隨機變量的個數(shù)，lambda1和lambda2表示兩個負(fù)指數(shù)分布的 λ 參數(shù)，這里指定進入分支1的概率為α1，進入分支2的概率為α2。由于2階超指數(shù)分布是2個 λ 參數(shù)分別為lambda1和lambda2的負(fù)指數(shù)分布的并聯(lián)，且以一定概率進入各分支，因此可以根據(jù)概率隨機的從兩個 λ 參數(shù)不同的負(fù)指數(shù)分布中抽取一個隨機變量作為一個超指數(shù)分布的隨機變量，如下面程序所示。通過該函數(shù)生成之后，可以將這些隨機值保存在文件中，以備分析和驗證，比如保存在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個 參數(shù)分別為0.2和0.5，α1=0.3、α2=0.7的超指數(shù)分布隨機變量，然后保存到hyper_exp.txt中，使用下面的matlab程序?qū)@些隨機數(shù)的性質(zhì)進行驗證，驗證的方法與上面相同，如果生成的隨機數(shù)正確，均值應(yīng)=α1/λ1+α2/λ2=0.3/0.2+0.7/0.5=2.9，方差應(yīng)=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階超指數(shù)分布 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%randomValues = load ('d:/hyper_exp.txt');%從文件導(dǎo)入生成的隨機數(shù) X = 1:1:80;%以長度為1的區(qū)間為統(tǒng)計單位，統(tǒng)計1~80內(nèi)的隨機數(shù)頻數(shù)m = mean(randomValues);%計算平均值，如果生成的隨機數(shù)正確，均值應(yīng)=α1/λ1+α2/λ2=0.3/0.2+0.7/0.5=2.9 d = var(randomValues);%計算方差，方差應(yīng)=2*(α1/λ1^2+α2/λ2^2)-(α1/λ1+α2/λ2)^2=12.19all_count = length(randomValues);%隨機數(shù)個數(shù)，方面后面將頻數(shù)轉(zhuǎn)成概率 [f,xout] = hist(randomValues, X);%按區(qū)間統(tǒng)計頻數(shù) for i=1:length(X)f(i) = f(i)/all_count;%頻數(shù)轉(zhuǎn)概率end;plot(X,f); grid on; % 顯示格線 title_str = sprintf('均值:%d 方差:%d', m, d); title(title_str);

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。