Regresi kurva eksponensial dan logistik untuk jumlah kasus di Indonesia beserta grafiknya

regresi-jumlahkasus-Indonesia.m

% Penghitung model kurva eksponensial dan logistik untuk jumlah kasus di
% Indonesia beserta grafiknya (khusus GNU Octave)
% @author Muhammad Rifqi Priyo Susanto
% @license MIT License

% Jalankan perintah berikut untuk memasang paket optim.
% pkg install -forge optim

% Jalankan perintah berikut untuk memuat paket optim.
% pkg load optim

clear; clc;

% jumlah kasus COVID-19 di Indonesia (hingga 24 Maret 2020)
% sumber data: covid19.go.id
tanggal = '24 Maret 2020';
kasus = [  2,   2,   2,   2,   4, ...
           4,   6,  19,  27,  34, ...
          34,  69,  96, 117, 134, ...
         172, 227, 309, 369, 450, ...
         514, 579, 686];

% variabel bantuan
jumlah = size(kasus)(2);
rerata = sum(kasus) / jumlah;
total_kuadrat = sum((kasus - rerata).^2);
kasus_terakhir = kasus(jumlah);
hari = 1:jumlah;

% model kurva eksponensial
f = @ (p, x) p(1) * (p(2).^x);
fp_awal = [.8; 1.4];

[fp, fy, f_cvg, f_klr] = nonlin_curvefit(f, fp_awal, hari, kasus);
f_selisih = [];
for i = hari f_selisih(i) = f(fp, i) - rerata; end
f_r2 = sum(f_selisih.^2) / total_kuadrat;
fprintf('eksponensial: y = %.2f * (%.2f^x)\n', fp(1), fp(2));
fprintf('    R^2 = %.1f%%\n', f_r2 * 100);

% model kurva logistik
g = @ (p, x) p(1) ./ (1 + p(2) * exp(p(3) * x));
gp_awal = [kasus_terakhir; kasus_terakhir / 2; .3];

[gp, gy, g_cvg, g_klr] = nonlin_curvefit(g, gp_awal, hari, kasus);
g_selisih = [];
for i = hari g_selisih(i) = g(gp, i) - rerata; end
g_r2 = sum(g_selisih.^2) / total_kuadrat;
fprintf('logistik: y = %.2f / (1 +  %.2f exp(%.2fx))\n', gp(1), gp(2), gp(3));
fprintf('    R^2 = %.1f%%\n', g_r2 * 100);

% perkiraan tiga hari ke depan
th_eksp = [round(f(fp, jumlah + 1)),
           round(f(fp, jumlah + 1)) / kasus_terakhir - 1];
th_logi = [round(g(gp, jumlah + 1)),
           round(g(gp, jumlah + 1)) / kasus_terakhir - 1];
for i = 2:3
  th_eksp(2 * i - 1) = round(f(fp, jumlah + i));
  th_eksp(2 * i) = round(f(fp, jumlah + i)) / th_eksp(2 * i - 3) - 1;
  th_logi(2 * i - 1) = round(g(gp, jumlah + i));
  th_logi(2 * i) = round(g(gp, jumlah + i)) / th_logi(2 * i - 3) - 1;
end
fprintf('perkiraan tiga hari ke depan\n');
fprintf('eksponensial\t%d\t+%.1f%%\t%d\t+%.1f%%\t%d\t+%.1f%%\n',
    th_eksp(1), th_eksp(2) * 100,
    th_eksp(3), th_eksp(4) * 100,
    th_eksp(5), th_eksp(6) * 100);
fprintf('logistik    \t%d\t+%.1f%%\t%d\t+%.1f%%\t%d\t+%.1f%%\n',
    th_logi(1), th_logi(2) * 100,
    th_logi(3), th_logi(4) * 100,
    th_logi(5), th_logi(6) * 100);

% perkiraan jumlah kasus maksimal
jumlah_maks = floor(gp(1));
fprintf('perkiraan jumlah kasus maksimal: %d kasus\n', jumlah_maks);

% perkiraan waktu kasus terakhir
waktu_akhir = ceil(log((gp(1) - floor(gp(1))) / (gp(2) * floor(gp(1)))) / gp(3));
fprintf('perkiraan waktu kasus terakhir: hari ke-%d\n', waktu_akhir);

% serba-serbi grafik
hari_trans = (jumlah):(jumlah + 1);
hari_eksp = (jumlah + 1):(jumlah * 2);
hari_logi = (jumlah + 1):(jumlah * 3);

% grafik 1: keseluruhan
grafik_tinggi1 = gp(1) * 1.5;
grafik_batasx1 = [0 (jumlah * 3)];
grafik_batasy1 = [(grafik_tinggi1 * -.05) (grafik_tinggi1 * 1.05)];
figure;
plot1 = plot(
  hari, kasus, '-og',
  hari_eksp, round(f(fp, hari_eksp)), '-or',
  hari_logi, round(g(gp, hari_logi)), '-ob',
  grafik_batasx1, [jumlah_maks jumlah_maks], '--b',
  hari_trans, [kasus(jumlah) round(f(fp, jumlah + 1))], '-r',
  hari_trans, [kasus(jumlah) round(g(gp, jumlah + 1))], '-b',
  [waktu_akhir waktu_akhir], grafik_batasy1, '--k'
);
set(plot1(1:3), 'linewidth', 2);
set(plot1(5:6), 'linewidth', 2);
xlim(grafik_batasx1);
ylim(grafik_batasy1);
jdl1 = title(sprintf(
    'Perkembangan Jumlah Kasus COVID-19 di Indonesia hingga %s', tanggal));
set(jdl1, 'fontsize', 14);
sbx1 = xlabel('hari ke-');
sby1 = ylabel('jumlah kasus');
set(sbx1, 'fontsize', 14);
set(sby1, 'fontsize', 14);
legenda1 = legend(
  'data saat ini',
  '(perkiraan) tanpa penanganan',
  '(perkiraan) dengan penanganan',
  '(perkiraan) jumlah kasus maks.'
);
set(legenda1, 'fontsize', 14);
set(legenda1, 'xlim', [0 .5]);
set(legenda1, 'position', [0.67 0.7 0.22 0.18]);
grid on;

% grafik 2: ketepatan model
grafik_tinggi2 = kasus_terakhir * 1.1;
grafik_batasx2 = [-1 (jumlah + 2)];
grafik_batasy2 = [(grafik_tinggi2 * -.05) (grafik_tinggi2 * 1.05)];
hari_khusus = 0:(jumlah + 1);
figure;
plot2 = plot(
  hari_khusus, round(f(fp, hari_khusus)), '--r',
  hari_khusus, round(g(gp, hari_khusus)), '--b',
  hari, kasus, 'og'
);
set(plot2, 'linewidth', 2);
set(plot2(3), 'linewidth', 3);
xlim(grafik_batasx2);
ylim(grafik_batasy2);
jdl2 = title(sprintf('Perbandingan Kurva hingga %s', tanggal));
set(jdl2, 'fontsize', 14);
sbx2 = xlabel('hari ke-');
sby2 = ylabel('jumlah kasus');
set(sbx2, 'fontsize', 14);
set(sby2, 'fontsize', 14);
legenda2 = legend(
  'kurva eksponensial (tanpa penanganan)',
  'kurva logistik (dengan penanganan)'
);
set(legenda2, 'fontsize', 14);
set(legenda2, 'xlim', [0 .5]);
set(legenda2, 'position', [0.15 0.8 0.25 0.09]);
grid on;