Using a “black box” likelihood function

In [1]:
%matplotlib inline
%load_ext Cython

%matplotlib inline

import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-darkgrid')
print('Running on PyMC3 v{}'.format(pm.__version__))

# for reproducibility here's some version info for modules used in this notebook
import theano
import theano.tensor as tt
import platform
import cython
import IPython
import matplotlib
import emcee
import corner
import os
print("Python version:     {}".format(platform.python_version()))
print("IPython version:    {}".format(IPython.__version__))
print("Cython version:     {}".format(cython.__version__))
print("GSL version:        {}".format(os.popen('gsl-config --version').read().strip()))
print("Numpy version:      {}".format(np.__version__))
print("Theano version:     {}".format(theano.__version__))
print("Matplotlib version: {}".format(matplotlib.__version__))
print("emcee version:      {}".format(emcee.__version__))
print("corner version:     {}".format(corner.__version__))
Running on PyMC3 v3.5
Python version:     3.6.6
IPython version:    6.5.0
Cython version:     0.28.5
GSL version:        2.1
Numpy version:      1.15.1
Theano version:     1.0.2
Matplotlib version: 2.2.3
emcee version:      2.2.1
corner version:     2.0.1

PyMC3 is a great tool for doing Bayesian inference and parameter estimation. It has a load of in-built probability distributions that you can use to set up priors and likelihood functions for your particular model. You can even create your own custom distributions.

However, this is not necessarily that simple if you have a model function, or probability distribution, that, for example, relies on an external code that you have little/no control over (and may even be, for example, wrapped C code rather than Python). This can be problematic went you need to pass parameters set as PyMC3 distributions to these external functions; your external function probably wants you to pass it floating point numbers rather than PyMC3 distributions!

import pymc3 as pm:
from external_module import my_external_func  # your external function!

# set up your model
with pm.Model():
    # your external function takes two parameters, a and b, with Uniform priors
    a = pm.Uniform('a', lower=0., upper=1.)
    b = pm.Uniform('b', lower=0., upper=1.)

    m = my_external_func(a, b)  # <--- this is not going to work!

Another issue is that if you want to be able to use the gradient-based step samplers like NUTS and Hamiltonian Monte Carlo (HMC), then your model/likelihood needs a gradient to be defined. If you have a model that is defined as a set of Theano operators then this is no problem - internally it will be able to do automatic differentiation - but if your model is essentially a “black box” then you won’t necessarily know what the gradients are.

Defining a model/likelihood that PyMC3 can use and that calls your “black box” function is possible, but it relies on creating a custom Theano Op. This is, hopefully, a clear description of how to do this, including one way of writing a gradient function that could be generally applicable.

In the examples below, we create a very simple model and log-likelihood function in Cython. Cython is used just as an example to show what you might need to do if calling external C codes, but you could in fact be using pure Python codes. The log-likelihood function used is actually just a Normal distribution, so defining this yourself is obviously overkill (and I’ll compare it to doing the same thing purely with the pre-defined PyMC3 Normal distribution), but it should provide a simple to follow demonstration.

First, let’s define a super-complicated™ model (a straight line!), which is parameterised by two variables (a gradient m and a y-intercept c) and calculated at a vector of points x. Here the model is defined in Cython and calls GSL functions. This is just to show that you could be calling some other C library that you need. In this example, the model parameters are all packed into a list/array/tuple called theta.

Let’s also define a really-complicated™ log-likelihood function (a Normal log-likelihood that ignores the normalisation), which takes in the list/array/tuple of model parameter values theta, the points at which to calculate the model x, the vector of “observed” data points data, and the standard deviation of the noise in the data sigma. This log-likelihood function calls the super-complicated™ model function.

In [2]:
%%cython -I/usr/include -L/usr/lib/x86_64-linux-gnu -lgsl -lgslcblas -lm

import cython
cimport cython

import numpy as np
cimport numpy as np

### STUFF FOR USING GSL (FEEL FREE TO IGNORE!) ###

# declare GSL vector structure and functions
cdef extern from "gsl/gsl_block.h":
    cdef struct gsl_block:
        size_t size
        double * data

cdef extern from "gsl/gsl_vector.h":
    cdef struct gsl_vector:
        size_t size
        size_t stride
        double * data
        gsl_block * block
        int owner

    ctypedef struct gsl_vector_view:
        gsl_vector vector

    int gsl_vector_scale (gsl_vector * a, const double x) nogil
    int gsl_vector_add_constant (gsl_vector * a, const double x) nogil
    gsl_vector_view gsl_vector_view_array (double * base, size_t n) nogil

###################################################


