Custom theano Op to do numerical integration

There is a small mistake in the gradient part:

...
for grad in grads:
    integrate = Integrate(grad, self._var, *self._extra_vars)
    ...

Also, the correct usage should be:

y_obs = 8.3

start = theano.shared(1.)
stop = theano.shared(2.)
with pm.Model() as basic_model:
    a = pm.Uniform('a', 1.5, 3.5)
    b = pm.Uniform('b', 4., 6.)

    # Define the function to integrate in plain theano
    t = tt.dscalar('t')
    t.tag.test_value = np.zeros(())
    a_ = tt.dscalar('a_')
    a_.tag.test_value = np.ones(())*2.
    b_ = tt.dscalar('b_')
    b_.tag.test_value = np.ones(())*5.
    func = t**a_ + b_
    integrate = Integrate(func, t, a_, b_)

    # Now we plug in the values from the model.
    # The `a_` and `b_` from above corresponds to the `a` and `b` here.
    mu = integrate(start, stop, a, b)
    y = pm.Normal('y', mu=mu, sd=0.4, observed=y_obs)
    trace = pm.sample(1500, tune=500, cores=2, chains=2)