ODE с неаналитическими параметрами, зависящими от времени, в PyMC3

Я работаю над решением следующего ODE с PyMC3:

def production( y, t, p ):
    return p[0]*getBeam( t ) - p[1]*y[0]

getBeam( t ) - это мой коэффициент, зависящий от времени. Эти коэффициенты задаются массивом данных, доступ к которым осуществляется по временному индексу следующим образом:

def getBeam( t ):
    nBeam = I[int(t/10)]*pow( 10, -6 )/q_e
    return nBeam

Я успешно реализовал это с помощью scipy.integrate.odeint, но мне трудно сделать это с pymc3.ode. На самом деле, используя следующее:

ode_model = DifferentialEquation(func=production, times=x, n_states=1, n_theta=3, t0=0)
with pm.Model() as model:
    a = pm.Uniform( "S-Factor", lower=0.01, upper=100 )
    ode_solution = ode_model(y0=[0], theta=[a, Yield, lambd])

Очевидно, я получаю ошибку TypeError: __trunc__ returned non-Integral (type TensorVariable), так как t является TensorVariable, поэтому его нельзя использовать для доступа к массиву, в котором хранятся коэффициенты.

Есть ли способ преодолеть эту трудность? Я думал об использовании theano.function, но не могу заставить его работать, так как, к сожалению, коэффициенты не могут быть выражены какой-либо аналитической функцией: они просто хранятся внутри массива I, индекс которого представляет временную переменную t.

Спасибо


python bayesian ode theano pymc3
person Jakub Skowroński    schedule 21.04.2021    source источник
comment
Почему вы не даете числовую константу pow( 10, -6 ) непосредственно как 1e-6?   -  person Lutz Lehmann    schedule 21.04.2021
comment
Да, вы правы, я могу так написать. Спасибо, но проблема осталась.   -  person Jakub Skowroński    schedule 21.04.2021

Ответы (1)


arrow_upward
0
arrow_downward

Поскольку у вас уже есть работающая реализация с scipy.integrate.odeint, вы можете использовать theano.compile.ops.as_op, хотя это сопряжено с некоторыми неудобствами (см. принадлежность-экземпляру-с-pymc3"> как подобрать метод, принадлежащий экземпляру, с помощью pymc3? и Как написать пользовательский детерминистический или стохастический анализ в pymc3 с помощью theano.op?)

Используя ваши точные определения для production и getBeam, у меня работает следующий код:

from scipy.integrate import odeint
from theano.compile.ops import as_op
import theano.tensor as tt
import pymc3 as pm

def ScipySolveODE(a):
    return odeint(production, y0=[0], t=x, args=([a, Yield, lambd],)).flatten()

@as_op(itypes=[tt.dscalar], otypes=[tt.dvector])
def TheanoSolveODE(a):
    return ScipySolveODE(a)

with pm.Model() as model:
    a = pm.Uniform( "S-Factor", lower=0.01, upper=100 )
    ode_solution = TheanoSolveODE(a)

Извините, я знаю, что это скорее обходной путь, чем фактическое решение...

person Alexander S. Brunmayr    schedule 18.05.2021