Blackbox likelihood example

Question

I'm trying to understand how to use a black box likelihood function in pymc. Basically, this is explained here. I have tried implementing this on my own with a very simple Python model (a double logistic function), and no gradient. In addition to the code in the link I mentioned, that sets up the theano wrappers around the loglikelihood function, I have the following code

# Define the model
def my_model(p, t, m_m, m_M):
    """An simple double logistic model"""
    return m_m + (m_M - m_m) * (
        1.0 / (1 + np.exp(-p[0] * (t - p[1] * 365)))
        + 1.0 / (1 + np.exp(p[2] * (t - p[3] * 365)))
        - 1
    )

# The loglikelihood function
def log_lklhood(theta, t, data, sigma):
    """Normal log-likelihoood"""
    y_pred = my_model(theta, t, data.max(), data.min())
    retval = -(0.5 / sigma ** 2) * np.sum((data - y_pred) ** 2)
    return retval

## Experimental data

data = np.array([0.104, 0.133, 0.131, 0.329, 0.669, 0.822, 0.847, 0.807, 0.638, 0.486, 0.162, 0.125])
t = np.array([1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336])
sigma = 0.02 # Uncertainty in the observations
# Straight outta internet...
logl = LogLike(log_lklhood, data, t, sigma)
ndraws = 500  # number of draws from the distribution
nburn = 100  # number of "burn-in points" (which we'll discard)

# use PyMC3 to sampler from log-likelihood
with pm.Model():
    theta1 = pm.Normal("theta_1", mu=0, sigma=1, testval=0.1)
    theta3 = pm.Normal("theta_3", mu=0, sigma=1, testval=0.1)
    theta2 = pm.Uniform("theta_2", lower=0.01, upper=0.5, testval=120.0 / 365)
    theta4 = pm.Uniform("theta_4", lower=0.51, upper=0.99, testval=250.0 / 365)

    # convert theta1...theta4 to a tensor vector
    theta = tt.as_tensor_variable([theta1, theta2, theta3, theta4])
    # use a DensityDist (use a lamdba function to "call" the Op)
    pm.DensityDist("likelihood", lambda v: logl(v), observed={"v": theta})
    # Use simple metropolis
    step = pm.Metropolis()
    trace = pm.sample(ndraws, tune=nburn, step=step, discard_tuned_samples=True)

This goes over the sampling but then fails with rather cryptic error:

Multiprocess sampling (2 chains in 2 jobs)
CompoundStep
>Metropolis: [theta_4]
>Metropolis: [theta_2]
>Metropolis: [theta_3]
>Metropolis: [theta_1]

100.00% [1010/1010 00:01<00:00 Sampling 2 chains, 0 divergences]

Sampling 2 chains for 5 tune and 500 draw iterations (10 + 1_000 draws total) took 2 seconds.

---------------------------------------------------------------------------
MissingInputError                         Traceback (most recent call last)
<ipython-input-14-3e5ab1ff0971> in <module>
     17     pm.DensityDist('likelihood', lambda v: logl(v), observed={'v': theta})
     18     step = pm.Metropolis()
---> 19     trace = pm.sample(ndraws, tune=nburn, step=step, discard_tuned_samples=True)

[...]


MissingInputError: Input 0 of the graph (indices start from 0), used to compute sigmoid(theta_4_interval__), was not provided and not given a value. Use the Theano flag exception_verbosity='high', for more information on this error.

