如何加速PyMC马尔可夫模型?

13

有没有办法加速这个简单的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。谢谢。


你使用的是哪个版本的pymc?在2.3.6的pymc文档中,我找不到Bernoulli函数,只有Bernoulli_like 文档 - CodeMonkey
我对优化也有类似的担忧(https://dev59.com/Z5_ha4cB1Zd3GeqPtRmW) - Stéphane
1个回答

6

马尔科夫链蒙特卡罗是一个已知的顺序问题。

这意味着它的运行时间与步数和您的适应度函数的运行时间成比例。

然而,有一些技巧可以使用:

  • 使用PyPy(需要重写,不支持pymc)
  • 使用吉布斯采样来改进下一步
  • 使用多个起点(并行)
  • 使用多个分支(并行)
  • 使用启发式方法提前停止链
  • 对于接近已经计算过的点,使用近似值

更难的方法:

  • 使用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)


作为对一个非常普遍的问题的回答,我尽量保持简洁。请针对具体情况进行评论... - Dima Tisnek
在pymc3中,默认采样是NUTS;Gibbs更快吗? - Itamar Mushkin
1
我不是专家,可以参考 https://stats.stackexchange.com/questions/311813/can-somebody-explain-to-me-nuts-in-english 。我的观点是:这可能取决于问题;https://gsejournal.biomedcentral.com/articles/10.1186/s12711-019-0515-1 表明总体性能是“可比的”。 - Dima Tisnek

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接