SEIR model for ebola

3 次查看（过去 30 天）
Arjun
Arjun 2024-6-18
回答： Kautuk Raj 2024-6-27
Ran in:
%% SEIR model function
function dydt = seir_model(t, y, params, N, fixed_gamma, fixed_delta)
beta = params(1);
q = params(2);
alpha = params(3); % Now alpha is part of the parameters to be estimated
S = y(1);
E = y(2);
I = y(3);
R = y(4);
dS = -beta * S * (q * E + I) / N;
dE = beta * S * (q * E + I) / N - fixed_delta * E;
dI = fixed_delta * E - fixed_gamma * I;
dR = fixed_gamma * I;
dydt = [dS; dE; dI; dR];
end
%% GMMP scheme to solve FDEs
function y = gmmp_scheme(f, y0, t, params, N, fixed_gamma, fixed_delta)
n = length(t);
m = length(y0);
y = zeros(m, n);
y(:, 1) = y0;
% Initial values
for i = 2:n
ti = t(i);
yi = y(:, i-1);
% Fractional derivative approximation using GMMP
yi_next = yi + (ti^(params(3)-1) - (i-1)^(params(3)-1)) * f(ti, yi, params, N, fixed_gamma, fixed_delta);
y(:, i) = yi_next;
end
end
%% Objective function for MGAM
function error = objective_function(params, f, y0, t, data, N, fixed_gamma, fixed_delta)
beta = params(1);
q = params(2);
alpha = params(3); % Now alpha is part of the parameters to be estimated
y = gmmp_scheme(f, y0, t, params, N, fixed_gamma, fixed_delta);
infected_model = y(3, :);
error = sum((infected_model - data).^2);
end
%% MGAM function for parameter estimation
function estimated_params = mgam(f, y0, t, data, N, fixed_gamma, fixed_delta)
param_guess = [0.5, 0.5, 0.5]; % Guess for beta, q, and alpha
obj_func = @(params) objective_function(params, f, y0, t, data, N, fixed_gamma, fixed_delta); % Pass data here
options = optimoptions('ga', 'Display', 'iter', 'PopulationSize', 50);
[estimated_params, ~] = ga(obj_func, 3, [], [], [], [], [0, 0, 0], [1, 1, 1], [], options);
end
%% Main script
N = 18805278;
m = 1; % <-- User input value
S0 = N * m/100;
E0 = 0;
I0 = 15;
R0 = 0;
y0 = [S0; E0; I0; R0];
t0 = 0;
tf = 250;
h = 0.1;
t = t0:h:tf;
data = y0(3) * exp(0.05 * (t - t0));
%% Define fixed gamma and delta here (replace with your desired values)
fixed_gamma = 1/7;
fixed_delta = 1/12;
estimated_params = mgam(@seir_model, y0, t, data, N, fixed_gamma, fixed_delta);
Single objective optimization: 3 Variables Options: CreationFcn: @gacreationuniform CrossoverFcn: @crossoverscattered SelectionFcn: @selectionstochunif MutationFcn: @mutationadaptfeasible Best Mean Stall Generation Func-count f(x) f(x) Generations 1 100 1.628e+15 1.628e+15 0 2 147 1.628e+15 1.628e+15 0 3 194 1.628e+15 1.628e+15 0 4 241 1.628e+15 1.628e+15 1 5 288 1.628e+15 1.628e+15 0 6 335 1.628e+15 1.628e+15 1 7 382 1.628e+15 1.628e+15 2 8 429 1.628e+15 1.628e+15 0 9 476 1.628e+15 1.628e+15 1 10 523 1.628e+15 1.628e+15 2 11 570 1.628e+15 1.628e+15 0 12 617 1.628e+15 1.628e+15 0 13 664 1.628e+15 1.628e+15 1 14 711 1.628e+15 1.628e+15 0 15 758 1.628e+15 1.628e+15 1 16 805 1.628e+15 1.628e+15 2 17 852 1.628e+15 1.628e+15 3 18 899 1.628e+15 1.628e+15 4 19 946 1.628e+15 1.628e+15 5 20 993 1.628e+15 1.628e+15 0 21 1040 1.628e+15 1.628e+15 1 22 1087 1.628e+15 1.628e+15 0 23 1134 1.628e+15 1.628e+15 0 24 1181 1.628e+15 1.628e+15 0 25 1228 1.628e+15 1.628e+15 0 26 1275 1.628e+15 1.628e+15 1 27 1322 1.628e+15 1.628e+15 2 28 1369 1.628e+15 1.628e+15 3 29 1416 1.628e+15 1.628e+15 4 30 1463 1.628e+15 1.628e+15 5 Best Mean Stall Generation Func-count f(x) f(x) Generations 31 1510 1.628e+15 1.628e+15 0 32 1557 1.628e+15 1.628e+15 0 33 1604 1.628e+15 1.628e+15 1 34 1651 1.628e+15 1.628e+15 0 35 1698 1.628e+15 1.628e+15 1 36 1745 1.628e+15 1.628e+15 0 37 1792 1.628e+15 1.628e+15 1 38 1839 1.628e+15 1.628e+15 2 39 1886 1.628e+15 1.628e+15 0 40 1933 1.628e+15 1.628e+15 1 41 1980 1.628e+15 1.628e+15 0 42 2027 1.628e+15 1.628e+15 0 43 2074 1.628e+15 1.628e+15 1 44 2121 1.628e+15 1.628e+15 0 45 2168 1.628e+15 1.628e+15 0 46 2215 1.628e+15 1.628e+15 1 47 2262 1.628e+15 1.628e+15 0 48 2309 1.628e+15 1.628e+15 0 49 2356 1.628e+15 1.628e+15 1 50 2403 1.628e+15 1.628e+15 2 51 2450 1.628e+15 1.628e+15 3 ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
beta = estimated_params(1);
q = estimated_params(2);
alpha = estimated_params(3);
fprintf('Estimated Parameters:\n');
Estimated Parameters:
fprintf('Beta: %.2f\n', beta);
Beta: 1.00
fprintf('Alpha: %.2f\n', alpha);
Alpha: 1.00
fprintf('Q: %.2f\n', q);
Q: 1.00
fprintf('\n');
y = gmmp_scheme(@seir_model, y0, t, [beta, q, alpha], N, fixed_gamma, fixed_delta);
figure;
plot(t, y(3, :), 'r-', 'LineWidth', 2);
hold on;
% Set the plot limits
xlim([0 250]); % X-axis limit from 0 to 250
ylim([0 12000]); % Y-axis limit from 0 to 12000
% Optionally, add labels and title for better visualization
xlabel('Time');
ylabel('Infections');
title('SEIR Model Simulation');
hold off;
this is the code which i used for plotting. but while estimating the parameters, each time I run the code i get different parameter values also the parameter values estimated are wrong. what would possibly be the problem?

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

回答（1 个）

Kautuk Raj
Kautuk Raj 2024-6-27
It looks like you are encountering issues with parameter estimation in your SEIR model using a genetic algorithm (GA). Here are a few suggestions which you can try:
1. Initial Guess and Bounds:
Expand the bounds to realistic ranges. For example:
lower_bounds = [0, 0, 0];
upper_bounds = [5, 5, 2];
2. Objective Function Scaling:
Normalize the error in your objective_function to improve optimization:
error = sum(((infected_model - data) ./ data).^2);
3. Optimization Options:
Adjust GA options to improve convergence:
options = optimoptions('ga', 'PopulationSize', 100, 'MaxGenerations', 200, 'FunctionTolerance', 1e-6);
4. Stochastic Nature of GA:
Set a random seed for reproducibility:
rng(1);
I hope these adjustments help improve the parameter estimation process in your SEIR model.

类别

在 Help Center 和 File Exchange 中查找有关 Genetic Algorithm 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by