# define your super-complicated model that uses loads of external codes
cpdef my_model(theta, np.ndarray[np.float64_t, ndim=1] x):
    """
    A straight line!

    Note:
        This function could simply be:

            m, c = thetha
            return m*x + x

        but I've made it more complicated for demonstration purposes
    """
    m, c = theta  # unpack line gradient and y-intercept

    cdef size_t length = len(x)  # length of x

    cdef np.ndarray line = np.copy(x)  # make copy of x vector
    cdef gsl_vector_view lineview      # create a view of the vector
    lineview = gsl_vector_view_array(<double *>line.data, length)

    # multiply x by m
    gsl_vector_scale(&lineview.vector, <double>m)

    # add c
    gsl_vector_add_constant(&lineview.vector, <double>c)

    # return the numpy array
    return line


# define your really-complicated likelihood function that uses loads of external codes
cpdef my_loglike(theta, np.ndarray[np.float64_t, ndim=1] x,
                 np.ndarray[np.float64_t, ndim=1] data, sigma):
    """
    A Gaussian log-likelihood function for a model with parameters given in theta
    """

    model = my_model(theta, x)

    return -(0.5/sigma**2)*np.sum((data - model)**2)

Now, as things are, if we wanted to sample from this log-likelihood function, using certain prior distributions for the model parameters (gradient and y-intercept) using PyMC3, we might try something like this (using a PyMC3 DensityDist):

import pymc3 as pm

# create/read in our "data" (I'll show this in the real example below)
x = ...
sigma = ...
data = ...

with pm.Model():
    # set priors on model gradient and y-intercept
    m = pm.Uniform('m', lower=-10., upper=10.)
    c = pm.Uniform('c', lower=-10., upper=10.)

    # create custom distribution
    pm.DensityDist('likelihood', my_loglike,
                   observed={'theta': (m, c), 'x': x, 'data': data, 'sigma': sigma})

    # sample from the distribution
    trace = pm.sample(1000)

But, this will give an error like:

ValueError: setting an array element with a sequence.

This is because m and c are Theano tensor-type objects.

So, what we actually need to do is create a Theano Op. This will be a new class that wraps our log-likelihood function (or just our model function, if that is all that is required) into something that can take in Theano tensor objects, but internally can cast them as floating point values that can be passed to our log-likelihood function. We will do this below, initially without defining a grad() method for the Op.

In [3]:
# define a theano Op for our likelihood function
class LogLike(tt.Op):

    """
    Specify what type of object will be passed and returned to the Op when it is
    called. In our case we will be passing it a vector of values (the parameters
    that define our model) and returning a single "scalar" value (the
    log-likelihood)
    """
    itypes = [tt.dvector] # expects a vector of parameter values when called
    otypes = [tt.dscalar] # outputs a single scalar value (the log likelihood)

    def __init__(self, loglike, data, x, sigma):
        """
        Initialise the Op with various things that our log-likelihood function
        requires. Below are the things that are needed in this particular
        example.

        Parameters
        ----------
        loglike:
            The log-likelihood (or whatever) function we've defined
        data:
            The "observed" data that our log-likelihood function takes in
        x:
            The dependent variable (aka 'x') that our model requires
        sigma:
            The noise standard deviation that our function requires.
        """

        # add inputs as class attributes
        self.likelihood = loglike
        self.data = data
        self.x = x
        self.sigma = sigma

    def perform(self, node, inputs, outputs):
        # the method that is used when calling the Op
        theta, = inputs  # this will contain my variables

        # call the log-likelihood function
        logl = self.likelihood(theta, self.x, self.data, self.sigma)

        outputs[0][0] = np.array(logl) # output the log-likelihood

Now, let’s use this Op to repeat the example shown above. To do this let’s create some data containing a straight line with additive Gaussian noise (with a mean of zero and a standard deviation of sigma). For simplicity we set uniform prior distributions on the gradient and y-intercept. As we’ve not set the grad() method of the Op PyMC3 will not be able to use the gradient-based samplers, so will fall back to using the Slice sampler.

In [4]:
# set up our data
N = 10  # number of data points
sigma = 1.  # standard deviation of noise
x = np.linspace(0., 9., N)

mtrue = 0.4  # true gradient
ctrue = 3.   # true y-intercept

truemodel = my_model([mtrue, ctrue], x)

# make data
np.random.seed(716742)  # set random seed, so the data is reproducible each time
data = sigma*np.random.randn(N) + truemodel

ndraws = 3000  # number of draws from the distribution
nburn = 1000   # number of "burn-in points" (which we'll discard)

# create our Op
logl = LogLike(my_loglike, data, x, sigma)

# use PyMC3 to sampler from log-likelihood
with pm.Model():
    # uniform priors on m and c
    m = pm.Uniform('m', lower=-10., upper=10.)
    c = pm.Uniform('c', lower=-10., upper=10.)

    # convert m and c to a tensor vector
    theta = tt.as_tensor_variable([m, c])

    # use a DensityDist (use a lamdba function to "call" the Op)
    pm.DensityDist('likelihood', lambda v: logl(v), observed={'v': theta})

    trace = pm.sample(ndraws, tune=nburn, discard_tuned_samples=True)

