新型コロナウィルスの第６波を予測して見た


# はじめに
今回は、新型コロナウィルスの第６波の予測を、SIRモデルを用いながらしたいと思います。

# 必要なモジュール

プログラムコードは、python

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline
import pandas as pd
plt.rcParams['font.family'] = 'Hiragino Sans'
import datetime
```

# データセット
[厚生労働省](https://www.mhlw.go.jp/stf/covid-19/open-data.html)のサイトの「新規陽性者数の推移（日別）」よりダウンロードし、`第６波生データ.csv`とファイル名を変更した。

```python
df = pd.read_csv('第６波生データ.csv', header=0)
df.head()
```

![Screen Shot 2022-01-28 at 21.45.46.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/c1bcbf88-01c0-5f0f-4a9e-dab4beb5c0a2.png)

# 表の可視化
## 大阪府

```python
osaka = df[['Date','Osaka']]
osaka.head()
```

![Screen Shot 2022-01-28 at 21.46.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/319675b4-5565-4c6a-e0fc-7386ba030814.png)

## 愛知県

```python
aichi = df[['Date','Aichi']]
aichi.head()
```

![Screen Shot 2022-01-28 at 21.46.31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/b71fbac8-da28-b7ba-0a57-c3e582592632.png)

## 愛知県、京都府、大阪府

```python
aichikyotoosaka = df[['Date','Aichi','Kyoto','Osaka']]
aichikyotoosaka.head()
```

![Screen Shot 2022-01-28 at 21.46.51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/38da8c4d-b063-e666-bc18-799396eccc9a.png)

# グラフによる可視化

## 大阪府

```python
plt.bar(np.arange(osaka['Date'].size), osaka['Osaka'].values)
plt.xlabel("コロナの日数")
plt.ylabel("新規感染者数")
plt.title('大阪府新規感染者数推移')
```

![Screen Shot 2022-01-28 at 21.43.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/3387cf13-8702-5508-801d-f46b7cc569c6.png)

## 愛知県

```python
plt.bar(np.arange(aichi['Date'].size), aichi['Aichi'].values)
plt.xlabel('コロナの日数')
plt.ylabel('新規感染者数')
plt.title ('愛知県新規感染者数推移')
```

![Screen Shot 2022-01-28 at 21.47.31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/83614d02-49e9-add6-fa68-d38a686d4cde.png)

# 第6波の抽出
図より700日後からを第6波とした。

```python
infected=osaka[700:]["Osaka"].values
plt.bar(np.arange(infected.size), infected)
```

![Screen Shot 2022-01-28 at 21.49.04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/88a3e748-e641-b185-3dfe-b1030ff4ab69.png)

```python
infected=aichi[700:]["Aichi"].values
plt.bar(np.arange(infected.size), infected)
```

![Screen Shot 2022-01-28 at 21.49.30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/5c628189-ff5e-3f5e-ba05-dda35649f73f.png)

# 予測モデル
今回は、SIRモデルを用いた。
図は[このサイト](https://club.informatix.co.jp/?p=140)より引用。

![Screen Shot 2022-01-30 at 10.39.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/fbad3f47-7261-2ea9-fc39-4d0238d8dc59.png)

![Screen Shot 2022-01-28 at 21.51.44.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/61ef0db6-0bbb-6680-e51c-3b6ad963ded9.png)

```python
def SIR_function (sir, t, beta, gamma):
    s,i,r=sir
    ds_dt=-beta*s*i/N
    di_dt=beta*s*i/N-gamma*i
    dr_dt=gamma*i
    return np.array([ds_dt,di_dt,dr_dt])
```
## 大阪府の予測
$N$に大阪府の人口、$\beta,\gamma$は適当な値で予測。

```python
N=8850000
beta0,gamma0=0.056,0.045
predict_times=np.arange(200)
times=np.arange(infected.size)

r0=osaka[:700]["Osaka"].values.sum()
i0=infected[0]
sir0=(N-r0+i0,i0,r0)

from scipy.integrate import odeint

predict=odeint(SIR_function,sir0,predict_times,args=(beta0,gamma0))
i=predict[:, 1]

plt.bar(times,infected,label="actual")
plt.plot(predict_times,i,label="predicted", color='orange')
plt.legend()
```

当然ながらうまく予測できない。

![Screen Shot 2022-01-28 at 21.52.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/78cf870b-e15d-7dfc-a8f8-7345b921b3bf.png)

そこで、`scipy`の`curve_fit`を用いて、$\beta,\gamma$の最適値を求める。

```python
from scipy.optimize import curve_fit

def return_infected (t, beta, gamma):
    predict=odeint(SIR_function,sir0,t,args=(beta,gamma))
    i=predict[:, 1]
    return i

best_params, cov=curve_fit(return_infected, times, infected, (beta0, gamma0), bounds=((0, 0,), (5, 5)))
```

$\beta=4.99999999,\gamma=4.68203601$が最適値であると算出された。

```python
best_params = array([4.99999999, 4.68203601])

best_beta, best_gamma = best_params
predict=odeint(SIR_function,sir0,predict_times,args=(best_beta,best_gamma))
predict_infected=predict[:, 1] 

plt.bar(times,infected,label="actual")
plt.plot(predict_times,predict_infected,label="predicted", color='orange')
plt.legend()
```

良い当てはまり具合だと言える。

![Screen Shot 2022-01-28 at 21.53.52.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/5f1d3d7f-51d8-244b-89bb-bfa10501ea8d.png)

## 横軸を日付に変更

```python
startdate = osaka[('Date')][700]

startdate
```
 
![Screen Shot 2022-01-28 at 21.54.32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/54f5a3dd-174e-8dca-f032-3fb3cfb40cc5.png)

```python
startdatetime = datetime.datetime.strptime(startdate,'%Y/%m/%d').date()
predict_datetimes = [startdatetime + datetime.timedelta(days=int(day)) for day in predict_times]
predict_datetimes
```


```python
best_beta, best_gamma = best_params
predict=odeint(SIR_function,sir0,predict_times,args=(best_beta,best_gamma))
predict_infected=predict[:, 1] 

plt.bar(predict_datetimes[:times.size],infected,label="actual")
plt.plot(predict_datetimes,predict_infected,label="predicted", color='orange')
plt.legend()
plt.xticks(rotation=90)
```

![Screen Shot 2022-01-29 at 9.24.45.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2518583/00dd80d0-e5ed-9d53-cdb8-01e131756667.png)

# 終わりに

筆者が、予測をした際は1月23日であった。その時の、大阪府の新規感染者数は６２１７人だった。このSIRモデルによると、大阪府の新規感染者数は、2月にピークを迎える。どんどん増え続けている新規感染者数であるが、2月には減少傾向を迎えるだろう。よって、SIRモデルがコロナの新規感染者数を予測するのに、重要なモデルであることを熟知した。これからは、日本国民の３回目のワクチン接種が進めば、第７波は防ぐことができるであろう。