Модель Pymc3 с ODE Solver с использованием Theano

Question

Я использую модель, где средний ответ зависит от решения ODE. Я пытаюсь приспособить эту модель, используя pymc3, но у меня возникают проблемы (связанные с отсутствием тестовых значений) при присоединении решателя ODE к модели.

модель

y_t Логнормально распределено со средним значением mu_t и сигмой стандартного отклонения.

mu_t - это решение для набора ODE в момент времени t.

проблема

Theano/ pymc3 дает ошибку, потому что теано-тензорные переменные, используемые при решении ODE, не имеют тестовых значений. Смотрите ниже копию ошибок. Я пытался установить

th.config.compute_test_value = 'ignore'

но я думаю, что pymc3 изменяет его обратно, чтобы требовать тестовые значения. Я довольно новичок в theano и pymc3, поэтому прошу прощения, если мне не хватает чего-то очевидного

Код

импорт

import pymc3 as pm
import theano as th
import theano.tensor as tt
from FormatData import *
import pandas as pd

Функции для решения ODE

# Runge Kutta integrator
def rungekuttastep(h, y, fprime, *args):
    k1 = h*fprime(y, *args)
    k2 = h*fprime(y + 0.5*k1, *args)
    k3 = h*fprime(y + 0.5*k2, *args)
    k4 = h*fprime(y + k3, *args)
    y_np1 = y + (1./6.)*k1 + (1./3.)*k2 + (1./3.)*k3 + (1./6.)*k4
    return y_np1

# ODE equations for my model
def ODE(y, *args):
    alpha = args
    yprime = tt.zeros_like(y)
    yprime = tt.set_subtensor(yprime[0],  alpha[0]*y[1] - alpha[1]*y[0])
    yprime = tt.set_subtensor(yprime[1], -alpha[2]*y[0]*y[1])
    return yprime

# Function to return ODE values given parameters
def calcFittedTitreVals(alpha, niter, hsize, inits):
    y0 = tt.vector('y0')
    h = tt.scalar('h')
    i = tt.iscalar('i')
    alpha0 = tt.scalar('alpha0')
    alpha1 = tt.scalar('alpha1')
    alpha2 = tt.scalar('alpha2')

    result, updates = th.scan(fn=lambda y0, h: rungekuttastep(h, y0, ODE, alpha0, alpha1, alpha2), 
                              outputs_info=[{'initial': y0}], non_sequences=h, n_steps=i)
    odeint = th.function(inputs=[h, y0, i, alpha0, alpha1, alpha2], outputs=result, updates=updates)
    z1 = odeint(h=hsize,   # size of the interval
                y0=inits,  # starting values
                i=niter,   # number of iterations
                alpha0=alpha[0], alpha1=alpha[1], alpha2=alpha[2])
    C = z1[:, 0]
    A = z1[:, 1]
    return C, A

Генерировать данные

t = np.arange(0, 45, 0.1)  # times at which to solve ODE
alpha = np.array([2, 0.4, 0.0001])  # true paramter values ODE
C, A = calcFittedTitreVals(alpha, niter=450, hsize=0.1, inits=[0, 1200])
td = np.arange(0, 45, 1) # times at which I observe data
sigma = 0.1
indices = np.array(np.searchsorted(t, td)).flatten()
DATA = pd.DataFrame(
    data={'observed': np.random.lognormal(np.log(C[indices]), sigma),
          'true': C[indices], 'time': td})

модельная функция pymc3

def titreLogNormal(Y, hsize, inits, times):
    Y = th.shared(Y)
    inits = th.shared(inits)
    timesG = np.arange(0, 45, step=hsize) 
    indices = np.array(np.searchsorted(timesG, times)).flatten()
    nTsteps = th.shared(timesG.shape[0])
    hsize = th.shared(hsize)

    y0 = tt.vector('y0')
    h = tt.scalar('h')
    i = tt.iscalar('i')
    alpha0 = tt.scalar('alpha0')
    alpha1 = tt.scalar('alpha1')
    alpha2 = tt.scalar('alpha2')

    result, updates = th.scan(fn=lambda y0, h: rungekuttastep(h, y0, ODE, alpha0, alpha1, alpha2),
                              outputs_info=[{'initial': y0}], non_sequences=h, n_steps=i)
    odeint = th.function(inputs=[h, y0, i, alpha0, alpha1, alpha2], outputs=result, updates=updates)

    model = pm.Model()

    with model:
        alpha = pm.Gamma('alpha', 0., 10., shape=3, testval=[2, 0.4, 0.001])
        sigma = pm.Gamma('sigma', 0.1, 0.1, testval=0.1)

        res = odeint(h=hsize, y=inits, i=nTsteps, alpha0=alpha[0], alpha1=alpha[1], alpha2=alpha[2])

        mu = pm.Deterministic("mu", res[indices, 0])
        y = pm.Lognormal('y', mu, sigma, observed=Y)
        return model