# plot the traces
_ = pm.traceplot(trace, lines={'m': mtrue, 'c': ctrue})

# put the chains in an array (for later!)
samples_pymc3 = np.vstack((trace['m'], trace['c'])).T
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Initializing NUTS failed. Falling back to elementwise auto-assignment.
Multiprocess sampling (2 chains in 2 jobs)
CompoundStep
>Slice: [c]
>Slice: [m]
Sampling 2 chains: 100%|██████████| 8000/8000 [00:04<00:00, 1928.75draws/s]
The number of effective samples is smaller than 25% for some parameters.
../_images/notebooks_blackbox_external_likelihood_7_1.png

What if we wanted to use NUTS or HMC? If we knew the analytical derivatives of the model/likelihood function then we could add a grad() method to the Op using that analytical form.

But, what if we don’t know the analytical form. If our model/likelihood is purely Python and made up of standard maths operators and Numpy functions, then the autograd module could potentially be used to find gradients (also, see here for a nice Python example of automatic differentiation). But, if our model/likelihood truely is a “black box” then we can just use the good-old-fashioned finite difference to find the gradients - this can be slow, especially if there are a large number of variables, or the model takes a long time to evaluate. Below, a function to find gradients has been defined that uses the finite difference (the central difference) - it uses an iterative method with successively smaller interval sizes to check that the gradient converges. But, you could do something far simpler and just use, for example, the SciPy approx_fprime function. Here, the gradient function is defined in Cython for speed, but if the function it evaluates to find the gradients is the performance bottle neck then having this as a pure Python function may not make a significant speed difference.

In [5]:
%%cython

import cython
cimport cython

import numpy as np
cimport numpy as np

import warnings

def gradients(vals, func, releps=1e-3, abseps=None, mineps=1e-9, reltol=1e-3,
              epsscale=0.5):
    """
    Calculate the partial derivatives of a function at a set of values. The
    derivatives are calculated using the central difference, using an iterative
    method to check that the values converge as step size decreases.

    Parameters
    ----------
    vals: array_like
        A set of values, that are passed to a function, at which to calculate
        the gradient of that function
    func:
        A function that takes in an array of values.
    releps: float, array_like, 1e-3
        The initial relative step size for calculating the derivative.
    abseps: float, array_like, None
        The initial absolute step size for calculating the derivative.
        This overrides `releps` if set.
        `releps` is set then that is used.
    mineps: float, 1e-9
        The minimum relative step size at which to stop iterations if no
        convergence is achieved.
    epsscale: float, 0.5
        The factor by which releps if scaled in each iteration.

    Returns
    -------
    grads: array_like
        An array of gradients for each non-fixed value.
    """

    grads = np.zeros(len(vals))

    # maximum number of times the gradient can change sign
    flipflopmax = 10.

    # set steps
    if abseps is None:
        if isinstance(releps, float):
            eps = np.abs(vals)*releps
            eps[eps == 0.] = releps  # if any values are zero set eps to releps
            teps = releps*np.ones(len(vals))
        elif isinstance(releps, (list, np.ndarray)):
            if len(releps) != len(vals):
                raise ValueError("Problem with input relative step sizes")
            eps = np.multiply(np.abs(vals), releps)
            eps[eps == 0.] = np.array(releps)[eps == 0.]
            teps = releps
        else:
            raise RuntimeError("Relative step sizes are not a recognised type!")
    else:
        if isinstance(abseps, float):
            eps = abseps*np.ones(len(vals))
        elif isinstance(abseps, (list, np.ndarray)):
            if len(abseps) != len(vals):
                raise ValueError("Problem with input absolute step sizes")
            eps = np.array(abseps)
        else:
            raise RuntimeError("Absolute step sizes are not a recognised type!")
        teps = eps

    # for each value in vals calculate the gradient
    count = 0
    for i in range(len(vals)):
        # initial parameter diffs
        leps = eps[i]
        cureps = teps[i]

        flipflop = 0

        # get central finite difference
        fvals = np.copy(vals)
        bvals = np.copy(vals)

        # central difference
        fvals[i] += 0.5*leps  # change forwards distance to half eps
        bvals[i] -= 0.5*leps  # change backwards distance to half eps
        cdiff = (func(fvals)-func(bvals))/leps

        while 1:
            fvals[i] -= 0.5*leps  # remove old step
            bvals[i] += 0.5*leps

            # change the difference by a factor of two
            cureps *= epsscale
            if cureps < mineps or flipflop > flipflopmax:
                # if no convergence set flat derivative (TODO: check if there is a better thing to do instead)
                warnings.warn("Derivative calculation did not converge: setting flat derivative.")
                grads[count] = 0.
                break
            leps *= epsscale

            # central difference
            fvals[i] += 0.5*leps  # change forwards distance to half eps
            bvals[i] -= 0.5*leps  # change backwards distance to half eps
            cdiffnew = (func(fvals)-func(bvals))/leps

            if cdiffnew == cdiff:
                grads[count] = cdiff
                break

            # check whether previous diff and current diff are the same within reltol
            rat = (cdiff/cdiffnew)
            if np.isfinite(rat) and rat > 0.:
                # gradient has not changed sign
                if np.abs(1.-rat) < reltol:
                    grads[count] = cdiffnew
                    break
                else:
                    cdiff = cdiffnew
                    continue
            else:
                cdiff = cdiffnew
                flipflop += 1
                continue

        count += 1

    return grads

