Matlab實現。
分為主函數 MyLine 和被調用函數 Func。
主函數 MyLine 實現在 Func 函數的基礎上實現序貫法， 將平均等待隊長作為每次模擬的 X，求出置信區間。
Func 函數實現單次模擬排隊系統，返回向量 vector。向量vector分別包括仿真平均排隊時間Cus_Queue_avg，仿真平均等待時間Cus_Wait_avg，仿真系統中平均等待隊長 Cus_Wait_Queue_avg，仿真系統中平均顧客數 Cus_Wait_CusNum_avg。
主函數 MyLine：
1 clear; 2 clc; 3 Alpha=0.1; %置信水平 4 Gama=0.1; %相對精度 5 Beta=0.1; 6 Lambda=0.2; %到達率 Lambda 7 Mu=0.25; %服務率 Mu 8 time=50; %單回模擬次數 9 %序貫法實現 10 All_vector=Func( Lambda,Mu );%函數返回的向量 11 for i=2:time 12 All_vector=[All_vector;Func( Lambda,Mu )]; 13 end 14 vect=sum(All_vector)/time;%各列求和取平均 15 %下面系統中平均等待隊長作為每次模擬的 X，用來評估運行結果 16 S=sum((All_vector(3)-vect(3)).*(All_vector(3)-vect(3)))/(ti 17 me-1);%后面假設 S 不變 18 Betan=sqrt(S/time)*tinv(1-Alpha/2,time-1); 19 Gaman=Betan/vect(3); 20 while(Gaman>=Gama) %Betan>=Beta|| 21 time=time+1; 22 All_vector=[All_vector;Func( Lambda,Mu )]; 23 vect=sum(All_vector)/time; 24 S=sum((All_vector(3)-vect(3)).*(All_vector(3)-vect(3)))/(ti 25 me-1);%后面假設 S 不變 26 Betan=sqrt(S/time)*tinv(1-Alpha/2,time-1); 27 Gaman=Betan/vect(3); 28 end 29 time 30 disp(['系統中平均等待隊長的置信區間下界 31 =',num2str(vect(3)-Betan)]); 32 disp(['系統中平均等待隊長的置信區間上界 33 =',num2str(vect(3)+Betan)]);
Func 函數：
1 function [ vector ] = Func( Lambda,Mu ) 2 %單次的排隊模擬，樣本數 CusTotal 3 CusTotal=10000; %仿真顧客總數%=input('請輸入仿真顧客 4 總數 CusTotal='); 5 Cus_Arrive=zeros(1,CusTotal);%到達時刻 6 Cus_Leave=zeros(1,CusTotal);%離開時刻 7 IntervaCus_Arrive=-log(rand(1,CusTotal))/Lambda;%到達時間間 8 隔 9 Cus_Arrive=cumsum(IntervaCus_Arrive); %每列 10 累加，形成到達初始時刻；如果只有一行，該行向后疊加 11 Interval_Serve=-log(rand(1,CusTotal))/Mu; %服務時間間隔 12 %為事件調度法做准備 13 Cus_Leave(1)=Cus_Arrive(1)+Interval_Serve(1);%顧客離開時間 14 for i=2:CusTotal 15 if Cus_Leave(i-1)<Cus_Arrive(i) 16 Cus_Leave(i)=Cus_Arrive(i)+Interval_Serve(i); 17 else 18 Cus_Leave(i)=Cus_Leave(i-1)+Interval_Serve(i); 19 end 20 end 21 Cus_Wait=Cus_Leave-Cus_Arrive; %各顧客在系統中的等待時間 22 %mean：如果是 n*m 的矩陣，mean 對各列分別進行求平均；當 n=1 時， 23 對一行求平均 24 Cus_Wait_avg=mean(Cus_Wait); %平均等待時間 25 Cus_Queue=Cus_Wait-Interval_Serve;%各顧客在系統中的排隊時間 26 Cus_Queue_avg=mean(Cus_Queue); %平均排隊時間 27 %TimePoint 系統調度時間 28 TimePoint=[Cus_Arrive,Cus_Leave];%系統中顧客數隨時間的變化 29 TimePoint=sort(TimePoint); 30 CusNum=zeros(size(TimePoint)); 31 temp=2; %指向事件表 32 CusNum(1)=1; 33 %統計 dt 時間內的排隊人數——事件調度法 34 %截止到 i，還有多少個事件 35 for i=2:length(TimePoint) 36 %后一事件的結束時間晚於前一事件的結束時間 37 %所以只要第 temp-1 個事件沒有結束， 后面的事件到了發生時間 38 (事件仿真鍾<=系統仿真鍾)，都要加入 CusNum 中計數 39 if 40 (temp<=length(Cus_Arrive))&&(TimePoint(i)==Cus_Arrive(temp) 41 ) 42 CusNum(i)=CusNum(i-1)+1; 43 temp=temp+1; 44 else 45 CusNum(i)=CusNum(i-1)-1; 46 end 47 end 48 %系統中平均顧客數計算 49 Time_interval=zeros(size(TimePoint)); 50 Time_interval(1)=Cus_Arrive(1); 51 for i=2:length(TimePoint) 52 Time_interval(i)=TimePoint(i)-TimePoint(i-1); 53 end 54 %這里畫一下圖即可.i 時間段*(i-1)時刻統計的事件數量。類似於向 55 后積分 56 CusNum_fromStart=[0 CusNum]; 57 Cus_Wait_CusNum_avg=sum(CusNum_fromStart.*[Time_interval 58 0] )/TimePoint(end); 59 %系統平均等待隊長 60 QueLength=zeros(size(CusNum)); 61 for i=1:length(CusNum) 62 if CusNum(i)>=2 63 %還有等待事件（滿足事件仿真鍾<=系統仿真鍾，但沒有發 64 生） 65 QueLength(i)=CusNum(i)-1; 66 else 67 QueLength(i)=0; 68 end 69 end 70 Cus_Wait_Queue_avg=sum([0 QueLength].*[Time_interval 71 0] )/TimePoint(end); 72 %仿真值與理論值比較 73 row=Lambda/Mu; 74 QueLength_avg=row*row/(1-row);%Q 75 CusNum_avg=row/(1-row);%L 76 Queue_avg=QueLength_avg/Lambda;%d 77 Wait_avg=CusNum_avg/Lambda;%w 78 %返回的向量 79 vector=[Cus_Queue_avg,Cus_Wait_avg,Cus_Wait_Queue_avg,Cus_W 80 ait_CusNum_avg]; 81 end
基本代碼都在了，報告自己寫吧。