Parameter estimation of SEIHRDBP model

조회 수: 5 (최근 30일)
University 2024년 3월 2일
댓글: Star Strider 2024년 3월 6일
% I want simulate "parameter estimation of SEIHRDBP Model", where
% S=susceptible, E=exposed, I=infected, H=hospitalised, R=recovered,
% D=death, B=buried, P=pathogen. I tried Matlab codes from the Matlab
% community but couldn't get right. The parameters here are assumed values
% Parameter value
% Delta = params(1); mu = params(2); delta = params(3); gamma = params(4); gamma1 = params(5);
% eta = params(6); xi = params(7); d = params(8); b = params(9); phip = params(10);
% d1 = params(11); b1 = params(12); alpha = params(13); sigma = params(14);
% eta1 = params(15); beta = params(16); beta1 = params(17); beta2 = params(18); beta3 = params(19);
% betap = params(20);
% params= [10, 0.5, 0.5, 0.005, 0.12, 0.12, 0.6, 0.8, 0.5, 0, 0.04, 0.04, 0.03, 0.04, 0.04, 0.006, 0.002, 0.012, 0.0012; 0.01];
% dxdt = zeros(8, 1); % initializing the state variables
% params = zeros(20, 1); % initializing the parameters
%
% dxdt(1) = params(1) - (params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) + params(20)*x(8))*x(1)-params(2)*x(1);
% dxdt(2) =(params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) + params(20)*x(8))*x(1) - (params(2) + params(3))*x(2);
% dxdt(3) = params(3)*x(2) - (params(2) + params(4) + params(5))*x(1);
% dxdt(4) = params(5)*x(3) - (params(6) + params(7) + params(2))*x(4);
% dxdt(5) = params(7)*x(4) - params(2)*x(5);
% dxdt(6) = params(4)*x(3) + params(6)*x(4) -params(8)*x(6);
% dxdt(7) = params(8)*x(6) - params(9)*x(7);
% dxdt(8) = params(10) + params(14)*x(3) + params(15)*x(4) + params(11)*x(6) + params(12)*x(7) -params(13)*x(8);
% This is code I tried. I basically used this help site: https://uk.mathworks.com/matlabcentral/answers/1710260-need-help-with-ode-parameter-estimation
% The c1---c8 are the state variables, while the theta1---theta20 are the parameter values.
t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];
c=[0.902 0.06997 0.02463 0.00218 0.902 0.06997 0.02463 0.00218
0.8072 0.1353 0.0482 0.008192 0.8072 0.1353 0.0482 0.008192
0.6757 0.2123 0.0864 0.0289 0.6757 0.2123 0.0864 0.0289
0.5569 0.2789 0.1063 0.06233 0.5569 0.2789 0.1063 0.06233
0.4297 0.3292 0.1476 0.09756 0.4297 0.3292 0.1476 0.09756
0.3774 0.3457 0.1485 0.1255 0.3774 0.3457 0.1485 0.1255
0.2149 0.3486 0.1821 0.2526 0.2149 0.3486 0.1821 0.2526
0.141 0.3254 0.194 0.3401 0.141 0.3254 0.194 0.3401
0.04921 0.2445 0.1742 0.5277 0.04921 0.2445 0.1742 0.5277
0.0178 0.1728 0.1732 0.6323 0.0178 0.1728 0.1732 0.6323
0.006431 0.1091 0.1137 0.7702 0.006431 0.1091 0.1137 0.7702
0.002595 0.08301 0.08224 0.835 0.002595 0.08301 0.08224 0.835];
theta0 = rand(20,1);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@SEIHRDBP,theta0,t,c,zeros(size(theta0)));
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = SEIHRDBP(theta, tv);
% figure
% hd = plot(t, c, 'p');
% for k1 = 1:size(c,2)
% CV(k1,:) = hd(k1).Color;
% hd(k1).MarkerFaceColor = CV(k1,:);
% end
% hold on
% hlp = plot(tv, Cfit);
% for k1 = 1:size(c,2)
% hlp(k1).Color = CV(k1,:);
% end
% hold off
% grid
% xlabel('Time')
% ylabel('Population')
% legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')
%
function C=SEIHRDBP(theta,t)
c0=theta(end-3:end);
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(8,1);
dcdt(1) = theta(1) - (theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) + theta(20)*c(8))*c(1)-theta(2)*c(1);
dcdt(2) =(theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) + theta(20)*c(8))*c(1) - (theta(2) + theta(3))*c(2);
dcdt(3) = theta(3)*c(2) - (theta(2) + theta(4) + theta(5))*c(1);
dcdt(4) = theta(5)*c(3) - (theta(6) + theta(7) + theta(2))*c(4);
dcdt(5) = theta(7)*c(4) - theta(2)*c(5);
dcdt(6) = theta(4)*c(3) + theta(6)*c(4) -theta(8)*c(6);
dcdt(7) = theta(8)*c(6) - theta(9)*c(7);
dcdt(8) = theta(10) + theta(14)*c(3) + theta(15)*c(4) + theta(11)*c(6) + theta(12)*c(7) -theta(13)*c(8);
dC=dcdt;
end
C=Cv;
end
% Please can some help check it for me. I have attached the observed caese of all the state variables(SEIHRDBP)

