在Python中加速Beta Pert分布的计算

Question

我正在计算每次循环迭代的beta PERT分布（还有其他一些事情，但计算这个分布是最耗时的）。

最开始我是在R语言中写的这个代码，但花了太长时间，所以我尝试用一个更快的工具。

我的一些数据集可能会很大，比如我刚刚处理了一个有153413个案例的情况，结果在Python中运行了大约8个小时（比R快，但还是有点慢）。

我对Python还很陌生，想知道有没有什么方法可以加快这样的计算？

示例代码：

af = lambda pmu, pmin, pmode, pmax: (pmu-pmin)*(2*pmode-pmin-pmax)/((pmode-pmu)*(pmax-pmin))
bf = lambda pmu, pmin, pmode, pmax: (pmax-pmu)/(pmu-pmin)*((pmu-pmin)*(2*pmode-pmin-pmax)/((pmode-pmu)*(pmax-pmin)))

e=5.
shape=4.
max=10.
mu_d = np.arange(0, 10, 0.05)                
d = np.arange(0.025, 60.025, 0.05)
nlocs=153413  # number of rows in dataset


f0_dist = np.zeros(len(mu_d))
f1_dist = np.zeros(len(mu_d))
f2_dist = np.zeros(len(mu_d))

f0 = st.norm.cdf(d, 0.9/2., 0.9/6.)
f1 = st.uniform.cdf(d, 0.001, 0.9)

tic = time.clock()     
    for i in xrange(nlocs):
       for j in xrange(len(mu_d)): # mu_d has 121 values
            Rp_min = mu_d[j] - 1.96*e
            Rp_mode = mu_d[j] - 0.75*e
            Rp_max = max
            Rp_mu=(Rp_min+Rp_max+shape*Rp_mode)/(shape+2)
   dist = st.beta.cdf(d, a=af(Rp_mu, Rp_min, Rp_mode, Rp_max), b=bf(Rp_mu, Rp_min, Rp_mode, Rp_max), loc=Rp_min, scale=1-Rp_min)

    f0_dist[j] = 1 - np.sum(dist*f0*0.05)
    f1_dist[j] = 1- np.sum(dist*f1*0.05)
    f2_dist[j] = 1 - np.sum(dist*0.05)
    temp = 0.4*f0_dist + 0.5*f1_dist + 0.1*f1_dist
    aggr_dist = aggr_dist + temp

toc = time.clock() - tic
print '\nTime elapsed: %.3f seconds\n' % toc

Answer 1

这里有一段稍微修改过的代码：

af = lambda pmu, pmin, pmode, pmax: (pmu-pmin)*(2*pmode-pmin-pmax)/((pmode-pmu)*(pmax-pmin))
bf = lambda pmu, pmin, pmode, pmax: (pmax-pmu)/(pmu-pmin)*((pmu-pmin)*(2*pmode-pmin-pmax)/((pmode-pmu)*(pmax-pmin)))

e=5.
shape=4.
max=10.
mu_d = np.arange(0, 10, 0.05)                
d = np.arange(0.025, 60.025, 0.05)

Rp_max = max
e1_96 = 1.96 * e
e0_75 = 0.75 * e
for i in xrange(nlocs): # e.g 153413
   for mu_d_j in mu_d: # mu_d has 121 values
        Rp_min = mu_d_j - e1_96
        Rp_mode = mu_d_j - e0_75
        Rp_mu=(Rp_min+Rp_max+shape*Rp_mode)/(shape+2)

   dist = st.beta.cdf(d, a=af(Rp_mu, Rp_min, Rp_mode, Rp_max), b=bf(Rp_mu, Rp_min, Rp_mode, Rp_max), loc=Rp_min, scale=1-Rp_min)

接下来是解释：

在循环内保存每个指令

• 把 Rp_max = max 移到循环外面
• 把常量（比如 e1_96 和 e0_75）的计算放到循环外面

避免深层引用

• 只执行一次 mu_d[j]，然后用一个局部变量来保存这个值，深层次的取值会浪费时间

用 for 循环来获取值，而不是用 lst[i]

下面的内容：

for j in xrange(len(mu_d)): # mu_d has 121 values
    mu_d_j = mu_d[j]

应该改成更高效（而且符合Python风格）的：

for mu_d_j in mu_d: # mu_d has 121 values
    #now use mu_d_j

测量时间

这是基本原则，每次修改都要评估一下效果。如果你设定了期望的速度（处理时间），那么你可能会很快停止优化。

免责声明

因为我无法运行这段代码，所以不能保证所有的修改都是正确的。有几行代码我不太确定它们的作用：

最后一行 dist =

   dist = st.beta.cdf(d, a=af(Rp_mu, Rp_min, Rp_mode, Rp_max), b=bf(Rp_mu, Rp_min, Rp_mode, Rp_max), loc=Rp_min, scale=1-Rp_min)

这个缩进正确吗？现在的情况是，它在每次 nloc 循环中都会执行一次。

结果的 dist 值在哪里使用呢？

如果它是最深层循环的一部分，那么还可以做更多的优化（比如减少变量名，合并一些代码）。