Fitting data to numerical solution of an ode in python

Question

I have a system of two first order ODEs, which are nonlinear, and hence difficult to solve analytically in a closed form. I want to fit the numerical solution to this system of ODEs to a data set. My data set is for only one of the two variables that are part of the ODE system. How do I go about this? This didn't help because there's only one variable there.

My code which is currently leading to an error is:

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import curve_fit

def f(y, t, a, b, g):
    S, I = y # S, I are supposed to be my variables
    Sdot = -a * S * I
    Idot = (a - b) * S * I + (b - g - b * I) * I
    dydt = [Sdot, Idot]
    return dydt

def y(t, a, b, g, y0):
    y = odeint(f, y0, t, args=(a, b, g))
    return y.ravel()

I_data =[] # I have data only for I, not for S
file = open('./ratings_showdown.csv')
for e_raw in file.read().split('\r\n'):
    try:
        e=float(e_raw); I_data.append(e)
    except ValueError:
        continue

data_t = range(len(I_data))
popt, cov = curve_fit(y, data_t, I_data, [.05, 0.02, 0.01, [0.99,0.01]]) 
#want to fit I part of solution to data for variable I
#ERROR here, ValueError: setting an array element with a sequence
a_opt, b_opt, g_opt, y0_opt = popt

print("a = %g" % a_opt)
print("b = %g" % b_opt)
print("g = %g" % g_opt)
print("y0 = %g" % y0_opt)

import matplotlib.pyplot as plt
t = np.linspace(0, len(data_y), 2000)
plt.plot(data_t, data_y, '.',
         t, y(t, a_opt, b_opt, g_opt, y0_opt), '-')
plt.gcf().set_size_inches(6, 4)
#plt.savefig('out.png', dpi=96) #to save the fit result
plt.show()

Answer 1

This type of ODE fitting becomes a lot easier in symfit, which I wrote specifically as a user friendly wrapper to scipy. I think it will be very useful for your situation because the decreased amount of boiler-plate code simplifies things a lot.

From the docs and applied roughly to your problem:

from symfit import variables, parameters, Fit, D, ODEModel

S, I, t = variables('S, I, t')
a, b, g = parameters('a, b, g')

model_dict = {
    D(S, t): -a * S * I,
    D(I, t): (a - b) * S * I + (b - g - b * I) * I
}

ode_model = ODEModel(model_dict, initial={t: 0.0, S: 0.99, I: 0.01})

fit = Fit(ode_model, t=tdata, I=I_data, S=None)
fit_result = fit.execute()

Check out the docs for more :)

Answer 2

So I figured out the problem. The curve_fit() function apparently returns a list as it's second return value. So, instead of passing the initial conditions as a list [0.99,0.01], I passed them separately as 0.99 and 0.01.