问题描述

我有一个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的和，我也得不到想要的结果。我想知道我的代码哪里错了。这是我用代码得到的图。

推荐答案

首先，您当前将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'])

除了一些标准化问题外，这开始看起来与您正在寻找的曲线大致相似。当然，平滑曲线对于估计绘图的形状是很好的，但这可能不是实际可视化模拟的最佳方法。您可能会注意到一个特别的问题，即曲线左端的值因求平均值而严重扭曲。您可以通过更改解释窗口的方式或使用不同的卷积内核来略微缓解这一问题，但总会有一些某些边缘效果。

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

优化#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就地集成这些步骤。

使用掩码可以很容易地找到人行道首先穿过-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)

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