在我的模型中，我需要使用复杂的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中定义一般确定性函数

更新

回答

相关问题