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