So, now we can just redefine our Op with a grad() method, right?

It’s not quite so simple! The grad() method itself requires that its inputs are Theano tensor variables, whereas our gradients function above, like our my_loglike function, wants a list of floating point values. So, we need to define another Op that calculates the gradients. Below, I define a new version of the LogLike Op, called LogLikeWithGrad this time, that has a grad() method. This is followed by anothor Op called LogLikeGrad that, when called with a vector of Theano tensor variables, returns another vector of values that are the gradients (i.e., the Jacobian) of our log-likelihood function at those values. Note that the grad() method itself does not return the gradients directly, but instead returns the Jacobian-vector product (you can hopefully just copy what I’ve done and not worry about what this means too much!).

In [6]:
# define a theano Op for our likelihood function
class LogLikeWithGrad(tt.Op):

    itypes = [tt.dvector] # expects a vector of parameter values when called
    otypes = [tt.dscalar] # outputs a single scalar value (the log likelihood)

    def __init__(self, loglike, data, x, sigma):
        """
        Initialise with various things that the function requires. Below
        are the things that are needed in this particular example.

        Parameters
        ----------
        loglike:
            The log-likelihood (or whatever) function we've defined
        data:
            The "observed" data that our log-likelihood function takes in
        x:
            The dependent variable (aka 'x') that our model requires
        sigma:
            The noise standard deviation that out function requires.
        """

        # add inputs as class attributes
        self.likelihood = loglike
        self.data = data
        self.x = x
        self.sigma = sigma

        # initialise the gradient Op (below)
        self.logpgrad = LogLikeGrad(self.likelihood, self.data, self.x, self.sigma)

    def perform(self, node, inputs, outputs):
        # the method that is used when calling the Op
        theta, = inputs  # this will contain my variables

        # call the log-likelihood function
        logl = self.likelihood(theta, self.x, self.data, self.sigma)

        outputs[0][0] = np.array(logl) # output the log-likelihood

    def grad(self, inputs, g):
        # the method that calculates the gradients - it actually returns the
        # vector-Jacobian product - g[0] is a vector of parameter values
        theta, = inputs  # our parameters
        return [g[0]*self.logpgrad(theta)]


class LogLikeGrad(tt.Op):

    """
    This Op will be called with a vector of values and also return a vector of
    values - the gradients in each dimension.
    """
    itypes = [tt.dvector]
    otypes = [tt.dvector]

    def __init__(self, loglike, data, x, sigma):
        """
        Initialise with various things that the function requires. Below
        are the things that are needed in this particular example.

        Parameters
        ----------
        loglike:
            The log-likelihood (or whatever) function we've defined
        data:
            The "observed" data that our log-likelihood function takes in
        x:
            The dependent variable (aka 'x') that our model requires
        sigma:
            The noise standard deviation that out function requires.
        """

        # add inputs as class attributes
        self.likelihood = loglike
        self.data = data
        self.x = x
        self.sigma = sigma

    def perform(self, node, inputs, outputs):
        theta, = inputs

        # define version of likelihood function to pass to derivative function
        def lnlike(values):
            return self.likelihood(values, self.x, self.data, self.sigma)

        # calculate gradients
        grads = gradients(theta, lnlike)

        outputs[0][0] = grads

Now, let’s re-run PyMC3 with our new “grad”-ed Op. This time it will be able to automatically use NUTS.

In [7]:
# create our Op
logl = LogLikeWithGrad(my_loglike, data, x, sigma)

# use PyMC3 to sampler from log-likelihood
with pm.Model() as opmodel:
    # uniform priors on m and c
    m = pm.Uniform('m', lower=-10., upper=10.)
    c = pm.Uniform('c', lower=-10., upper=10.)

    # convert m and c to a tensor vector
    theta = tt.as_tensor_variable([m, c])

    # use a DensityDist
    pm.DensityDist('likelihood', lambda v: logl(v), observed={'v': theta})

    trace = pm.sample(ndraws, tune=nburn, discard_tuned_samples=True)

