Error using nlinfit (Check FunVals)
5 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi there.. i am trying to fit my odel equations to some data but i keep on getiing error.. please help
below is my code
% Data fitting for the covid-19 outbreak in Lagos
covid19_nlinfit
function covid19_nlinfit
clear all
clc
global a x1 beta_c ep_g c_g sigma nu theta psi d_o d_d gamma_i alpha ...
delta ep_s ep_i ep_a E0 A0 I0
ep_g = 0.9; c_g = 0.9; delta = 0; ep = 0; ep_s = ep; ep_a=ep; ep_i=ep;
sigma = 1/5.2; nu = 0.5;
gamma_i = 1/15;
gamma_a = 0.13978; gamma_o = 0.13978;
d_o = 0.015; d_d = 0.015;
alpha = 0.5;
%==========================================================================
% Lagos data, starting from March 16, 2020 till May 3, 2020.
day = [1 3 4 6 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 ...
43 44 45 46 47 48];
Totalcases = [1 5 9 19 22 25 29 32 44 59 68 81 82 91 98 109 109 120 120 130 145 158 166 174 176 189 214 232 251 283 ...
306 376 376 430 504 582 657 689 731 764 844 931 976 1006 1068];
Activecases = [1 5 8 19 22 25 29 32 43 58 67 75 76 85 81 86 86 91 87 96 111 117 118 119 116 123 139 140 154 182 200 ...
266 266 309 383 460 525 548 584 602 682 718 756 753 791];
%==========================================================================
%Initial guesses for parameters to be estimated
beta_c=0; theta=0;psi=0;
% E0 = 16; A0 = 20; I0 = 10;
%======================================================================
% Using NLINFIT function
a = [beta_c theta psi]; % E0 A0 I0;
x1 = day;
options = statset('Display', 'iter');
[values] = nlinfit(day, Totalcases, @covidinfection, [a(1) a(2) a(3)], options);
%[values] = nlinfit(day, Activecases, @covidinfection, [a(1) a(2) a(3) a(4) a(5) a(6)], options);
beta_c = values(1)
theta = values(2)
psi = values(3)
% E0 = values(4)
% A0 = values(5)
% I0 = values(6)
%=====================================================================
%Control Reproduction Number
R_c = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s)*((alpha*nu*(1-ep_a)/...
(theta+gamma_a)) + ((1-nu)*(1-ep_i)/(psi+d_o+gamma_o)));
[t,y] = ode15s(@covidmodel1, [0,50],[14368332,16,20,10,1,0,0]); %E0 A0 I0;
plot(t,y(:,7),'-b',day,Totalcases,'rs')
xlabel('Time (days)'),ylabel('Cumulative number of reported cases')
legend('Model','Lagos COVID-19 Data')
%plot(t,y(:,5),'-b',day,Activecases,'rs')
%========================================================================
% Using NLINFIT function for data fitting
function y = covidinfection(a,x1)
beta_c = a(1);
theta = a(2);
psi = a(3);
% E0 = a(4);
% A0 = a(5);
% I0 = a(6);
[t,y] = ode15s(@covidmodel1, [0,x1(end)],[14368332,16,20,10,1,0,0]); %E0 A0 I0;
y = interp1(t,y(:,7),x1); % cumulative reported cases
%y = interp1(t,y(:,5),x1); % active cases
end
%=========================================================================
function dx = covidmodel1(t,x)
dx = [0;0;0;0;0;0;0];
N = x(1)+x(2)+x(3)+x(4)+x(5)+x(6);
beta = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s);
dx(1) = -beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/N;
dx(2) = beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/N - sigma*x(2);
dx(3) = nu*sigma*x(2) - theta*x(3) - gamma_a*x(3);
dx(4) = (1-nu)*sigma*x(2) - psi*x(4) - d_o*x(4) - gamma_o*x(4);
dx(5) = theta*x(3) + psi*x(4) - gamma_i*x(5) -d_d*x(5); % Active cases, detected (asymptomatic and symptomatic)
dx(6) = gamma_i*x(5) + gamma_a*x(3) + gamma_o*x(4);
dx(7) = theta*x(3) + psi*x(4); % Cumulative number of new cases reported (asymptomatic and symptomatic)
end
end
0 commentaires
Réponses (1)
Matt J
le 27 Août 2024
Modifié(e) : Matt J
le 27 Août 2024
For one thing,
y = interp1(t,y(:,7),x1,'spline','extrap'); % cumulative reported cases
For another, some sort of bounds or constraints on the variables are needed to prevent the ODE solver from considering ill-behaved systmes of equations. nlinfit does not let you impose such constraints, so you may have to look to another solver, e.g. lsqcurvefit.
0 commentaires
Voir également
Catégories
En savoir plus sur Biological and Health Sciences dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!