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)]);

 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