Wie schreibt man eine benutzerdefinierte Deterministische oder Stochastische in pymc3 mit theano.op?

Question

Ich mache etwas pymc3 und ich möchte benutzerdefinierte Stochastik erstellen, aber es scheint nicht viel Dokumentation darüber zu sein, wie es gemacht wird. Ich weiß, wie man die as_op way verwendet, aber anscheinend macht es das unmöglich, den NUTS-Sampler zu verwenden, in welchem ​​Fall ich nicht den Vorteil von pymc3 gegenüber pymc sehe.

Das Tutorial erwähnt, dass es getan werden kann, indem er von theano.Op erbt. Aber kann mir jemand zeigen, wie das funktionieren würde (ich fange immer noch an auf theano)? Ich habe zwei Stochastiken, die ich definieren möchte.

Die erste einfacher sein sollte, ist es ein N Dimension Vectors F, die nur konstante Mutter Variablen hat:

with myModel: 
    F = DensityDist('F', lambda value: pymc.skew_normal_like(value, F_mu_array, F_std_array, F_a_array), shape = N) 

ich eine Skew-Normalverteilung wollen, die noch in pymc3 umgesetzt zu werden scheint nicht, Ich habe gerade die pymc2-Version importiert. Leider sind F_mu_array, F_std_array, F_a_array and F alle N-dimensionale Vektoren und das Lambda-Ding scheint nicht mit einer N-Dimension-Liste zu arbeiten value.

Erstens, gibt es eine Möglichkeit, den Lambda-Eingang zu einem N-dimensionalen Array zu machen? Wenn nicht, würde ich den Stochastic F direkt definieren müssen, und hier nehme ich an, dass ich theano.Op brauche, damit es funktioniert.

---

Das zweite Beispiel ist eine kompliziertere Funktion anderer Stochastik. Hier ist, wie ich will es definieren (fälschlicherweise im Moment):

with myModel: 
    ln2_var = Uniform('ln2_var', lower=-10, upper=4) 
    sigma = Deterministic('sigma', exp(0.5*ln2_var))   
    A = Uniform('A', lower=-10, upper=10, shape=5) 
    C = Uniform('C', lower=0.0, upper=2.0, shape=5) 
    sw = Normal('sw', mu=5.5, sd=0.5, shape=5) 

    # F from before 
    F = DensityDist('F', lambda value: skew_normal_like(value, F_mu_array, F_std_array, F_a_array), shape = N) 
    M = Normal('M', mu=M_obs_array, sd=M_stdev, shape=N) 

    # Radius forward-model (THIS IS THE STOCHASTIC IN QUESTION) 
    R = Normal('R', mu = R_forward(F, M, A, C, sw, N), sd=sigma, shape=N) 

Wo die Funktion R_forward(F,M,A,C,sw,N) ist naiv wie folgt definiert:

from theano.tensor import lt, le, eq, gt, ge 

def R_forward(Flux, Mass, A, C, sw, num): 
    for i in range(num): 
     if lt(Mass[i], 0.2): 
      if lt(Flux[i], sw[0]): 
       muR = C[0] 
      else: 
       muR = A[0]*log10(Flux[i]) + C[0] - A[0]*log10(sw[0]) 
     elif (le(0.2, Mass[i]) or le(Mass[i], 0.5)): 
      if lt(Flux[i], sw[1]): 
       muR = C[1] 
      else: 
       muR = A[1]*log10(Flux[i]) + C[1] - A[1]*log10(sw[1]) 
     elif (le(0.5, Mass[i]) or le(Mass[i], 1.5)): 
      if lt(Flux[i], sw[2]): 
       muR = C[2] 
      else: 
       muR = A[2]*log10(Flux[i]) + C[2] - A[2]*log10(sw[2]) 
     elif (le(1.5, Mass[i]) or le(Mass[i], 3.5)): 
      if lt(Flux[i], sw[3]): 
       muR = C[3] 
      else: 
       muR = A[3]*log10(Flux[i]) + C[3] - A[3]*log10(sw[3]) 
     else: 
      if lt(Flux[i], sw[4]): 
       muR = C[4] 
      else: 
       muR = A[4]*log10(Flux[i]) + C[4] - A[4]*log10(sw[4]) 
    return muR 

Dies vermutlich nicht natürlich funktionieren. Ich kann sehen, wie ich as_op verwenden würde, aber ich möchte die NUTS-Sampling beibehalten.

Answer 1

Ich weiß, das ist jetzt ein bisschen spät, aber ich dachte, ich würde die Frage (ziemlich vage) sowieso beantworten.

Wenn Sie eine stochastische Funktion (zB eine Wahrscheinlichkeitsverteilung) definieren wollen, dann müssen Sie ein paar Dinge tun:

Zunächst definieren eine Unterklasse von entweder Discrete (pymc3.distributions.Discrete) oder Continuous , die mindestens die Methode logp hat, die die Log-Likelihood Ihres Stochastic zurückgibt. Wenn Sie dies als eine einfache symbolische Gleichung (x + 1) definieren, glaube ich, dass Sie sich nicht um irgendwelche Gradienten kümmern müssen (aber zitieren Sie mich nicht dazu; see the documentation about this). Ich komme unten zu komplizierteren Fällen. In dem unglücklichen Fall, dass Sie etwas komplizierteres tun müssen, wie in Ihrem zweiten Beispiel (pymc3 hat nun übrigens eine Skew-Normalverteilung implementiert), müssen Sie die dafür erforderlichen Operationen (in der LogP-Methode verwendet) als definieren ein Theano Op. Wenn Sie keine Derivate benötigen, dann erledigt der as_op die Aufgabe, aber wie Sie sagten, sind Gradienten die Art von pymc3.