# plot the traces
_ = pm.traceplot(trace, lines={'m': mtrue, 'c': ctrue})

# put the chains in an array (for later!)
samples_pymc3_2 = np.vstack((trace['m'], trace['c'])).T
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [c, m]
Sampling 2 chains: 100%|██████████| 8000/8000 [00:08<00:00, 898.88draws/s]
../_images/notebooks_blackbox_external_likelihood_13_1.png

Now, finally, just to check things actually worked as we might expect, let’s do the same thing purely using PyMC3 distributions (because in this simple example we can!)

In [8]:
with pm.Model() as pymodel:
    # uniform priors on m and c
    m = pm.Uniform('m', lower=-10., upper=10.)
    c = pm.Uniform('c', lower=-10., upper=10.)

    # convert m and c to a tensor vector
    theta = tt.as_tensor_variable([m, c])

    # use a Normal distribution
    pm.Normal('likelihood', mu=(m*x + c), sd = sigma, observed=data)

    trace = pm.sample(ndraws, tune=nburn, discard_tuned_samples=True)

# plot the traces
_ = pm.traceplot(trace, lines={'m': mtrue, 'c': ctrue})

# put the chains in an array (for later!)
samples_pymc3_3 = np.vstack((trace['m'], trace['c'])).T
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [c, m]
Sampling 2 chains: 100%|██████████| 8000/8000 [00:03<00:00, 2049.77draws/s]
../_images/notebooks_blackbox_external_likelihood_15_1.png

To check that they match let’s plot all the examples together and also find the autocorrelation lengths.

In [9]:
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)  # supress emcee autocorr FutureWarning

matplotlib.rcParams['font.size'] = 18

hist2dkwargs = {'plot_datapoints': False,
                'plot_density': False,
                'levels': 1.0 - np.exp(-0.5 * np.arange(1.5, 2.1, 0.5) ** 2)} # roughly 1 and 2 sigma

colors = ['r', 'g', 'b']
labels = ['Theanp Op (no grad)', 'Theano Op (with grad)', 'Pure PyMC3']

for i, samples in enumerate([samples_pymc3, samples_pymc3_2, samples_pymc3_3]):
    # get maximum chain autocorrelartion length
    autocorrlen = int(np.max(emcee.autocorr.integrated_time(samples, c=3)));
    print('Auto-correlation length ({}): {}'.format(labels[i], autocorrlen))

    if i == 0:
        fig = corner.corner(samples, labels=[r"$m$", r"$c$"], color=colors[i],
                            hist_kwargs={'density': True}, **hist2dkwargs,
                            truths=[mtrue, ctrue])
    else:
        corner.corner(samples, color=colors[i], hist_kwargs={'density': True},
                      fig=fig, **hist2dkwargs)

fig.set_size_inches(8, 8)
Auto-correlation length (Theanp Op (no grad)): 6
Auto-correlation length (Theano Op (with grad)): 3
Auto-correlation length (Pure PyMC3): 3
../_images/notebooks_blackbox_external_likelihood_17_1.png

We can now check that the gradient Op works was we expect it to. First, just create and call the LogLikeGrad class, which should return the gradient directly (note that we have to create a Theano function to convert the output of the Op to an array). Secondly, we call the gradient from LogLikeWithGrad by using the Theano tensor gradient function. Finally, we will check the gradient returned by the PyMC3 model for a Normal distribution, which should be the same as the log-likelihood function we defined. In all cases we evaluate the gradients at the true values of the model function (the straight line) that was created.

In [10]:
# test the gradient Op by direct call
theano.config.compute_test_value = "ignore"
theano.config.exception_verbosity = "high"

var = tt.dvector()
test_grad_op = LogLikeGrad(my_loglike, data, x, sigma)
test_grad_op_func = theano.function([var], test_grad_op(var))
grad_vals = test_grad_op_func([mtrue, ctrue])

print('Gradient returned by "LogLikeGrad": {}'.format(grad_vals))

# test the gradient called through LogLikeWithGrad
test_gradded_op = LogLikeWithGrad(my_loglike, data, x, sigma)
test_gradded_op_grad = tt.grad(test_gradded_op(var), var)
test_gradded_op_grad_func = theano.function([var], test_gradded_op_grad)
grad_vals_2 = test_gradded_op_grad_func([mtrue, ctrue])

print('Gradient returned by "LogLikeWithGrad": {}'.format(grad_vals_2))

# test the gradient that PyMC3 uses for the Normal log likelihood
test_model = pm.Model()
with test_model:
    m = pm.Uniform('m', lower=-10., upper=10.)
    c = pm.Uniform('c', lower=-10., upper=10.)

    pm.Normal('likelihood', mu=(m*x + c), sd=sigma, observed=data)

    gradfunc = test_model.logp_dlogp_function([m, c], dtype=None)
    gradfunc.set_extra_values({'m_interval__': mtrue, 'c_interval__': ctrue})
    grad_vals_pymc3 = gradfunc(np.array([mtrue, ctrue]))[1]  # get dlogp values

