Parameter estimation of SEIHRDBP model
5 views (last 30 days)
Show older comments
% 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)
0 Comments
Accepted Answer
Star Strider
on 2 Mar 2024
Edited: Star Strider
on 2 Mar 2024
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')
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)));
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')
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 Comments
More Answers (0)
See Also
Categories
Find more on Geometry and Mesh in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!