如何在PyMC中定义一般确定性函数

Question

在我的模型中，我需要通过一段复杂的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的阶段出现了错误，因为在FindFromGrid(x,w,z)函数内部的int()函数需要的是整数，而不是FreeRV。

从预先计算的网格中找到y_hat是很重要的，因为我真实的y_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采样器，因为它需要梯度。

而Metropolis方法似乎没有收敛。

在这种情况下，我应该使用哪种步进方法呢？

暂无标签

如何在PyMC中定义一般确定性函数

更新

1 个回答

撰写回答