print('Gradient returned by PyMC3 "Normal" distribution: {}'.format(grad_vals_pymc3))
Gradient returned by "LogLikeGrad": [-7.17656625 -1.39486358]
Gradient returned by "LogLikeWithGrad": [-7.17656625 -1.39486358]
Gradient returned by PyMC3 "Normal" distribution: [-7.17656625 -1.39486358]

We can also do some profiling of the Op, as used within a PyMC3 Model, to check performance. First, we’ll profile using the LogLikeWithGrad Op, and then doing the same thing purely using PyMC3 distributions.

In [11]:
# profile logpt using our Op
opmodel.profile(opmodel.logpt).summary()
Function profiling
==================
  Message: /home/matthew/.conda/envs/testing/lib/python3.6/site-packages/pymc3/model.py:909
  Time in 1000 calls to Function.__call__: 4.123712e-02s
  Time in Function.fn.__call__: 2.523708e-02s (61.200%)
  Time in thunks: 2.199244e-02s (53.332%)
  Total compile time: 7.463312e-02s
    Number of Apply nodes: 8
    Theano Optimizer time: 5.898285e-02s
       Theano validate time: 5.080700e-04s
    Theano Linker time (includes C, CUDA code generation/compiling): 4.571915e-03s
       Import time 0.000000e+00s
       Node make_thunk time 4.115105e-03s
           Node Elemwise{Composite{((i0 + Switch(Cast{int8}((GE(i1, i2) * LE(i1, i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 1)](TensorConstant{2.995732273553991}, c, TensorConstant{-10.0}, TensorConstant{10.0}, TensorConstant{-2.995732273553991}, TensorConstant{-inf}, TensorConstant{2.0}, c_interval__) time 8.506775e-04s
           Node Elemwise{Composite{((i0 + Switch(Cast{int8}((GE(i1, i2) * LE(i1, i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 1)](TensorConstant{2.995732273553991}, m, TensorConstant{-10.0}, TensorConstant{10.0}, TensorConstant{-2.995732273553991}, TensorConstant{-inf}, TensorConstant{2.0}, m_interval__) time 8.454323e-04s
           Node Elemwise{Composite{(i0 + (i1 * scalar_sigmoid(i2)))}}(TensorConstant{-10.0}, TensorConstant{20.0}, c_interval__) time 6.031990e-04s
           Node Elemwise{Composite{(i0 + (i1 * scalar_sigmoid(i2)))}}(TensorConstant{-10.0}, TensorConstant{20.0}, m_interval__) time 5.083084e-04s
           Node <__main__.LogLikeWithGrad object at 0x7f1d606e8fd0>(MakeVector{dtype='float64'}.0) time 3.864765e-04s

Time in all call to theano.grad() 6.874464e-01s
Time since theano import 43.624s
Class
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>
  91.8%    91.8%       0.020s       2.02e-05s     Py    1000       1   __main__.LogLikeWithGrad
   4.5%    96.3%       0.001s       2.48e-07s     C     4000       4   theano.tensor.elemwise.Elemwise
   2.8%    99.1%       0.001s       3.04e-07s     C     2000       2   theano.tensor.opt.MakeVector
   0.9%   100.0%       0.000s       2.06e-07s     C     1000       1   theano.tensor.elemwise.Sum
   ... (remaining 0 Classes account for   0.00%(0.00s) of the runtime)

Ops
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>
  91.8%    91.8%       0.020s       2.02e-05s     Py    1000        1   <__main__.LogLikeWithGrad object at 0x7f1d606e8fd0>
   2.8%    94.6%       0.001s       3.04e-07s     C     2000        2   MakeVector{dtype='float64'}
   2.5%    97.1%       0.001s       2.79e-07s     C     2000        2   Elemwise{Composite{(i0 + (i1 * scalar_sigmoid(i2)))}}
   2.0%    99.1%       0.000s       2.16e-07s     C     2000        2   Elemwise{Composite{((i0 + Switch(Cast{int8}((GE(i1, i2) * LE(i1, i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 1)]
   0.9%   100.0%       0.000s       2.06e-07s     C     1000        1   Sum{acc_dtype=float64}
   ... (remaining 0 Ops account for   0.00%(0.00s) of the runtime)

Apply
------
<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>
  91.8%    91.8%       0.020s       2.02e-05s   1000     5   <__main__.LogLikeWithGrad object at 0x7f1d606e8fd0>(MakeVector{dtype='float64'}.0)
   2.0%    93.8%       0.000s       4.42e-07s   1000     1   Elemwise{Composite{(i0 + (i1 * scalar_sigmoid(i2)))}}(TensorConstant{-10.0}, TensorConstant{20.0}, m_interval__)
   1.7%    95.5%       0.000s       3.73e-07s   1000     6   MakeVector{dtype='float64'}(__logp_m_interval__, __logp_c_interval__, __logp_likelihood)
   1.5%    97.0%       0.000s       3.30e-07s   1000     4   Elemwise{Composite{((i0 + Switch(Cast{int8}((GE(i1, i2) * LE(i1, i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 1)](TensorConstant{2.995732273553991}, m, TensorConstant{-10.0}, TensorConstant{10.0}, TensorConstant{-2.995732273553991}, TensorConstant{-inf}, TensorConstant{2.0}, m_interval__)
   1.1%    98.1%       0.000s       2.35e-07s   1000     2   MakeVector{dtype='float64'}(m, c)
   0.9%    99.0%       0.000s       2.06e-07s   1000     7   Sum{acc_dtype=float64}(MakeVector{dtype='float64'}.0)
   0.5%    99.5%       0.000s       1.16e-07s   1000     0   Elemwise{Composite{(i0 + (i1 * scalar_sigmoid(i2)))}}(TensorConstant{-10.0}, TensorConstant{20.0}, c_interval__)
   0.5%   100.0%       0.000s       1.03e-07s   1000     3   Elemwise{Composite{((i0 + Switch(Cast{int8}((GE(i1, i2) * LE(i1, i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 1)](TensorConstant{2.995732273553991}, c, TensorConstant{-10.0}, TensorConstant{10.0}, TensorConstant{-2.995732273553991}, TensorConstant{-inf}, TensorConstant{2.0}, c_interval__)
   ... (remaining 0 Apply instances account for 0.00%(0.00s) of the runtime)

Here are tips to potentially make your code run faster
                 (if you think of new ones, suggest them on the mailing list).
                 Test them first, as they are not guaranteed to always provide a speedup.
  - Try the Theano flag floatX=float32
  - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.
In [12]:
# profile using our PyMC3 distribution
pymodel.profile(pymodel.logpt).summary()
Function profiling
==================
  Message: /home/matthew/.conda/envs/testing/lib/python3.6/site-packages/pymc3/model.py:909
  Time in 1000 calls to Function.__call__: 1.982212e-02s
  Time in Function.fn.__call__: 5.885363e-03s (29.691%)
  Time in thunks: 2.546072e-03s (12.845%)
  Total compile time: 1.308653e-01s
    Number of Apply nodes: 11
    Theano Optimizer time: 1.083095e-01s
       Theano validate time: 8.940697e-04s
    Theano Linker time (includes C, CUDA code generation/compiling): 8.928537e-03s
       Import time 1.179457e-03s
       Node make_thunk time 8.301497e-03s
           Node Elemwise{Composite{(i0 * scalar_sigmoid(i1))}}(TensorConstant{20.0}, c_interval__) time 1.806974e-03s
           Node Elemwise{Composite{((i0 + Switch(Cast{int8}((GE((i1 + i2), i1) * LE((i1 + i2), i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 2)](TensorConstant{2.995732273553991}, TensorConstant{-10.0}, Elemwise{Composite{(i0 * scalar_sigmoid(i1))}}.0, TensorConstant{10.0}, TensorConstant{-2.995732273553991}, TensorConstant{-inf}, TensorConstant{2.0}, c_interval__) time 1.499176e-03s
           Node InplaceDimShuffle{x}(Elemwise{Composite{(i0 * scalar_sigmoid(i1))}}.0) time 8.523464e-04s
           Node Elemwise{Composite{((i0 + Switch(Cast{int8}((GE(i1, i2) * LE(i1, i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 1)](TensorConstant{2.995732273553991}, m, TensorConstant{-10.0}, TensorConstant{10.0}, TensorConstant{-2.995732273553991}, TensorConstant{-inf}, TensorConstant{2.0}, m_interval__) time 7.934570e-04s
           Node InplaceDimShuffle{x}(m) time 7.789135e-04s

Time in all call to theano.grad() 6.874464e-01s
Time since theano import 43.954s
Class
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>
  53.3%    53.3%       0.001s       2.26e-07s     C     6000       6   theano.tensor.elemwise.Elemwise
  22.9%    76.2%       0.001s       2.91e-07s     C     2000       2   theano.tensor.elemwise.DimShuffle
  14.1%    90.3%       0.000s       3.58e-07s     C     1000       1   theano.tensor.opt.MakeVector
   9.7%   100.0%       0.000s       1.23e-07s     C     2000       2   theano.tensor.elemwise.Sum
   ... (remaining 0 Classes account for   0.00%(0.00s) of the runtime)

Ops
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>
  22.9%    22.9%       0.001s       2.91e-07s     C     2000        2   InplaceDimShuffle{x}
  14.1%    37.0%       0.000s       3.58e-07s     C     1000        1   MakeVector{dtype='float64'}
  11.5%    48.5%       0.000s       2.93e-07s     C     1000        1   Elemwise{Composite{(i0 + (-sqr((i1 - ((i2 * i3) + i4)))))}}
  11.3%    59.8%       0.000s       2.88e-07s     C     1000        1   Elemwise{Composite{(i0 + (i1 * scalar_sigmoid(i2)))}}
  10.8%    70.6%       0.000s       2.75e-07s     C     1000        1   Elemwise{Composite{((i0 + Switch(Cast{int8}((GE((i1 + i2), i1) * LE((i1 + i2), i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 2)]
   9.7%    80.3%       0.000s       2.47e-07s     C     1000        1   Elemwise{Composite{((i0 + Switch(Cast{int8}((GE(i1, i2) * LE(i1, i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 1)]
   9.7%    90.0%       0.000s       1.23e-07s     C     2000        2   Sum{acc_dtype=float64}
   6.0%    96.0%       0.000s       1.52e-07s     C     1000        1   Elemwise{Composite{(i0 * scalar_sigmoid(i1))}}
   4.0%   100.0%       0.000s       1.03e-07s     C     1000        1   Elemwise{Mul}[(0, 1)]
   ... (remaining 0 Ops account for   0.00%(0.00s) of the runtime)

Apply
------
<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>
  14.1%    14.1%       0.000s       3.58e-07s   1000     9   MakeVector{dtype='float64'}(__logp_m_interval__, __logp_c_interval__, __logp_likelihood)
  13.5%    27.5%       0.000s       3.43e-07s   1000     3   InplaceDimShuffle{x}(m)
  11.5%    39.0%       0.000s       2.93e-07s   1000     4   Elemwise{Composite{(i0 + (-sqr((i1 - ((i2 * i3) + i4)))))}}(TensorConstant{(1,) of -1..0664093453}, TensorConstant{[13.059668...46571405]}, InplaceDimShuffle{x}.0, TensorConstant{[0. 1. 2. .. 7. 8. 9.]}, InplaceDimShuffle{x}.0)
  11.3%    50.4%       0.000s       2.88e-07s   1000     1   Elemwise{Composite{(i0 + (i1 * scalar_sigmoid(i2)))}}(TensorConstant{-10.0}, TensorConstant{20.0}, m_interval__)
  10.8%    61.2%       0.000s       2.75e-07s   1000     5   Elemwise{Composite{((i0 + Switch(Cast{int8}((GE((i1 + i2), i1) * LE((i1 + i2), i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 2)](TensorConstant{2.995732273553991}, TensorConstant{-10.0}, Elemwise{Composite{(i0 * scalar_sigmoid(i1))}}.0, TensorConstant{10.0}, TensorConstant{-2.995732273553991}, TensorConstant{-inf}, TensorConstant{2.0}, c_interval__)
   9.7%    70.9%       0.000s       2.47e-07s   1000     6   Elemwise{Composite{((i0 + Switch(Cast{int8}((GE(i1, i2) * LE(i1, i3))), i4, i5)) - ((i6 * scalar_softplus((-i7))) + i7))}}[(0, 1)](TensorConstant{2.995732273553991}, m, TensorConstant{-10.0}, TensorConstant{10.0}, TensorConstant{-2.995732273553991}, TensorConstant{-inf}, TensorConstant{2.0}, m_interval__)
   9.4%    80.3%       0.000s       2.40e-07s   1000     2   InplaceDimShuffle{x}(Elemwise{Composite{(i0 * scalar_sigmoid(i1))}}.0)
   6.0%    86.3%       0.000s       1.52e-07s   1000     0   Elemwise{Composite{(i0 * scalar_sigmoid(i1))}}(TensorConstant{20.0}, c_interval__)
   5.2%    91.4%       0.000s       1.31e-07s   1000    10   Sum{acc_dtype=float64}(MakeVector{dtype='float64'}.0)
   4.5%    96.0%       0.000s       1.15e-07s   1000     7   Sum{acc_dtype=float64}(Elemwise{Composite{(i0 + (-sqr((i1 - ((i2 * i3) + i4)))))}}.0)
   4.0%   100.0%       0.000s       1.03e-07s   1000     8   Elemwise{Mul}[(0, 1)](TensorConstant{0.5}, Sum{acc_dtype=float64}.0)
   ... (remaining 0 Apply instances account for 0.00%(0.00s) of the runtime)

Here are tips to potentially make your code run faster
                 (if you think of new ones, suggest them on the mailing list).
                 Test them first, as they are not guaranteed to always provide a speedup.
  - Try the Theano flag floatX=float32
  - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.

Authors