python - 自定义 Theano Op 进行数值积分

Question

我正在尝试编写一个自定义 Theano Op，它对两个值之间的函数进行数值积分。 Op 是 PyMC3 的自定义可能性，它涉及一些积分的数值评估。我不能简单地使用 @as_op 装饰器，因为我需要使用 HMC 来执行 MCMC 步骤。任何帮助将不胜感激，因为这个问题似乎已经出现了好几次但从未得到解决(例如 https://stackoverflow.com/questions/36853015/using-theano-with-numerical-integration ， Theano: implementing an integral function )。

很明显，一种解决方案是在 Theano 中编写一个数值积分器，但是当已经有非常好的积分器可用时，这似乎是一种浪费，例如通过 scipy.integrate。

为了将此作为最小示例，让我们尝试在 Op 中集成一个介于 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.

我对此感到惊讶，尤其是 TypeError，因为我认为我已经将 output_storage 变量转换为张量，但它似乎在这里相信它仍然是一个 ndarray。

Answer 1

我发现了你的问题，因为我试图在 PyMC3 中构建一个随机变量，它表示一般点过程(Hawkes、Cox、Poisson 等)并且似然函数具有积分。我真的希望能够使用哈密顿蒙特卡罗或 NUTS 采样器，所以我需要关于时间的积分是可微的。

从您的尝试开始，我制作了一个 integrateOut theano Op，它似乎可以正确处理我需要的行为。我已经在几个不同的输入上对其进行了测试(目前还没有在我的统计模型上进行测试，但它看起来很有希望!)。我完全是 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)