Comment écrire un déterministe ou stochastique personnalisé dans pymc3 avec theano.op?

Question

Je fais du pymc3 et je voudrais créer des stochastiques personnalisés, mais il ne semble pas y avoir beaucoup de documentation sur la façon dont cela est fait. Je sais comment utiliser le as_op way, mais apparemment cela rend impossible l'utilisation de l'échantillonneur NUTS, auquel cas je ne vois pas l'avantage de pymc3 sur pymc.

Le tutoriel mentionne que cela peut être fait en héritant de theano.Op. Mais est-ce que quelqu'un peut me montrer comment cela fonctionnerait (je commence toujours à le faire)? J'ai deux stochastiques que je veux définir.

Le premier devrait être plus facile, il est un vecteur de dimension N F qui n'a que des variables parentales constantes:

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

Je veux un biais distribution normale, ce qui ne semble pas être encore mis en œuvre dans pymc3, Je viens d'importer la version pymc2. Malheureusement, F_mu_array, F_std_array, F_a_array and F sont tous des vecteurs N-dimensionnels, et la chose lambda ne semble pas fonctionner avec une liste N-dimension value.

Premièrement, existe-t-il un moyen de faire de l'entrée lambda un réseau de dimension N? Si non, je suppose que je devrais définir le stochastique F directement, et c'est là où je suppose que j'ai besoin de theano.Op pour le faire fonctionner.

---

Le deuxième exemple est une fonction plus compliquée d'autres stochastiques. Voici comment je veux le définir (à tort pour le 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) 

Lorsque la fonction est R_forward(F,M,A,C,sw,N) définie comme naïvement:

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 

Ce probablement ne fonctionnera pas bien sûr. Je peux voir comment j'utiliserais as_op, mais je veux préserver l'échantillonnage NUTS.

Answer 1

Je me rends compte que c'est un peu tard maintenant, mais je pensais que je répondrais à la question (plutôt vaguement) de toute façon.

Si vous souhaitez définir une fonction stochastique (par exemple, une distribution de probabilité), alors vous devez faire deux choses:

D'abord, définir une sous-classe de l'une discrète (pymc3.distributions.Discrete) ou continu , qui a au moins la méthode logp, qui renvoie la log-vraisemblance de votre stochastique. Si vous définissez cela comme une simple équation symbolique (x + 1), je crois que vous n'avez pas besoin de prendre en compte les gradients (mais ne me citez pas sur ceci: see the documentation about this). Je vais passer à des cas plus compliqués ci-dessous. Dans le cas malheureux où vous avez besoin de faire quelque chose de plus complexe, comme dans votre deuxième exemple (pymc3 a maintenant une distribution normale de skew implémentée, d'ailleurs), vous devez définir les opérations nécessaires (utilisées dans la méthode logp) comme un Theano Op. Si vous n'avez besoin d'aucun dérivé, alors le as_op fait le travail, mais comme vous l'avez dit, les dégradés sont une sorte d'idée de pymc3.

C'est là que ça devient compliqué. Si vous voulez utiliser NUTS (ou avoir besoin de gradients pour une raison quelconque), alors vous devez implémenter votre opération utilisée dans logp en tant que sous-classe de theano.gof.Op. Votre nouvelle classe op (appelons-la juste Op à partir de maintenant) aura besoin de deux ou trois méthodes au moins. Le premier définit les entrées/sorties de l'op (check the Op documentation). La méthode perform() (ou les variantes que vous pouvez choisir) est celle qui effectue l'opération que vous voulez (votre fonction R_forward, par exemple). Cela peut être fait en python pur, si vous le souhaitez.La troisième méthode, grad(), est l'endroit où vous définissez le gradient de la sortie de votre performance() par rapport aux entrées. La sortie réelle de grad() est un peu différente, mais pas un gros problème.

Et c'est dans l'application grad() que l'utilisation de Theano est payante. Si vous définissez votre performance entière() dans Theano, alors il se peut que vous puissiez facilement utiliser la différenciation automatique (theano.tensor.grad ou theano.tensor.jacobian) pour faire le travail pour vous (voir l'exemple ci-dessous). Cependant, cela ne va pas nécessairement être facile. Dans votre deuxième exemple, cela impliquerait l'implémentation de votre fonction R_forward dans Theano, ce qui pourrait être compliqué.

Ici, j'inclus un exemple un peu minimal d'une Op que j'ai créé en apprenant à faire ces choses.

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], 

Cette op peut ensuite être utilisé dans la méthode LogP() d'un stochastique pymc3:

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] 

J'espère que cela aide tout (sans espoir) débutant pymc3/Théano là-bas.