Python Itertools蒙特卡洛模拟结果异常问题求助
slothrop评论后编辑：
看起来takewhile会消耗迭代器中最后一个被检查的元素，这意味着无法判断下一个元素是否也满足条件，导致计算出现问题。在无限随机实现的迭代器上使用takewhile似乎不是好主意（如果需要保留这些实现的话）。
问题
我编写了两次蒙特卡洛模拟，一次用itertools实现，一次用传统方式。前者出现异常行为，结果不正确，而后者运行正常，想请教原因，抱歉附上了较多代码。
猜想
我猜测takewhile函数结合随机变量（RVs）时存在异常行为。注意simulate_freqs_iter函数中有一些被注释的写法（如zip(takewhile(lambda x: x < t, honest_arrival_times), attack_arrival_times)），它们的结果差异很大，我无法解释。
说明：我知道n=1e5太小，原本计划用n=1e8或1e9，但不会让您运行20小时。n=1e5只需约10秒，且结果差异已很明显！
非常感谢您的解答，我非常喜欢itertools（即便这个实现既不快也无其他效率优势）。
代码
from random import expovariate, seed from itertools import accumulate, takewhile, count, islice, chain seed("31.05.2023") def simulate_freqs_iter(k:int, mu:float, lam:float, n:int, conditional_on_finite:bool=True, verbose:bool=False) -> dict: """Simulates n monte-carlo runs of a slow poisson process catching up with a fast one returns a dictionary containing the empirical cdf for t={10, 11, ..., 50}""" abs_freqs = {t: 0 for t in range(10, 51)} for i in range(n): # Waiting times for each poisson process honest_arrival_times = accumulate(((expovariate(mu)) for _ in count(start=0, step=1))) attack_arrival_times = islice( # attack_arrival_times is k blocks behind, but only needs to catch up, not surpass accumulate(((expovariate((lam))) for _ in count(start=0, step=1))), k-1, None) # Optional printing function if verbose and (i%50000)==0: first_10_honest = list(islice(honest_arrival_times, 10)) first_10_attack = list(islice(attack_arrival_times, 10)) print(list(zip(first_10_honest, first_10_attack))) honest_arrival_times = chain(first_10_honest, honest_arrival_times) attack_arrival_times = chain(first_10_attack, attack_arrival_times) # Simulates poisson waiting times and counts number of hits where 1 catches up with 2 for t in abs_freqs.keys(): # Every time N_2 surpasses N_1 we have a hit hits = list((t, honest, attack) for honest, attack in # Take all jump times of N_1 and see if N_2 was there before N_1 within t # takewhile(lambda x: min(x[0], x[1]) < t, zip(honest_arrival_times, attack_arrival_times)) # zip(takewhile(lambda x: x < t, honest_arrival_times), attack_arrival_times) zip(honest_arrival_times, takewhile(lambda x: x < t, attack_arrival_times)) # zip(takewhile(lambda x: x < t, honest_arrival_times), takewhile(lambda x: x < t, attack_arrival_times)) if honest > attack) if hits: if verbose: print(hits[0]) abs_freqs[t]+=1 # Don't continue this path once found a hit break # Condition on the probability that it ever catches up prob_finite = min(((lam/mu)**k), 1) if conditional_on_finite else 1 # Accumulate relative frequencies to get cdf results = {key:((value/n) / prob_finite) for key, value in dict( accumulate(abs_freqs.items(), lambda prev,curr: (curr[0], (curr[1]+prev[1])))).items()} return results def simulate_freqs_orthodox(k:int, mu:float, lam:float, n:int, conditional_on_finite:bool=True, verbose:bool=False) -> dict: """Simulates n monte-carlo runs of a slow poisson process catching up with a fast one returns a dictionary containing the empirical cdf for t={10, 11, ..., 50}""" abs_freqs = {t: 0 for t in range(10, 51)} # We run the monte carlo n times for _ in range(n): # Initialise first block arrival next_honest_arrival = expovariate(mu) next_attack_arrival = expovariate(lam) # the attacker chain is behind k blocks for _ in range(k-1): next_attack_arrival += expovariate(lam) for t in abs_freqs.keys(): while next_attack_arrival < t: # count when a double spend attack is successful if next_honest_arrival > next_attack_arrival: abs_freqs[t] += 1 if verbose: print(next_honest_arrival, next_attack_arrival) # stop searching for all t on this path once found a hit break # Add another block to each chain next_honest_arrival += expovariate(mu) next_attack_arrival += expovariate(lam) else: # If while-loop terminates without a hit, continue with greater t continue # Terminate search for this path after a successful hit break # Condition on the probability that it ever catches up prob_finite = min(((lam/mu)**k), 1) if conditional_on_finite else 1 # Accumulate relative frequencies to get cdf results = {key:((value/n) / prob_finite) for key, value in dict( accumulate(abs_freqs.items(), lambda prev,curr: (curr[0], (curr[1]+prev[1])))).items()} return results # Monte Carlo simulation of probabilities for T_k k=10 mu=7/10 lam=3/10 n=int(1e5) print(simulate_freqs_iter(k=k, mu=mu, lam=lam, n=n, conditional_on_finite=True, verbose=False)) print(simulate_freqs_orthodox(k=k, mu=mu, lam=lam, n=n, conditional_on_finite=True, verbose=False))
解答
核心问题出在takewhile对迭代器的消耗方式上，和随机变量本身无关，是itertools版本的逻辑和正统版本完全不一致：
迭代器的一次性消耗特性：
honest_arrival_times和attack_arrival_times都是无限迭代器，每次取元素都会消耗掉该元素，无法复用。你在循环每个t时，每次用takewhile截取迭代器，会导致前一个t处理时已经消耗了部分元素，后一个t只能拿到剩余元素，完全违背了正统版本中“每个t都从初始状态开始模拟”的逻辑。takewhile的截断逻辑：以你当前使用的
zip(honest_arrival_times, takewhile(lambda x: x < t, attack_arrival_times))为例，takewhile会一直取attack_arrival_times的元素直到第一个不满足x < t的元素，并且这个不满足条件的元素也会被消耗掉，不会放回迭代器。下一次循环更大的t时，attack_arrival_times已经跳过了这个元素，模拟状态直接断片，和正统版本逐步累积到t的逻辑完全不符。版本逻辑差异：正统版本中，每个蒙特卡洛样本从初始状态开始，对每个
t逐步模拟到该时间点；而itertools版本中，每个样本的迭代器是共享的，处理完t=10后，迭代器已经走到t=10之后的状态，处理t=11时是从中间状态继续，而非重新从初始状态模拟到t=11，这直接导致统计结果错误。
修复建议
如果想保留itertools的写法，需要为每个t重新生成独立的迭代器（但无限迭代器无法保存状态）。更合理的做法是，对每个蒙特卡洛样本，先生成到足够大的t的所有到达时间，再针对每个t统计：
# 在simulate_freqs_iter的每个样本循环内 honest_times = [] attack_times = [] # 先生成到t=50之后的足够多的点 current_honest = 0.0 current_attack = 0.0 for _ in range(k-1): current_attack += expovariate(lam) while current_honest < 50 or current_attack < 50: current_honest += expovariate(mu) honest_times.append(current_honest) current_attack += expovariate(lam) attack_times.append(current_attack) # 然后对每个t统计 for t in abs_freqs.keys(): hits = [(t, h, a) for h, a in zip(honest_times, attack_times) if a < t and h > a] if hits: abs_freqs[t] +=1 break
这样既保留了迭代器的思想，又保证每个t的统计基于完整的初始状态模拟。
内容的提问来源于stack exchange，提问作者WonderfulWonder
