PyMC3 : Using a parameter as a limit of integral

Questions
1

I have been using PyMC for a while, and recently started using PyMC3.
I am currently using a model where I my parameter is a normal distribution whose centre is given by an array. With PyMC, I define it as :
xp=pymc.Normal('xp',mu=x,tau=0.1)
However, I get a few errors when I try doing this in PyMC3. Thus, I considered the error term to be a parameter instead :
e=pymc3.Normal('e',mu=0,tau=0.1)
and used (x+e) in my calculations.

In my model, the likelihood function has an integral whose limit is (x+e).
In PyMC, I defined the centre (mu) of the likelihood as :

@pymc.deterministic
def mu(om=om,xp=xp):
    f=lambda y:(1+om*((1+y)**float(3)-1))**float(-0.5)
    numpy.vectorize(f)
    def f1(xp):
        return integrate.quad(f,0,xp)
    f1=numpy.vectorize(f1)
    f2=f1(xp)
    return (f2[0]*(1+xp))

However, this didn’t really work with PyMC3. I also tried using a Custom theano Op (from https://stackoverflow.com/questions/42678490/custom-theano-op-to-do-numerical-integration/43154053). My code currently looks like this:

basic_model = pymc3.Model()

with basic_model:
    om=pymc3.Uniform('omega_m',0,1)

    a=pymc3.Uniform('a',0.01,0.1)

    b=pymc3.Uniform('b',0.001,0.01)

    e_r=pymc3.Normal("e_r",mu=0,sd=a+b*x)
    t = T.dscalar('t')
    z= (1+om*((1+t)**float(3)-1))**float(-0.5)
    numpy.vectorize(z)
    def f1(xp):
        intZ = integrateOut(z,t,0,xp)()

        funcIntZ = theano.function([],intZ)
        return funcIntZ()
    f1=numpy.vectorize(f1)
    f2=f1(x+e)
    mu=f2[0]*(1+x+e)
    y=pymc3.Normal('y',mu=mu,sd=0.5,observed=y_obs)

This is full traceback of the error:
Traceback (most recent call last):

  File "<ipython-input-9-4264bc40650f>", line 12, in <module>
    z= (1+om*((1+t)**float(3)-1))**float(-0.5)

  File "/anaconda3/lib/python3.6/site-packages/theano/tensor/var.py", line 227, in __radd__
    return theano.tensor.basic.add(other, self)

  File "/anaconda3/lib/python3.6/site-packages/theano/gof/op.py", line 639, in __call__
    (i, ins, node, detailed_err_msg))

ValueError: Cannot compute test value: input 1 (t) of Op Elemwise{add,no_inplace}(TensorConstant{1}, t) missing default value.  
Backtrace when that variable is created:

  File "/anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py", line 233, in dispatch_shell
    handler(stream, idents, msg)
  File "/anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py", line 399, in execute_request
    user_expressions, allow_stdin)
  File "/anaconda3/lib/python3.6/site-packages/ipykernel/ipkernel.py", line 208, in do_execute
    res = shell.run_cell(code, store_history=store_history, silent=silent)
  File "/anaconda3/lib/python3.6/site-packages/ipykernel/zmqshell.py", line 537, in run_cell
    return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)
  File "/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2728, in run_cell
    interactivity=interactivity, compiler=compiler, result=result)
  File "/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2850, in run_ast_nodes
    if self.run_code(code, result):
  File "/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2910, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-9-4264bc40650f>", line 11, in <module>
    t = T.dscalar('t')

Is there any other way of achieving this? I have tried various methods, but they don’t seem to work when I include a parameter as a limit of integration.

2

Try rewrite the function f into a theano function, remove the numpy.vectorize, and use the integration theano.Ops from below.

3

@junpenglao, thanks for your reply. I modified my code as per your suggestions. However, I got this error when I tried doing so :

AsTensorError: ('Cannot convert <function mu at 0x1c2ab42ea0> to TensorType', <class 'function'>)

I tried modifying the code further by converting the return value of the function to a tensor. This is the modified code :

basic_model = pymc3.Model()

with basic_model:
    start=pymc3.Uniform('start',0,0)
    a=pymc3.Uniform('a',0.01,0.1)

    b=pymc3.Uniform('b',0.001,0.01)

    e=pymc3.Normal("e",mu=0,sd=a+b*x)
    def mu(e):
        t = tt.dscalar('t')
        t.tag.test_value = np.zeros(())
        z=(1+0.2*((1+t)**float(3)-1))**float(-0.5)
        def f1(xp):
            integrate = Integrate(z,t)

            val=integrate(start, x+e)
            return val
        f2=f1(x+e)
        return tf.convert_to_tensor((f2*(1+x+e)))
    y=pymc3.Normal('y',mu=mu,sd=0.5,observed=y_obs)

But I still get the same error. This is the full traceback :

Traceback (most recent call last):

  File "<ipython-input-10-b47f07234c39>", line 21, in <module>
    y=pymc3.Normal('y',mu=mu,sd=0.5,observed=y_obs)

  File "/anaconda3/lib/python3.6/site-packages/pymc3/distributions/distribution.py", line 38, in __new__
    dist = cls.dist(*args, **kwargs)

  File "/anaconda3/lib/python3.6/site-packages/pymc3/distributions/distribution.py", line 49, in dist
    dist.__init__(*args, **kwargs)

  File "/anaconda3/lib/python3.6/site-packages/pymc3/distributions/continuous.py", line 386, in __init__
    self.mean = self.median = self.mode = self.mu = mu = tt.as_tensor_variable(mu)

  File "/anaconda3/lib/python3.6/site-packages/theano/tensor/basic.py", line 200, in as_tensor_variable
    raise AsTensorError("Cannot convert %s to TensorType" % str_x, type(x))

AsTensorError: ('Cannot convert <function mu at 0x1c2ab42ea0> to TensorType', <class 'function'>)

Is there any alternate way to define this function? I also tried using tf.py_func(), but this too doesn’t seem to work. Thanks in advance.

4

Oh I see, the end point of the Integral is an array, that’s why it’s not working.
You will need to rewrite the IntegrateOps into something that accept vector input for start point and end point, alternatively, try wrapping it in a for-loop:

x = np.arange(100)
y_obs = np.random.randn(100)

start = theano.shared(0.)
with pm.Model() as basic_model:
    a = pm.Uniform('a', 0.01, 0.1)
    b = pm.Uniform('b', 0.001, 0.01)
    xp = pm.Normal('xp', mu=x, sd=a + b * x, shape=100)

    t = tt.dscalar('t')
    t.tag.test_value = np.zeros(())
    func = (1 + 0.2 * ((1 + t)**3. - 1))**(-0.5)
    integrate = Integrate(func, t)
    
    mu = tt.stack([integrate(start, xp[i]) for i in range(100)])

    y = pm.Normal('y', mu=mu, sd=0.5, observed=y_obs)
5

@junpenglao, thanks a lot for your help! The code runs fine if the number of elements in x/y_obs is small. However, when I increase the number to more than 200, I get the following message :

 fatal error: bracket nesting level exceeded maximum of 256.         ]

I found this problem discussed already (https://github.com/pymc-devs/pymc3/issues/624), and increased the bracket nesting level. However, this has slowed down the program by a huge factor.

Should I try rewriting the IntegrateOps to accept vector inputs, would this make the program faster? Or are there any alternatives to bypass this problem? Thank you again for your help.

6

I would imagine that’s the best solution. make sure you validate the gradient like in @aseyboldt’s reply.

7

I’ll try doing it sometime this week, and will post it here if I’m successful. Thank you for your help!

1 Like