AssertionError in PyMC3-Code mit deterministischen Distributionen

Question

Ich habe in den letzten Tagen versucht, sich an PyMC zu gewöhnen, einige MCMC Sampling von Verteilungen von einigen Modellen zu machen, für die ich direkte Codes habe (ich bin letztlich interessiert an Parameterschätzung).

Nach meinem Wissen gibt es nicht so viele Beispiele von Leuten, die ihre Codes zeigen (wie ein externer C- oder FORTRAN-Code), dass sie erfolgreich mit PyMC3 arbeiten. Ich habe here oder here Bestechungsgelder bisher gefunden.

Um mit einfachen Problemen zu beginnen, habe ich versucht, mit PyMC3 einige Ergebnisse von bestehenden Python-Codes zu reproduzieren, die MCMC mit "komplizierten" Distributionen ausführen, die nicht PyMC verwenden, nämlich diese einfache result, die in läuft 2 Sekunden für 10000 Proben auf meinem Laptop.

deterministische Variablen in PyMC3 zu definieren, sollte man den @theano.compile.ops.as_op Dekorateur verwenden (es war @deterministic in PyMC, die jetzt in PyMC3 veraltet ist), und das ist, was ich tat.

Hier ist mein Code (ich benutze Spyder und so IPython), durch die PyMC inspiriert documentation, in dem ich ein AssertionError nach der zweiten Zelle (im Schätzungsprozess, so vor der Probenahme) zu begegnen. Ich habe versucht, dies für die letzten zwei Tage zu lösen, aber ich verstehe wirklich nicht, was der Fehler sein kann. Ich glaube, es sollte eine Art PyMC oder Theano Trick sein, den ich noch nicht verstanden habe, da ich glaube, dass ich sehr nahe bin.

#%% Define model 
import numpy,math 
import matplotlib.pyplot as plt 
import random as random 

import theano.tensor as t 
import theano 

random.seed(1) # set random seed 

# copy-pasted function from the specific model used in the source and adapted with as_op 
@theano.compile.ops.as_op(itypes=[t.iscalar, t.dscalar, t.fscalar, t.fscalar],otypes=[t.dvector]) 
def sampleFromSalpeter(N, alpha, M_min, M_max):   
    log_M_Min = math.log(M_min) 
    log_M_Max = math.log(M_max) 
    maxlik = math.pow(M_min, 1.0 - alpha) 
    Masses = [] 
    while (len(Masses) < N):    
     logM = random.uniform(log_M_Min,log_M_Max) 
     M = math.exp(logM)    
     likelihood = math.pow(M, 1.0 - alpha)    
     u = random.uniform(0.0,maxlik) 
     if (u < likelihood): 
      Masses.append(M) 
    return Masses 

# SAME function as above, used to make test data (so no Theano here) 
def sampleFromSalpeterDATA(N, alpha, M_min, M_max):   
    log_M_Min = math.log(M_min) 
    log_M_Max = math.log(M_max)   
    maxlik = math.pow(M_min, 1.0 - alpha)   
    Masses = []   
    while (len(Masses) < N):    
     logM = random.uniform(log_M_Min,log_M_Max) 
     M = math.exp(logM)    
     likelihood = math.pow(M, 1.0 - alpha)    
     u = random.uniform(0.0,maxlik) 
     if (u < likelihood): 
      Masses.append(M) 
    return Masses 

# Generate toy data. 
N  = 1000000 # Draw 1 Million stellar masses. 
alpha = 2.35 
M_min = 1.0 
M_max = 100.0 
Masses = sampleFromSalpeterDATA(N, alpha, M_min, M_max) 

#%% Estimation process 
import pymc3 as pm 
basic_model = pm.Model() 
with basic_model: 

    # Priors for unknown model parameters 
    alpha2 = pm.Normal('alpha2', mu=3, sd=10) 

    N2=t.constant(1000000) 
    M_min2 = t.constant(1.0) 
    M_max2 = t.constant(100.0) 

    # Expected value of outcome 
    m = sampleFromSalpeter(N2, alpha2, M_min2, M_max2) 

    # Likelihood (sampling distribution) of observations 
    Y_obs = pm.Normal('Y_obs', mu=m, sd=10, observed=Masses) 

#%% Sample 
with basic_model: 
    step = pm.Metropolis() 
    trace = pm.sample(10000, step=step) 

#%% Plot posteriors 
_ = pm.traceplot(trace) 
pm.summary(trace) 

Answer 1

Ich bekam eine Antwort auf die Discourse of PyMC. Ich war in der Tat ziemlich nah, da ich nur die letzte Zeile der Funktionsdefinition ändern musste, um [numpy.array (Masses)] zurückzugeben.

Nun scheint es, dass dieser Code etwa 300h benötigt (statt 2 Sekunden in der ursprünglichen Implementierung ohne PyMC), und ich fragte mich, ob jemand eine Idee hatte, warum es so viel langsamer ist? Ich hatte tatsächlich wirklich einfache Distributionen früher und mit as_op versucht, und bemerkte eine erhebliche Verlangsamung, aber nicht einmal um so viel. Ich denke, der ursprüngliche Code funktioniert direkt auf dem LogP, ist es, was macht den ganzen Unterschied?