计算随机游动的python库？

Question

我试图评估随机游动的终点位置的概率，但我的程序速度有点问题。基本上，我要做的是把包含随机游走概率的字典作为输入（例如p={0:0.5，1:0.2）。-1： 0.3}意味着50%的概率X保持在0，20%的概率X增加1，30%的概率X降低1），然后计算n次迭代后所有可能的未来状态的概率。在

例如p={0:0.5，1:0.2。-1： 然后返回{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

多处理代码

^{pr2}$

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的状态将是：

^{pr2}$

用A和{}相乘，可以预测系统的一步演化。在

所以当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