Pymc3 likelihood math with non-theano function

Questions
1

I’m new to doing Bayesian inference. I’m trying to adapt a grid search code I wrote to Bayesian Monte Carlo Markov Chain approach, and I’m using PyMC3. My problem is that the code has to call a function that can’t be rewritten in theano syntax. (The function relies on a piece of Fortran code in an f2py wrapper.)

Here’s the code I’m working with:

with pm.Model() as model:
    # Independent parameters
    x = pm.Normal('x', sx, sd=incx*float(nrangex)/2.0)
    y = pm.Normal('y', sy, sd=incy*float(nrangey)/2.0)
    depth = pm.Normal('depth', sdepth, sd=incdepth*float(nrangedepth)/2.0)
    length = pm.Normal('length', slength, sd=inclength*float(nrangelength)/2.0)
    width = pm.Normal('width', swidth, sd=incwidth*float(nrangewidth)/2.0)
    strike = pm.Normal('strike', sstrike, sd=incstrike*float(nrangestrike)/2.0)
    dip = pm.Normal('dip', sdip, sd=incdip*float(nrangedip)/2.0)
    rake = pm.Normal('rake', srake, sd=incrake*float(nrangerake)/2.0)

    # Model (This is the part that doesn't work.)
    los_disp =  zne2los(getdisp(lon2km(dsidata['lon'], x), lat2km(dsidata['lat'], y), depth, length, width, strike, dip, rake))

    # Likelihood
    odisp = pm.Normal('los_disp', mu = los_disp, sd = 0.5, observed = dsidata['disp'])

    # Sampling
    trace = pm.sample(100)

What this code is trying to do is invert for earthquake source parameters from ground displacement data. The dsidata data frame contains ground displacement data as a function of latitude and longitude. I’m trying to solve for the eight earthquake source parameters that produce the best fit to this ground surface displacement.

The getdisp function simply cannot be rewritten for theano because it calls a piece of Fortran that forward models ground surface displacement from earthquake source parameters. Is there a way to compile non-theano code into a theano form? Is there another way do this?

Since I’m new to this, and I can’t find many great examples to look at, there may well be other errors in the code. Are there other mistakes I’m making?

2

Did u try to sample with metropolis?

3

My problem is not sampling. It’s building the likelihood function. Metropolis sampling would solve the problem of not having a gradient. I feel like the answer has something to do with op decorators, but I can’t find enough documentation online to figure them out.

4

I see, you can wrap it as a theano operation, you can see in the Doc and Example.

5

I implemented exactly that problem here in beat:

Let me know if you want to try to get it running. If you want to continue using the Okada half-space you can look into the theanof.py code there to follow how to wrap an external fortran code.

In case it is not too late by now :wink: and you solved it anyways…
Cheers!