定製Theano運算數值積分

Question

我試圖寫一個自定義Theano歐普其數字集成了兩個值之間的函數。 Op是一種PyMC3的定製可能性，涉及一些積分的數值計算。我不能簡單地使用@as_op裝飾器，因爲我需要使用HMC來執行MCMC步驟。任何幫助將不勝感激，因爲這個問題似乎已經走到了好幾次，但一直沒有得到解決（例如https://stackoverflow.com/questions/36853015/using-theano-with-numerical-integration，Theano: implementing an integral function）。

顯然，一個解決辦法是寫中Theano數值積分，但是這似乎是白費力氣的時候很好的集成已經上市，例如通過scipy.integrate。

爲了保持作爲一個最小的例子，讓我們只是嘗試和0和1之間的整合運算內部函數。以下內容集成了Op之外的Theano函數，並且在我的測試結束後生成正確的結果。

import theano 
import theano.tensor as tt 
from scipy.integrate import quad 

x = tt.dscalar('x') 
y = x**4 # integrand 
f = theano.function([x], y) 

print f(0) 
print f(1) 

ans = integrate.quad(f, 0, 1)[0] 

print ans 

但是，嘗試在Op中進行整合看起來要困難得多。我現在最大的努力爲：

import numpy as np 
import theano 
import theano.tensor as tt 
from scipy import integrate 

class IntOp(theano.Op): 
    __props__ =() 

    def make_node(self, x): 
     x = tt.as_tensor_variable(x) 
     return theano.Apply(self, [x], [x.type()]) 

    def perform(self, node, inputs, output_storage): 
     x = inputs[0] 
     z = output_storage[0] 

     f_to_int = theano.function([x], x) 
     z[0] = tt.as_tensor_variable(integrate.quad(f_to_int, 0, 1)[0]) 

    def infer_shape(self, node, i0_shapes): 
     return i0_shapes 

    def grad(self, inputs, output_grads): 
     ans = integrate.quad(output_grads[0], 0, 1)[0] 
     return [ans] 

intOp = IntOp() 

x = tt.dmatrix('x') 
y = intOp(x) 

f = theano.function([x], y) 

inp = np.asarray([[2, 4], [6, 8]], dtype=theano.config.floatX) 
out = f(inp) 

print inp 
print out 

其中提供了以下錯誤：

Traceback (most recent call last): 
    File "stackoverflow.py", line 35, in <module> 
    out = f(inp) 
    File "/usr/local/lib/python2.7/dist-packages/theano/compile/function_module.py", line 871, in __call__ 
    storage_map=getattr(self.fn, 'storage_map', None)) 
    File "/usr/local/lib/python2.7/dist-packages/theano/gof/link.py", line 314, in raise_with_op 
    reraise(exc_type, exc_value, exc_trace) 
    File "/usr/local/lib/python2.7/dist-packages/theano/compile/function_module.py", line 859, in __call__ 
    outputs = self.fn() 
    File "/usr/local/lib/python2.7/dist-packages/theano/gof/op.py", line 912, in rval 
    r = p(n, [x[0] for x in i], o) 
    File "stackoverflow.py", line 17, in perform 
    f_to_int = theano.function([x], x) 
    File "/usr/local/lib/python2.7/dist-packages/theano/compile/function.py", line 320, in function 
    output_keys=output_keys) 
    File "/usr/local/lib/python2.7/dist-packages/theano/compile/pfunc.py", line 390, in pfunc 
    for p in params] 
    File "/usr/local/lib/python2.7/dist-packages/theano/compile/pfunc.py", line 489, in _pfunc_param_to_in 
    raise TypeError('Unknown parameter type: %s' % type(param)) 
TypeError: Unknown parameter type: <type 'numpy.ndarray'> 
Apply node that caused the error: IntOp(x) 
Toposort index: 0 
Inputs types: [TensorType(float64, matrix)] 
Inputs shapes: [(2, 2)] 
Inputs strides: [(16, 8)] 
Inputs values: [array([[ 2., 4.], 
     [ 6., 8.]])] 
Outputs clients: [['output']] 

Backtrace when the node is created(use Theano flag traceback.limit=N to make it longer): 
    File "stackoverflow.py", line 30, in <module> 
    y = intOp(x) 
    File "/usr/local/lib/python2.7/dist-packages/theano/gof/op.py", line 611, in __call__ 
    node = self.make_node(*inputs, **kwargs) 
    File "stackoverflow.py", line 11, in make_node 
    return theano.Apply(self, [x], [x.type()]) 

HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node. 

我感到意外，特別是類型錯誤，因爲我以爲我已經轉換的output_storage變量爲張量，但它似乎相信它仍然是一個ndarray。

Answer 1

我發現你的問題，因爲我試圖建立PyMC3表示一般的點過程（霍克斯，考克斯，泊松等）的隨機變量和似然函數有一個整體。我真的希望能夠使用哈密爾頓蒙特卡羅或NUTS採樣器，所以我需要時間積分來區分。