댓글을 달려면 로그인하십시오.

채택된 답변

Star Strider 2024년 3월 2일
편집: Star Strider 2024년 3월 2일
Thank you for quoting my code!
The problem was here:
c0=theta(end-3:end);
use this instead:
c0=theta(end-7:end);
Since my code uses the last elements of the paramter vector for the initial conditions of the differential equations (it fits them as parameters, as well as the rate constants), that needs to be accounted for. Here there are 8 differential equations, so there must be 8 initial conditions.
Try this —
% I want simulate "parameter estimation of SEIHRDBP Model", where
% S=susceptible, E=exposed, I=infected, H=hospitalised, R=recovered,
% D=death, B=buried, P=pathogen. I tried Matlab codes from the Matlab
% community but couldn't get right. The parameters here are assumed values
% Parameter value
% Delta = params(1); mu = params(2); delta = params(3); gamma = params(4); gamma1 = params(5);
% eta = params(6); xi = params(7); d = params(8); b = params(9); phip = params(10);
% d1 = params(11); b1 = params(12); alpha = params(13); sigma = params(14);
% eta1 = params(15); beta = params(16); beta1 = params(17); beta2 = params(18); beta3 = params(19);
% betap = params(20);
% params= [10, 0.5, 0.5, 0.005, 0.12, 0.12, 0.6, 0.8, 0.5, 0, 0.04, 0.04, 0.03, 0.04, 0.04, 0.006, 0.002, 0.012, 0.0012; 0.01];
% dxdt = zeros(8, 1); % initializing the state variables
% params = zeros(20, 1); % initializing the parameters
%
% dxdt(1) = params(1) - (params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) + params(20)*x(8))*x(1)-params(2)*x(1);
% dxdt(2) =(params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) + params(20)*x(8))*x(1) - (params(2) + params(3))*x(2);
% dxdt(3) = params(3)*x(2) - (params(2) + params(4) + params(5))*x(1);
% dxdt(4) = params(5)*x(3) - (params(6) + params(7) + params(2))*x(4);
% dxdt(5) = params(7)*x(4) - params(2)*x(5);
% dxdt(6) = params(4)*x(3) + params(6)*x(4) -params(8)*x(6);
% dxdt(7) = params(8)*x(6) - params(9)*x(7);
% dxdt(8) = params(10) + params(14)*x(3) + params(15)*x(4) + params(11)*x(6) + params(12)*x(7) -params(13)*x(8);
T1 = readtable('observed.csv')
T1 = 35×8 table
susceptible exposed infected Hospitalized Recovered deaths Buried Pathogens ___________ _______ ________ ____________ _________ ______ ______ _________ 3444 323 164 164 0 55 55 0 3444 426 1 1 0 1 1 0 3444 300 5 5 0 5 5 0 3444 335 11 11 0 9 9 0 3444 386 12 12 0 6 6 0 3444 350 23 23 0 12 12 0 3444 404 130 130 0 55 55 0 3444 416 3470 3470 0 2287 2287 0 3444 449 54 54 0 33 33 0 3444 446 8 8 0 4 4 0 3444 374 66 66 0 49 49 0 34440 38646 28646 28646 0 11323 11323 0 3444 377 6 6 0 3 3 0 3444 491 36 36 0 13 13 0 3444 491 11 11 0 4 4 0 3444 437 1 1 0 1 1 0
VN = T1.Properties.VariableNames;
L = size(T1,1);
t = linspace(0, L-1, L).'; % Time Vector (NECESSARY) Units = ?
% % This is code I tried. I basically used this help site: https://uk.mathworks.com/matlabcentral/answers/1710260-need-help-with-ode-parameter-estimation
% % The c1---c8 are the state variables, while the theta1---theta20 are the parameter values.
% t=[0.1
% 0.2
% 0.4
% 0.6
% 0.8
% 1
% 1.5
% 2
% 3
% 4
% 5
% 6];
% c=[0.902 0.06997 0.02463 0.00218 0.902 0.06997 0.02463 0.00218
% 0.8072 0.1353 0.0482 0.008192 0.8072 0.1353 0.0482 0.008192
% 0.6757 0.2123 0.0864 0.0289 0.6757 0.2123 0.0864 0.0289
% 0.5569 0.2789 0.1063 0.06233 0.5569 0.2789 0.1063 0.06233
% 0.4297 0.3292 0.1476 0.09756 0.4297 0.3292 0.1476 0.09756
% 0.3774 0.3457 0.1485 0.1255 0.3774 0.3457 0.1485 0.1255
% 0.2149 0.3486 0.1821 0.2526 0.2149 0.3486 0.1821 0.2526
% 0.141 0.3254 0.194 0.3401 0.141 0.3254 0.194 0.3401
% 0.04921 0.2445 0.1742 0.5277 0.04921 0.2445 0.1742 0.5277
% 0.0178 0.1728 0.1732 0.6323 0.0178 0.1728 0.1732 0.6323
% 0.006431 0.1091 0.1137 0.7702 0.006431 0.1091 0.1137 0.7702
% 0.002595 0.08301 0.08224 0.835 0.002595 0.08301 0.08224 0.835];
c = table2array(T1);
theta0 = rand(28,1);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@SEIHRDBP,theta0,t,c,zeros(size(theta0)));
Local minimum possible. lsqcurvefit stopped because the size of the current step is less than the value of the step size tolerance.
fprintf(1,'\tRate Constants:\n')
Rate Constants:
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
Theta(1) = 2.40444 Theta(2) = 0.03125 Theta(3) = 1.41359 Theta(4) = 1.11936 Theta(5) = 1.11908 Theta(6) = 0.30996 Theta(7) = 0.06624 Theta(8) = 0.60836 Theta(9) = 0.41146 Theta(10) = 0.29405 Theta(11) = 0.29017 Theta(12) = 0.13046 Theta(13) = 1.72754 Theta(14) = 0.69723 Theta(15) = 0.38384 Theta(16) = 1.29961 Theta(17) = 0.99986 Theta(18) = 0.78319 Theta(19) = 0.54491 Theta(20) = 0.44378 Theta(21) = 0.41763 Theta(22) = 1.15486 Theta(23) = 1.41118 Theta(24) = 0.89700 Theta(25) = 0.14331 Theta(26) = 0.62479 Theta(27) = 0.24113 Theta(28) = 0.66698
tv = linspace(min(t), max(t));
Cfit = SEIHRDBP(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Population')
% legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')
legend(hlp, VN, 'Location','best')
%
function C=SEIHRDBP(theta,t)
c0=theta(end-7:end);
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(8,1);
dcdt(1) = theta(1) - (theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) + theta(20)*c(8))*c(1)-theta(2)*c(1);
dcdt(2) =(theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) + theta(20)*c(8))*c(1) - (theta(2) + theta(3))*c(2);
dcdt(3) = theta(3)*c(2) - (theta(2) + theta(4) + theta(5))*c(1);
dcdt(4) = theta(5)*c(3) - (theta(6) + theta(7) + theta(2))*c(4);
dcdt(5) = theta(7)*c(4) - theta(2)*c(5);
dcdt(6) = theta(4)*c(3) + theta(6)*c(4) -theta(8)*c(6);
dcdt(7) = theta(8)*c(6) - theta(9)*c(7);
dcdt(8) = theta(10) + theta(14)*c(3) + theta(15)*c(4) + theta(11)*c(6) + theta(12)*c(7) -theta(13)*c(8);
dC=dcdt;
end
C=Cv;
end
% Please can some help check it for me. I have attached the observed caese of all the state variables(SEIHRDBP)
That works, however it is not an escellent fit. I would keep funning it with random initial parameter vectyors until you get a good fit. Using the ga (genetic algorithm) function might produce better results.
.
댓글 수: 10이전 댓글 8개 표시이전 댓글 8개 숨기기
University 2024년 3월 6일
Thank you so much. I was about to reply you. I have gotten it.
Star Strider 2024년 3월 6일
As always, my pleasure!

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Earth, Ocean, and Atmospheric Sciences에 대해 자세히 알아보기

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by