Как написать пользовательский детерминированный или стохастический в pymc3 с theano.op?

Question

Я делаю несколько pymc3, и я хотел бы создать пользовательский Stochastics, однако, похоже, нет много документации о том, как это делается. Я знаю, как использовать as_op way, однако, по-видимому, это делает невозможным использование пробоотборника NUTS, и в этом случае я не вижу преимущества pymc3 над pymc.

В учебнике упоминается, что это можно сделать, наследуя от theano.Op. Но может ли кто-нибудь показать мне, как это будет работать (я все еще начинаю на анано)? У меня есть два Стохастика, которые я хочу определить.

Первый должен быть легче, что это вектор размерности N F, что имеет только постоянные родительские переменные:

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

хочу перекоса нормальное распределение, которое, кажется, не будет реализован в pymc3 еще, Я только что импортировал версию pymc2. К сожалению, F_mu_array, F_std_array, F_a_array and F - все N-мерные векторы, и предмет лямбда, похоже, не работает с N-мерным списком value.

Во-первых, есть ли способ сделать лямбда-вход N-мерным массивом? Если нет, я думаю, мне нужно было бы определить Stochastic F напрямую, и именно здесь я предполагаю, что мне нужен theano.Op, чтобы он работал.

---

Второй пример - более сложная функция других стохастик. Вот как я хочу, чтобы определить его (неправильно в данный момент):

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) 

Если функция R_forward(F,M,A,C,sw,N) наивно определяется как:

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 

Предположительно, это не будет работать, конечно. Я вижу, как я буду использовать as_op, но я хочу сохранить выборку NUTS.

Answer 1

Я понимаю, что это немного поздно, но я думал, что отвечу на вопрос (скорее смутно) в любом случае.

Если вы хотите определить вероятностную функцию (например, распределение вероятностей), то вам нужно сделать несколько вещей:

Первое, определить подкласс либо дискретной (pymc3.distributions.Discrete) или Continuous , который имеет хотя бы метод logp, который возвращает логарифмическую вероятность вашего стохастика. Если вы определяете это как простое символическое уравнение (x + 1), я считаю, что вам не нужно заботиться о каких-либо градиентах (но не цитируйте меня на этом: see the documentation about this). Ниже я расскажу о более сложных случаях. В неудачном случае, когда вам нужно сделать что-то более сложное, как в вашем втором примере (кстати, у pymc3 теперь есть косое нормальное распределение), вам нужно определить необходимые ему операции (используемые в методе logp), так как Theano Op. Если вам не нужны производные инструменты, то as_op выполняет эту работу, но, как вы сказали, градиенты - это своего рода идея pymc3.

Здесь все усложняется. Если вы хотите использовать NUTS (или нужны градиенты по любой причине), тогда вам нужно реализовать свою операцию, используемую в logp, в качестве подкласса theano.gof.Op. Ваш новый класс op (назовем его просто Op с этого момента) потребует, по крайней мере, два или три метода. Первый определяет входы/выходы для Op (check the Op documentation). Метод perform() (или варианты, которые вы можете выбрать) - это тот, который выполняет требуемую операцию (например, ваша функция R_forward). Это можно сделать в чистом питоне, если хотите.Третий метод grad() - это то, где вы определяете градиент вашего вывода perform() по входам. Фактический вывод в grad() немного отличается, но не большой.

И именно в grad() используется использование Anano. Если вы определите всю свою работу() в Theano, возможно, вы можете легко использовать автоматическое дифференцирование (theano.tensor.grad или theano.tensor.jacobian), чтобы выполнить эту работу для вас (см. Пример ниже). Однако это не обязательно будет легко.

В вашем втором примере это означало бы реализацию вашей функции R_forward в Theano, которая может быть сложной.

Здесь я включаю несколько минимальный пример Op, который я создал, учась делать это.

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

Этот Op может быть использован внутри метода LogP() стохастического в 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] 

Я надеюсь, что это помогает любому (безнадежной) pymc3/Theano новичку там.