I'm very new to pymc, so don't really know what this error means, or what I am doing wrong. There are a number of SO questions(eg0 on this that suggest I am calling a python function with a theano array as the basis of the problem. However, I don't see where this is happening.

Answer 1

So it turns out that there's an issue with the blackbox likelihood example: Don't use pm.DensityDist, but rather pm.Potential (see this arviz issue). The example now works correctly, even using scipy.optimize.approx_fprime to approximate the gradient of the log-likelihood:

import numpy as np
import pymc3 as pm
import theano.tensor as tt
import arviz as az
import scipy.optimize # Can't be bothered calculating derivatives by hand

## Experimental data

data = np.array([0.104, 0.133, 0.131, 0.329, 0.669, 0.822, 0.847, 0.807, 0.638, 0.486, 0.162, 0.125])
t = np.array([1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336])

## LogLikelihood and gradient of the LogLikelihood functions
def log_likelihood(theta, data, t):
    if type(theta) == list:
        theta = theta[0]
    (
        theta1,
        theta2,
        theta3,
        theta4,
        sigma,
    ) = theta

    m_m = data.min()
    m_M = data.max()
    y_pred = m_m + (m_M - m_m) * (
        1.0 / (1 + np.exp(-theta1 * (t - theta2)))
        + 1.0 / (1 + np.exp(theta3 * (t - theta4)))
        - 1
    )

    logp = -len(data) * np.log(np.sqrt(2.0 * np.pi) * sigma)
    logp += -np.sum((data - y_pred) ** 2.0) / (2.0 * sigma ** 2.0)
    return logp


def der_log_likelihood(theta, data, t):
    def lnlike(values):
        return log_likelihood(values, data, t)

    eps = np.sqrt(np.finfo(float).eps)
    grads = scipy.optimize.approx_fprime(theta[0], lnlike, eps * np.ones(len(theta)))
    return grads

## Wrapper classes to theano-ize LogLklhood and gradient...
class Loglike(tt.Op):
    itypes = [tt.dvector]
    otypes = [tt.dscalar]

    def __init__(self, data, t):
        self.data = data
        self.t = t
        self.loglike_grad = LoglikeGrad(self.data, self.t)

    def perform(self, node, inputs, outputs):
        logp = log_likelihood(inputs, self.data, self.t)
        outputs[0][0] = np.array(logp)

    def grad(self, inputs, grad_outputs):
        (theta,) = inputs
        grads = self.loglike_grad(theta)
        return [grad_outputs[0] * grads]


class LoglikeGrad(tt.Op):
    itypes = [tt.dvector]
    otypes = [tt.dvector]

    def __init__(self, data, t):
        self.der_likelihood = der_log_likelihood
        self.data = data
        self.t = t

    def perform(self, node, inputs, outputs):
        (theta,) = inputs
        grads = self.der_likelihood(inputs, self.data, self.t)
        outputs[0][0] = grads

## Sample!
loglike = Loglike(data, t)
with pm.Model() as model:
    theta1 = pm.Normal("theta_1", mu=0, sigma=1, testval=0.1)
    theta3 = pm.Normal("theta_3", mu=0, sigma=1, testval=0.1)
    theta2 = pm.DiscreteUniform("theta_2", lower=1, upper=180, testval=120)
    theta4 = pm.DiscreteUniform("theta_4", lower=181, upper=365, testval=250)
    sigma = pm.HalfNormal("sigma", sigma=0.6, testval=0.01)
    # convert parameters to a tensor vector
    theta = tt.as_tensor_variable([theta1, theta2, theta3, theta4, sigma])
    # Do not use a DensityDist!
    pm.Potential("like", loglike(theta))

with model:
    trace = pm.sample(step=pm.NUTS())
    print(pm.summary(trace).to_string())

Happily, the above code results in

Multiprocess sampling (2 chains in 2 jobs)
CompoundStep
>NUTS: [sigma, theta_3, theta_1]
>CompoundStep
>>Metropolis: [theta_4]
>>Metropolis: [theta_2]

100.00% [4000/4000 00:16<00:00 Sampling 2 chains, 0 divergences]

Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 17 seconds.
The number of effective samples is smaller than 25% for some parameters.

            mean     sd   hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_mean  ess_sd  ess_bulk  ess_tail  r_hat
theta_1    0.073  0.019    0.049    0.095      0.001    0.001     258.0   230.0     641.0     305.0   1.01
theta_3    0.052  0.010    0.037    0.070      0.000    0.000     404.0   361.0     681.0     379.0   1.01
theta_2  104.616  3.004  100.000  110.000      0.178    0.126     285.0   283.0     307.0     312.0   1.01
theta_4  270.178  3.139  264.000  275.000      0.147    0.104     455.0   455.0     454.0     606.0   1.01
sigma      0.040  0.013    0.020    0.063      0.001    0.001     324.0   298.0     410.0     422.0   1.01