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

Question

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

我想绘制第一次撞到位于x = -10的线性障碍的概率F(t)（粒子撞到这个点后会消失）。我开始计算每次撞击陷阱的模拟的粒子数 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'])

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

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

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

优化 #2 纯粹是为了速度。无需运行多个 for 循环，您可以以纯矢量化的方式完成所有操作。我将展示 new RNG API as well, since I've generally found it to be much faster than the 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.

就地集成这些步骤

walk第一个穿过-10的地方用masking很容易找到：

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

argmax returns 第一次出现最大值。数组 (random_walk == -10) 由布尔值组成，因此它将 return 为每列中第一次出现的 -10 的索引。在 simulation_count 步内从未穿过 -10 的粒子将在其列中包含所有 False 值，因此 argmax 将 return 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 秒。我怀疑使用您拥有的循环实现会花费很多的时间。

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

How to simulate first passage time probability in python for a random walk?

python

random

simulation