Python 中带边界的向量化随机游走

我正在尝试在 python 中模拟带有边界的二维随机游走（粒子/对象将无法跨越边界并且必须返回）。然而我的版本没有矢量化并且非常慢。如何在不使用（或最小化）循环的情况下实现它。

这是我的方法

def bound_walk():

    Origin = [0, 0] #Starting Point

    #All Possible directiom
    directions = ((-1, 1), (0, 1), (1, 1),
                 (-1, 0)        , (1, 0),
                 (-1, -1), (0, -1), (1, -1))


    #Directions allowed when x-coordinate reaches boundary
    refelectionsx = ((-1, 1), (0, 1),
                 (-1, 0),(-1, -1), (0, -1))

    #Directions allowed when y-coordinate reaches boundary
    refelectionsy = ((-1, 0)        , (1, 0),
                 (-1, -1), (0, -1), (1, -1))

    points = [(0, 0)]
    for i in range(20000):
        direction = choice(directions)
        reflection1 = choice(refelectionsx)
        reflection2 = choice(refelectionsy)
        if Origin[0]>50: #Boundary==50
            Origin[0] += reflection1[0]
        elif Origin[0]<-50:
            Origin[0] -= reflection1[0]
        else:
            Origin[0] += direction[0]

        if Origin[1]>50:
            Origin[1] += reflection2[1]
        elif Origin[1] < -50:
            Origin[1] -= reflection2[1]
        else:
            Origin[1] += direction[1]
        points.append(Origin[:])
    return points

这是一种快速的方法，但不能 100% 等同于您的实现。不同之处在于，在我的实现中，在边界处，沿着与边界平行的方向之一前进的机会是从边界后退的方向的一半。如果您将方向视为对连续方向进行分箱的结果，那么这可以说是更好的模型，因为边界将相关分箱切成两半。

如果您尝试一下，您会发现它可以立即或多或少执行 1000 万步。

诀窍是我们简单地“展开”空间，这样我们就可以模拟无约束的随机游走，这是很便宜的，然后最后我们将它折叠回边界矩形。

# parameters
>>> directions = np.delete(np.indices((3, 3)).reshape(2, -1), 4, axis=1).T - 1
>>> boundaries = np.array([(-50, 50), (-50, 50)])
>>> start = np.array([0, 0])>>> steps = 10**7
>>>
# "simulation" 
>>> size = np.diff(boundaries, axis=1).ravel()
>>> 
>>> trajectory = np.cumsum(directions[np.random.randint(0, 8, (steps,))], axis=0)
>>> trajectory = np.abs((trajectory + start - boundaries[:, 0] + size) % (2 * size) - size) + boundaries[:, 0]
>>> 
# some sanity checks
# boundaries are respected
>>> print(trajectory.min(axis=0))
[-50 -50]
>>> print(trajectory.max(axis=0))
[50 50]
# step size looks ok
>>> print(np.diff(trajectory, axis=0).min(axis=0))
[-1 -1]
>>> print(np.diff(trajectory, axis=0).max(axis=0))
[1 1]
# histograms of time spent at coordinates looks flat
>>> print(np.bincount(trajectory[:, 0] - boundaries[0, 0]))
[ 50276 100134 100395 100969 101218 101388 101708 100688 101460 102667
 103613 103652 103540 103296 102676 102105 102766 102855 101786 101246
 101442 101152 101020 100498 100637 100588 100100  99745 100034  99878
  99120  98076  98193  98126  97715  98317  98343  97693  97391  96854
  96576  96906  96423  96445  96779  96672  96376  95747  95732  95881
  96833  97149  98490  99692  99519  98800  99497 100070 100065  99816
  99838 100470 100466 100887 100461 100033  99405  99425 100537 100227
 100796 101668 101218 101413 101559 101258 101416 101292 100567 100022
 100266 100770 100882 100519 100326 100795 101066 101293 101667 101666
 101040 101221 101019 100868 101681 100778 100121  98500  98174  98308
  49254]
>>> print(np.bincount(trajectory[:, 1] - boundaries[1, 0]))
[ 52316 104725 104235 103801 102936 102269 102604 102557 102514 103063
 102130 101805 101699 102285 102456 102464 102590 104010 103502 103105
 102784 102927 103430 104750 104671 104836 104547 103280 102131 101548
 101173 101806 101345 101959 101525 101061 101260 100774 100126  98806
  99209 100105  99686 100418 101056 101434 101078 101680 103042 103732
 103003 102047 100832 100489 100809 100429 101325 102420 102282 102205
 101341 100644  99827  99482  98931  98588  97911  97981  97053  96794
  96818  97364  97025  97093  97807  98594  98280  98406  98474  98516
  98555  98713  98381  98296  97600  97374  97423  97092  96238  95771
  95547  95325  94710  94115  93332  92219  91309  91780  92399  92345
  45461]