Trying to fit parameters for an ODE model to real data using lsqcurvefit

Hi! I'm trying to fit parameters for an ODE model to real data using lsqcurvefit. The issue I'm having is that lsqcurvefit always returns the initial guess. I've looked at previous forum posts similar to this problem and found a function to use as an example (Igor_Moura.m) and I can't seem to find any real differences between the code I've constructed and that example. I think the issue may be a result of the parameter not being able to change once fed into the ODE function, but I'm not sure how to fix it if that's the case. (My file is called SIRpvectorattempt.m)

function SIRpvectorattempt





%define function so solve for data
function G = SIRsolve(pvector,t)
x0=[0,5.704*10^6];
%ODE solver
[T,x]= ode45(@ODESIRlsq,t,x0);

     function dx = ODESIRlsq(t,x)
        %x(1)=i(t), x(2)=s(t), pvector(1)=lambda, pvector(2)=mu
        %define ODE's
        dx=zeros(2,1);
        dx(1)=pvector(1).*x(2).*x(1)-pvector(2).*x(1);
        dx(2)=-pvector(1).*x(2).*x(1); 
     end
%infected and susceptible populations over timeframe
G=x;  %  (:,1:2);
end


t=[1;2;3;4;5;6;7;8;9;10];

%real data from xcel file
ydata=zeros(10,2);
ydata(:,1) = xlsread('Covid Case Data Sorted SIR with columns for singapore','Singapore','H2:H11'); %453 itreal
ydata(:,2) = xlsread('Covid Case Data Sorted SIR with columns for singapore','Singapore','I2:I11'); %streal


pvecguess = [1;1];


pvector = lsqcurvefit(@SIRsolve,pvecguess,t,ydata)

end

Answer 1

[Star Strider]

‘... I can't seem to find any real differences between the code I've constructed and that example.’

I agree that the initial parameter estimates are likely the problem. If you have the Global Optimization Toolbox, this version of that same code (attached), using the ga function to estimate the parameters, will likely do what you want. (This will require a few tweaks for your current code to work with ga.) It may take a few runs until ga discovers a suitable set of parameters, however it will likly succeed.

Note that if you want to use lsqcurvefit to produce results that could be used with nlparci to estimate the parameter confidence intervals (and to slightly tweak the paramters themselves), use the parameter estimates that ga produces as the initial parameter estimates for lsqcurvefit.

[Brian Voke]

Hi Star Strider!

I managed to format the code to use ga and as I was testing the error between the result from that and my data I realized what the problem was. My initial condition for the ode system was actually a steady state solution so no matter what parameter was passed into the system the result was the same, that's why the initial condition was being returned. I made the correction and tried to run the lsqcurvefit script, but I waited for over an hour and it still hadn't finished running (that was only with 11 real data points out of 452 available). Do you know why the script might be taking so long or if there's anyway to fix it. It's currently the same as the one I posted in my question I only changed the first input to the ode intial condition to be a 1 instead of a zero.

[Star Strider]

Usually situations such as that are related to the differential equation system being ‘stiff’, with parameters of the differential equations being separated by several orders-of-magnitude. The usual solution to the problem is to use a stiff solver, such as ode15s, rather than ode45. (The ode45 function will eventually integrate it over several hours, however ode15s will integrate it much more quickly, because of the difference in the algorithms.) The calling syntax is the same (in most instances), so simply change the name from ode45 to ode15s, and that problem should resolve itself.

If that does not solve the problem, the ga solution will likely work. One advantage of ga is that if the differential equation system begins to encounter singularities and the integration stops, while lsqcurvefit would throw errors related to the output matrix not having the same dimensions as the data matrix, and stop, ga will simply start over with a new set of parameters, and continue until it converges on some likely vector of parameters.

[Brian Voke]

Changing the ode solver to ode15s worked, thanks so much!

[Star Strider]

As always, my pleasure!