{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
},
"colab": {
"name": "2020Kermack_McKendrick_model.ipynb",
"provenance": [],
"include_colab_link": true
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "nIAyV7yU40Fj",
"colab_type": "text"
},
"source": [
"# 疫学の数理モデルによる導入\n",
"\n",
"- Date: 2020-0320\n",
"- Author: 浅川伸一 asakawa@ieee.org\n",
"\n",
"CoVID-19 の理解に向けて\n",
"\n",
"Kermack と McKendrick は 1905 年から 1906 年にかけてインドで発生したコレラの大流行についての疫学的流行モデルを提案した。\n",
"\n",
"- $y$: population infected 感染者数\n",
"- $x$: number of individuals still unaffected 未感染者数\n",
"- $z$: number who have been removed by recover and death 地域社会から隔離された人口(理由は問わない,隔離,死去,免疫獲得など)\n",
"- $N$: population density 地域社会の総人口 $N=x+y+z$\n",
"\n",
"$\\kappa$ を感染率,$\\ell$ を隔離係数として Kermack=McKendrick モデルは次式で表される:\n",
"\n",
"$$\n",
"\\begin{cases}\n",
"\\frac{dx}{dt} &= - \\kappa x y\\\\\n",
"\\frac{dy}{dt} &= \\kappa x y - \\ell y\\\\\n",
"\\frac{dz}{dt} &= \\ell y\n",
"\\end{cases}\n",
"$$\n",
"\n",
"- Reference:\n",
"Kermack and McKendrick (1927) A Contribution to the Mathematical Theory of Epidemics.\n"
]
},
{
"cell_type": "code",
"metadata": {
"id": "YZCYF2jZ40Fm",
"colab_type": "code",
"colab": {}
},
"source": [
"import numpy as np\n",
"from IPython.display import Image\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "oD497lby40Fp",
"colab_type": "code",
"colab": {}
},
"source": [
"def KM(x0=99, y0=1, length=100, kappa=1.0, ell=1.0, dt=0.005):\n",
" \"\"\"return Kermack and McKendrick models output\"\"\"\n",
"\n",
" X, Y = [], []\n",
" x, y = x0, y0\n",
" for _ in range(length):\n",
" dx = (-kappa * x * y) * dt\n",
" dy = (kappa * x * y - ell * y) * dt\n",
" x += dx\n",
" y += dy\n",
" X.append(x)\n",
" Y.append(y)\n",
"\n",
" return X, Y\n",
"\n",
"kappa = 1.0\n",
"l1 = 10.\n",
"l2 = 50.\n",
"\n",
"_, Y1 = KM(kappa=1.0, ell=l1, length=100)\n",
"X, Y2 = KM(kappa=1.0, ell=l2, length=100)\n",
"T = range(len(X))\n",
"\n",
"plt.figure(figsize=(8,6)) # width, height inches\n",
"plt.plot(T, X, c='g', label='non infected') # X軸を T, Y軸を X(非感染者) 色(c)を緑(g)で描画\n",
"plt.plot(T, Y1, c='y', label='Y1 infected') # X軸を T, Y軸を Y1(感染者, l=l1) 色(c)を黄色(y)で描画\n",
"plt.plot(T, Y2, c='r', label='Y2 infected') # X軸を T, Y軸を X2(感染者, l=l2) 色(c)を緑(r)で描画\n",
"plt.legend()\n",
"# plt.savefig('KM_model_output.png') グラフの保存"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "ic1bFlJN40Ft",
"colab_type": "text"
},
"source": [
"$\\sinh(x) = \\frac{\\exp(x)-\\exp(-x)}{2}$\n",
"\n",
"$\\cosh(x) = \\frac{\\exp(x)+\\exp(-x)}{2}$\n",
"\n",
"$\\tanh(x) = \\frac{\\exp(x)-\\exp(-x)}{\\exp(x)+\\exp(-x)}$\n",
"\n",
"$\\text{sech}(x) = \\frac{1}{\\cosh(x)}$\n",
"\n",
"$\\text{cosech}(x) = \\frac{1}{\\sinh(x)}$\n",
"\n",
"$\\coth(x) = \\frac{\\cosh(x)}{\\sinh(x)}$\n"
]
},
{
"cell_type": "code",
"metadata": {
"id": "as8i3XY540Fu",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 198
},
"outputId": "db2521b3-5b73-41a4-a392-24ac1ac3b753"
},
"source": [
"x = np.linspace(-3,3)\n",
"plt.figure(figsize=(8,8))\n",
"plt.plot(x, 1/np.cosh(x) **2)"
],
"execution_count": 1,
"outputs": [
{
"output_type": "error",
"ename": "NameError",
"evalue": "ignored",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcosh\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mNameError\u001b[0m: name 'np' is not defined"
]
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "t3gVOhGn40F0",
"colab_type": "code",
"colab": {}
},
"source": [
"!wget https://raw.githubusercontent.com/ShinAsakawa/ShinAsakawa.github.io/master/assets/1927Kermack_McKendrick_fig.png"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "osfTvpzXz2QS",
"colab_type": "code",
"colab": {}
},
"source": [
"import matplotlib.pyplot as plt\n",
"import PIL.Image \n",
"img = PIL.Image.open('1927Kermack_McKendrick_fig.png')\n",
"\n",
"plt.figure(figsize=(10, 16))\n",
"plt.axis(False); plt.imshow(img)\n"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "CiTgpCHh0ESm",
"colab_type": "code",
"colab": {}
},
"source": [
""
],
"execution_count": null,
"outputs": []
}
]
}