これは、最近のコロナ禍に関してこの先の知識を深めていきたいと感じた、いちプログラマが
一つずつ調べながらコロナのシミュレーションを書いてみている記録の様なものです。
素人が一から調べながら書き続けているものなので、間違いを含んでいたり今後も修正を加えられていくところがあるかもしれません。
書かれるコードも本業でPythonを書いているわけではないので、甘い部分が多々あると思いますのでご了承願います。
疫学について簡単に調べよう
まずは、コロナ専門の知識というのはネットに多々ありますが、実際に感染症を学ぶ時の最初の知識ってなんだろうかを調べて、基礎っぽい本をAmazonで調べてみました
これを読んでみた感想としては、自分が知りたかったコロナの影響でこれから患者はどう増えていくのか、社会はどうなって、どういう対策を行えばいいのかの基礎的な知識は得られませんでした。
分かったのは以下の事です。
- 疫学のノウハウは今は感染症ではなく、生活習慣病などの研究に活用の主軸を移している
- 本で書かれていることの主体は、調査のためのサンプルの取り方とちょっとだけ結果のまとめ方
なのでもっとシミュレーション寄りの本を買ってみようと思ったのですが、本自体も3000円程度のものでななく、もっと高価なものでしたので、試し読みもせずに買う事を躊躇してしまいました。
一応Amazonのリンクだけ貼っておきます
それに、この時制で外出も難しく読めそうな図書館なども近場にないのでネットを調べていると以下のpdfファイルを見つけました。
非常に簡潔にまとまっている資料なので、そのままなぞる形での紹介になってしまうのですが
分かることをまとめると以下の認識で正しい様です
シミュレーションのモデルは代表的なものは以下の３つ。
- SIRモデル
- SEIRモデル
- MASモデル
そして書かれている、SIRモデル、SEIRモデルの両方とも、計算のモデルはそんなに難しくない微分方程式ということでした。
プログラムを書いてみよう
上記のpdfファイルに書かれているものをあてにするなら、人を４つのグループに分けて状態を遷移させている様です。
- 無免疫者(Suspectiable)
- 感染しているが発症していない者(Expoed)
- 発症者(Infected)
- 回復者(Recoverd)
コロナの場合、感染してから即、人に症状が出で人に感染させるわけではないので。基本的にはSEIRモデルとSIRモデルでは、SIRモデルをベースに調整をかける方が正しい様ですので、まずはSIRモデルのプログラムを書いてみました。
SIRモデルだと、人間は、未感染→発症→回復、と状態を遷移させていきますがそれそれの状態の発生するSIRモデルの数式は、以下の様なシンプルな微分方程式でした。
\frac{dS(t)}{dt} = -bS(t)I(t) \\
\frac{dI(t)}{dt} = bS(t)I(t) - gI(t) \\
\frac{R(t)}{dt} = gI(t) \\
b:感染率、g:回復率
これくらいの数式でしたら、Pythonで簡単にプログラムに落とせそうなので、自分なりに書いて見ました
今回は勉強のためにモデルを理解してみようということで
ネットから得た幾つかの数値から、ざっくりパラメーターを以下の様に設定してpythonでプログラムを書いて見ました。
- 国内の感染者の初期値：100人
- 2週間で2.5人に感染させる
- 2週間で抗体を得て他の人に感染させなくなる
書いたプログラムは、numpyとmatplitlibさえ入っていれば実行可能な、簡素な数値計算のコードです。
これで実行をしてみた結果、毎日発生する新規の感染者の推移を一月分グラフにプロットしてみました。
新規の感染者が指数関数的に増えているのが良くわかります。
これをベースに色々とパラメータを弄ってみることで、長期的な患者の推移や、感染率をどれくらい減少させればいいのかを実際に計算して進めることができそうです。
試しにまずはこのパラメーターで１年分の新規感染者の推移を見てた結果です。
最高で１日２５０万人を超える感染者が現れて、半年ほどの地獄を見てから集団免疫を獲得することになるようです。
あとは感染率を倍に引き上げた結果を見てみます。
３０日後の感染者数が、あっという間に341人から５万２千まで増加しました。
ちょっとした感染率の変化がイタリアやニューヨークなどの地獄を引き起こしていたという事が良く見てとれました。
シミュレーションをもっと細かくしてみよう
上記でなんとなくモデルは理解が出来た気がするので、あとはもっと現実に即した形でシミュレーションを行いたいので、検索をし直したら、以下の様なページを見つけました。
この記事を読む限り、先ほどの様にSIRモデルを選んだのは正しかった様です。
ただ、感染してから発症するまでに潜伏期間があるので、発症者の数を計算するためには発症期間をずらすためのプログラムの変更が必要になりました。
これに合わせて、これまで見つかっている感染者数の推移から、パラメータをを以下の様に修正しました。
感染率: 0.488
回復率: 0.0455
感染から14日後に発症する。
こちらの、プログラムは以下のURLにあります。
出力されたグラフが以下の通り、これは日本の2/28〜3/12までの推移をベースに３月末日の新規感染者数を予想したものです。
3500人を超える患者に増加する計算結果になりましたが、オリジナルのプログラムの結果のグラフが下のグラフになっていました。
3500人と4000人で３月末日時点での新規患者数にズレがありますが、オリジナルのソースコードが公開されていない中では、それなりに再現が出来たと思われます。
というより、計算結果が近くなるまで何度も記事を読みなおして何が抜けているか１日かけて修正を繰り返したので、このあたりで心が折れてきました。
なんとなく理解できたので、このまま次に進ませて頂きます。
あとは、３月前半までの患者数の推移で行った場合、最終的にどんな結末になったのかシミュレーションの続きを見ていきます。
まずは、新規の感染者数の推移が以下のグラフでした。
3/1時点では242人だった新規の感染者は、132日後に一日の発生数が385万3377人になります。
その後減少して182日目に一日の発生数が1000人を切って、このあたりで集団免疫を獲得したかなといえる数値になってきました。
イギリスが最初にやろうとしていたり、現在スウェーデンの行っている作戦がこれにあたるのでしょう。
新規の感染者が発生するという事はその人たちが回復するまでは、患者として積みあがっていきますから、病院なり自宅なりで治療を受ける人の数がどれくらいになるかというのが本当の問題なので
それを計算しなおしてみました。
結果としては、157日目に患者数はピークを迎えて、5576万6923人が感染中になりました。
ピーク時には人口の４割以上が感染している状態になりますね。
ニュースでは感染している人の８割は無症状者というので、ピーク時には人口の１割以上がなんらかの症状が出ている事になるようです。
このレベルまで感染や発症が進むと、社会の必要なインフラに関わる人に対しても一定数の欠員が発生するので、ある程度の期間は、医療以外にも電力や警察などのサービスにも問題が発生すると思われます。
仮に感染者数１０００万人を、その基準とすると、その期間はシミュレーションを開始した2/28をスタートとして123〜199日の間。およそ２か月半程度インフラが不安定になると計算から推測出来ます。
ここまでで得られた数値から何を読み取るのかは人によって様々なのですが、個人的な感想としては、スウェーデンの様に集団免疫を獲得する作戦も、集団免疫を獲得するまでの期間もそれなりに長いですし、インフラを含めた社会全体を不安定にさせる期間も短くないので、ワクチンの開発と量産の期間が１年２年で収まるのであれば、ワクチンの完成まで我慢をした方が良さそうに見えてきました。
減少させるにはどうすれば良いのか？
これで、感染率と回復率が分かれば今後の推移を計算出来る様になりましたが
感染率と回復率自体は、色々な値を入れてシミュレーションを繰り返して、現実の数値に合う値に寄せていくしか無いそうです。
ここはヒューリスティックに答えを求めていくアルゴリズムは幾つかありますが、実際にはどれを使っているかは分かりませんでした。
今回は３月半ばまでの数値を元にシミュレーションを行いましたが、その後のシミュレーションも感染率を、都知事の注意喚起や総理の緊急事態宣言などの政治イベントのタイミングで変化させる様な式に変えれば実現出来そうです。
あとは、新規の患者数が何所かで横ばいか減少に転じないと最終的な感染者の総数は変わらないので、どこかで減少に転ずるにはどうすればいいかを考えれば良いのですが、実際に患者数を減少させるためには、以下の式を満たせば良い様です。
感染率 < 回復率
効果的な治療薬を、患者の数だけ量産でもしない限り回復率を下げるのは難しいので、現実的には感染率を下げる作戦を立てる事になるのですが、感染率と行動抑制を繋げる式を、見つけることができませんでした。
このあたりがもう少し分かれば記事を更新していく予定です。