sim.py

# 基于《统计模拟和R实现》一书例子
import numpy as np
from numpy import random as rm
import pandas as pd

def simulation(T=4200):
    t = 0;nA =0;nD=0;n=0;A=[];D=[];N=[];S=[]
    tA = rm.exponential(10,1) # 客来时间
    tD = inf  # 客走时间
    while True:
        # print tA,tD,n,S
        if  (tA <= tD)  & (tA <= T): # 客来时点早于客走时点，且在营业
            t = tA # 更新当前时间
            nA = nA + 1 # 更新来客人数
            n = n+1 # n保存服务系统中人数
            tA = t + rm.exponential(10,1) # 产生新的客来时间
            if n==1: # 系统中只有一个客人时才更新tD
                tS = rm.uniform(5,10,1)
                tD = t + tS# 客走时间等于客走时间加服务时间
                S.append(tS) 
            A.append(t) # 保存客来时间序列        
        elif (tD <= tA) & (tD <=T): #如果客来时点晚于客走时点
            t = tD; n= n-1
            nD = nD + 1
            if n==0:  tD = inf #如果无人等待
            else:  
                tS = rm.uniform(5,10,1)
                tD = t + tS
                S.append(tS) # 保存每人服务时间长度
            D.append(t) #保存客走时间序列
            N.append(n) #保存系统中人数序列
        elif (tA>T) & (tD>T):  break # 过点关门
    while True:
        if n <=0:  break
        t = tD;n=n-1;nD=nD+1 #还有最后几个人在服务
        D.append(t) 
        N.append(n) 
        if n>0:
            tS = rm.uniform(5,10,1)
            tD = t +  tS
            S.append(tS) 
    Tp = max(t-T,0) #关门时间
    # A表示客来时点，D表示客走时点，N表示客走时系统中还有几人，S表示此人服务时长
    raw = {'A':A,'D':D,'S':S,'N':N}
    data = pd.DataFrame(raw)
    return {'count': data.N.count(),'wcount':sum(data.N>0),'avgwait':float(mean(data.D-data.A-data.S))}

# 模拟100次，将结果存为dataframe
res = [simulation() for i in range(100)]
res = pd.DataFrame(res)
res

# 画出平均等待时间的直方图
import matplotlib.pyplot  as plt
plt.hist(res.avgwait)