用theano解决pymc3中的简单化学网络方程

Question

我正在尝试使用贝叶斯推理求解简单的化学网络 A->B（反应率 k1）和 A1->B（反应率 k2）。我希望得到 k1 和 k2 的敏感性分析。如果 A、A1 和 B 是我的常量变量，唯一合乎逻辑的事情是，如果例如 k1 减少，k2 应该按比例增加，反之亦然。但是我在 pymc3 中遇到了一些 ODE 的问题。所以这是我的尝试：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp
import seaborn
import pymc3 as pm
import theano.tensor as T
from theano.compile.ops import as_op
from sys import exit
time = 10
Nt = 11
tt = np.linspace(0,time, Nt)
y0 = [1,2,0]
k1, k2 = 1, 1
#Actual Solution of the Differential Equation(Used to generate data)

def real(t,c):
        da_dt = -k1*c[0]
        da1_dt = -k2*c[1]
        db_dt = k1*c[0] + k2*c[1]
        return da_dt, da1_dt, db_dt

c_est = solve_ivp(real, t_span = [0,time], t_eval = tt, y0 = y0)


#Method For Solving the ODE
def lv(xdata, k1=1, k2=1):
    def equat(c,t):
        da_dt = -k1*c[0]
        da1_dt = -k2*c[1]
        db_dt = k1*c[0] + k2*c[1]
        return da_dt, da1_dt, db_dt
    Y, dict  = odeint(equat,y0,xdata,full_output=True)
    return Y

#Generating Data for Bayesian Inference
k1, k2 = 1, 1
ydata = c_est.y

# Adding some error to the ydata points
yerror = 10*np.random.rand(Nt)
ydata += np.random.normal(0.0, np.sqrt(yerror))
ydata = np.ravel(ydata)

@as_op(itypes=[T.dscalar, T.dscalar], otypes=[T.dvector])

def func(al,be):
    Q = lv(tt, k1=al, k2=be)
    return np.ravel(Q)


# Number of Samples and Initial Conditions
nsample = 5000
y0 = 1.0
sd = 0.2

# Model for Bayesian Inference
model = pm.Model()
with model:
    # Priors for unknown model parameters
    k1 = pm.HalfNormal('k1', sd = sd)
    k2 = pm.HalfNormal('k2', sd = sd)

    # Expected value of outcome
    mu = func(k1,k2)


    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=mu, sd=yerror, observed=y_data)

    trace = pm.sample(nsample, nchains=1)


pm.traceplot(trace)
plt.show()

但它不会 "loop" 通过 equat 函数。输出错误：

Traceback (most recent call last):

  File "<ipython-input-16-14ca425a8735>", line 1, in <module>
    runfile('/folder/code.py', wdir='/folder')

  File "/anaconda3/lib/python3.7/site-packages/spyder_kernels/customize/spydercustomize.py", line 786, in runfile
    execfile(filename, namespace)

  File "/anaconda3/lib/python3.7/site-packages/spyder_kernels/customize/spydercustomize.py", line 110, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/code.py", line 77, in <module>
    mu = func(k1,k2)

  File "/anaconda3/lib/python3.7/site-packages/theano/gof/op.py", line 674, in __call__
    required = thunk()

  File "/anaconda3/lib/python3.7/site-packages/theano/gof/op.py", line 892, in rval
    r = p(n, [x[0] for x in i], o)

  File "/anaconda3/lib/python3.7/site-packages/theano/compile/ops.py", line 555, in perform
    outs = self.__fn(*inputs)

  File "/code.py", line 60, in func
    Q = lv(tt, k1=al, k2=be)

  File "/code.py", line 42, in lv
    Y, dict  = odeint(equat,y0,xdata,full_output=True)

  File "/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py", line 233, in odeint
    int(bool(tfirst)))

  File "/code.py", line 39, in equat
    da1_dt = -k2*c[1]

IndexError: index 1 is out of bounds for axis 0 with size 1

我快疯了。 :( 我什至不知道我是否走在正确的道路上。编辑，更正，但现在它显示另一个错误。

Answer 1

如果其他人在这里有困难我解决了！

from scipy.integrate import odeint, solve_ivp
import numpy as np
import matplotlib.pyplot as plt
from theano.compile.ops import as_op
import theano.tensor as T
import pymc3 as pm
import copy
from sys import exit
time = 10
Nt = 11
tt = np.linspace(0,time, Nt+1)
y0 = [1,2,0]
k1, k2 = 1, 1
def real_equat(t,c):
    da_dt = -k1*c[0]
    da1_dt = -k2*c[1]
    db_dt = k1*c[0] + k2*c[1]
    return da_dt, da1_dt, db_dt
z = solve_ivp(real_equat, t_span=[0,time], t_eval= tt, y0 = y0)
def lv(xdata, k1=k1, k2=k2):
    def equat(c,t):
        da_dt = -k1*c[0]
        da1_dt = -k2*c[1]
        db_dt = k1*c[0] + k2*c[1]
        return da_dt, da1_dt, db_dt
    Y, dict = odeint(equat,y0,tt,full_output=True)
    return Y
a = z.y
ydata = copy.copy(a)

yerror = 10*np.random.rand(Nt+1)
ydata += np.random.normal(0.0, np.sqrt(yerror))
ydata = np.ravel(ydata)

@as_op(itypes=[T.dscalar, T.dscalar], otypes=[T.dvector])

def func(al,be):
    Q = lv(tt, k1 = al, k2 = be)
    return np.ravel(Q)
niter = 10
model = pm.Model()
with model:
    # Priors for unknown model parameters
    k1 = pm.Uniform('k1', upper = 1.2, lower = 0.8)
    k2 = pm.Uniform('k2', upper = 1.2, lower = 0.8)

    # Expected value of outcome
    mu = func(k1,k2)


    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=mu, sd=0.2, observed=ydata)

    trace = pm.sample(niter = niter, nchains=4)

用theano解决pymc3中的简单化学网络方程

Solving simple chemical network odes in pymc3 with theano

bayesian

python-3.x

theano

pymc3