2021-04-18

新型コロナウイルス感染者数予測をSIRモデルで表現（Python）

$\frac{dy(t)}{dt} = \beta x(t)y(t)-\gamma y(t)$

新型の病気のため、人々は免疫を持っていない
モデルとする人間は、平均的に不特定多数の人と接触している
人口の流入、流出はない（ロックダウンされた国や都市を想定したモデル）
接触によって、一定の確率 $p$ で感染する
感染者は快復後に免疫を獲得し、再び罹患者と接触しても再発しない、未感染者に罹患させることもない。

ここでは、S,I,Rを用いて
S：(Susceptable)未感染者
I : (Infected)感染者
R : (Recovered)快復者
を表すこととする。

総人口を $N$ 、未感染者 $(S)$ が $x$ 人、感染者 $(I)$ が $y$ 人、快復した人 $(R)$ が $z$ 人と置く。また、経過時間を変数 $t$ （日）で表すことにする。集団の中で、一人当たり平均 $m$ 人と接触するとする。そのうち感染している人の割合は $y(t)/N$ であるから、未感染者 $x(t)$ 人が一日あたりに接触する感染者の人数は $m\frac{y(t)}{N} x(t)$ 人と表せる。さらに、未感染者が感染者と接触した際に感染する確率を $p$ とすれば、 $\Delta t$ 日の間に生じる新たな感染者は、 $m\frac{y(t)}{N} x(t)p\Delta t$ となる。

$\beta = \frac{mp}{N}$ と置くと、 $t$ 日目から $t+\Delta t$ 日目の間に、新たな感染者によって非感染者の数が

$x(t+\Delta t)-x(t) = -\beta x(t)y(t)\Delta t$
だけ変化することになる。この式の $\Delta t$ を左辺の分母に移行し、
$\frac{x(t+\Delta t)-x(t)}{\Delta t} = -\beta x(t)y(t)$
と書き直して、 $\Delta t \rightarrow 0$ の極限をとると
$\frac{dx(t)}{dt} = -\beta x(t)y(t)$
を得る。次に、感染者は一日当たり一定の率 $\gamma$ で快復してゆくとすると、回復者の増加分は
$\frac{dz(t)}{dt} = \gamma y(t)$
人口は一定であるため、(非感染者)+(感染者)+(快復者)＝ $N$
$x(t)+y(t)+z(t) = N$
両辺を $t$ で微分すると
$\begin{align} \frac{dy}{dt} &= -\frac{dx(t)}{dt}-\frac{dz(t)}{dt}\\ &= \beta x(t)y(t)-\gamma y(t)\cdots (1) \end{align}$

これをPythonで数値シミュレーションする。4次のルンゲクッタ法を用いて微分方程式からグラフ描写を行う。
$\begin{align} &k_1 \leftarrow f(x, y(x)) \\ &k_2 \leftarrow f(x+h/2, y(x)+k_1 h/2)\\ &k_3 \leftarrow f(x+h/2, y(x)+k_2 h/2)\\ &k_4 \leftarrow f(x+h, y(x)+k_3 h)\\ &y(x+h) \leftarrow y(x)+(k_1 +2k_2+2k_3+k_4)h/6 \end{align}$

import math
import matplotlib.pyplot as plt

N = 1000000
x = N-10
y = 10
z = 0

taxis = []
xaxis = []
yaxis = []
zaxis = []

beta = 0.2/N
gamma = 0.07
dt = 0.001

t = 0
cnt = 0
while t < 180:
    if cnt%100 == 0:
        taxis.append(t)
        xaxis.append(x)
        yaxis.append(y)
        zaxis.append(z)
# step 1
    kx1 = -beta*x*y
    ky1 = beta*x*y - gamma*y
# step 2
    x2 = x+kx1*dt/2
    y2 = y+ky1*dt/2
    kx2 = -beta*x2*y2
    ky2 = beta*x2*y2 - gamma*y2
# step3
    x3 = x+kx2*dt/2
    y3 = y+ky2*dt/2
    kx3 = -beta*x3*y3
    ky3 = beta*x3*y3 - gamma*y3
# step4
    x4 = x+kx3*dt
    y4 = y+ky3*dt
    kx4 = -beta*x4*y4
    ky4 = beta*x4*y4 - gamma*y4

