How to get scipy.special functions to take pymc3 RV's as input

Question

I'm trying to use pymc3 to fit a model involving the voigt function (coming from scipy.special). The inputs to the voigt function are supposed to be arrays, whereas a,b are pymc3 classes. How do I get scipy.special functions to take pymc3 RV's as input? Running the code attached below produces an error:

import pymc3 as pm
from scipy.special import voigt_profile
import numpy as np

with pm.Model() as linear_model:
    a = pm.Lognormal('a',mu=0, sigma=2.)
    b = pm.Lognormal('b',mu=0, sigma=2.)
    x = np.linspace(-1,1)
    c = voigt_profile(x,a,b)

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

Answer 1

For better or for worse, you'll need to (re)implement the function using theano. Here is one naive version that works: notice that you can not use erfc, because theano errors out.

import theano.tensor as tt

def faddeeva(z):
    m = tt.exp(-tt.square(z))
    return (m - m * tt.erf(1.j * z))

def voigt_profile(x, sigma, gamma):
    z = (x + 1.j * gamma) / (tt.sqrt(2.) * sigma)
    return faddeeva(z).real / (sigma * tt.sqrt(2 * np.pi))