Решение простых химических сетевых одов в pymc3 с помощью theano - PullRequest
Новая
       85

Решение простых химических сетевых одов в pymc3 с помощью theano

0 голосов
/ 12 июня 2019

Я пытаюсь решить простую химическую сеть A-> B (скорость реакции k1) и A1-> B (скорость реакции k2) с байесовским выводом. Я надеюсь получить анализ чувствительности k1 и k2. Если A, A1 и B являются моими постоянными переменными, то вполне логично, что если, например, k1 уменьшается, k2 должно увеличиваться на некоторую пропорциональную величину и наоборот. Но у меня проблемы с ODE в pymc3. Итак, вот моя попытка: 

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()

Но это не "цикл" через функцию экватора. Ошибка вывода: 

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

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

1 Ответ

0 голосов
/ 13 июня 2019

Если у кого-то еще есть трудности, я решил их! 

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)
...