Как ускорить модель марки PyMC?

Есть ли способ ускорить эту простую модель 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, если это необходимо. Спасибо.

Ответы

Ответ 1

Марковская цепь Монте-Карло является известной последовательной проблемой.

Это означает, что его время выполнения пропорционально количеству шагов и времени выполнения вашей фитнес-функции.

Есть некоторые трюки, которые вы можете сделать, однако:

  • Использовать 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)