用于评估随机游走的python库？

Question

我正在尝试评估随机游走结束位置的概率，但我的程序速度有些问题。基本上我要做的是将包含随机游走概率的字典作为输入（例如 p = {0:0.5, 1:0.2. -1:0.3} 意味着 X 有 50% 的概率停留在0，20%概率X增加1，30%概率X减少1)，然后计算n次迭代后所有可能的未来状态的概率。

例如，如果 p = {0:0.5, 1:0.2。 -1:0.3} 和 n = 2 那么它将返回 {0:0.37, 1:0.2, -1:0.3, 2:0.04, -2:0.09} 如果 p = {0:0.5, 1:0.2。 -1:0.3} 和 n = 1 那么它将返回 {0:0.5, 1:0.2。 -1:0.3}

我有工作代码，如果 n 很低并且 p 字典很小，它运行相对较快，但是当 n > 500 并且字典有大约 50 个值时，计算需要超过 5 分钟。我猜这是因为它只在一个处理器上执行，所以我继续修改它，以便它使用 python 的多处理模块（因为我读到多线程不会因为 GIL 而提高并行计算性能）。

我的问题是，多处理并没有太大的改进，现在我不确定是因为我实现错误还是因为 python 中的多处理开销。我只是想知道当 n > 500 并行时是否有一个库可以评估随机游走的所有可能性的所有概率？如果我找不到任何东西，我的下一步是用 C 编写我自己的函数作为扩展，但这将是我第一次这样做，尽管我已经用 C 编码了一段时间。

原始的非多处理代码

def random_walk_predictor(probabilities_tree, period):
    ret = probabilities_tree
    probabilities_leaves = ret.copy()
    for x in range(period):
        tmp = {}
        for leaf in ret.keys():
            for tree_leaf in probabilities_leaves.keys():
                try:
                    tmp[leaf + tree_leaf] = (ret[leaf] * probabilities_leaves[tree_leaf]) + tmp[leaf + tree_leaf]
                except:
                    tmp[leaf + tree_leaf] = ret[leaf] * probabilities_leaves[tree_leaf]
        ret = tmp
    return ret

多处理代码

from multiprocessing import Manager,Pool
from functools import partial

def probability_calculator(origin, probability, outp, reference):
    for leaf in probability.keys():
        try:
            outp[origin + leaf] = outp[origin + leaf] + (reference[origin] * probability[leaf])
        except KeyError:
            outp[origin + leaf] = reference[origin] * probability[leaf]

def random_walk_predictor(probabilities_leaves, period):
    probabilities_leaves = tree_developer(probabilities_leaves)
    manager = Manager()
    prob_leaves = manager.dict(probabilities_leaves)
    ret = manager.dict({0:1})
    p = Pool()

    for x in range(period):
        out = manager.dict()
        partial_probability_calculator = partial(probability_calculator, probability = prob_leaves, outp = out, reference = ret.copy())

        p.map(partial_probability_calculator, ret.keys())
        ret = out

    return ret.copy()

Answer 1

往往有分析解决方案可以精确解决这种看起来类似于二项分布的problem，但我假设您确实在为更一般的问题类别寻求计算解决方案。

与使用 Python 字典相比，从基本数学问题的角度来考虑这一点更容易。构建一个矩阵A，描述从一种状态到另一种状态的概率。建立一个状态x，描述某个时间在给定位置的概率。

因为在n 转换之后，您最多可以从原点（在任一方向）步进n - 您的状态需要有 2n+1 行，A 需要是正方形且大小为 2n+ 1 乘 2n+1。

对于两个时间步长的问题，您的转换矩阵将是 5x5，如下所示：

[[ 0.5  0.2  0.   0.   0. ]
 [ 0.3  0.5  0.2  0.   0. ]
 [ 0.   0.3  0.5  0.2  0. ]
 [ 0.   0.   0.3  0.5  0.2]
 [ 0.   0.   0.   0.3  0.5]]

你在时间 0 的状态将是：

[[ 0.]
 [ 0.]
 [ 1.]
 [ 0.]
 [ 0.]]

系统的一步演化可以通过将A 和x 相乘来预测。

所以在 t = 1 时，

 x.T = [[ 0.   0.2  0.5  0.3  0. ]]

在 t = 2 时，

x.T = [[ 0.04  0.2   0.37  0.3   0.09]]

因为即使是适度数量的时间步长，这也可能会占用相当多的存储空间（A 需要 n^2 存储空间），但是非常稀疏，我们可以使用稀疏矩阵来减少存储空间（并加快我们的计算）。这样做意味着A 需要大约 3n 个元素。

import scipy.sparse as sp
import numpy as np

def random_walk_transition_probability(n, left = 0.3, centre = 0.5, right = 0.2):
    m = 2*n+1
    A  = sp.csr_matrix((m, m))
    A += sp.diags(centre*np.ones(m), 0)
    A += sp.diags(left*np.ones(m-1), -1)
    A += sp.diags(right*np.ones(m-1),  1)
    x = np.zeros((m,1))
    x[n] = 1.0
    for i in xrange(n):
        x = A.dot(x)
    return x

print random_walk_transition_probability(4)

时间

%timeit random_walk_transition_probability(500)
100 loops, best of 3: 7.12 ms per loop

%timeit random_walk_transition_probability(10000)
1 loops, best of 3: 1.06 s per loop