在pymc3之外使用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函数返回的元组的一部分。它现在是模型的确定性。