2015-12-08 19 views
12

Gibt es eine Möglichkeit, dieses einfache PyMC-Modell zu beschleunigen? Bei 20-40 Datenpunkten dauert es ~ 5-11 Sekunden, um zu passen.Wie beschleunigt man PyMC Markov Modell?

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) 

Mit nur 40 Datenpunkte, dauert es 11 Sekunden auszuführen:

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

(mit 80 Punkten Es dauert 20 Sekunden). Dies ist ein Spielzeugbeispiel. Die Ausdrücke innerhalb von Lambda(), die Übergänge bestimmen, sind in der Praxis komplexer. Diese grundlegende Codestruktur ist flexibel (während die Codierung des Modells mit Übergangsmatrizen weniger flexibel ist). Gibt es eine Möglichkeit, eine ähnliche Codestruktur beizubehalten, aber eine bessere Leistung zu erzielen? Ich freue mich, wenn nötig, zu PyMC3 wechseln zu können. Vielen Dank.

+0

Welche Version von pymc verwenden Sie? In der Pymc-Dokumentation für 2.3.6 kann ich die Bernoulli-Funktion nicht finden, nur Bernoulli_like [Doc] (https://pymc-devs.github.io/pymc/). – CodeMonkey

+0

es existiert in 2.2 – slushy

+0

Ich habe eine ähnliche Sorge über die Optimierung (https://stackoverflow.com/questions/42205123/how-to-fit-a-method-belonging-to-an-instance-with-pymc3) –

Antwort

3

Markov-Kette Monte Carlo ist ein bekanntes sequenzielles Problem.

Das bedeutet, dass die Laufzeit proportional zur Anzahl der Schritte und der Laufzeit Ihrer Fitnessfunktion ist.

Es gibt einige Tricks, die Sie tun können, aber:

  • Verwendung PyPy (erfordert umschreiben, pymc nicht unterstützt)
  • Verwenden Gibbs Sampling nächste Schritt
  • Verwenden Sie mehrere Startpunkte zu verbessern (in parallel)
  • Verwenden mehrere Zweige (parallel)
  • Verwenden heuristischer die Kette zu stoppen früheren
  • Verwendung Approximation f oder Punkte, die nahe an berechnet bereits

Härtere Ansätze:

  • Verwenden Numba (kompiliert zu Maschinencode Fitness-Funktion)
  • umschreiben Ihre Fitness-Funktion in C (oder ähnlich)
  • Verwendung nativer MCMC-Code (nicht Python, erfordert das oben)

Endlich gibt es eine Menge Forschung dort draußen:

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)

+0

Being eine Antwort auf eine sehr allgemeine Frage, ich versuchte, knapp zu bleiben. Bitte kommentieren für Details ... –