如何使用theano.op編寫定製的pymc3中的確定性或隨機性？

Question

我正在做一些pymc3，我想創建自定義的隨機指標，但似乎並沒有很多關於如何完成的文檔。我知道如何使用as_op way，但顯然這使得不可能使用NUTS採樣器，在這種情況下，我沒有看到pymc3相對於pymc的優勢。

該教程提到可以通過從theano.Op繼承來完成。但任何人都可以告訴我如何工作（我仍然開始對theano）？我有兩個我想定義的隨機指標。

第一個應該更容易，它是具有唯一不變的父變量的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維矢量，並且lambda似乎不適用於N維列表value。

首先，有沒有辦法讓lambda輸入成爲一個N維數組？如果沒有，我想我需要直接定義隨機指標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）的子類或持續，它至少有方法logp，它返回隨機的對數似然值。如果你將其定義爲一個簡單的符號方程（x + 1），我相信你不需要照顧任何漸變（但不要在此引用我; see the documentation about this）。我將在下面討論更復雜的案例。在不幸的情況下，你需要做更復雜的事情，就像在第二個例子中一樣（順便說一句pymc3現在實現了偏斜正態分佈），你需要定義它所需的操作（在logp方法中使用）爲一個Theano操作。如果你不需要衍生物，那麼as_op就可以完成這項工作，但正如你所說的那樣，漸變是pymc3的想法。

這就是它變得複雜的地方。如果你想使用NUTS（或者出於任何原因需要漸變），那麼你需要將你在logp中使用的操作作爲theano.gof.Op的子類來實現。您的新操作類（從現在開始稱它爲「操作」）至少需要兩個或三個方法。第一個定義了Op的輸入/輸出（check the Op documentation）。 perform（）方法（或可能選擇的變體）是執行所需操作（例如，R_forward函數）的方法。如果你願意，這可以在純Python中完成。第三種方法grad（）是您定義perform（）的輸出和輸入的梯度的地方。到grad（）的實際輸出有點不同，但不是什麼大不了的。

而且在使用Theano的grad（）中也有回報。如果您在Theano中定義完整的perform（），那麼可能很容易使用自動區分（theano.tensor.grad或theano.tensor.jacobian）爲您完成工作（請參見下面的示例）。但是，這並不一定容易。

在第二個例子中，這意味着在Theano中實現您的R_forward函數，這可能很複雜。

這裏我包含了一個在學習做這些事情時創建的作業的一個有點簡單的例子。

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

這個作品就可以隨機的的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新手在那裏。