如何在python中模拟随机游走的首次通过时间概率？

Question

我有一个 2D 随机游走，其中粒子具有相同的概率向左、向右、向上、向下移动或保持在同一位置。我生成一个从 1 到 5 的随机数来决定粒子将向哪个方向移动。粒子会执行n 步骤，我重复模拟几次。

我想绘制F(t) 第一次击中位于x = -10 的线性屏障的概率（粒子在击中该点后会消失）。我开始计算每个模拟命中陷阱的粒子数量fp，每次在x = -10 位置有一个粒子时添加值1。在此之后，我绘制了fp，第一次撞击陷阱的粒子数，与t，时间步长。

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

n = 1000
n_simulations=1000

x = numpy.zeros((n_simulations, n))
y = numpy.zeros((n_simulations, n))
steps = np.arange(0, n, 1)
for i in range (n_simulations):
    for j in range (1, n):
        val=random.randint(1, 5)
        if val == 1:
            x[i, j] = x[i, j - 1] + 1
            y[i, j] = y[i, j - 1] 
        elif val == 2:
            x[i, j] = x[i, j - 1] - 1
            y[i, j] = y[i, j - 1] 
        elif val == 3:
            x[i, j] = x[i, j - 1]
            y[i, j] = y[i, j - 1] + 1
        elif val == 4:
            x[i, j] = x[i, j - 1]
            y[i, j] = y[i, j - 1] - 1
        else:
            x[i, j] = x[i, j - 1]
            y[i, j] = y[i, j - 1]
        if x[i, j] == -10:
            break

fp = np.zeros((n_simulations, n)) # number of paricles that hit the trap for each simulation. 
for i in range(n_simulations):
    for j in range (1, n):
        if x[i, j] == -10:
            fp[i, j] = fp[i, j - 1] + 1
        else:
            fp[i, j] = fp[i, j - 1]
s = [sum(x) for x in zip(*fp)]
plt.xlim(0, 1000)
plt.plot(steps, s)
plt.show()

我应该有以下情节：

但是我得到的图是不同的，因为曲线总是在增加，并且对于大的t 它应该减少（对于大的t，大多数粒子已经击中目标并且概率降低）。即使不使用fp 的总和，我也没有想要的结果。我想知道我的代码哪里错了。这是我用我的代码得到的情节。

Answer 1

首先，您当前将fp 计算为穿过陷阱的所有粒子的累积总和。这个数字必然与n 渐近。您正在寻找的是累积和的导数，即每单位时间穿过陷阱的粒子数。

在第二个循环中需要进行非常简单的更改。更改以下条件

if x[i, j] == -10:
    fp[i, j] = fp[i, j - 1] + 1
else:
    fp[i, j] = fp[i, j - 1]

到

fp[i, j] = int(x[i, j] == -10)

这是可行的，因为布尔值已经是 int 的子类，并且您希望在每一步都存储 1 或 0。这相当于从if 语句的两个分支中的赋值的RHS 中删除fp[i, j - 1]。

你得到的情节是

这看起来很奇怪，但希望你已经可以看到你想要的情节的微光。奇怪的原因是穿过陷阱的粒子密度低。您可以通过增加粒子密度或平滑曲线来修复外观，例如具有运行平均值。

首先，让我们尝试使用np.convolve 的平滑方法：

x1 = np.convolve(fp.sum(0), np.full(11, 1/11), 'same')
x2 = np.convolve(fp.sum(1), np.full(101, 1/101), 'same')

plt.plot(s, x1)
plt.plot(s, x2)
plt.legend(['Raw', 'Window Size 11', 'Window Size 101'])

这开始看起来与您正在寻找的曲线大致相似，除非存在一些标准化问题。当然，平滑曲线有利于估计绘图的形状，但它可能不是实际可视化模拟的最佳方法。您可能注意到的一个特殊问题是曲线左端的值因平均而变得高度失真。您可以通过更改窗口的解释方式或使用不同的卷积核来稍微缓解这种情况，但总会有一些边缘效应。

要真正提高结果质量，您需要增加样本数量。在这样做之前，我建议先优化一下你的代码。

如 cmets 中所述，优化 #1 是您不需要为这个特定问题同时生成 x 和 y 坐标，因为陷阱的形状允许您解耦这两个方向。相反，你有 1/5 的概率踩到 -x 和 1/5 的概率踩到 +x。

优化#2 纯粹是为了速度。无需运行多个 for 循环，您可以以纯矢量化的方式执行所有操作。我还将展示new RNG API 的示例，因为我通常发现它比legacy API 快得多。

优化#3 是为了提高可读性。如果没有完整的文档，像 n_simulations、n 和 fp 这样的名称信息量不大。我将在下面的示例中重命名一些东西以使代码自记录：

particle_count = 1000000
step_count = 1000

# -1 always floor divides to -1, +3 floor divides to +1, the rest zero
random_walk = np.random.default_rng().integers(-1, 3, endpoint=True, size=(step_count, particle_count), dtype=np.int16)
random_walk //= 3  # Do the division in-place for efficiency
random_walk.cumsum(axis=0, out=random_walk)

这个 sn-p 计算 random_walk 作为一系列步骤，首先使用聪明的地板除法技巧来确保每个步骤的比率正好是 1/5。然后使用cumsum 就地集成这些步骤。

使用遮罩很容易找到步行首先穿过 -10 的位置：

steps = (random_walk == -10).argmax(axis=0)

argmax 返回最大值的第一次出现。数组(random_walk == -10) 由布尔值组成，因此它将返回每列中第一次出现-10 的索引。在simulation_count 步骤中从不跨越-10 的粒子将在其列中包含所有False 值，因此argmax 将返回0。由于0 永远不是有效的步数，因此很容易过滤掉。

步数的直方图将为您提供所需的准确信息。对于整数数据，np.bincount 是计算直方图的最快方法：

histogram = np.bincount(steps)
plt.plot(np.arange(2, histogram.size + 1), hist[1:] / particle_count)

histogram 的第一个元素是在step_count 步内从未到达-10 的粒子数。 histogram 的前 9 个元素应该始终为零，除非 argmax 是如何工作的。因为histogram[0]名义上代表一步后的计数，所以显示范围移动了一位。

在我的中等功率机器上，生成 10 亿个样本并将它们求和需要不到 30 秒。我怀疑使用您拥有的循环实现会花费更多更长的时间。