Wie ein Behälter in pymc3 erstellen

Question

Ich versuche, ein Modell für die Wahrscheinlichkeitsfunktion eines bestimmten Ergebnisses einer Langevin Gleichung (Brownsche Teilchen in einem harmonischen Potential) zu bauen:

Hier mein Modell in pymc2 ist, das scheint zu arbeiten: https://github.com/hstrey/BayesianAnalysis/blob/master/Langevin%20simulation.ipynb

#define the model/function to be fitted. 
def model(x): 
    t = pm.Uniform('t', 0.1, 20, value=2.0) 
    A = pm.Uniform('A', 0.1, 10, value=1.0) 

    @pm.deterministic(plot=False) 
    def S(t=t): 
     return 1-np.exp(-4*delta_t/t) 

    @pm.deterministic(plot=False) 
    def s(t=t): 
     return np.exp(-2*delta_t/t) 

    path = np.empty(N, dtype=object) 

    path[0]=pm.Normal('path_0',mu=0, tau=1/A, value=x[0], observed=True) 
    for i in range(1,N): 
     path[i] = pm.Normal('path_%i' % i, 
         mu=path[i-1]*s, 
         tau=1/A/S, 
         value=x[i], 
         observed=True) 
     return locals() 

mcmc = pm.MCMC(model(x)) 
mcmc.sample(20000, 2000, 10) 

die Grundidee ist, dass jeder Punkt in der Kette (Markow-Kette) auf dem vorherigen Punkt abhängt. Btw, x ist ein Array von Daten, N ist seine Länge, delta_t ist der Zeitschritt = 0,01. Irgendeine Idee, wie man das in pymc3 implementiert? Ich habe versucht:

# define the model/function for diffusion in a harmonic potential 
DHP_model = pm.Model() 
with DHP_model: 
    t = pm.Uniform('t', 0.1, 20) 
    A = pm.Uniform('A', 0.1, 10) 

    S=1-pm.exp(-4*delta_t/t) 

    s=pm.exp(-2*delta_t/t) 

    path = np.empty(N, dtype=object) 

    path[0]=pm.Normal('path_0',mu=0, tau=1/A, observed=x[0]) 
    for i in range(1,N): 
     path[i] = pm.Normal('path_%i' % i, 
         mu=path[i-1]*s, 
         tau=1/A/S, 
         observed=x[i]) 

Leider stürzt das Modell, sobald ich versuche, es zu starten. Ich habe einige pymc3-Beispiele (Tutorial) auf meinem Rechner ausprobiert und das funktioniert.

Vielen Dank im Voraus. Ich hoffe wirklich, dass die neuen Sampler in pymc3 mir bei diesem Modell helfen werden. Ich versuche Bayessche Methoden auf Einzelmolekül-Experimente anzuwenden.

Answer 1

Statt viele einzelne normalverteilte 1-D-Variablen in einer Schleife zu erstellen, können Sie eine benutzerdefinierte Verteilung (durch Erweitern Continuous) erstellen, die die Formel zum Berechnen der Log-Wahrscheinlichkeit Ihres gesamten Pfads kennt. Sie können diese Likelihood-Formel von der Normal-Likelihood-Formel, die pymc3 bereits kennt, aus laden. Ein Beispiel finden Sie in the built-in AR1 class.

Da Ihr Teilchen die Markov-Eigenschaft folgt, sieht Ihre Wahrscheinlichkeit wie

import theano.tensor as T 

def logp(path): 
    now = path[1:] 
    prev = path[:-1] 

    loglik_first = pm.Normal.dist(mu=0., tau=1./A).logp(path[0]) 
    loglik_rest = T.sum(pm.Normal.dist(mu=prev*ss, tau=1./A/S).logp(now)) 
    loglik_final = loglik_first + loglik_rest 

    return loglik_final 

Ich vermute, dass Sie in diesem Fall einen Wert für ss bei jedem Zeitschritt ziehen möchten, sollten Sie sicherstellen, angeben ss = pm.exp(..., shape=len(x)-1), so dass prev*ss im obigen Block als elementweise Multiplikation interpretiert wird.

Dann können Sie nur geben Sie Ihre Beobachtungen mit

path = MyLangevin('path', ..., observed=x) 

Dieses viel schneller laufen soll.

Answer 2

Da ich keine Antwort auf meine Frage sah, lassen Sie mich selbst antworten. Ich kam mit der folgenden Lösung:

# now lets model this data using pymc 
# define the model/function for diffusion in a harmonic potential 
DHP_model = pm.Model() 
with DHP_model: 
    D = pm.Gamma('D',mu=mu_D,sd=sd_D) 
    A = pm.Gamma('A',mu=mu_A,sd=sd_A) 

    S=1.0-pm.exp(-2.0*delta_t*D/A) 

    ss=pm.exp(-delta_t*D/A) 

    path=pm.Normal('path_0',mu=0.0, tau=1/A, observed=x[0]) 
    for i in range(1,N): 
     path = pm.Normal('path_%i' % i, 
         mu=path*ss, 
         tau=1.0/A/S, 
         observed=x[i]) 

    start = pm.find_MAP() 
    print(start) 
    trace = pm.sample(100000, start=start) 

leider dieser Code nimmt an N = 50 irgendwo zwischen 6 Stunden bis zu 2 Tagen zu kompilieren. Ich laufe auf einem ziemlich schnellen PC (24Gb RAM) mit Ubuntu. Ich habe versucht, die GPU zu verwenden, aber das läuft etwas langsamer. Ich vermute Speicherprobleme, da es beim Ausführen 99,8% des Speichers belegt. Ich habe die gleiche Berechnung mit Stan versucht und es dauert nur 2 Minuten zu laufen.