-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsiroutput.asv
49 lines (38 loc) · 1.52 KB
/
siroutput.asv
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
%% This function takes three inputs
% x - a set of parameters
% t - the number of time-steps you wish to simulate
% data - actual data that you are attempting to fit
function f = siroutput(x,t,data)
% set up transmission constants
k_infections = x(1);
k_fatality = x(2);
k_recover = x(3);
% set up initial conditions
ic_susc = x(4);
ic_inf = x(5);
ic_rec = x(6);
ic_fatality = x(7);
% Set up SIRD within-population transmission matrix
A = [(1-k_infections) 0 0 0;
k_infections (1-k_recover-k_fatality) 0 0;
0 k_recover 1 0;
0 k_fatality 0 1];
B = zeros(4,1);
% Set up the vector of initial conditions
x0 = [ic_susc, ic_inf, ic_rec, ic_fatality];
% simulate the SIRD model for t time-steps
sys_sir_base = ss(A,B,eye(4),zeros(4,1),1)
y = lsim(sys_sir_base,zeros(t,1),linspace(0,t-1,t),x0);
% return a "cost". This is the quantitity that you want your model to
% minimize. Basically, this should encapsulate the difference between your
% modeled data and the true data. Norms and distances will be useful here.
% Hint: This is a central part of this case study! choices here will have
% a big impact!
%stl = data(data{:, 5} == 2, :);
%population = ic_susc-stl{:,4};
population = 1-data(:, 2);
casesModel = y(:, 2) + y(:, 3) + y(:, 4);
%1: Suscpetible, 2: infected, 3:
% f =sum( (population-y(:, 1)).^2+ (casesModel-stl{:, 3}).^2 +(stl{:, 4}-y(:, 2)).^2);
f =sum( (population-y(:, 1)-y(:, 1)).^2+ (casesModel-data(:, 1)).^2 +(data(:, 2)-y(:, 2)).^2);
end