Fitting model to data set by estimating parameters

I am trying to estimate the parameters in a model which is based on a system of ordinary differential equations. I have tried the following, but there is something wrong with my method. Here is some of the code for context

load('data.mat');
t = time;
T_exp = data(1,:);
I_exp = data(2,:);

p = 2.23 * 10 ^ - 2;
c = 0.1;
beta
T0 = 10000;
I0 = 0;
V0 = 0.05;
initial_values = [T0; I0; V0];

beta_initial_guess = 10 ^- 5;
delta_initial_guess = 0.04;

estimated_params = fminsearch(@(params) objective_function(params, t, initial_values, T_exp, I_exp), [beta_initial_guess, delta_initial_guess])

function error = objective_function(params, t, initial_values, T_exp, I_exp)
    beta = params(1);
    delta = params(2);
    p = 2.23 * 10 ^ - 2;
    c = 0.1;
    
    [tSol, YSol] = ode45(@(t, Y) virus_solver(t, Y, beta, delta, p, c), t, initial_values);
    
    % Model predictions
    T_model = YSol(:, 1);
    I_model = YSol(:, 2);
    
    % Calculate the objective function (sum of squared differences)
    error = sum((T_model - T_exp).^2) + sum((I_model - I_exp).^2);
end

function dYdt = virus_solver(t, Y, beta, delta, p, c) 
    T = Y(1);
    I = Y(2);
    V = Y(3);
    dTdt = -beta * V * T;
    dIdt = beta * T * V - delta * I;
    dVdt = p * I - c * V;
    dYdt = [dTdt; dIdt; dVdt];
end

Unrecognized function or variable 'beta_true'.

The error i am getting with this code is :

Unable to perform assignment because the size of the left
side is 1-by-1 and the size of the right side is 1-by-28.

Error in fminsearch (line 209)
fv(:,1) = funfcn(x,varargin{:});

Error in exercise_3_4_4 (line 24)
estimated_params = fminsearch(@(params) objective_function(params, t, initial_values, T_exp, I_exp), [beta_true, delta_true]);
 

Is this a good approach to estimate the parameters? Or am thinking about this all wrong.

Answer 1

[Star Strider]

The calculated ‘YSol’ should return a (Nx3) matrix where ‘N’ equals ‘numel(t)’ so if it does not, you need to determine the reason. Common problems include ode45 (in this instance) stopping because it has encountered a singularity or other significant discontinuity, and returning a matrix that does not have the same row size as the data. I would need ‘data.mat’ to see what the problem is. (I have written a number of these sorts of parameter estimation solutions.)

Determining what the problem is requires running the code, and that requires ‘data.mat’.

[Oskar]

I have tried to figure out what it is based on your reponse, but I have not figured it out yet. I have attatched the data in this comment. Is there anything more you can tell me?

load('data_es4.mat');
t = time

t = 1×28

0 2 5 9 13 20 29 38 47 55 64 72 80 88 95 102 108 113 118 123 127 133 139 146 157 171 186 200

T_exp = data(1,:)

T_exp = 1×28

1.0e+04 * 1.0000 1.0488 0.9502 0.9351 1.0088 0.9704 1.0907 0.9439 0.9990 0.9614 0.9620 0.9321 0.9002 0.6588 0.2836 0.1216 0.1809 0.0632 -0.0331 -0.0066 0.0011 0.0566 0.0338 0.0024 -0.0159 0.0401 0.1050 0.0523

I_exp = data(2,:)

I_exp = 1×28

1.0e+03 * 0 0.1687 -0.2616 0.4588 0.3777 0.9225 -0.0607 -0.3365 -0.0167 0.1950 -0.5081 0.3725 0.7172 2.6404 4.3372 4.8280 5.4770 3.9718 3.6379 2.2146 1.9121 1.6446 1.1767 1.1226 1.7251 0.9629 0.3814 -0.5832

[Star Strider]

The data are all row vectors. They need to be column vectors, then it works and gives a very good fit to the data —

% load('data.mat');

LD = load('data_es4.mat');

t = LD.time(:);                                 % Force Column Vector
T_exp = LD.data(1,:).';                         % Transpose To Column Vector
I_exp = LD.data(2,:).';                         % Transpose To Column Vector

p = 2.23 * 10 ^ - 2;
c = 0.1;
% beta
T0 = 10000;
I0 = 0;
V0 = 0.05;
initial_values = [T0; I0; V0];

beta_initial_guess = 10 ^- 5;
delta_initial_guess = 0.04;

estimated_params = fminsearch(@(params) objective_function(params, t, initial_values, T_exp, I_exp), [beta_initial_guess, delta_initial_guess])

estimated_params = 1×2

0.0001 0.0412

tv = linspace(min(t), max(t)).';
beta = estimated_params(1);
delta = estimated_params(2);

[tSol, YSol] = ode45(@(t, Y) virus_solver(t, Y, beta, delta, p, c), tv, initial_values);

figure
hp11 = plot(tSol, YSol(:,1), 'DisplayName','T(t) Fit');
hold on
hp12 = plot(tSol, YSol(:,2), 'DisplayName','I(t) Fit');
hp13 = plot(tSol, YSol(:,3), 'DisplayName','V(t) Estimate');
hp2 = plot(t, T_exp, 'p', 'DisplayName','T(t) Data', 'Color',hp11.Color);
hp3 = plot(t, I_exp, 'p', 'DisplayName','I(t) Data', 'Color',hp12.Color);
hold off
grid
legend('Location','best')
xlabel('Time')
ylabel('Value')

function error = objective_function(params, t, initial_values, T_exp, I_exp)
    beta = params(1);
    delta = params(2);
    p = 2.23 * 10 ^ - 2;
    c = 0.1;
    
    [tSol, YSol] = ode45(@(t, Y) virus_solver(t, Y, beta, delta, p, c), t, initial_values);
    
    % Model predictions
    T_model = YSol(:, 1);
    I_model = YSol(:, 2);
    
    % Calculate the objective function (sum of squared differences)
    error = sum((T_model - T_exp).^2) + sum((I_model - I_exp).^2);
end

function dYdt = virus_solver(t, Y, beta, delta, p, c) 
    T = Y(1);
    I = Y(2);
    V = Y(3);
    dTdt = -beta * V * T;
    dIdt = beta * T * V - delta * I;
    dVdt = p * I - c * V;
    dYdt = [dTdt; dIdt; dVdt];
end

.

[Oskar]

Thank you very much for that explaination. It solved the problem.

[Star Strider]

As always, my pleasure!