В моей модели мне нужно получить значение моей детерминированной переменной из набора родительских переменных, используя сложную функцию 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_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, так как ему нужен градиент.
Мегаполис, похоже, не сходится.
Какой пошаговый метод следует использовать в подобном сценарии?
