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

Question

我有一个 2D 随机游走，其中粒子向左、向右、向上、向下移动或保持在相同位置的概率相等。我生成一个从 1 到 5 的随机数来决定粒子移动的方向。粒子将执行 n 步骤，我重复多次模拟。
我想绘制第一次击中位于 F(t) 处的线性障碍的概率 x = -10 (粒子将在击中该点后消失)。我开始计算每次击中陷阱的粒子 fp 的数量，每次在 1 位置有粒子时添加值 x = -10 。在此之后，我绘制了 fp ，第一次击中陷阱的粒子数，vs 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)

这是有效的，因为 bool 值已经是 int 的子类，并且您希望在每一步都存储 1 或 0。它相当于从 fp[i, j - 1] 语句的两个分支中的赋值的 RHS 中删除 if。
你得到的情节是

这看起来很奇怪，但希望你能看到你想要的情节的一丝曙光。奇怪的原因是穿过陷阱的粒子密度低。您可以通过增加粒子密度或平滑曲线来修复外观，例如与运行平均值。
首先，让我们尝试使用 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'])

除了一些规范化问题外，这开始看起来与您正在寻找的曲线大致相似。当然，平滑曲线有利于估计绘图的形状，但它可能不是实际可视化模拟的最佳方法。您可能会注意到的一个特殊问题是，曲线左端的值因求平均值而严重失真。您可以通过更改窗口的解释方式或使用不同的卷积核来稍微缓解这种情况，但总会有一些边缘效应。
为了真正提高结果的质量，您需要增加样本数量。在这样做之前，我建议先优化您的代码。
如评论中所述，优化 #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)

此代码段首先使用巧妙的楼层划分技巧将 random_walk 计算为一系列步骤，以确保每个步骤的比率恰好为 1/5。然后使用 cumsum 就地集成这些步骤。
使用 mask 很容易找到步行第一次穿过 -10 的地方:

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

argmax 返回最大值的第一次出现。数组 (random_walk == -10) 由 bool 值组成，因此它将返回每列中第一次出现 -10 的索引。在 -10 步骤中从未穿过 simulation_count 的粒子将在其列中包含所有 False 值，因此 argmax 将返回 0 。由于 0 从来都不是有效的步数，因此很容易过滤掉。
步数的直方图将为您提供您想要的确切信息。对于整数数据， np.bincount 是计算直方图的最快方法:

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

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

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