How to implement MultiStart to find a good fit of two curves at the same time with real data such as our curves are governed by an ODE system?

6 次查看（过去 30 天）

Khadija 2024-8-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2142581-how-to-implement-multistart-to-find-a-good-fit-of-two-curves-at-the-same-time-with-real-data-such-as

评论： Khadija 2024-8-4

% Set up the problem for MultiStart
problem = createOptimProblem('lsqcurvefit','x0',k0,'objective',@simulatedhs,...
    'lb',zeros(size(k0)),'xdata',tforward ,'ydata',[Hdata,HSdata] );
ms = MultiStart;
[k,fval,Exitflag,Output,Solutions] = run(ms,problem,50);
simulated_data = simulatedhs(k,tforward);
X=simulated_data(:,1);
Z=simulated_data(:,2);
%plot the result
figure(5)
plot(tforward,Hdata ,tforward,X)
plot(tforward,Hdata,tforward,Y)
legend('Data','Fitted result')

13 个评论
显示 11更早的评论隐藏 11更早的评论

Torsten 2024-8-3

在 MATLAB Online 中打开

There are problems with ode15s, but MultiStart seems to work in principle. Do you use an older MATLAB version ?

% 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 ],odeset('RelTol',1e-10,'AbsTol',1e-10));

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]);

% Set up the problem for MultiStart

problem = createOptimProblem('lsqcurvefit','x0',k0,'objective',@simulatedhs,...

'lb',zeros(size(k0)),'xdata',tforward ,'ydata',[Hdata,HSdata] );

ms = MultiStart;

[k,fval,Exitflag,Output,Solutions] = run(ms,problem,5);

Warning: Failure at t=2.009000e+03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.637979e-12) at time t.

MultiStart completed some of the runs from the start points. 1 out of 5 local solver runs converged with a positive local solver exitflag.

simulated_data = simulatedhs(k,tforward);

X=simulated_data(:,1);

Z=simulated_data(:,2);

%plot the result

figure(5)

plot(tforward,Hdata ,tforward,X)

plot(tforward,HSdata,tforward,Z)

legend('Data','Fitted result')

%kstart = k0.*(1+0.01*rand(size(k0))) % Change k by a little amount

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

%k(:)-k0(:)

%simulated_data = simulatedhs(k,tforward);

%H = simulated_data(:,1);

%HS = simulated_data(:,2);

%max(H-Hh)

%max(HS-HShs)

%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 ],odeset('RelTol',1e-10,'AbsTol',1e-10));

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

Torsten 2024-8-4

在 MATLAB Online 中打开

lb = [0 0 1 1 1 1 0 0 0 1 1];
ub = [1 1 inf inf inf inf 1 1 1 inf inf];
options = optimset('MaxFunEvals',10000,'MaxIter',10000);
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@simulatedhs,k0,tforward,[Hdata,HSdata],lb,ub,options);

Khadija 2024-8-4

for now, MultiStart does not give a better result, I have not been able to solve the problem. I keep the estimate without multistart which gives a result close to the suitable one. Thank you again for your patience!!

请先登录，再进行评论。

请先登录，再回答此问题。