0で除算しました。

しばらくコロナ禍専用、時々生活と美しい物。

2019-nCov(新型コロナウイルス) のExcelシミュレーション

[tex: ]

ある動画より

昨日自宅でこの動画を見ていた。

www.youtube.com

前半で、中国・武漢市内の悲惨な動画がいくつか紹介される。手術をする医師や会議室での女性がバタバタ倒れていく。こちらもなかなかショッキングだか、後半で動画主がExcelで計算している。武漢の状況を日本にあてはめるとどうなるか?を計算したら、そんなバカな!という状況になっている。

試算の準備

例によってジョンズ・ホプキンズ大学の2019-nCovのページを参照する。トータルで感染者数1万人を上回るのも時間の問題である。

それでは、上期の動画主はどのようにしてExcelで感染者数を試算していたのか?わからないので、予測モデルとして、ここでは以下の微分方程式を使う。

ロジスティック方程式 - Wikipedia

N:感染者の数(時刻 t の関数)、r:個体あたり増加率、K:環境収容力

          f:id:divisionby0:20200131125651p:plain

 口ジスティック方程式の解析解は以下の通り。これをロジスティック曲線と呼ぶ。  

           f:id:divisionby0:20200131125952p:plain

※まだはてなブログTeX記法を習得していないため、とりあえずこれでお茶を濁す

Nをロジスティック方程式に当てはめるには、いくつかの仮定が必要である。以下上記WiKipediaから引用する。

1.環境内には単一の種か、あるいは同等とみなせる種のみが存在する。

2.対象の生物の各世代(親子)は連続的に重なっている。すなわち、連続的に子が生まれ、親と子が共存する期間が存在する。

3.個体は一定の大きさの環境内に常に存在する。すなわち、環境から移出したり、外部から移入が無い。(用語としては閉じた個体群とも呼ばれる)

4, 環境の大きさは変わらず、一定状態が保たれる。

5. 個体群のために、食糧や資源が一定して供給される。

  武漢市の場合は、5.は知らないけれど他は条件を満たしているとする。

Excelでl計算してみる

実際にExcelで計算してみるには、tを変数として動かしてみれば良いが、その前に2つのパラメータKとrを決める必要がある。

まずKであるが、t→∞でN→Kであることから

K=(コロナウイルスが流行している)都市の人口

として良いだろう。あとはrであるが、ロジスティック曲線に色々代入して計算するよりも、関数を定義して強引にExcelで出したほうが早い。

  1. Excelを新規に立ち上げる。
  2. [Alt] + [F11] でVBA開発画面へ
  3. 以下をコード画面にコピペする。

Function ncov(t As Double, n0 As Double, n1 As Double) As Double
    r = 0.17
    ncov = n1 / (1 + (n1 / n0 - 1) * Exp(-r * t))
End Function 

  上で定義した関数ncov(t,n0,n1)の使い方:

  • おなじt : 最初の患者が発生してから何日目か(t=1,2,3...
  • n0: t=1 のときの感染者数(0以上の数)
  • n1: 都市の人口

武漢の場合はn1=1150万を入れる。tを増やしながらセルにプロットし、グラフを書いてみよう。なお、rの値は、12/8に最初の患者が発生し、本日1/31時点で湖北省の患者が(公式発表の)9658人に近くなるように調整した。汎用的にr=0.17を推奨する。

f:id:divisionby0:20200201093425p:plain

*もしかしてこの領域ではほぼ指数関数でないかと思い、単純にN=exp(0.17t)をプロットしてみたら、12/8〜1/31ではほぼ同じでした。短期的にはdN/dt=rN のマルサスモデルで十分だったんですね…

このモデルで国内都市を推定

全く同じモデルを、1/30 に初めて患者が発生した京都市に当てはめてシミュレーションしてみる。人口は150万とした。

f:id:divisionby0:20200201133041p:plain

あまり面白くない結果が出た。人口は100万だろうが1,000万だろうがあまり関係なく、ほぼ同じN=exp(0.17t)の曲線と同じ。そして3月末には感染者3万人を突破する。