SIMULASI COVID DIAMOND PRINCESS

clc
clear all
S(1) = 3699;
I(1) = 1;
R(1) = 0;
N = 3700;

ttes = [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
];
Ites = [0	0	0	0	0	0	61	61	64	135	135	175	175	218	285	355	454	542	621	634	634	634	691	691	691	705	705	705	705	705	705	706	706	706	696
];

t(1) = 0;
dt = 1;
T = 50;

for j=1:10000;
    betha(j) = 0.6 + 0.2*rand(1);
    gamma(j) = 0.2 + 0.1*rand(1);
    %------------SIR-------------
    for i=1:T;
        S(i+1) = S(i) - (betha(j)*S(i)*I(i)/N)*dt;
        I(i+1) = I(i) + (betha(j)*S(i)*I(i)/N - gamma(j)*I(i))*dt;
        R(i+1) = R(i) + gamma(j)*I(i)*dt;
        t(i+1) = t(i) + dt;
    end

    Y = I(1:length(Ites));
    e = (Ites-Y).^2;
    LS(j) = sum(e);
end
min(LS) % 1.0846e+06
x = find(LS == min(LS));
b = betha(x)
g = gamma(x)

%% ==== PLOT ====
b = 0.6287;
g = 0.2532;
for i=1:T;
    S(i+1) = S(i) - (b*S(i)*I(i)/N)*dt;
    I(i+1) = I(i) + (b*S(i)*I(i)/N - g*I(i))*dt;
    R(i+1) = R(i) + g*I(i)*dt;
    t(i+1) = t(i) + dt;
end
plot(t,S,t,I,t,R)
hold on
plot(ttes, Ites,'ok');
title('SIR Model')
legend('S','I','R','I real')
xlabel('Time (days)')
ylabel('Number of people')