Using pymc3 / posterior probability outside of pymc3: how?

For comparison, I want to use the back density function outside of PyMC3.

In my research project, I want to know how well PyMC3 performs compared to my own code. So I need to compare it with my own probes and likelihood functions.

I think I figured out how to call the inner PyMC3 from behind, but it feels very awkward and I want to know if there is a better way. Right now I am converting the variables manually, whereas I should just pass the parameter dictionary to pymc and get the back density. Is this possible straightforwardly?

Thank you so much!

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)))

      

+3


source to share





All Articles