Hier wird es kompliziert. Wenn Sie NUTS verwenden möchten (oder aus irgendeinem Grund Gradienten benötigen), müssen Sie Ihre in logp als eine Unterklasse von theano.gof.Op verwendete Operation implementieren. Ihre neue Op-Klasse (von nun an einfach Op genannt) benötigt mindestens zwei oder drei Methoden. Die erste definiert Eingänge/Ausgänge für den Op (check the Op documentation). Die perform() - Methode (oder die Varianten, die Sie auswählen können) ist diejenige, die die gewünschte Operation ausführt (z. B. Ihre R_forward-Funktion). Dies kann in reinem Python gemacht werden, wenn Sie dies wünschen.Mit der dritten Methode grad() definieren Sie den Farbverlauf Ihrer perform() -Ausgabe für die Eingaben. Die tatsächliche Ausgabe von grad() ist ein bisschen anders, aber keine große Sache.

Und es ist in grad(), dass die Verwendung von Theano auszahlt. Wenn Sie Ihre gesamte perform() in Theano definieren, dann können Sie die automatische Differenzierung (theano.tensor.grad oder theano.tensor.jacobian) einfach verwenden, um die Arbeit für Sie zu erledigen (siehe Beispiel unten). Dies wird jedoch nicht unbedingt einfach sein.

In Ihrem zweiten Beispiel würde es bedeuten, Ihre R_forward-Funktion in Theano zu implementieren, was kompliziert sein könnte.

Hier schließe ich ein etwas minimales Beispiel eines Op ein, das ich schuf, während ich lernte, diese Dinge zu tun.

def my_th_fun(): 
    """ Some needed auxiliary functions. 
    """ 
    X = th.tensor.vector('X') 
    SCALE = th.tensor.scalar('SCALE') 

    X.tag.test_value = np.array([1,2,3,4]) 
    SCALE.tag.test_value = 5. 

    Scale, upd_sm_X = th.scan(lambda x, scale: scale*(scale+ x), 
           sequences=[X], 
           outputs_info=[SCALE]) 
    fun_Scale = th.function(inputs=[X, SCALE], outputs=Scale) 
    D_out_d_scale = th.tensor.grad(Scale[-1], SCALE) 
    fun_d_out_d_scale = th.function([X, SCALE], D_out_d_scale) 
    return Scale, fun_Scale, D_out_d_scale, fun_d_out_d_scale 

class myOp(th.gof.Op): 
    """ Op subclass with a somewhat silly computation. It uses 
    th.scan and th.tensor.grad is used to calculate the gradient 
    automagically in the grad() method. 
    """ 
    __props__ =() 
    itypes = [th.tensor.dscalar] 
    otypes = [th.tensor.dvector] 
    def __init__(self, *args, **kwargs): 
     super(myOp, self).__init__(*args, **kwargs) 
     self.base_dist = np.arange(1,5) 
     (self.UPD_scale, self.fun_scale, 
     self.D_out_d_scale, self.fun_d_out_d_scale)= my_th_fun() 

    def perform(self, node, inputs, outputs): 
     scale = inputs[0] 
     updated_scale = self.fun_scale(self.base_dist, scale) 
     out1 = self.base_dist[0:2].sum() 
     out2 = self.base_dist[2:4].sum() 
     maxout = np.max([out1, out2]) 
     exp_out1 = np.exp(updated_scale[-1]*(out1-maxout)) 
     exp_out2 = np.exp(updated_scale[-1]*(out2-maxout)) 
     norm_const = exp_out1 + exp_out2 
     outputs[0][0] = np.array([exp_out1/norm_const, exp_out2/norm_const]) 

    def grad(self, inputs, output_gradients): #working! 
     """ Calculates the gradient of the output of the Op wrt 
     to the input. As a simple example, the input is scalar. 

     Notice how the output is actually the gradient multiplied 
     by the output_gradients, which is an input provided by 
     theano when calculating gradients. 
     """ 
     scale = inputs[0] 
     X = th.tensor.as_tensor(self.base_dist) 
     # Do I need to recalculate all this or can I assume that perform() has 
     # always been called before grad() and thus can take it from there? 
     # In any case, this is a small enough example to recalculate quickly: 
     all_scale, _ = th.scan(lambda x, scale_1: scale_1*(scale_1+ x), 
           sequences=[X], 
           outputs_info=[scale]) 
     updated_scale = all_scale[-1] 

     out1 = self.base_dist[0:1].sum() 
     out2 = self.base_dist[2:3].sum() 
     maxout = np.max([out1, out2]) 

     exp_out1 = th.tensor.exp(updated_scale*(out1 - maxout)) 
     exp_out2 = th.tensor.exp(updated_scale*(out2 - maxout)) 
     norm_const = exp_out1 + exp_out2 

     d_S_d_scale = th.theano.grad(all_scale[-1], scale) 
     Jac1 = (-(out1-out2)*d_S_d_scale* 
       th.tensor.exp(updated_scale*(out1+out2 - 2*maxout))/(norm_const**2)) 
     Jac2 = -Jac1 
     return Jac1*output_gradients[0][0]+ Jac2*output_gradients[0][1], 

Diese Op kann dann innerhalb des logP() -Methode eines stochastischen in pymc3 verwendet werden:

import pymc3 as pm 

class myDist(pm.distributions.Discrete): 
    def __init__(self, invT, *args, **kwargs): 
     super(myDist, self).__init__(*args, **kwargs) 
     self.invT = invT 
     self.myOp = myOp() 
    def logp(self, value): 
     return self.myOp(self.invT)[value] 

Ich hoffe, es hilft jeden (hoffnungslos) pymc3/Theanos Neuling da draußen.