Poisson分布光子抽样Python代码重复波束结构异常排查

你的代码出现光子数量异常、耗时剧增的核心原因是抽样逻辑错误，同时还有循环范围和索引处理的问题，下面逐一拆解：

1. 光子数量爆炸的核心原因 #

你在循环里写了这行代码：

spkt = np.random.rand(i+1)<fr*dt

这意味着在遍历第i个bin时，你生成了i+1个随机数（比如遍历到第935个bin时，一次性生成936个随机数），然后把所有满足条件的结果都算作该bin的光子。这完全违背了“每个bin独立抽样一次”的需求——原本每个有效bin只需要做1次伯努利抽样（因为fr*dt极小，Poisson分布可以用伯努利近似），你却随着循环次数增加，抽样次数越来越多，自然会产生远超预期的光子。

另外，当你替换为y1/y2/y3时，你并没有修改循环范围（依然只遍历原始y的936个bin），但你的意图是遍历重复后的所有bin，这也会导致逻辑偏差。

2. 耗时剧增的原因 #

随着循环到后面的bin，i+1的数值越来越大（比如模拟10000次重复时，最后一个bin的i是9359999，要生成近1000万个随机数），这种指数级增长的计算量直接导致耗时飙升。

修正后的代码 #

下面是调整后的代码，完全贴合你的需求：

import numpy as np
from matplotlib import pyplot as plt

# 定义单个波束结构
num_bins = 936
charge_per_bin = 62e-9  # 0.62nC
dt = 2e-9  # 每个bin时长2ns
cycle_time = num_bins * dt  # 单圈周期1.872us

# 构建波束结构：前900个bin有电荷，后36个为0
y = np.zeros(num_bins)
y[:900] = charge_per_bin

# 设置模拟参数
repeat_times = 10000  # 重复10000次波束结构
fr = 10000  # 输入计数率
lambda_p = fr * dt  # 每个bin的Poisson均值（伯努利概率）

# 生成重复后的波束结构
y_repeated = np.tile(y, repeat_times)
total_time = len(y_repeated) * dt

# 抽样生成光子
spks_t = []
for idx, charge in enumerate(y_repeated):
    if charge != 0:
        # 每个有效bin做1次伯努利抽样
        if np.random.rand() < lambda_p:
            # 记录光子时间：这里用bin的起始时间，也可以改成中间时间或随机时间
            photon_time = idx * dt
            spks_t.append(photon_time)

# 剔除堆积光子（时间间隔小于80ns的）
corrected_times = []
prev_time = -np.inf
for t in spks_t:
    if t - prev_time >= 80e-9:
        corrected_times.append(t)
        prev_time = t

# 计算结果
corrected_photons = len(corrected_times)
firing_rate = corrected_photons / total_time

print(f"有效光子数：{corrected_photons}")
print(f"输出计数率：{firing_rate:.2f} Hz")

关键修正点说明 #

• 抽样逻辑：每个有效bin仅做1次随机抽样，符合Poisson分布（小λ下的伯努利近似）的要求，不会产生多余的光子。
• 循环遍历：用enumerate遍历重复后的波束结构所有bin，确保每个bin都被处理。
• 堆积光子处理：通过遍历光子时间列表，跳过与前一个光子间隔小于80ns的样本，避免索引越界问题。
• 性能优化：避免了随循环次数增长的随机数生成，耗时与模拟的bin总数线性相关，10000次重复的模拟也能高效运行。

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