Scipy与Matplotlib绘制CCDF图的代码排查及正确性验证

我尝试使用Scipy和Matplotlib绘制CCDF图，参考相关教程实现的代码能生成图表，但不确定结果是否正确，请帮忙检查代码中的问题。

完整代码 #

# compare the completion times and service times and check that they are different
# do this for multiple runs 
# 

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

class Job:
    def __init__(self, arrivalTime, size):
        self.arrivalTime = arrivalTime
        self.size = size
        self.departureTime = None

departures = 0
maxDepartures = 100000

# Function that generates the job size
def generateJobSize():
    return np.random.exponential(1 / 3)

# Function that generates an interarrival time
def generateInterarrivalTime():
    return np.random.exponential(1 / 2)

jobQueue = [] 
serverEmpty = True
nextArrTime = generateInterarrivalTime()
nextDepTime = float('inf')
completionTimes = []
# initiate pointer for the job currently being serviced
servingJob = None
clock = 0.0
loop = 0

def handleArr():
  # declare these variables as global so that they can be edited outside the function
  global job, serverEmpty, servingJob, nextArrTime, nextDepTime, completionTimes
  # generate the size of the job
  size = generateJobSize()
  # create job object with those attriubtes
  job = Job(arrivalTime = clock, size = size)
  # if the server is empty upon this job's arrival into the system, immediately start servicing it 
  if serverEmpty == True:
    # create pointer that points to the job that is currently being serviced
    servingJob = job
    # determine the job's departure time based on the size and the current time
    job.departureTime = clock + job.size
    nextDepTime = job.departureTime
    # update the status of the server to busy
    serverEmpty = False
  else:
    # if a job is already being serviced, append the new arrival to the queue
    jobQueue.append(job)
  # generate the arrival time of the next job
  interArrTime = generateInterarrivalTime ()
  # set the next arrival time based upon that
  nextArrTime = clock + interArrTime
 
 

# fcn to handle departures
def handleDep():
  global job, departures, serverEmpty, servingJob, nextDepTime, completionTimes
  # determine the job's departure time based on the size and the current time
  job.departureTime = clock + job.size
  completionTimes.append(job.departureTime - job.arrivalTime)
  # update the time of the next departure based on this new information
  nextDepTime = clock + job.size
# when all the jobs have finished
# inc the counter by 1
  departures += 1
  # if the queue is not empty, pop the next job in line
  if len(jobQueue) != 0:
    job = jobQueue.pop(0)
    servingJob = job
  else: 
    serverEmpty = True
    nextDepTime = float("inf")
    servingJob = None

while departures < maxDepartures: 
  # determine the next event
  if nextArrTime <= nextDepTime:
    clock = nextArrTime
    handleArr()
  else: 
    loop += 1
    clock = nextDepTime
    handleDep()

print(len(completionTimes))

# --------- GRAPHS A CCDF PLOT -----------------

# # Sort completion times in ascending order
# completionTimes.sort()

# # Calculate probabilities for completion times
# probabilities_completion = np.arange(1, len(completionTimes) + 1) / len(completionTimes)

# # Calculate the CCDF for completion times
# ccdf_completion = 1 - probabilities_completion

# # Interpolate the data for a smooth line
# interp_func_completion = interp1d(completionTimes, ccdf_completion, kind='linear')

# # Generate interpolated data points for completion times
# completion_times_interp = np.linspace(min(completionTimes), max(completionTimes), num=1000)
# ccdf_completion_interp = interp_func_completion(completion_times_interp)

# # Plot completion times vs CCDF
# plt.plot(completion_times_interp, ccdf_completion_interp)
# plt.xlabel('Completion Time')
# plt.ylabel('CCDF')
# plt.title('Completion Time vs CCDF')
# plt.grid(True)
# plt.show()

x = np.sort(completionTimes)

#calculate CDF values
y = 1. * np.arange(len(completionTimes)) / (len(completionTimes) - 1)

#plot CDF
plt.plot(x, y)

代码问题分析 #

1. 事件处理逻辑错误（核心问题） #

handleDep函数存在两处关键错误，直接导致completionTimes数据失真：

• 错误覆盖departureTime：函数开头的job.departureTime = clock + job.size完全多余，且会破坏正确值。因为clock此时已经等于nextDepTime，也就是当前服务job的预设离开时间，重新计算会覆盖之前在handleArr中设置的正确departureTime，导致completionTime（离开时间-到达时间）计算错误。
• 未设置新服务job的离开时间：当队列中有job时，弹出后仅赋值servingJob = job，但未计算该job的departureTime和更新nextDepTime，后续事件调度会完全混乱，生成的completionTimes数据毫无意义。

修正后的handleDep函数：

def handleDep():
  global departures, serverEmpty, servingJob, nextDepTime, completionTimes
  # 使用预先设置好的departureTime计算完成时间
  completionTimes.append(servingJob.departureTime - servingJob.arrivalTime)
  departures += 1
  
  # 处理下一个job
  if len(jobQueue) != 0:
    servingJob = jobQueue.pop(0)
    servingJob.departureTime = clock + servingJob.size
    nextDepTime = servingJob.departureTime
  else: 
    serverEmpty = True
    nextDepTime = float("inf")
    servingJob = None

注意：这里移除了不必要的global job，直接使用servingJob来操作当前服务的job，避免变量混淆。

2. CCDF/CDF绘图的正确性问题 #

• 被注释的CCDF代码：逻辑基本正确，但需要确保completionTimes数据正确（先修复上述事件处理错误）。排序后，probabilities_completion = np.arange(1, len(completionTimes)+1)/len(completionTimes)是标准的经验CDF计算（表示P(X ≤ x)），因此ccdf_completion = 1 - probabilities_completion即为CCDF（P(X > x)），插值绘图的逻辑没问题。
• 未注释的CDF代码：y = 1. * np.arange(len(completionTimes)) / (len(completionTimes)-1)存在问题，当样本量maxDepartures=1e5时，数值差异不大，但这不是标准的经验CDF计算方式。标准计算应该是y = np.arange(1, len(completionTimes)+1)/len(completionTimes)，这样第k个排序后的值对应概率k/n，符合经验CDF的定义。

修正后绘图示例 #

修复事件处理逻辑后，使用以下代码绘制CCDF即可得到正确结果：

# 排序完成时间
completionTimes.sort()

# 计算经验CDF
cdf = np.arange(1, len(completionTimes)+1) / len(completionTimes)
# 计算CCDF
ccdf = 1 - cdf

# 插值平滑（可选，大数据量下也可以直接绘制散点）
interp_func = interp1d(completionTimes, ccdf, kind='linear')
x_interp = np.linspace(min(completionTimes), max(completionTimes), 1000)
y_interp = interp_func(x_interp)

plt.plot(x_interp, y_interp)
plt.xlabel('完成时间')
plt.ylabel('CCDF')
plt.title('完成时间CCDF分布图')
plt.grid(True)
plt.show()

内容的提问来源于stack exchange，提问作者Sabrina Hatch