如何用Python模拟二维随机游走的首次通过时间概率并排查代码错误？

问题根源 #

你混淆了累计触达粒子数和**t时刻新增首次触达粒子数（概率密度）**两个统计量：

1. 你当前统计的s是到t步为止总共触达的粒子数，这个值本身就是单调非降的，所有粒子要么触达了算进去，要么没触达不算，当然不可能出现下降段，这是你曲线不符合预期的核心原因。
2. 代码细节bug：你导入numpy时取了别名np，但后续初始化数组时直接写numpy.zeros会直接报错，运行前需要先修正命名一致性问题。
3. 粒子触达屏障break后，数组x中该粒子后续所有步的位置都保留为初始0，你统计fp时会把该粒子的触达计数一直携带到后续所有时间步，进一步强化了累计求和的单调特性。

修正方案 #

你需要统计的是每个时间步新发生首次触达的粒子占比，也就是先记录每个粒子的首次触达时间，再按时间步统计频数即可，修正后代码如下：

import matplotlib.pyplot as plt 
import numpy as np
import random

n = 1000
n_simulations = 1000

# 记录每个模拟的首次触达时间，初始设为n表示未触达
first_hit = np.full(n_simulations, n)

for i in range(n_simulations):
    x_pos = 0
    for j in range(1, n):
        val = random.randint(1, 5)
        if val == 1:
            x_pos += 1
        elif val == 2:
            x_pos -= 1
        # 上下移动和原地不动都不改变x坐标，不用额外处理
        if x_pos == -10:
            first_hit[i] = j
            break

# 统计每个时间步的首次触达频数
fp_count = np.zeros(n)
for t in first_hit:
    if t < n:
        fp_count[t] += 1

# 转换为概率，如果你需要累计触达概率可以替换为np.cumsum(fp_count)/n_simulations
F_t = fp_count / n_simulations
steps = np.arange(n)

plt.xlim(0, 1000)
plt.plot(steps, F_t)
plt.xlabel('时间步t')
plt.ylabel('首次触达概率F(t)')
plt.show()

效果说明 #

运行上述代码后得到的曲线就是你预期的先升后降形态：前10步不可能有粒子触达（初始x=0，最少需要10步连续左移才能到x=-10），所以前10步值为0，随后新增触达粒子数逐渐上升到峰值，再随着剩余未触达粒子减少逐渐下降，最后趋近于0。

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