1

私のモデルでは、複雑な Python 関数を使用して一連の親変数から決定論的変数の値を取得する必要があります。

それは可能ですか？

以下は、単純化されたケースで何をしようとしているのかを示す pyMC3 コードです。

import numpy as np
import pymc as pm


#Predefine values on two parameter Grid (x,w) for a set of i values (1,2,3)
idata = np.array([1,2,3])
size= 20
gridlength = size*size
Grid = np.empty((gridlength,2+len(idata)))
for x in range(size):
    for w in range(size):
        # A silly version of my real model evaluated on grid.
        Grid[x*size+w,:]= np.array([x,w]+[(x**i + w**i) for i in idata])

# A function to find the nearest value in Grid and return its product with third variable z
def FindFromGrid(x,w,z):
    return Grid[int(x)*size+int(w),2:] * z

#Generate fake Y data with error
yerror = np.random.normal(loc=0.0, scale=9.0, size=len(idata))
ydata = Grid[16*size+12,2:]*3.6 + yerror  # ie. True x= 16, w= 12 and z= 3.6


with pm.Model() as model:

    #Priors
    x = pm.Uniform('x',lower=0,upper= size)
    w = pm.Uniform('w',lower=0,upper =size)
    z = pm.Uniform('z',lower=-5,upper =10)

    #Expected value
    y_hat = pm.Deterministic('y_hat',FindFromGrid(x,w,z))

    #Data likelihood
    ysigmas = np.ones(len(idata))*9.0 
    y_like = pm.Normal('y_like',mu= y_hat, sd=ysigmas, observed=ydata)

    # Inference...
    start = pm.find_MAP() # Find starting value by optimization
    step = pm.NUTS(state=start) # Instantiate MCMC sampling algorithm
    trace = pm.sample(1000, step, start=start, progressbar=False) # draw 1000 posterior samples using NUTS sampling


print('The trace plot')
fig = pm.traceplot(trace, lines={'x': 16, 'w': 12, 'z':3.6})
fig.show()

このコードを実行すると、y_hat の段階でエラーが発生します。これは、int()関数内のFindFromGrid(x,w,z)関数が FreeRV ではなく整数を必要とするためです。

y_haty_hatの実際のモデルには表現する分析形式がないため、事前に計算されたグリッドからの検索は重要です。

以前に OpenBUGS を使用しようとしましたが、ここで OpenBUGS でこれを行うことができないことがわかりました。PyMC で可能ですか?

アップデート

pyMC github ページの例に基づいて、関数に次のデコレータを追加する必要があることがわかりましたFindFromGrid(x,w,z)。

@pm.theano.compile.ops.as_op(itypes=[t.dscalar, t.dscalar, t.dscalar],otypes=[t.dvector])

これにより、上記の問題が解決されるようです。しかし、グラデーションが必要なため、NUTSサンプラーを使用できなくなりました。

メトロポリスは収束していないようです。

このようなシナリオでは、どのステップ メソッドを使用すればよいですか?

4

1 に答える 1

1

で正しい解が見つかりましたas_op。

収斂について：ひょっとしてpm.Metropolis()代わりに使ってるpm.NUTS()？これが収束できなかった理由の 1 つはMetropolis()、デフォルトではジョイント スペースのサンプルであるのに対し、Metropolis 内の Gibbs の方が効果的であることが多いためです (これは pymc2 のデフォルトでした)。そうは言っても、これをマージしただけです: https://github.com/pymc-devs/pymc/pull/587Metropolisこれにより、サンプラーのデフォルトの動作がデフォルトSliceでブロックされないように変更されます (Gibbs 内)。そのような他のサンプラーNUTSは、主にジョイント スペースをサンプリングするように設計されていますが、既定ではブロックされています。これは kwarg でいつでも明示的に設定できますblocked=True。

とにかく、pymc を最新のマスターで更新し、収束が改善するかどうかを確認します。Sliceそうでない場合は、サンプラーを試してください。

于 2014-10-18T09:23:08.787 に答える