# update
    x = x+(kx1+2*kx2+2*kx3+kx4)*dt/6
    y = y+(ky1+2*ky2+2*ky3+ky4)*dt/6
    z = N-x-y
    t = t+dt
    cnt+=1

plt.title("SIR MODEL")
plt.plot(taxis,xaxis, color=(0.0,1,0.0), linewidth=1.0, label='S')
plt.plot(taxis,yaxis, color=(1.0,0,0.0), linewidth=1.0, label='I')
plt.plot(taxis,zaxis, color=(0.0,0,1.0), linewidth=1.0, label='R')
plt.xlim(0,180)
plt.legend()
plt.xlabel('DAY')
plt.ylabel('X, Y, Z')
plt.grid(True)

f:id:nya__nya:20210418005530p:plain

参考文献

情報基礎「Pythonプログラミング」（ステップ６：SIRモデル）
情報基礎「Pythonプログラミング」（ステップ６：ルンゲクッタ法）

2021-04-17

Nim言語動的配列

var x: seq[int]
x = @[1, 2, 3, 4, 5, 6]
for index, value in x:
    echo index, ":", value 

var v = newseq[int](5)
v.fill(1)
v.add(2)
for e in v:
    echo e

2021-04-15

Nim言語　プロシージャの使い方

int型の整数を引数に持ち、string型を返す関数f
nが偶数の時に、"Even", 奇数の時に"Odd"を返す。
modは割った時の余りを求める演算子、C++で言うところの%になります。

proc f(n: int): string =
    if n mod 2 == 0:
        return "Even"
    else:
        return "Odd"

var t = stdin.readLine.parseInt()
echo f(t)

次に戻り値がない関数

proc f(n: int) =
    if n mod 2 == 0:
        echo "Even"
    else:
        echo "Odd"

var t = stdin.readLine.parseInt()
f(t)

参考文献
flat-leon.hatenablog.com

2021-04-02

Google Apps Script でwebスクレイピングした地震情報をmail送信で知らせるbotをつくる

今回実装した、Google Apps Scriptの全てのコードがこちら。

function myFunction() {
  const url = "https://earthquake.tenki.jp/bousai/earthquake/entries/";
  let html = UrlFetchApp.fetch(url).getContentText('UTF-8');
  
  var data = Parser.data(html).from('<table class="earthquake-entries-table">').to('<table class="earthquake-entries-table">').build();
  
  var data_list = Parser.data(html).from('<tr').to('</tr>').iterate();

  var sheet = SpreadsheetApp.getActiveSheet()

  for(let i=0; i<data_list.length; ++i){
    var s = data_list[i].replace(/<("[^"]*"|'[^']*'|[^'">])*>|\r\n|\n|\r/g,'');
    s = s.replace(/>.*?2/g, "2");
    s = s.replace(/2.*?年/g, "");
    sheet.getRange(i+1, 1, 1, 1).setValue(s);
  }

  var text = sheet.getRange(2, 1, 4, 1).getValues();

  const to = '@mailaddres';
  const subject = '地震情報';
  var body = "日本の地震情報をお届けします。\n\n"+text[0]+"\n"+text[1]+"\n"+text[2]+"\n"+text[3]+"\n\n";
  body+="以上が日本で発生した直近4つまでの地震となります。\n";
  body+="詳しくは気象庁のホームページをご覧ください。\n";
  body+="https://earthquake.tenki.jp/bousai/earthquake/entries/";
  
  const options = {
    cc: '',
    bcc:''
  }
  GmailApp.sendEmail(to, subject, body, options); 
}

1.環境セッティング

まず最初に、google driveからスプレッドシートを開く、好きな名前に変更したら、ツール欄からスクリプトエディタを開く。エディタの名前を変更する。
次に、Parserというライブラリを追加して、webスクレイピングをおこなう。
ライブラリ追加のidの欄に以下の文字列を入力してParserをライブラリに追加する。
2021年4月2日現在ではParserの8が最新バージョンとなっている。以下のIdもそれを追加するためのものであるため、適宜最新バージョンを選んで追加してください。

1Mc8BthYthXx6CoIz90-JiSzSafVnT6U3t0z_W3hLTAX5ek4w0G_EIrNw

2. webスクレイピングをする

