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
基本代码都在了，报告自己写吧。