我发现了您的问题,因为我试图在PyMC3中构建一个表示一般点过程(Hawkes,Cox,Poisson等)的随机变量,并且似然函数具有一个积分。我真的很想能够使用Hamiltonian
Monte Carlo或NUTS采样器,因此我需要相对于时间的积分是可区分的。
从您的尝试开始,我制作了一个IntegratedOut theano
Op,它似乎可以正确处理我需要的行为。我已经在一些不同的输入上对其进行了测试(尚未在我的统计模型中进行过测试,但是看起来很有希望!)。我总是theano
n00b,所以请原谅任何愚蠢。如果有人有任何反馈,我将不胜感激。不确定它是否正是您要的东西,但这是我的解决方案(示例位于底部和doc字符串中)。*编辑:简化了一些残留的拧干方式。
import theanoimport theano.tensor as Tfrom scipy.integrate import quadclass 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)*tintZ = 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)


