how to get a good estimate of the positive parameters that will give a good fit of the curves to real data?
clear
close all
clc
% Données spécifiques
specific_data = [
2009 2 8;
2010 10 22;
2011 30 45;
2012 111 75;
2013 125 96;
2014 255 192;
2015 379 227;
2016 384 238
2017 360 279;
2018 399 229;
2019 235 128
];
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
HSdata = specific_data(:, 3);
tforward = 2009:1:2019;
%tmeasure = [ 1:100:1001]’;
% initial values
gamma = 1.5;
phi_S = 0.0006; % transmission prob
phi_H = 0.000051; % trans proba
c1=3;
c2=1.5;
theta1 = 100; % djustment parameters for syph
theta2 = 4; %djustment parameters for hi
alpha = 0.6; % progression rate
beta = 0.2; %Complications rate
rho1 = 0.4; % adjustment parameters
rho2 = 1.5; % adjustment parameters
rho3 = 1.5; % adjustment parameters
k0 = [phi_S phi_H c1 c2 theta1 theta2 alpha beta rho1 rho2 rho3 ];
% solve equ with initial value of parameters
[t, Y] = ode23s(@(t, y)modelhs(t, y, k0), tforward, [ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ]);
yintM = Y(:,1);
yintS = Y(:,2);
yintH1 = Y(:,3);
yintH2 = Y(:,4);
yintH1S=Y(:,5);
yintH2S=Y(:,6);
Hh=phi_H * ( yintH1 + c1 * yintH1S ).*yintM + rho1 * gamma * yintH1S + alpha* yintH2; % new cases from monoH
HShs=theta1*phi_S * ( yintS + c2 * yintH1S ).*yintH1 + theta2*phi_H * ( yintH1 + c1 * yintH1S ).*yintS+ rho2*alpha*yintH1S; %new case of co h-s
%H2q = Y(:,4);% assignts the y-coordinates of .
% Plotting specific data and solutions
% Display the results
figure(1)
%subplot(1,2,1);
plot(tdata, Hdata, ‘r*’);
hold on
plot(tdata, Hh, ‘b-‘);
xlabel(‘time in days’);
ylabel(‘Number of monhiv cases’);
axis([2009 2019 0 500]);
figure(2)
%subplot(1,2,1);
plot(tdata, HSdata, ‘r*’);
hold on
plot(tdata, HShs, ‘b-‘);
xlabel(‘time in days’);
ylabel(‘Number of Coinfection cases’);
axis([2009 2019 0 500]);
% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval
% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval
[k,fval] = fminsearch(@moderHS,k0)
%print final values of alpha and beta
disp(k);
%Draw the data with the final ODE
[T, Y] = ode45(@(t,y)(modelhs(t,y,k)),tforward,[ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ]);
yintM = Y(:,1);
yintS = Y(:,2);
yintH1 = Y(:,3);
yintH2 = Y(:,4);
yintH1S=Y(:,5);
yintH2S=Y(:,6);
Hh=phi_H * ( yintH1 + c1 * yintH1S ).*yintM + rho1 * gamma * yintH1S + alpha* yintH2; % new cases from mono-HIV
HShs=theta1*phi_S * ( yintS + c2 * yintH1S ).*yintH1 + theta2*phi_H * ( yintH1 + c1 * yintH1S ).*yintS+ rho2*alpha*yintH1S; %new case of coinfection hiv+syphilis
%H2q = Y(:,4);% assignts the y-coordinates of …
residuals = (Hdata+HSdata – Hh-HShs)./2;
%subplot(1,2,2);
figure(3)
plot(tdata,Hdata,’r*’);
hold on
plot(tdata,Hh,’b-‘);
xlabel(‘Time in days’);
ylabel(‘Number of mono-HIV cases’);
axis([2009 2019 0 1000]);
legend(‘Data’, ‘Model estimation’);
figure(4)
plot(tdata,HSdata,’r*’);
hold on
plot(tdata,HShs,’b-‘);
xlabel(‘Time in days’);
ylabel(‘Number of Coinfection cases’);
axis([2009 2019 0 1000]);
legend(‘Data’, ‘Model estimation’);
function dy=modelhs(~,y,k)
delta = 0.01; % Taux de mortalité
delta_S = 0.05; % Taux de mort de Syphilis.
delta_H = 0.4;
Lambda =4.04 *100;
gamma=1.5;
phi_S =k(1);
phi_H =k(2);
c1 = k(3);
c2=k(4) ;
theta1 =k(5) ;
theta2 = k(6);
alpha = k(7);
beta = k(8);
rho1 = k(9);
rho2= k(10);
rho3=k(11);
dy = zeros(7,1);
%lambda_s=phi_S * ( y(2) + c2 * y(5))
%lambda_H= phi_H * ( y(3) + c1 * y(5) )
dy(1) = Lambda + gamma * y(2) – (phi_S * ( y(2) + c2 * y(5)) + phi_H * (y(3) + c1 * y(5)) + delta ) * y(1) ;%M
dy(2)= phi_S * ( y(2) + c2 * y(5) ) * y(1) – ( gamma + theta2 * phi_H * ( y(3) + c1 * y(5) ) + delta_S ) * y(2) ;%S
dy(3) = phi_H * ( y(3) + c1 * y(5) ) * y(1) + rho1 * gamma * y(5) – (theta1 * phi_S * ( y(2) + c2 * y(5)) + delta + delta_H + alpha) * y(3) ;%H1
dy(4) = alpha * y(3) – (beta + delta + delta_H) * y(4) ;%H2
dy(5) = theta2 * phi_H * ( y(3) + c1 * y(5) ) * y(2) + ( theta1 * phi_S * ( y(2) + c2 * y(5) ) ) * y(3) – ( rho1 * gamma + rho2 * alpha + delta_H + delta_S + delta ) * y(5) ;%H1S
dy(6)= rho2 * alpha * y(5) – ( rho3 * beta + delta_S + delta_H + delta ) * y(6) ;%H2S
dy(7)= beta * y(4) + rho3 * beta * y(6) – ( delta + delta_H ) * y(7) ;%C
end
function error_in_data = moderHS(k)
% Données spécifiques
specific_data = [
2009 2 8;
2010 10 22;
2011 30 45;
2012 111 75;
2013 125 96;
2014 255 192;
2015 379 227;
2016 384 238
2017 360 279;
2018 399 229;
2019 235 128
];
gamma = 1.5;
phi_S =k(1);
phi_H =k(2);
c1 = k(3);
c2=k(4) ;
theta1 =k(5) ;
theta2 = k(6);
alpha = k(7);
beta = k(8);
rho1 = k(9);
rho2= k(10);
rho3=k(11);
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
HSdata = specific_data(:, 3);
tforward = 2009:1:2019;
[T, Y] = ode23s(@(t,y)(modelhs(t,y,k)),tforward,[ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ]);
M = Y(:,1);
S = Y(:,2);
H1 = Y(:,3);
H2 = Y(:,4);
H1S=Y(:,5);
H2S=Y(:,6);
H=phi_H * ( H1 + c1 * H1S ).*M + rho1 * gamma*H1S + alpha* H2; % new cases from mono-HIV
HS=theta1*phi_S * ( S + c2 * H1S ).*H1 + theta2*phi_H * ( H1 + c1 * H1S ).*S+ rho2*alpha*H1S; %new case of coinfection hiv+syphilis
%H2q = Y(:,4);% assignts the y-coordinates of …
%the solution at
D=mean(Hdata).^2;
D1=mean(HSdata).^2;
A=(H – Hdata).^2;
B=(HS – HSdata).^2;
error_in_data =sum(A)./(11*D)+ sum(B)./(11*D1);
%%
endclear
close all
clc
% Données spécifiques
specific_data = [
2009 2 8;
2010 10 22;
2011 30 45;
2012 111 75;
2013 125 96;
2014 255 192;
2015 379 227;
2016 384 238
2017 360 279;
2018 399 229;
2019 235 128
];
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
HSdata = specific_data(:, 3);
tforward = 2009:1:2019;
%tmeasure = [ 1:100:1001]’;
% initial values
gamma = 1.5;
phi_S = 0.0006; % transmission prob
phi_H = 0.000051; % trans proba
c1=3;
c2=1.5;
theta1 = 100; % djustment parameters for syph
theta2 = 4; %djustment parameters for hi
alpha = 0.6; % progression rate
beta = 0.2; %Complications rate
rho1 = 0.4; % adjustment parameters
rho2 = 1.5; % adjustment parameters
rho3 = 1.5; % adjustment parameters
k0 = [phi_S phi_H c1 c2 theta1 theta2 alpha beta rho1 rho2 rho3 ];
% solve equ with initial value of parameters
[t, Y] = ode23s(@(t, y)modelhs(t, y, k0), tforward, [ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ]);
yintM = Y(:,1);
yintS = Y(:,2);
yintH1 = Y(:,3);
yintH2 = Y(:,4);
yintH1S=Y(:,5);
yintH2S=Y(:,6);
Hh=phi_H * ( yintH1 + c1 * yintH1S ).*yintM + rho1 * gamma * yintH1S + alpha* yintH2; % new cases from monoH
HShs=theta1*phi_S * ( yintS + c2 * yintH1S ).*yintH1 + theta2*phi_H * ( yintH1 + c1 * yintH1S ).*yintS+ rho2*alpha*yintH1S; %new case of co h-s
%H2q = Y(:,4);% assignts the y-coordinates of .
% Plotting specific data and solutions
% Display the results
figure(1)
%subplot(1,2,1);
plot(tdata, Hdata, ‘r*’);
hold on
plot(tdata, Hh, ‘b-‘);
xlabel(‘time in days’);
ylabel(‘Number of monhiv cases’);
axis([2009 2019 0 500]);
figure(2)
%subplot(1,2,1);
plot(tdata, HSdata, ‘r*’);
hold on
plot(tdata, HShs, ‘b-‘);
xlabel(‘time in days’);
ylabel(‘Number of Coinfection cases’);
axis([2009 2019 0 500]);
% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval
% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval
[k,fval] = fminsearch(@moderHS,k0)
%print final values of alpha and beta
disp(k);
%Draw the data with the final ODE
[T, Y] = ode45(@(t,y)(modelhs(t,y,k)),tforward,[ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ]);
yintM = Y(:,1);
yintS = Y(:,2);
yintH1 = Y(:,3);
yintH2 = Y(:,4);
yintH1S=Y(:,5);
yintH2S=Y(:,6);
Hh=phi_H * ( yintH1 + c1 * yintH1S ).*yintM + rho1 * gamma * yintH1S + alpha* yintH2; % new cases from mono-HIV
HShs=theta1*phi_S * ( yintS + c2 * yintH1S ).*yintH1 + theta2*phi_H * ( yintH1 + c1 * yintH1S ).*yintS+ rho2*alpha*yintH1S; %new case of coinfection hiv+syphilis
%H2q = Y(:,4);% assignts the y-coordinates of …
residuals = (Hdata+HSdata – Hh-HShs)./2;
%subplot(1,2,2);
figure(3)
plot(tdata,Hdata,’r*’);
hold on
plot(tdata,Hh,’b-‘);
xlabel(‘Time in days’);
ylabel(‘Number of mono-HIV cases’);
axis([2009 2019 0 1000]);
legend(‘Data’, ‘Model estimation’);
figure(4)
plot(tdata,HSdata,’r*’);
hold on
plot(tdata,HShs,’b-‘);
xlabel(‘Time in days’);
ylabel(‘Number of Coinfection cases’);
axis([2009 2019 0 1000]);
legend(‘Data’, ‘Model estimation’);
function dy=modelhs(~,y,k)
delta = 0.01; % Taux de mortalité
delta_S = 0.05; % Taux de mort de Syphilis.
delta_H = 0.4;
Lambda =4.04 *100;
gamma=1.5;
phi_S =k(1);
phi_H =k(2);
c1 = k(3);
c2=k(4) ;
theta1 =k(5) ;
theta2 = k(6);
alpha = k(7);
beta = k(8);
rho1 = k(9);
rho2= k(10);
rho3=k(11);
dy = zeros(7,1);
%lambda_s=phi_S * ( y(2) + c2 * y(5))
%lambda_H= phi_H * ( y(3) + c1 * y(5) )
dy(1) = Lambda + gamma * y(2) – (phi_S * ( y(2) + c2 * y(5)) + phi_H * (y(3) + c1 * y(5)) + delta ) * y(1) ;%M
dy(2)= phi_S * ( y(2) + c2 * y(5) ) * y(1) – ( gamma + theta2 * phi_H * ( y(3) + c1 * y(5) ) + delta_S ) * y(2) ;%S
dy(3) = phi_H * ( y(3) + c1 * y(5) ) * y(1) + rho1 * gamma * y(5) – (theta1 * phi_S * ( y(2) + c2 * y(5)) + delta + delta_H + alpha) * y(3) ;%H1
dy(4) = alpha * y(3) – (beta + delta + delta_H) * y(4) ;%H2
dy(5) = theta2 * phi_H * ( y(3) + c1 * y(5) ) * y(2) + ( theta1 * phi_S * ( y(2) + c2 * y(5) ) ) * y(3) – ( rho1 * gamma + rho2 * alpha + delta_H + delta_S + delta ) * y(5) ;%H1S
dy(6)= rho2 * alpha * y(5) – ( rho3 * beta + delta_S + delta_H + delta ) * y(6) ;%H2S
dy(7)= beta * y(4) + rho3 * beta * y(6) – ( delta + delta_H ) * y(7) ;%C
end
function error_in_data = moderHS(k)
% Données spécifiques
specific_data = [
2009 2 8;
2010 10 22;
2011 30 45;
2012 111 75;
2013 125 96;
2014 255 192;
2015 379 227;
2016 384 238
2017 360 279;
2018 399 229;
2019 235 128
];
gamma = 1.5;
phi_S =k(1);
phi_H =k(2);
c1 = k(3);
c2=k(4) ;
theta1 =k(5) ;
theta2 = k(6);
alpha = k(7);
beta = k(8);
rho1 = k(9);
rho2= k(10);
rho3=k(11);
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
HSdata = specific_data(:, 3);
tforward = 2009:1:2019;
[T, Y] = ode23s(@(t,y)(modelhs(t,y,k)),tforward,[ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ]);
M = Y(:,1);
S = Y(:,2);
H1 = Y(:,3);
H2 = Y(:,4);
H1S=Y(:,5);
H2S=Y(:,6);
H=phi_H * ( H1 + c1 * H1S ).*M + rho1 * gamma*H1S + alpha* H2; % new cases from mono-HIV
HS=theta1*phi_S * ( S + c2 * H1S ).*H1 + theta2*phi_H * ( H1 + c1 * H1S ).*S+ rho2*alpha*H1S; %new case of coinfection hiv+syphilis
%H2q = Y(:,4);% assignts the y-coordinates of …
%the solution at
D=mean(Hdata).^2;
D1=mean(HSdata).^2;
A=(H – Hdata).^2;
B=(HS – HSdata).^2;
error_in_data =sum(A)./(11*D)+ sum(B)./(11*D1);
%%
end clear
close all
clc
% Données spécifiques
specific_data = [
2009 2 8;
2010 10 22;
2011 30 45;
2012 111 75;
2013 125 96;
2014 255 192;
2015 379 227;
2016 384 238
2017 360 279;
2018 399 229;
2019 235 128
];
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
HSdata = specific_data(:, 3);
tforward = 2009:1:2019;
%tmeasure = [ 1:100:1001]’;
% initial values
gamma = 1.5;
phi_S = 0.0006; % transmission prob
phi_H = 0.000051; % trans proba
c1=3;
c2=1.5;
theta1 = 100; % djustment parameters for syph
theta2 = 4; %djustment parameters for hi
alpha = 0.6; % progression rate
beta = 0.2; %Complications rate
rho1 = 0.4; % adjustment parameters
rho2 = 1.5; % adjustment parameters
rho3 = 1.5; % adjustment parameters
k0 = [phi_S phi_H c1 c2 theta1 theta2 alpha beta rho1 rho2 rho3 ];
% solve equ with initial value of parameters
[t, Y] = ode23s(@(t, y)modelhs(t, y, k0), tforward, [ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ]);
yintM = Y(:,1);
yintS = Y(:,2);
yintH1 = Y(:,3);
yintH2 = Y(:,4);
yintH1S=Y(:,5);
yintH2S=Y(:,6);
Hh=phi_H * ( yintH1 + c1 * yintH1S ).*yintM + rho1 * gamma * yintH1S + alpha* yintH2; % new cases from monoH
HShs=theta1*phi_S * ( yintS + c2 * yintH1S ).*yintH1 + theta2*phi_H * ( yintH1 + c1 * yintH1S ).*yintS+ rho2*alpha*yintH1S; %new case of co h-s
%H2q = Y(:,4);% assignts the y-coordinates of .
% Plotting specific data and solutions
% Display the results
figure(1)
%subplot(1,2,1);
plot(tdata, Hdata, ‘r*’);
hold on
plot(tdata, Hh, ‘b-‘);
xlabel(‘time in days’);
ylabel(‘Number of monhiv cases’);
axis([2009 2019 0 500]);
figure(2)
%subplot(1,2,1);
plot(tdata, HSdata, ‘r*’);
hold on
plot(tdata, HShs, ‘b-‘);
xlabel(‘time in days’);
ylabel(‘Number of Coinfection cases’);
axis([2009 2019 0 500]);
% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval
% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval
[k,fval] = fminsearch(@moderHS,k0)
%print final values of alpha and beta
disp(k);
%Draw the data with the final ODE
[T, Y] = ode45(@(t,y)(modelhs(t,y,k)),tforward,[ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ]);
yintM = Y(:,1);
yintS = Y(:,2);
yintH1 = Y(:,3);
yintH2 = Y(:,4);
yintH1S=Y(:,5);
yintH2S=Y(:,6);
Hh=phi_H * ( yintH1 + c1 * yintH1S ).*yintM + rho1 * gamma * yintH1S + alpha* yintH2; % new cases from mono-HIV
HShs=theta1*phi_S * ( yintS + c2 * yintH1S ).*yintH1 + theta2*phi_H * ( yintH1 + c1 * yintH1S ).*yintS+ rho2*alpha*yintH1S; %new case of coinfection hiv+syphilis
%H2q = Y(:,4);% assignts the y-coordinates of …
residuals = (Hdata+HSdata – Hh-HShs)./2;
%subplot(1,2,2);
figure(3)
plot(tdata,Hdata,’r*’);
hold on
plot(tdata,Hh,’b-‘);
xlabel(‘Time in days’);
ylabel(‘Number of mono-HIV cases’);
axis([2009 2019 0 1000]);
legend(‘Data’, ‘Model estimation’);
figure(4)
plot(tdata,HSdata,’r*’);
hold on
plot(tdata,HShs,’b-‘);
xlabel(‘Time in days’);
ylabel(‘Number of Coinfection cases’);
axis([2009 2019 0 1000]);
legend(‘Data’, ‘Model estimation’);
function dy=modelhs(~,y,k)
delta = 0.01; % Taux de mortalité
delta_S = 0.05; % Taux de mort de Syphilis.
delta_H = 0.4;
Lambda =4.04 *100;
gamma=1.5;
phi_S =k(1);
phi_H =k(2);
c1 = k(3);
c2=k(4) ;
theta1 =k(5) ;
theta2 = k(6);
alpha = k(7);
beta = k(8);
rho1 = k(9);
rho2= k(10);
rho3=k(11);
dy = zeros(7,1);
%lambda_s=phi_S * ( y(2) + c2 * y(5))
%lambda_H= phi_H * ( y(3) + c1 * y(5) )
dy(1) = Lambda + gamma * y(2) – (phi_S * ( y(2) + c2 * y(5)) + phi_H * (y(3) + c1 * y(5)) + delta ) * y(1) ;%M
dy(2)= phi_S * ( y(2) + c2 * y(5) ) * y(1) – ( gamma + theta2 * phi_H * ( y(3) + c1 * y(5) ) + delta_S ) * y(2) ;%S
dy(3) = phi_H * ( y(3) + c1 * y(5) ) * y(1) + rho1 * gamma * y(5) – (theta1 * phi_S * ( y(2) + c2 * y(5)) + delta + delta_H + alpha) * y(3) ;%H1
dy(4) = alpha * y(3) – (beta + delta + delta_H) * y(4) ;%H2
dy(5) = theta2 * phi_H * ( y(3) + c1 * y(5) ) * y(2) + ( theta1 * phi_S * ( y(2) + c2 * y(5) ) ) * y(3) – ( rho1 * gamma + rho2 * alpha + delta_H + delta_S + delta ) * y(5) ;%H1S
dy(6)= rho2 * alpha * y(5) – ( rho3 * beta + delta_S + delta_H + delta ) * y(6) ;%H2S
dy(7)= beta * y(4) + rho3 * beta * y(6) – ( delta + delta_H ) * y(7) ;%C
end
function error_in_data = moderHS(k)
% Données spécifiques
specific_data = [
2009 2 8;
2010 10 22;
2011 30 45;
2012 111 75;
2013 125 96;
2014 255 192;
2015 379 227;
2016 384 238
2017 360 279;
2018 399 229;
2019 235 128
];
gamma = 1.5;
phi_S =k(1);
phi_H =k(2);
c1 = k(3);
c2=k(4) ;
theta1 =k(5) ;
theta2 = k(6);
alpha = k(7);
beta = k(8);
rho1 = k(9);
rho2= k(10);
rho3=k(11);
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
HSdata = specific_data(:, 3);
tforward = 2009:1:2019;
[T, Y] = ode23s(@(t,y)(modelhs(t,y,k)),tforward,[ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ]);
M = Y(:,1);
S = Y(:,2);
H1 = Y(:,3);
H2 = Y(:,4);
H1S=Y(:,5);
H2S=Y(:,6);
H=phi_H * ( H1 + c1 * H1S ).*M + rho1 * gamma*H1S + alpha* H2; % new cases from mono-HIV
HS=theta1*phi_S * ( S + c2 * H1S ).*H1 + theta2*phi_H * ( H1 + c1 * H1S ).*S+ rho2*alpha*H1S; %new case of coinfection hiv+syphilis
%H2q = Y(:,4);% assignts the y-coordinates of …
%the solution at
D=mean(Hdata).^2;
D1=mean(HSdata).^2;
A=(H – Hdata).^2;
B=(HS – HSdata).^2;
error_in_data =sum(A)./(11*D)+ sum(B)./(11*D1);
%%
end curve fitting, parameters estimation, optimization, multiple curve fitting MATLAB Answers — New Questions