Hi there,
The demo has been very useful for me to understand how to use an external likelihood function. However, I’m a bit confused about the number of likelihood evaluations: if a Metropolis sampling is used, I expect the number of likelihood evaluations to be (almost) the same as the number of draws, as we need 1 evaluation per Markov step; However, in the attached minimal code, the number of likelihood evaluations is almost 4 times of draws. Could you please give a hint on why it’s like this?
Regards,
Charlie
The minimal (but complete) code is only slightly simplified from the demo so that the likelihood is evaluated within the Python environment, instead of Cython.
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 os
# global variable to save number of likelihood evalutions
no_of_eva = 0
def my_model(theta, x):
global no_of_eva
no_of_eva += 1
m, c = theta
return m*x + c
# define your really-complicated likelihood function that uses loads of external codes
def my_loglike(theta, x, data, sigma):
"""
A Gaussian log-likelihood function for a model with parameters given in theta
"""
# print("called")
model = my_model(theta, x)
return -(0.5/sigma**2)*np.sum((data - model)**2)
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
import warnings
# 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})
step = pm.Metropolis()
trace = pm.sample(ndraws, tune=nburn, step = step, discard_tuned_samples=True, chains=1, cores=1)
# print number of likelihood evalutions
print("number of Markov steps (ndraws+nburn): %d"%(ndraws+nburn))
print("number of likelihood evaluations: %d"%no_of_eva)
# plot the traces
_ = pm.traceplot(trace, lines={'m': mtrue, 'c': ctrue})
plt.show()
Screen output:
number of Markov steps (ndraws+nburn): 4000
number of likelihood evaluations: 16004