とっても暇なのです。
コードを書かないと頭が鈍ってしまう。
なので、今日は頭の体操に、最近よくニュースで発表されている「実効再生産数」を自分で計算してみたいと思います。
まず、計算アルゴリズム。
有名な西浦先生のやつみてみます。
http://statmodeling.hatenablog.com/entry/covid19-estimate-effective-reproduction-number
難しすぎる・・・。
予備知識がないと何を言っているのかわからない。
正確に計算しようとするとベイズ統計とか回帰解析とか最尤推定とかかなりめんどくさい。
おじさん、この「犬」みたいな漢字が何十年間も読めない・・・。
統計と推定なので結構ディープラーニングの知識と近い気がする。
しかし、いつものごとく、このブログはC言語限定。
プログラムコンテストみたいに数時間でできる簡単なコードでないとだめなのです。
ということで、別のアルゴリズムを探します。
過去の教授のひとりごと
http://www.med.u-toyama.ac.jp/biostat/hitori.html
実効再生産数の論文
https://www.medrxiv.org/content/10.1101/2020.01.27.20018952v1.full.pdf
こっちは簡単そう。
とりあえず上記のアルゴリズムでプログラムを組んでみます。
次にデータ。
おじさんの住んでいる東京都のデータを取ってきます。
https://stopcovid19.metro.tokyo.lg.jp/
最近は東京都のホームページからオープンデータとして感染者のデータが取れるのね。
すごいねぇ。
なんかこの前やったKaggleみたい。
これらのデータとアルゴリズムで実効再生産数を求めて、エクセルでグラフにてみました。
できたのかな？？？？
発症日じゃなくて発表日であるとか、配列の日にちをすこしずらさないといけないとか細かいところはたくさん違うけど、まぁしょうがない。
コードは以下の通りです。
-------------------------------------
#include <stdio.h>
#include <string.h>
#include <math.h>
#define MAX_LIST 200
struct infected {
char date[16];
int n;
double rt;
};
static int mon_day[] = { 31,29,31,30,31,30,31,31,30,31,30,31,0 };
struct infected infe[MAX_LIST];
static void create_infe_base()
{
int i, j,c=0;
for (i = 0; i < 100; i++) {
if (mon_day[i] == 0)break;
for (j = 0; j < mon_day[i]; j++) {
if (c >= MAX_LIST)break;
sprintf(infe[c].date, "2020-%02d-%02d", i+1, j+1);
c++;
}
}
}
char *mystrtok(char *s1, const char *s2) {
static char *str = 0;
register int i, j;
if (s1)str = s1;
else s1 = str;
if (!s1) { return(0); }
j = strlen(s2);
while (1) {
if (!*str) {
str = 0;
return(s1);
}
for (i = 0; i < j; i++) {
if (*str == s2[i]) {
*str++ = 0;
return(s1);
}
}
str++;
}
}
static int input_infe_data()
{
int ret = -1,id,ct,i;
FILE* fp;
char buf[256],date[256];
char* ptr;
fp = fopen("130001_tokyo_covid19_patients.csv","rb");
if (fp == NULL)goto end;
while (1) {
id = 0;ct = 0;buf[0] = 0;date[0] = 0;
fgets(buf, sizeof(buf), fp);
if (buf[0] == 0)break;
ptr = mystrtok(buf, ",");
if (ptr == NULL)break;
sscanf(ptr,"%d", &id);
if (id < 1)continue;
while (ptr != NULL) {
ct++;
ptr = mystrtok(NULL, ",");
if (ptr != NULL && ct==4) {
strcpy(date, ptr);
break;
}
}
if (date[0] == 0)continue;
for (i = 0; i < MAX_LIST; i++) {
if (strcmp(infe[i].date, date))continue;
infe[i].n++;
break;
}
}
ret = 0;
end:
if (fp)fclose(fp);
return ret;
}
static void print_infe()
{
int i;
for (i = 10; i < MAX_LIST; i++) {
//printf("%4d %s %5d %.2f\n",i, infe[i].date, infe[i].n, infe[i].rt);
printf("%d,%s,%d,%.2f\n", i, infe[i].date, infe[i].n, infe[i].rt);
}
}
static void calc_r()
{
int i,j;
double R, K, L = 7, D = 9, A, B;
for (i = 9; i < MAX_LIST;i++) {
A = 0,B=0;
for (j = 0; j < 5; j++) {
A += infe[i - j].n;
B += infe[i - 5-j].n;
}
if (A == 0 || B == 0)continue;
K = (log10(A) - log10(B))/5;
R = K * K*(L*D) + K * (L + D) + 1;
if (R < 0)continue;
infe[i].rt = R;
}
}
int main()
{
create_infe_base();
input_infe_data();
calc_r();
print_infe();
return 0;
}
-------------------------------------
今日は実効再生産数を求めるでした。
