4 
pymc.__version__ = '3.0'  
theano.__version__ = '0.6.0.dev-RELEASE'

複雑な尤度関数で PyMC3 を使用しようとしています:

最初の質問: これは可能ですか?

以下は、 Thomas Wiecki の投稿をガイドとして使用した私の試みです。

import numpy as np
import theano as th
import pymc as pm
import scipy as sp

# Actual data I'm trying to fit
x = np.array([52.08, 58.44, 60.0, 65.0, 65.10, 66.0, 70.0, 87.5, 110.0, 126.0])
y = np.array([0.522, 0.659, 0.462, 0.720, 0.609, 0.696, 0.667, 0.870, 0.889,  0.919])
yerr = np.array([0.104,  0.071,  0.138,  0.035,  0.102, 0.096,  0.136,  0.031, 0.024, 0.035])

th.config.compute_test_value = 'off'
a = th.tensor.dscalar('a')

with pm.Model() as model:
    # Priors
    alpha = pm.Normal('alpha', mu=0.3, sd=5)
    sig_alpha = pm.Normal('sig_alpha', mu=0.03, sd=5)
    t_double = pm.Normal('t_double', mu=4, sd=20)
    t_delay = pm.Normal('t_delay', mu=21, sd=20)
    nu = pm.Uniform('nu', lower=0, upper=20)

    # Some functions needed for calculation of the y estimator
    def T(eqd):
        doses = np.array([52.08, 58.44, 60.0, 65.0, 65.10, 
                          66.0, 70.0, 87.5, 110.0, 126.0])
        tmt_times = np.array([29,29,43,29,36,48,22,11,7,8])
        return np.interp(eqd, doses, tmt_times)

    def TCP(a):
        time = T(x)
        BCP = pm.exp(-1E7*pm.exp(-alpha*x*1.2 + 0.69315/t_delay(time-t_double)))
        return pm.prod(BCP)

    def normpdf(a, alpha, sig_alpha):
        return 1./(sig_alpha*pm.sqrt(2.*np.pi))*pm.exp(-pm.sqr(a-alpha)/(2*pm.sqr(sig_alpha)))

    def normcdf(a, alpha, sig_alpha):
        return 1./2.*(1+pm.erf((a-alpha)/(sig_alpha*pm.sqrt(2))))

    def integrand(a):
        return normpdf(a,alpha,sig_alpha)/(1.-normcdf(0,alpha,sig_alpha))*TCP(a)

    func = th.function([a,alpha,sig_alpha,t_double,t_delay], integrand(a))

    y_est = sp.integrate.quad(func(a, alpha, sig_alpha,t_double,t_delay), 0, np.inf)[0]

    likelihood = pm.T('TCP', mu=y_est, nu=nu, observed=y_tcp)

    start = pm.find_MAP()
    step = pm.NUTS(state=start)
    trace = pm.sample(2000, step, start=start, progressbar=True)

y_est の式に関する次のメッセージが生成されます。

TypeError: ('名前 ":42" のインデックス 0 (0 ベース) の theano 関数への入力引数が正しくありません', '配列のようなオブジェクトが必要でしたが、変数が見つかりました: おそらく、(おそらく共有)数値配列の代わりに変数?')

ここまで来るのに様々なハードルを乗り越えてきましたが、ここが行き詰っています。それで、私の最初の質問への答えが「はい」であるならば、私は正しい方向に進んでいますか? どんなガイダンスも役に立ちます！

注:これは、私が見つけた同様の質問と、別の質問です。

免責事項：私はこれに非常に慣れていません。私の唯一の以前の経験は、Thomas の投稿の線形回帰の例をうまく再現したことです。また、Theano テスト スイートを正常に実行したので、動作することがわかりました。

4

1 に答える 1

3

はい、複雑または恣意的な可能性で何かを作ることは可能です。それはあなたがここでやっていることのようには見えませんが。ある変数から別の変数への複雑な変換、統合ステップがあるようです。

あなたの特定の例外は、 integrate.quad が pymc ではなく、numpy 配列を期待していることですVariable。pymc 内で quad を実行する場合は、そのためのカスタム theanoOp (派生物を含む) を作成する必要があります。

于 2015-02-11T21:34:02.313 に答える