Error in Model Fitting

2 Ansichten (letzte 30 Tage)
Abdullahi
Abdullahi am 7 Aug. 2024
Bearbeitet: Torsten am 8 Aug. 2024
Hello! i have been to run this model fitting code but keeps on getting these error messages as:
  1. FUN must have two input arguments.
  2. userFcn_ME = initEvalErrorHandler(userFcn_ME,funfcn_x_xdata{3}, ...
  3. [params_fit, ~] = lsqcurvefit(obj_fun, params0, time, observed_I, [], [], options);
Find attached, the code,, Please help
function dydt = sir_model(t, y, beta, gamma)
S = y(1);
I = y(2);
R = y(3);
N = S + I + R; % Total population
dSdt = -beta * S * I / N;
dIdt = beta * S * I / N - gamma * I;
dRdt = gamma * I;
dydt = [dSdt; dIdt; dRdt];
end
function error = model_error(params, time, observed_I, y0)
beta = params(1);
gamma = params(2);
% Solve the differential equations
[~, y] = ode45(@(t, y) sir_model(t, y, beta, gamma), time, y0);
% Extract the infectious data from the solution
model_I = y(:, 2);
% Compute the error as the difference between model and observed data
error = model_I - observed_I';
end
% Define the time vector and observed data
time = [0, 1, 2, 3, 4, 5]; % Example time points
observed_I = [10, 15, 20, 25, 30, 35]; % Example observed infectious data
% Initial guess for parameters
beta_guess = 0.3;
gamma_guess = 0.1;
params0 = [beta_guess, gamma_guess];
% Initial conditions
S0 = 1000; % Initial number of susceptible individuals
I0 = observed_I(1); % Initial number of infectious individuals
R0 = 0; % Initial number of recovered individuals
y0 = [S0; I0; R0];
% Define the objective function for fitting
obj_fun = @(params) model_error(params, time, observed_I, y0);
% Perform the curve fitting
options = optimoptions('lsqcurvefit', 'Display', 'off');
[params_fit, ~] = lsqcurvefit(obj_fun, params0, time, observed_I, [], [], options);
% Display the fitted parameters
beta_fit = params_fit(1);
gamma_fit = params_fit(2);
disp(['Fitted beta: ', num2str(beta_fit)]);
disp(['Fitted gamma: ', num2str(gamma_fit)]);

Akzeptierte Antwort

Torsten
Torsten am 7 Aug. 2024
Bearbeitet: Torsten am 7 Aug. 2024
% Define the time vector and observed data
time = [0, 1, 2, 3, 4, 5]; % Example time points
observed_I = [10, 15, 20, 25, 30, 35]; % Example observed infectious data
% Initial guess for parameters
beta_guess = 0.3;
gamma_guess = 0.1;
params0 = [beta_guess, gamma_guess];
% Initial conditions
S0 = 1000; % Initial number of susceptible individuals
I0 = observed_I(1); % Initial number of infectious individuals
R0 = 0; % Initial number of recovered individuals
y0 = [S0; I0; R0];
% Define the objective function for fitting
obj_fun = @(params,time) model_error(params, time, y0);
% Perform the curve fitting
options = optimoptions('lsqcurvefit', 'Display', 'off');
[params_fit, ~] = lsqcurvefit(obj_fun, params0, time, observed_I, [], [], options);
% Display the fitted parameters
beta_fit = params_fit(1);
gamma_fit = params_fit(2);
disp(['Fitted beta: ', num2str(beta_fit)]);
Fitted beta: 1.6778
disp(['Fitted gamma: ', num2str(gamma_fit)]);
Fitted gamma: 1.2954
hold on
plot(time,observed_I,'o')
plot(time,model_error(params_fit,time,y0))
hold off
grid on
function error = model_error(params, time, y0)
beta = params(1);
gamma = params(2);
% Solve the differential equations
[~, y] = ode45(@(t, y) sir_model(t, y, beta, gamma), time, y0);
% Extract the infectious data from the solution
model_I = y(:, 2);
% Compute the error as the difference between model and observed data
error = model_I.';
end
function dydt = sir_model(t, y, beta, gamma)
S = y(1);
I = y(2);
R = y(3);
N = S + I + R; % Total population
dSdt = -beta * S * I / N;
dIdt = beta * S * I / N - gamma * I;
dRdt = gamma * I;
dydt = [dSdt; dIdt; dRdt];
end
  6 Kommentare
Abdullahi
Abdullahi am 8 Aug. 2024
I don't actually know the difference between the two..
what do you recommend?
Torsten
Torsten am 8 Aug. 2024
Bearbeitet: Torsten am 8 Aug. 2024
Both solvers are equivalent - only the user interfaces are different.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Gamma Functions finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by