過去の地震情報 (日付の新しい順) - 日本気象協会 tenki.jp
今回、上のwebページから日本で発生した直近4つの地震の情報ををGmailで送信するbotを作成する。
まず、webページを開いたら、F12キーを押して開発者ツールを開く。
開発者ツールを開くと384行目辺りから以下のようなhtmlのコードが存在する。これが、直近の地震情報が記述されている部分である。まず最初にこのブロック部分のhtmlを取得する。

  <div class="pager-index">1件～100件(全30103件)</div>

  <table class="earthquake-entries-table">
    <tr>
      <th class="image">&nbsp;</th>
      <th class="datetime"><a href="/bousai/earthquake/entries/asc/" class="text-link">▼発生時刻</a></th>
      <th class="center">震源地</th>
      <th class="magnitude"><a href="/bousai/earthquake/entries/magnitude/" class="text-link">マグニチュード</a></th>
      <th class="max-level"><a href="/bousai/earthquake/entries/max-level/" class="text-link">最大震度</a></th>
    </tr>
    <tr>
      <td class="image"><a href="/bousai/earthquake/detail/2021/04/02/2021-04-02-09-31-27.html"><img src="https://earthquake.tenki.jp/static-images/earthquake/detail/2021/04/02/2021-04-02-09-31-27-small.jpg" width="60" height="45"></a></td>
      <td class="datetime"><a href="/bousai/earthquake/detail/2021/04/02/2021-04-02-09-31-27.html" class="text-link">2021年04月02日09時31分頃</a></td>
      <td class="center">鳥取県中部</td>
      <td class="magnitude">M2.8</td>
      <td class="max-level"><img src="https://static.tenki.jp/images/icon/earthquake/level_1.png" alt="1" width="21" height="21"></td>
    </tr>
    
    <tr>
      <td class="image"><a href="/bousai/earthquake/detail/2021/04/02/2021-04-02-08-36-44.html"><img src="https://earthquake.tenki.jp/static-images/earthquake/detail/2021/04/02/2021-04-02-08-36-44-small.jpg" width="60" height="45"></a></td>
      <td class="datetime"><a href="/bousai/earthquake/detail/2021/04/02/2021-04-02-08-36-44.html" class="text-link">2021年04月02日08時36分頃</a></td>
      <td class="center">十勝地方中部</td>
      <td class="magnitude">M3.4</td>
      <td class="max-level"><img src="https://static.tenki.jp/images/icon/earthquake/level_2.png" alt="2" width="21" height="21"></td>
    </tr>

..................以下略


    <tr>
      <td class="image"><a href="/bousai/earthquake/detail/2021/03/15/2021-03-15-00-59-19.html"><img src="https://earthquake.tenki.jp/static-images/earthquake/detail/2021/03/15/2021-03-15-00-59-19-small.jpg" width="60" height="45"></a></td>
      <td class="datetime"><a href="/bousai/earthquake/detail/2021/03/15/2021-03-15-00-59-19.html" class="text-link">2021年03月15日00時59分頃</a></td>
      <td class="center">和歌山県北部</td>
      <td class="magnitude">M3.3</td>
      <td class="max-level"><img src="https://static.tenki.jp/images/icon/earthquake/level_3.png" alt="3" width="21" height="21"></td>
    </tr>
      </table><!-- /#earthquake-entries-table -->

まず最初にこのブロック部分のhtmlを取得する。

var data = Parser.data(html).from('<table class="earthquake-entries-table">').to('</table><!-- /#earthquake-entries-table -->').build();

次に地震1回ごとに分割してリスト構造としてdata_listに取り出す。

var data_list = Parser.data(html).from('<tr').to('</tr>').iterate();

2021-03-31

Nim言語で bit全探索

集合 $\{1, 2, 3, 4\}$ の部分集合を全て列挙するプログラム

import sequtils

let n = 4
for bit in 0..<(1 shl n):
    var vec = newSeq[int]()
    for i in 0..<n:
        if(bit and (1 shl i)) != 0:
            vec.add(i+1)
    echo bit, " : ", vec

実行結果は以下の通り

0 : @[]
1 : @[1]
2 : @[2]
3 : @[1, 2]
4 : @[3]
5 : @[1, 3]
6 : @[2, 3]
7 : @[1, 2, 3]
8 : @[4]
9 : @[1, 4]
10 : @[2, 4]
11 : @[1, 2, 4]
12 : @[3, 4]
13 : @[1, 3, 4]
14 : @[2, 3, 4]
15 : @[1, 2, 3, 4]

