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)

 
Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: 0.060905 

k = 1x11

0.0002 0.0000 6.2511 4.1461 449.9916 -2.9656 -0.2168 -0.2082 0.8523 0.3901 -0.1350

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval = 0.0609

%print final values of alpha and beta
disp(k);

    0.0002    0.0000    6.2511    4.1461  449.9916   -2.9656   -0.2168   -0.2082    0.8523    0.3901   -0.1350

%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

[Torsten]

What exactly is your question ?

If you want positive parameters, use a minimizer that accepts lower and upper bounds on the parameters, e.g. "lsqcurvefit" instead of "fminsearch".

If you want parameters that better fit your data, vary the initial guesses for the parameters (e.g. by using "MultiStart").

[Khadija]

Hi, I report that I am new to Matlab programming. I have two curves to fit at the same time. I can't code well to have a good fit for both?

In the input of the lsqcurvfit function, I would like to know how to store the real data and the values ​​of the variables generated by an ODE system, with positive obtained parameters?

I hope my question is clear?

Answer 1

[Torsten]

You should avoid defining the same data, parameters and variable definitions multiple times in your code (specific_data, gamma, phi_S,phi_H,H,HS,...). Pass them to functions where needed and define functions for repeating expressions.

% 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] = ode15s(@(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]);
kstart = k0.*(1+0.01*rand(size(k0)))   % Change k by a little amount

kstart = 1x11

0.0006 0.0001 3.0135 1.5061 100.6674 4.0112 0.6059 0.2012 0.4028 1.5099 1.5101

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[k,fval] = lsqcurvefit(@simulatedhs,kstart,tforward,[Hh,HShs],zeros(numel(k0),1),[],optimset('MaxFunEvals',10000,'MaxIter',10000));

Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the value of the step size tolerance.

k(:)-k0(:)

ans = 11x1

0.0000 -0.0000 0.0056 0.0101 1.4530 -0.0012 0.0736 0.0063 0.0039 -0.1543

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

simulated_data = simulatedhs(k,tforward);
H = simulated_data(:,1);
HS = simulated_data(:,2);
max(H-Hh)

ans = 0.5319

max(HS-HShs)

ans = 0.6126

figure(3)
plot(tforward.',[H,Hh])

figure(4)
plot(tforward.',[HS,HShs])

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 simulated_data = simulatedhs(k,tdata)
% 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);
[T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata,[ 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
simulated_data = [H,HS];
end

[Khadija]

Hi Mr Torsten,

Thanks for your comment, I'm trying to understand the changes you made. For the lsqcurvfit function: should the matrix [Hh, HShs] inserted as its argument match the data values?

[k,fval] = lsqcurvefit(@simulatedhs,kstart,tforward,[Hh,HShs],zeros(numel(k0),1),[],optimset('MaxFunEvals',10000,'MaxIter',10000));

[Torsten]

[Hh,HShs] are the data values that are to be matched best-possible by simulated_data = [H(k),HS(k)]; (where k is the parameter vector to be determined).

Best-possible means: sum((Hh-H(k)).^2) + sum((HShs-HS(k)).^2) should be minimal.

[Khadija]

Hi Mr Torsten,

I point out that [Hdata, HSdata] are the data values ​​for which I am trying to fit to the variables [H(k), HS(k)]; the data values ​​are represented by stars in figure 1 and 2.

I replaced the matrix [Hh, HShs] with [Hdata, HSdata], but it does not work!!

[Torsten]

What do you mean by "it does not work" ? If the found parameter vector k does not fit [Hdata, HSdata] well, you can try changing the initial guess for the k-vector or use MATLAB's "MultiStart". But this is your part of the work - I cannot help you with this.

If you have a technical problem with your code, you are of course invited to ask.

[Khadija]

Hi Mr Torsten,

I will follow your instructions to solve my adjustment problem. For the Mstart function, I will try to introduce it to my code to solve this problem. Thanks again for your patience!!

[Torsten]

It's usually hard to get a good fit for models with so many free parameters, especially if the curves are computed from differential equations.

The exactness of the solution of the differential equations is very important. So at last, I would suggest that you strengthen the tolerances for the ode integrator in function "simulatedhs":

[T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata,[ 5000.0 20.0 2.0  0.0  6.0 2.0  0.0 ],odeset('RelTol',1e-10,'AbsTol',1e-10));

instead of

[T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata,[ 5000.0 20.0 2.0  0.0  6.0 2.0  0.0 ]);

[Khadija]

Hi Mr Torsten, I followed your instruction; and it is starting to improve the fit of my curves. How do I judge if it is a good fit?

[Torsten]

Visual inspection of the solution and the data curve.

I don't think statistical methods will lead you somewhere here. If you want to try, you can change from lsqcurvefit to fitnlm.

[Khadija]

Thank you very much Mr Torsten, I will first try to understand and implement this MultiStart function well to be able to give an adjustment. All my respects!

[Khadija]

hi sir, the input arguments of fitnlm are not the same of lsqcurvefit?

[Torsten]

I have not yet used this solver. Whether you use nlinfit, lsqcurvefit, fmincon, fitnlm, ga - the changes to substitute one solver by another should be minimal.