python 2.7 - Using the pymc3 likelihood/posterior outside of pymc3: how? -


for comparison purposes, want utilize posterior density function outside of pymc3.

for research project, want find out how pymc3 performing compared own custom made code. such, need compare our own in-house samplers , likelihood functions.

i think figured out how call internal pymc3 posterior, feels awkward, , want know if there better way. right hand-transforming variables, whereas should able pass pymc parameter dictionary , posterior density. possible in straightforward manner?

thanks lot!

demo code:

import numpy np import pymc3 pm import scipy.stats st  # simple data, sigma = 4. want estimate sigma sigma_inject = 4.0 data = np.random.randn(10) * sigma_inject  # prior interval sigma a, b = 0.0, 20.0  # build pymc model pm.model() model:     sigma = pm.uniform('sigma', a, b)      # prior uniform between 0.0 , 20.0     likelihood = pm.normal('data', 0.0, sd=sigma, observed=data)  # write 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 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))) 


Comments

Popular posts from this blog

python - No exponential form of the z-axis in matplotlib-3D-plots -

php - Best Light server (Linux + Web server + Database) for Raspberry Pi -

c# - "Newtonsoft.Json.JsonSerializationException unable to find constructor to use for types" error when deserializing class -