微分方程式モデルを勉強する話
前回↓
今日もこの本
今回から「数値シミュレーション」編らしい。
で、とりあえずそれには"微分方程式"が重要らしい。いまいちイメージ湧かないけど、微分方程式で時間経過ごとのグラフの動きがわかれば数値がどう動くか予想できる？ってことかな。
- パラメータ設定
- 初期化
- 時間発展方程式
- グラフ描画
という構成でなっている
SIRモデル
数値シミュレーションの1つとして、感染症の予測をするための感染症モデルを考える。
今回は代表例としてSIRモデルを考える。
SIRモデルは
- Susceptible(感受性保持者)
- Infected(感染者)
- Recovered(免疫保持者)
の頭文字で命名されたもので、それぞれの経時的変化量は
ΔR = I * γ (感染者の一定割合が免疫保持者になる)
ΔI = -ΔS + ΔR (感染者は感受性保持者の減少分増え、免疫保持者の増加分減る)
ΔS = S * I * β (感受性保持者の一定割合が感染者になるが、それは感染者数に依存する)
で表せる。 ※γ、βは定数
これを任意の時間t-1→tの変化量として、tに対してforループを回せばグラフが書ける。
これを時間発展方程式というらしい。
実際には感染者の限界や、他のウイルスによる競合があるので、上記以外の微分方程式も考える必要がある。
ねずみ算
最も簡便な微分方程式はねずみ算で、制限がない状態での拡散の表現に使う。
時間tでのねずみの数をn[t]とすると、Δtでねずみが１匹あたり産む数を定数aとして
n[t] = Δn+n[t-1]
Δn = a*n[t-1]
として表せる。
しかし実際にはこの式通りにネズミは増える事はなく、頭打ちになる。
ロジスティック方程式
ねずみ算の課題であった、頭打ちを考慮した方程式
上記のネズミの例で考えると、頭打ちの数値＝capacityを設定して
Δn = a*n[t-1]*(1-n[t-1]/capacity)
で表せる。
n[t-1]の値がcapacityに近づくほどΔnは小さくなっていくので、最終的にはcapacityの数値で頭打ちになるようになっている。
逆にいうと、aとcapacityを設定すればいつ頭打ちになるかもわかる。
すごいぞ、ロジスティック方程式。
ロトカボルテラ方程式(競争系)
言いにくい。「ろとかぼるてらほうていしき」。
ねずみが増える際には餌が必要だが、他種と餌を取り合うようになることもある。
この「競争」という概念を組み込んだのがロトカボルテラ方程式
ねずみをn1、他種をn2とすると
Δn1 = a1*n1[t-1]*(1-(n1[t-1]+x*n2[t-1])/capacity1))
Δn2 = a2*n2[t-1]*(1-(n2[t-1]+y*n1[t-1])/capacity2))
x,yは他種への圧迫度合いを示す定数で、初期値とこの定数、capacityによってどちらが競争優位になるかが話かかる。
ロトカボルテラ方程式(捕食系)
捕食、被食という概念を取り込んだバージョン。
ねずみn1を食べる捕食者n2がいるとすると、
Δn1 =n1[t-1]*(α-β*n2[i-1])
Δn2 =-n2[t-1]*(γ-δ*n1[i-1])
と表せる
α〜δは定数でそれぞれ子孫を増やす度合いだったり捕食される度合いを表しているのだが、n1とn2が相互に関係しているので、αを増やせばn1>n2になるわけでもない。
これらを駆使することで、現状のデータがどう動いていくかを推定できる。
考えた人頭いいなぁ。何食ったら思いつくんだろう。
