12

有没有办法加速这个简单的 PyMC 模型？在 20-40 个数据点上，拟合大约需要 5-11 秒。

import pymc
import time
import numpy as np
from collections import OrderedDict

# prior probability of rain
p_rain = 0.5
variables = OrderedDict()
# rain observations
data = [True, True, True, True, True,
        False, False, False, False, False]*4
num_steps = len(data)
p_rain_given_rain = 0.9
p_rain_given_norain = 0.2
p_umbrella_given_rain = 0.8
p_umbrella_given_norain = 0.3
for n in range(num_steps):
    if n == 0:
        # Rain node at time t = 0
        rain = pymc.Bernoulli("rain_%d" %(n), p_rain)
    else:
        rain_trans = \
          pymc.Lambda("rain_trans",
                      lambda prev_rain=variables["rain_%d" %(n-1)]: \
                      prev_rain*p_rain_given_rain + (1-prev_rain)*p_rain_given_norain)
        rain = pymc.Bernoulli("rain_%d" %(n), p=rain_trans)
    umbrella_obs = \
      pymc.Lambda("umbrella_obs",
                  lambda rain=rain: \
                  rain*p_umbrella_given_rain + (1-rain)*p_umbrella_given_norain)
    umbrella = pymc.Bernoulli("umbrella_%d" %(n), p=umbrella_obs,
                              observed=True,
                              value=data[n])
    variables["rain_%d" %(n)] = rain
    variables["umbrella_%d" %(n)] = umbrella

print "running on %d points" %(len(data))
all_vars = variables.values()
t_start = time.time()
model = pymc.Model(all_vars)
m = pymc.MCMC(model)
m.sample(iter=2000)
t_end = time.time()
print "\n%.2f secs to run" %(t_end - t_start)

只有 40 个数据点，运行需要 11 秒：

running on 40 points
 [-----------------100%-----------------] 2000 of 2000 complete in 11.5 sec
11.54 secs to run

（80分需要20秒）。这是一个玩具示例。Lambda()确定转换的内部表达式实际上更复杂。这种基本的代码结构是灵活的（而使用转换矩阵对模型进行编码不太灵活）。有没有办法保持类似的代码结构但获得更好的性能？如有必要，很高兴切换到 PyMC3。谢谢。

4

1 回答 1

6

Markov-Chain Monte Carlo 是一个已知的序列问题。

这意味着它的运行时间与健身功能的步数和运行时间成正比。

但是，您可以使用一些技巧：

  • 使用 PyPy（需要重写，不支持 pymc）
  • 使用 Gibbs 抽样来改进您的下一步
  • 使用多个起点（并行）
  • 使用多个分支（并行）
  • 使用启发式更早地停止链条
  • 对接近已计算的点使用近似值

更难的方法：

  • 使用 Numba（将适应度函数编译为机器码）
  • 用 C（或类似语言）重写你的适应度函数
  • 使用本机 MCMC 代码（非 Python，需要以上内容）

最后有很多研究：

http://www.mas.ncl.ac.uk/~ndjw1/docs/pbc.pdf

https://sites.google.com/site/parallelmcmc/

http://pyinsci.blogspot.com/2010/12/efficcient-mcmc-in-python.html (pypy)

于 2015-12-18T11:04:04.257 回答