Using the pymc3 likelihood/posterior outside of pymc3: how?

Question

For comparison purposes, I want to utilize the posterior density function outside of PyMC3.

For my research project, I want to find out how well PyMC3 is performing compared to my own custom made code. As such, I need to compare it to our own in-house samplers and likelihood functions.

I think I figured out how to call the internal PyMC3 posterior, but it feels very awkward, and I want to know if there is a better way. Right now I am hand-transforming variables, whereas I should just be able to pass pymc a parameter dictionary and get the posterior density. Is this possible in a straightforward manner?

Thanks a lot!

Demo code:

import numpy as np
import pymc3 as pm
import scipy.stats as st

# Simple data, with sigma = 4. We want to estimate sigma
sigma_inject = 4.0
data = np.random.randn(10) * sigma_inject

# Prior interval for sigma
a, b = 0.0, 20.0

# Build PyMC model
with pm.Model() as model:
    sigma = pm.Uniform('sigma', a, b)      # Prior uniform between 0.0 and 20.0
    likelihood = pm.Normal('data', 0.0, sd=sigma, observed=data)

# Write my own likelihood
def logpost_self(sig, data):
    loglik = np.sum(st.norm(loc=0.0, scale=sig).logpdf(data))   # Gaussian
    logpr = np.log(1.0 / (b-a))                                 # Uniform prior
    return loglik + logpr

# Utilize PyMC likelihood (Have to hand-transform parameters)
def logpost_pymc(sig, model):
    sigma_interval = np.log((sig - a) / (b - sig))    # Parameter transformation
    ldrdx = np.log(1.0/(sig-a) + 1.0/(b-sig))         # Jacobian
    return model.logp({'sigma_interval':sigma_interval}) + ldrdx

print("Own posterior:   {0}".format(logpost_self(1.0, data)))
print("PyMC3 posterior: {0}".format(logpost_pymc(1.0, model)))

Answer 1

It's been over 5 years, but I figured this deserves an answer.

Firstly, regarding the transformations, you need to decide within the pymc3 definitions whether you want these parameters transformed. Here, sigma was being transformed with an interval transform to avoid hard boundaries. If you are interested in accessing the posterior as a function of sigma, then set transform=None. If you do transform, then the 'sigma' variable will be accessible as one of the deterministic parameters of the model.

Regarding accessing the posterior, there is a great description here. With the example given above, the code becomes:

import numpy as np
import pymc3 as pm
import theano as th
import scipy.stats as st

# Simple data, with sigma = 4. We want to estimate sigma
sigma_inject = 4.0
data = np.random.randn(10) * sigma_inject

# Prior interval for sigma
a, b = 0.1, 20.0

# Build PyMC model
with pm.Model() as model:
    sigma = pm.Uniform('sigma', a, b, transform=None)      # Prior uniform between 0.0 and 20.0
    likelihood = pm.Normal('data', mu=0.0, sigma=sigma, observed=data)

# Write my own likelihood
def logpost_self(sig, data):
    loglik = np.sum(st.norm(loc=0.0, scale=sig).logpdf(data))   # Gaussian
    logpr = np.log(1.0 / (b-a))                                 # Uniform prior
    return loglik + logpr

with model:
    # Compile model posterior into a theano function
    f = th.function(model.vars, [model.logpt] + model.deterministics)

    def logpost_pymc3(params):
        dct = model.bijection.rmap(params)
        args = (dct[k.name] for k in model.vars)
        results = f(*args)
        return tuple(results)

print("Own posterior:   {0}".format(logpost_self(1.0, data)))
print("PyMC3 posterior: {0}".format(logpost_pymc3([1.0])))

Note that if you remove the 'transform=None' part from the sigma prior, then the actual value of sigma becomes part of the tuple that is returned by the logpost_pymc3 function. It's now a deterministic of the model.