curve fitting to an ordinary differential equation

Question

I want to fit a curve to the ODE shown below-

dA/dt = k1*profit + k2

I have the observed time series of the variables A and profit and I would like to get the optimal values of k1 and k2 using a curve fitting technique in python. I am able to write the code below for that, but the solutions do not fit well, or maybe my approach is wrong.

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

def fitfunc(t, k1, k2):
    'Function that returns A computed from an ODE for k1 and k2'
    def myode(area, t):
        profit = get_profit(t)
        return k1*profit + k2

    A0 = 200000    #initial value of A
    out = odeint(myode, A0, t)
    return out[:,0]

k_fit, kcov = curve_fit(fitfunc, time_span, obs_A)  #time span from 1999-2019 and obs_A is the observed values of A

modeled_A = fitfunc(time_span, k_fit[0], k_fit[1])

The profit and obs_A data for the 20 year period is:

profit = [ 7.65976374e+06, -6.13172279e+06,  1.03946093e+07,  2.59937877e+06,
       -7.88358386e+06, -1.38918115e+04, -3.13403157e+06, -4.74348806e+06,
        1.87296164e+07,  4.13680709e+07, -1.77191198e+07,  2.39249499e+06,
        1.38521564e+07,  6.52548348e+07, -5.78102494e+07, -5.72469988e+07,
       -5.99056006e+06, -1.72424523e+07,  1.78509987e+07,  9.27860105e+06,
       -9.96709853e+06]

obs_A = [200000., 165000., 150000., 180000., 190000., 195000., 200000.,
       165000., 280000., 235000., 250000., 250000., 250000., 295000.,
       295000., 285000., 245000., 315000., 235000., 245000., 305000.]

time_span = np.arange(1999,2020)

Here get_profit is a function that outputs the value of profit at a given t, its created using interpolating the observed profit data, as below-

profit_fun = interp1d(t, profit.values, 1, fill_value="extrapolate")

def get_profit(t):
    return profit_fin(t)

I am not sure about how to use the profit variables here as it changes at each time step. Is my approach correct?

Could you provide data you're trying to fit, or at least some toy data that reproduces your behavior? — K.Cl, Jul 19 '22 at 19:49
I got the magnitude of the model closer to the data by supplying initial values to k1 and k2 of 1E-2 and 1E4, but the shape of the curve (looks like two joined peaks) is practically unaffected by whatever parameter I choose. That looked like a problem in the model, so I removed the interpolation of `t` (which I frankly don't know why it's useful, but I don't know ODEs), and the result was a bit better I guess, but it looked like a smooth curve. I don't know what to do from here, sorry. — K.Cl, Jul 20 '22 at 11:52
@K.Cl could you please post your code in the answer? Did you remove the interpolation or profit? — lsr729, Jul 20 '22 at 15:58
Ok, I posted it, but I don't know if it answers your question. Hope it helps guide you. — K.Cl, Jul 20 '22 at 16:22

score 1 · Answer 1 · answered Jul 20 '22 at 16:21

(As requested, here's the code)

First, setting things up. Only added fitfun2, which modifies fitfunc removing the call to get_profit (and consequently doesn't interpolate the data).

import numpy as np
from scipy.optimize import curve_fit
from scipy.integrate import odeint
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

def fitfunc(t, k1, k2):  # Original
    'Function that returns A computed from an ODE for k1 and k2'
    def myode(area, t):
        profit = get_profit(t)
        return k1*profit + k2

    A0 = 20000    #initial value of A
    out = odeint(myode, A0, t)
    return out[:,0]

def fitfunc2(t, k1, k2): # Modified
    'Modified fitfunc, removing the call to `profit_fun`'
    def myode(area, t):
        return k1*t+k2

    A0 = 20000    #initial value of A
    out = odeint(myode, A0, t)
    return out[:,0]

profit = np.array([ 7.65976374e+06, -6.13172279e+06,  1.03946093e+07,  2.59937877e+06,
       -7.88358386e+06, -1.38918115e+04, -3.13403157e+06, -4.74348806e+06,
        1.87296164e+07,  4.13680709e+07, -1.77191198e+07,  2.39249499e+06,
        1.38521564e+07,  6.52548348e+07, -5.78102494e+07, -5.72469988e+07,
       -5.99056006e+06, -1.72424523e+07,  1.78509987e+07,  9.27860105e+06,
       -9.96709853e+06])

obs_A = np.array([200000., 165000., 150000., 180000., 190000., 195000., 200000.,
       165000., 280000., 235000., 250000., 250000., 250000., 295000.,
       295000., 285000., 245000., 315000., 235000., 245000., 305000.])

time_span = np.arange(1999,2020)

profit_fun = interp1d(time_span, profit, 1, fill_value="extrapolate")

def get_profit(t):
    return profit_fun(t)

Now, fitting and plotting the results

p0 = (1E-2, 1E4)
k_fit, kcov = curve_fit(fitfunc, time_span, obs_A, p0=p0)
k_fit2, kcov2 = curve_fit(fitfunc2, time_span, obs_A, p0=p0)

modeled_A = fitfunc(time_span, *k_fit)
guess_A = fitfunc(time_span, *p0)

modeled_A2 = fitfunc2(time_span, *k_fit2)
guess_A2 = fitfunc2(time_span, *p0)

plt.plot(time_span, obs_A, marker='o', lw=0, label='data')
plt.plot(time_span, modeled_A, label='model A (original)')
plt.plot(time_span, modeled_A2, label='model A (modified)')
plt.plot(time_span, guess_A, label='initial guess (original)')
plt.plot(time_span, guess_A2, label='initial guess (modified)')

plt.legend()

This is the graph:

As I mentioned, modifying the parameters k does not affect the curve shape of the original model. It still looks kinda 'stepwise'. Removing the call to get_profit, the curve becomes much smoother, but I don't know if it's what @Isr729 expected.

Thanks, it doesn't answer my question, but its helpful. – lsr729 Jul 20 '22 at 16:43 — lsr729, Jul 20 '22 at 16:43

curve fitting to an ordinary differential equation

1 Answers1