This example may clarify the problem. Basically, I’ve data from some input or control features X(t) that interact with a system. The system produces an outcome y(t). For simplicity, I’m assuming that y'=X(t)\beta, but it could be a non-linear relationship like eqs 2 and 3 of this reference. Given that I have access to the observed data y_{obs}(t) (red line in plot below) and the input data X(t), I’d like to use PyMC3 to recover the vector of coefficients \beta = [2,3]. To emphasize, I’m modeling the rate of change in y, not y itself.

```
# For reproducibility
np.random.seed(20394)
def f(y, t, p):
return np.dot(p[0][int(t),:], np.array([2,3]))
# Times for observation
times = np.arange(0,20,1)
# Input matrix
X = np.random.randn(20,2)
y = odeint(f, t=times, y0=0, args=tuple([[X]]))
yobs = np.random.normal(y,1)
fig, ax = plt.subplots(dpi=120)
plt.plot(times, yobs, label='observed data', linestyle='dashed', marker='o', color='red')
plt.plot(times, y, label='true function', color='k', alpha=0.5)
plt.legend()
plt.xlabel('Time (Seconds)')
plt.ylabel(r'$y(t)$');
plt.show()
```