Создать модель с данными

model = titreLogNormal(
    Y=np.array(DATA[['observed']]).flatten(),
    hsize=0.1, inits={'a': 0, 'p': 1200},
    times=np.array(DATA[['time']]).flatten())

ошибки

Traceback (most recent call last):
  File "/home/millerp/.local/lib/python3.5/site-packages/theano/gof/op.py", line 625, in __call__
    storage_map[ins] = [self._get_test_value(ins)]
  File "/home/millerp/.local/lib/python3.5/site-packages/theano/gof/op.py", line 581, in _get_test_value
    raise AttributeError('%s has no test value %s' % (v, detailed_err_msg))
AttributeError: y0 has no test value  
Backtrace when that variable is created:


  File "/home/millerp/pycharm/pycharm-edu-3.5.1/helpers/pydev/_pydev_bundle/pydev_console_utils.py", line 251, in add_exec
more = self.do_add_exec(code_fragment)
  File "/home/millerp/pycharm/pycharm-edu-3.5.1/helpers/pydev/_pydev_bundle/pydev_ipython_console.py", line 41, in do_add_exec
    res = bool(self.interpreter.add_exec(codeFragment.text))
  File "/home/millerp/pycharm/pycharm-edu-3.5.1/helpers/pydev/_pydev_bundle/pydev_ipython_console_011.py", line 455, in add_exec
self.ipython.run_cell(line, store_history=True)
  File "/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py", line 2717, in run_cell
interactivity=interactivity, compiler=compiler, result=result)
  File "/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py", line 2821, in run_ast_nodes
if self.run_code(code, result):
  File "/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py", line 2881, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-54c976fefe1e>", line 99, in <module>
times=np.array(DATA[['time']]).flatten()
  File "<ipython-input-2-54c976fefe1e>", line 71, in titreLogNormal
y0 = tt.vector('y0')


During handling of the above exception, another exception occurred:
Traceback (most recent call last):
  File "/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py", line 2881, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-54c976fefe1e>", line 99, in <module>
times=np.array(DATA[['time']]).flatten()
  File "<ipython-input-2-54c976fefe1e>", line 86, in titreLogNormal
outputs_info=[{'initial': y0}], non_sequences=h, n_steps=i)
  File "/home/millerp/.local/lib/python3.5/site-packages/theano/scan_module/scan.py", line 660, in scan
tensor.shape_padleft(actual_arg), 0),
  File "/home/millerp/.local/lib/python3.5/site-packages/theano/tensor/basic.py", line 4429, in shape_padleft
return DimShuffle(_t.broadcastable, pattern)(_t)
  File "/home/millerp/.local/lib/python3.5/site-packages/theano/gof/op.py", line 639, in __call__
(i, ins, node, detailed_err_msg))
ValueError: Cannot compute test value: input 0 (y0) of Op InplaceDimShuffle{x,0}(y0) missing default value.
Backtrace when that variable is created:


  File "/home/millerp/pycharm/pycharm-edu-3.5.1/helpers/pydev/_pydev_bundle/pydev_console_utils.py", line 251, in add_exec
more = self.do_add_exec(code_fragment)
  File "/home/millerp/pycharm/pycharm-edu-3.5.1/helpers/pydev/_pydev_bundle/pydev_ipython_console.py", line 41, in do_add_exec
res = bool(self.interpreter.add_exec(codeFragment.text))
  File "/home/millerp/pycharm/pycharm-edu-3.5.1/helpers/pydev/_pydev_bundle/pydev_ipython_console_011.py", line 455, in add_exec
self.ipython.run_cell(line, store_history=True)
  File "/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py", line 2717, in run_cell
interactivity=interactivity, compiler=compiler, result=result)
  File "/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py", line 2821, in run_ast_nodes
if self.run_code(code, result):
  File "/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py", line 2881, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-54c976fefe1e>", line 99, in <module>
times=np.array(DATA[['time']]).flatten()
  File "<ipython-input-2-54c976fefe1e>", line 71, in titreLogNormal
y0 = tt.vector('y0')