2021-03-29

ABC197-D Opposite, 複素数の計算

atcoder.jp
$p_0 = (x_0, y_0), p_{n/2}=(x_{n/2}, y_{n/2})$
$p_h = (\frac{x_0+x_{n/2}}{2}, \frac{y_0 + y_{n/2}}{2})=(x_h, y_h)$
$\alpha = \frac{2\pi}{n}$
$z_0 = x_0+iy_0$
$z_h=x_h+iy_h$
$t = (cos\alpha+i sin\alpha)$
$z_0-z_h = r(cos\theta+i sin\theta )$
$z_1 = (z_0-z_h)\times t+z_h = r(cos(\theta+\alpha)+i sin(\theta+\alpha))+(x_h+iy_h)$

from math import cos, sin, atan2
n = int(input())
x0, y0 = map(int, input().split())
x2, y2 = map(int, input().split())
xh = (x0+x2)/2
yh = (y0+y2)/2

r = ((x0-x2)**2+(y0-y2)**2)**0.5
r/=2

A = atan2((y0-yh), (x0-xh))
x = r*cos(A+(2*pi)/n)+xh
y = r*sin(A+(2*pi)/n)+yh

print(x, y)

from cmath import rect
from math import pi
n = int(input())
x0, y0 = map(int, input().split())
x2, y2 = map(int, input().split())

p0 = complex(x0, y0)
p2 = complex(x2, y2)

o = (p0+p2)/2
r = rect(1, 2*pi/n)
a = (p0-o)*r+o

print("%.10f %.10f"%(a.real, a.imag))

2021-03-28

ABC197 A, B, C, D問題をNim言語で解く

atcoder.jp

let s = stdin.readLine
echo s[1..2] & s[0]

atcoder.jp

import sequtils, strutils, strformat, algorithm, math, sugar

var
  h, w, x, y, ans: int

(h, w, x, y) = stdin.readLine.split.map(parseInt)
let s = (0..<h).mapIt(stdin.readLine)

x-=1
y-=1
ans = 0

if s[x][y]=='.': ans+=1

for j in y+1..<w:
  if s[x][j]=='#': break
  else: ans+=1

for j in 1..y:
  if s[x][y-j]=='#': break
  else: ans+=1

for i in x+1..<h:
  if s[i][y]=='#': break
  else: ans+=1

for i in 1..x:
  if s[x-i][y]=='#': break
  else: ans+=1

echo ans

atcoder.jp

import sequtils, strutils, strformat, algorithm, math, sugar, complex

let 
    n = stdin.readLine.parseInt()
    A = stdin.readLine.split.map(parseInt)
var ans:int64 = 10000000000000000

for bit in 0..<(1 shl (n-1)):
    var xorSum:int64 = 0
    var orSum:int64 = A[0]
    for i in 0..<(n-1):
        if(bit and (1 shl i)) != 0:
            xorSum = (xorSum xor orSum)
            orSum = A[i+1]
            if i == n-2: xorSum=(xorSum xor A[n-1])
        else:
            orSum=(orSum or A[i+1])
            if i == n-2: xorSum=(xorsum xor orSum)
    if xorSum<ans: ans = xorSum

if n == 1: echo A[0]
else: echo ans

atcoder.jp

import complex, math, strutils, sequtils

var
    x0, y0, x2, y2: float64
var n: float64 = stdin.readLine.parseFloat()
(x0, y0) = stdin.readLine.split.map(parseFloat)
(x2, y2) = stdin.readLine.split.map(parseFloat)

let
    p0 = complex[float64](x0, y0)
    o = complex[float64]((x0+x2)/2, (y0+y2)/2)
    r = rect[float64](1.0, 2*PI/n)

var ans = (p0-o)*r+o
echo ans.re, " ", ans.im

メモ帳

楽しいアウトプットの場所

新型コロナウイルス感染者数予測をSIRモデルで表現（Python）

参考文献

Nim言語動的配列

Nim言語　プロシージャの使い方

Google Apps Script でwebスクレイピングした地震情報をmail送信で知らせるbotをつくる

1.環境セッティング

2. webスクレイピングをする

Nim言語で bit全探索

ABC197-D Opposite, 複素数の計算

ABC197 A, B, C, D問題をNim言語で解く