在我的模型中，我需要使用複雜的python函數從一組父變量中獲取我的確定性變量的值。如何在PyMC中定義一般確定性函數

有沒有可能這樣做？

以下是一個python代碼，它顯示了我正在嘗試做的簡化情況。

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階段，因爲FindFromGrid(x,w,z)函數內int()功能需要整數FreeRV。

從預先計算的網格中找到y_hat非常重要，因爲我的y_hat真實模型沒有分析形式來表示。

我已經嘗試過使用OpenBUGS，但是我發現here在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採樣器了，因爲它需要漸變。

大都會似乎並沒有收斂。

我應該在這種情況下使用哪種方法？

來源

2014-10-16 indiajoe

您找到正確的解決方案as_op。

關於收斂性：您是否正在使用pm.Metropolis()而不是pm.NUTS()？這個不能收斂的一個原因是Metropolis()默認樣本在聯合空間中，而Metropolis中的Gibbs通常更有效（這是pymc2中的默認值）。話雖如此，我只是合併了這個：https://github.com/pymc-devs/pymc/pull/587它改變了Metropolis和Slice採樣器的默認行爲默認情況下是非阻塞的（所以在Gibbs中）。其他採樣器，如NUTS，主要用於採樣關節空間，但仍默認爲阻塞。你總是可以用kwarg blocked=True明確地設置它。

無論如何，用最新的主人更新pymc，看看是否收斂改善。如果不是，請嘗試使用Slice採樣器。

來源

2014-10-18 09:23:08 twiecki

感謝您的幫助。我想，我上面的示例模型中發生的事情是變量高度相關。所以它在幾個MC鏈之後被卡住了。有沒有辦法在MC跳躍中增加步長？我還想知道爲什麼我爲我的函數獲取一個沒有grad（）屬性的錯誤，而使用pm.find_MAP（）時，如果它使用的是Nelder-Mead優化，而不需要派生。 – indiajoe 2014-10-20 08:21:00

是的，使用Nelder-Mead的'find_MAP（）'應該可以工作。你可以在github上用一些代碼重現這個問題嗎？ – twiecki 2014-10-27 07:20:56

關於卡住，這些是MCMC的痛苦。如果相關性是問題，您可以嘗試調整Metropolis-Hastings的提案分佈或者只運行更多樣本。 – twiecki 2014-10-27 07:21:39

如何在PyMC中定義一般確定性函數

更新

回答

相關問題