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

But it doesn’t “loop” through equat function. Output error:

```
Traceback (most recent call last):
File "<ipython-input-217-14ca425a8735>", line 1, in <module>
runfile('/codepy', 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 79, in <module>
mu = func(alpha,beta)
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 62, in func
Q = lv(tt, k1=al, k2=be)
File "/code.py", line 43, in lv
Y, dict = odeint(equat_2,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 40, in equat_2
da1_dt = -k2*c[1]
```

I tried something with this as well, I don’t know why it doesn’t vary k1 and k2.