PyMC3 中具有非分析时间相关参数的 ODE答案

【问题标题】：ODE with non-analytical time-dependent parameters in PyMC3PyMC3 中具有非分析时间相关参数的 ODE
【发布时间】：2021-04-21 10:14:43
【问题描述】：

我正在使用 PyMC3 解决以下 ODE：

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。

谢谢

【问题讨论】：

为什么不把数值常量pow( 10, -6 )直接给成1e-6？
是的，你是对的，我可以这样写。谢谢，但问题仍然存在。

标签： python theano bayesian ode pymc3

【解决方案1】：

由于您已经有了使用scipy.integrate.odeint 的有效实现，您可以使用theano.compile.ops.as_op，尽管它会带来一些不便（请参阅how to fit a method belonging to an instance with pymc3? 和How to write a custom Deterministic or Stochastic in pymc3 with 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)

对不起，我知道这更像是一种解决方法，而不是实际的解决方案......

【讨论】：