5

比較のために、PyMC3 の外部で事後密度関数を利用したいと思います。

私の研究プロジェクトでは、PyMC3 が独自のカスタム コードと比較してどの程度うまく機能しているかを調べたいと考えています。そのため、社内のサンプラーや尤度関数と比較する必要があります。

内部 PyMC3 posterior を呼び出す方法を理解したと思いますが、非常にぎこちなく感じます。より良い方法があれば知りたいです。現在、変数を手動で変換していますが、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)))
4

1 に答える 1

2

5年以上経ちましたが、これは答えに値すると思いました。

まず、変換に関して、pymc3 定義内で、これらのパラメーターを変換するかどうかを決定する必要があります。ここで、厳密な境界を避けるために、シグマは区間変換で変換されていました。シグマの関数として事後変数にアクセスすることに関心がある場合は、transform=None を設定します。変換を行う場合、「シグマ」変数は、モデルの決定論的パラメーターの 1 つとしてアクセスできます。

後部へのアクセスに関しては、ここに素晴らしい説明があります。上記の例では、コードは次のようになります。

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

'transform=None' 部分を sigma before から削除すると、sigma の実際の値が logpost_pymc3 関数によって返されるタプルの一部になることに注意してください。これは、モデルの決定論です。

于 2020-11-03T11:11:08.720 に答える