Use of scipy functions in pyMC3

Questions
1

Hi everyone.
I’ve struggle to use PyMC3 but got stuck in a problem.
A likelihood function I want to define involves the exponential integral function (expi) imported from scipy library.
But it doesn’t work when I defined the model in PyMC3 and the I see an error message as below

TypeError: ufunc ‘expi’ not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ‘‘safe’’

How can I use some special functions such as exponential integral and Bessel functions in PyMC3?

Please find the attached python script I am tackling with.
Thank you in advance
trialExpInt.py (751 Bytes)

2

You need to either find a implementation in theano, or wrap the scipy function in theano using @as_op.

3

Thank you for your answer.
Actually, I’ve already tried the decorator @as_op
However, It showed me an error message as below.

AttributeError: ‘int’ object has no attribute ‘type’

Is here something I am missing?

# In[1]:

import os
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
from scipy.special import expi
get_ipython().magic('matplotlib inline')


# In[2]:

Cs = 2.5E6
rb = 0.0825
Leff = 1.9
Rb = 0.162
q = 50
t = np.array([i*3600 for i in range(51)])
temp_mea = -(q/(4.0*Leff))*expi(-Cs*(rb**2)/(4.0*t*Leff)) + 6


# In[3]:

import theano
import theano.tensor as tt

from theano.compile.ops import as_op

@as_op(itypes=[tt.dscalar, tt.dvector, tt.dscalar, tt.dscalar, tt.dscalar],
       otypes=[tt.dvector])

def ils(q, t, Cs, rb, Leff):
    temp_bore = -(q/(4.0*Leff))*expi(-Cs*(rb**2)/(4.0*t*Leff))
    return temp_bore


#  In[4]:

basic_model = pm.Model()

with basic_model:

    Leff = pm.Normal('Leff', mu=1.8, sd=0.5)
    Rb = pm.Normal('Rb', mu=0.16, sd=0.1) 

    temp_cal = ils(q, t, Cs, rb, Leff) + Rb*q

    Y_obs = pm.Normal('Y_obs', mu=temp_cal, sd=0.4, observed=temp_mea)
4

The input to the function need to have the same type as you specified in the decorator, something like this should work:

# Make sure the input being shared are float
q = 50.
t = np.array([i*3600 for i in range(51)], dtype=np.float64)
...
q_shared, t_shared, Cs_shared, rb_shared = theano.shared(q), theano.shared(t), theano.shared(Cs), theano.shared(rb)
with basic_model:
    ...
    temp_cal = ils(q_shared, t_shared, Cs_shared, rb_shared, Leff) + Rb*q
    ...