python-2.7 - 在 pymc3 : how? 之外使用 pymc3 可能性/后验

Question

出于比较目的，我想利用 PyMC3 之外的后验密度函数。

对于我的研究项目，我想了解 PyMC3 与我自己定制的代码相比的性能如何。因此，我需要将它与我们自己的内部采样器和似然函数进行比较。

我想我想出了如何调用内部PyMC3后验，但感觉很尴尬，想知道是否有更好的方法。现在我正在手动转换变量，而我应该能够将参数字典传递给 pymc 并获得后验密度。这可能以直接的方式进行吗？

非常感谢!

演示代码:

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

已经超过 5 年了，但我认为这值得一个答案。
首先，关于转换，您需要在 pymc3 定义中决定是否要转换这些参数。在这里，sigma 使用区间变换进行变换以避免硬边界。如果您有兴趣访问作为 sigma 函数的后验，则设置 transform=None。如果您进行转换，则可以将“sigma”变量作为模型的确定性参数之一进行访问。
关于访问后部，有一个很好的描述here .使用上面给出的示例，代码变为:

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

请注意，如果先从 sigma 中删除“transform=None”部分，则 sigma 的实际值将成为 logpost_pymc3 函数返回的元组的一部分。它现在是模型的确定性。