出發的嘗試，我做了一個integrateOut theano作品，似乎與我所需要的行爲正常工作。我已經在幾個不同的輸入上進行了測試（不是我的統計模型，但看起來很有前途！）。我總是theano n00b，所以請原諒任何愚蠢。如果有人有任何問題，我將不勝感激。不確定這正是你在找什麼，但這裏是我的解決方案（例如在底部和文檔字符串中）。 *編輯：簡化了一些殘留的方式來做到這一點。

import theano 
import theano.tensor as T 
from scipy.integrate import quad 

class integrateOut(theano.Op): 
    """ 
    Integrate out a variable from an expression, computing 
    the definite integral w.r.t. the variable specified 
    !!! Only implemented in this for scalars !!! 


    Parameters 
    ---------- 
    f : scalar 
     input 'function' to integrate 
    t : scalar 
     the variable to integrate out 
    t0: float 
     lower integration limit 
    tf: float 
     upper integration limit 

    Returns 
    ------- 
    scalar 
     a new scalar with the 't' integrated out 

    Notes 
    ----- 

    usage of this looks like: 
    x = T.dscalar('x') 
    y = T.dscalar('y') 
    t = T.dscalar('t') 

    z = (x**2 + y**2)*t 

    # integrate z w.r.t. t as a function of (x,y) 
    intZ = integrateOut(z,t,0.0,5.0)(x,y) 
    gradIntZ = T.grad(intZ,[x,y]) 

    funcIntZ = theano.function([x,y],intZ) 
    funcGradIntZ = theano.function([x,y],gradIntZ) 

    """ 
    def __init__(self,f,t,t0,tf,*args,**kwargs): 
     super(integrateOut,self).__init__() 
     self.f = f 
     self.t = t 
     self.t0 = t0 
     self.tf = tf 

    def make_node(self,*inputs): 
     self.fvars=list(inputs) 
     # This will fail when taking the gradient... don't be concerned 
     try: 
      self.gradF = T.grad(self.f,self.fvars) 
     except: 
      self.gradF = None 
     return theano.Apply(self,self.fvars,[T.dscalar().type()]) 

    def perform(self,node, inputs, output_storage): 
     # Everything else is an argument to the quad function 
     args = tuple(inputs) 
     # create a function to evaluate the integral 
     f = theano.function([self.t]+self.fvars,self.f) 
     # actually compute the integral 
     output_storage[0][0] = quad(f,self.t0,self.tf,args=args)[0] 

    def grad(self,inputs,grads): 
     return [integrateOut(g,self.t,self.t0,self.tf)(*inputs)*grads[0] \ 
      for g in self.gradF] 

x = T.dscalar('x') 
y = T.dscalar('y') 
t = T.dscalar('t') 

z = (x**2+y**2)*t 

intZ = integrateOut(z,t,0,1)(x,y) 
gradIntZ = T.grad(intZ,[x,y]) 
funcIntZ = theano.function([x,y],intZ) 
funcGradIntZ = theano.function([x,y],gradIntZ) 
print funcIntZ(2,2) 
print funcGradIntZ(2,2) 

Answer 2

SymPy被證明比預期困難，但在任何情況下，同時是查找有用的話，我也將指出如何修改這個作品允許更改最後的時間點，而無需創建一個新的作品。如果你有一個點過程，或者你的時間測量有不確定性，這可能很有用。

class integrateOut2(theano.Op): 
    def __init__(self, f, int_var, *args,**kwargs): 
     super(integrateOut2,self).__init__() 
     self.f = f 
     self.int_var = int_var 

    def make_node(self, *inputs): 
     tmax = inputs[0] 
     self.fvars=list(inputs[1:]) 

     return theano.Apply(self, [tmax]+self.fvars, [T.dscalar().type()]) 

    def perform(self, node, inputs, output_storage): 
     # Everything else is an argument to the quad function 
     tmax = inputs[0] 
     args = tuple(inputs[1:]) 

     # create a function to evaluate the integral 
     f = theano.function([self.int_var]+self.fvars, self.f) 

     # actually compute the integral 
     output_storage[0][0] = quad(f, 0., tmax, args=args)[0] 

    def grad(self, inputs, grads): 
     tmax = inputs[0] 
     param_grads = T.grad(self.f, self.fvars) 

     ## Recall fundamental theorem of calculus 
     ## d/dt \int^{t}_{0}f(x)dx = f(t) 
     ## So sub in t_max to the graph 
     FTC_grad = theano.clone(self.f, {self.int_var: tmax}) 

     grad_list = [FTC_grad*grads[0]] + \ 
        [integrateOut2(grad_fn, self.int_var)(*inputs)*grads[0] \ 
        for grad_fn in param_grads] 

     return grad_list