よちよちpython

独習 python/Qpython/Pydroid3/termux/Linux

トップ > numpy > 【Numpyだけで単回帰分析】np.polyfit()とnp.poly1d()でコロナ検査数から陽性者数を予測する

【Numpyだけで単回帰分析】np.polyfit()とnp.poly1d()でコロナ検査数から陽性者数を予測する

numpy matplotlib Android JupyterNotebook pandas termux グラフ作成 機械学習

np.polyfit()とnp.poly1d()でコロナ陽性者数を単回帰分析予測

今回は、コロナの検査数と陽性者数のデータから、Numpyを使って単回帰分析を行ってみます。久しぶりでやり方忘れてるので復習。

独立変数xを検査件数、目的変数yを陽性者数として回帰分析をし、 できたモデルに任意の検査数を当てはめると、一体どれほどの陽性者数が出るかを予測します。



【実行環境】

  • Android
  • Termux
  • Python3.9
  • Jupyter Notebook
  • 外部ライブラリ
    • pandas, numpy, matplotlib, japanize-matplotlib(日本語表示用)



目次



データとライブラリのインポート

# データファイル(厚労省オープンデータからDL)
f1='pcr_tested_daily.csv'
f2='pcr_positive_daily.csv'

# インポート
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
from matplotlib import dates as mdates



データフレームに読み込み

# ファイル読み込み 時系列でデータフレーム化
df1=pd.read_csv(f1,index_col=0, parse_dates=True)
df2=pd.read_csv(f2,index_col=0, parse_dates=True)

# 結合
df = df1.join(df2)
# 表示
df
PCR 検査実施件数(単日) PCR 検査陽性者数(単日)
日付
2020-02-05 4 2
2020-02-06 19 0
2020-02-07 9 0
2020-02-08 4 0
2020-02-09 10 0
... ... ...
2021-06-29 61916 1375
2021-06-30 67576 1811
2021-07-01 77063 1741
2021-07-02 67044 1774
2021-07-03 38122 1864

512 rows × 2 columns

日付ごとの検査件数と陽性者数のデータフレームができました。



相関関係を見てみる

検査数と陽性者数の間には何らかの関係があるのだろうか。



相関関係を見ます。
相関関係を示す相関係数は、-1～1の範囲の値を取ります。
相関係数の意味は、

  • マイナス1に近い : 強い負の相関関係
  • 0 : 相関関係はない
  • プラス1に近い : 強い正の相関関係



  • 正の相関関係がある場合は、片一方が上がるともう片一方も上がるような関係、
  • 負の相関関係の場合は、片一方が上がるともう片一方は下がるような関係。
  • 相関係数が0に近いと、関係性はない。
# データフレームまるごと相関係数を出す
df.corr()
PCR 検査実施件数(単日) PCR 検査陽性者数(単日)
PCR 検査実施件数(単日) 1.000000 0.728319
PCR 検査陽性者数(単日) 0.728319 1.000000

リーグ表のように見ます。
1行目の2列目「陽性者数」の値は0.72なので、強い正の相関関係があるようです。

「検査数が上がると陽性者数も上がるような関係がある」ことを示しています。

陽性者数は当然ですけど検査をしないと判明しません。どれぐらい検査すればどれだけの陽性者数になるかが、回帰分析で大体予想できるかもしれない。



グラフで確かめる

  • 横軸 : 時間軸
  • 縦軸 : 「検査数」「陽性者数」の2本軸

でグラフを描いてみます。

# 縦2軸の折れ線グラフ
fig = plt.figure(facecolor='lightgray')
ax1 = fig.add_subplot(111)

# プロット
l1 = ax1.plot(df.iloc[:,0], 'm-', label='検査件数')
ax2 = ax1.twinx()
l2 = ax2.plot(df.iloc[:,1], 'c2', label='陽性者数')

# 凡例の設置
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2) 

ax1.set_xlabel('日付')    # 横軸 ラベル設置
ax1.set_ylabel('検査件数') # 縦軸1 ラベル設置
ax2.set_ylabel('陽性者数') # 縦軸2 ラベル設置

ax1.grid() # グリッド線設置
plt.title('検査件数-陽性者数2軸グラフ') # タイトル設置
plt.gcf().autofmt_xdate() # 横軸重なり解消
plt.savefig('検査件数と陽性者数の2軸グラフ.jpg') # グラフ画像保存
plt.show()

f:id:chayarokurokuro:20210719173821j:plain

紫の「検査件数」と、青の「陽性者数」は形が似ています。強い正の相関関係がある為です。





単回帰分析

上のグラフは横軸を時間軸に、縦軸を検査数と陽性者数の2軸で表していますが、では、

  • 横軸 : 検査数
  • 縦軸 : 陽性者数

にしたグラフを描くとどのような形になるでしょうか。散布図を描きます。

相関関係が正の強い相関でしたので、右肩上がりのグラフになりそうですが。

陽性者数/検査件数のグラフ

# 陽性者数/検査件数のグラフ
plt.figure(facecolor='w')

plt.scatter(df.iloc[:,0], df.iloc[:,1])

plt.title('陽性者数/検査件数グラフ')
plt.xlabel('検査件数')
plt.ylabel('陽性者数')
plt.grid()

plt.savefig('陽性者数-検査件数グラフ.jpg')
plt.show()

f:id:chayarokurokuro:20210719173937j:plain

バラけてはいますが、何となく右肩上がりの一次方程式で近似されそうな感じ？

冬はウィルス感染症が増えます。インフルエンザも他の風邪も大体冬に流行ります。一方で夏は少ない。現在コロナが少ないのはワクチンのおかげではなく、夏だからでしょう。

同じ検査数でも季節や気温、湿度などによって陽性者数が変わって来ますので、上のようにバラついたグラフになっているかと思われます。季節を絞ったデータでグラフを描けばもっとシャープな感じの散布図になるかも。

フィッティング

近似式を作ります。上のバラけた散布図を１つの方程式で表そうという強引な方法です。

その方程式に検査件数を代入すれば、陽性者数が予測できる。Ψ(｀∀´)Ψｹｹｹ



numpy.polyfit(x,y,deg=次数)を使うと、

  • y : 陽性者数
  • x : 検査件数



として、

  • y = a1(xn) + a2(xn-1) + a3(xn-2)・・・an(x) + b

みたいな関係式(多項回帰式)の係数a1,a2,a3,...,と、切片bがリスト型の戻り値で得られます。


polyfitの内部では最小二乗法というアルゴリズムが実装されているようです。

np.polyfit()の引数deg=1とすると、もっとも単純な一次方程式

  • y=ax+b

の傾きaと切片bを得ます。

やってみましょ。

# 列を分けて取り出す(最終行を除く)
n1 = df.iloc[:-2,0]
n2 = df.iloc[:-2,1]

最終行を除いたのは、回帰分析でできたモデルに最終行の検査件数を入れるとその陽性者数が出るかを確認する用に外しました。

# データフレームの最終行を表示
df.tail(1)
PCR 検査実施件数(単日) PCR 検査陽性者数(単日)
日付
2021-07-03 38122 1864

この検査数38122を回帰分析で出来上がったモデルに入れて出た値が陽性者数1864に近ければ、そのモデルは予測精度が高いと言える。



傾きと切片を出します。

# フィッティング (一次近似式)
fit = np.polyfit(n1, n2, deg=1)
fit

array([4.09474545e-02, 2.93552085e+02])

バラついたグラフを一次の式y=ax+bで近似した傾きaと切片bがリストで得られました。

で、y=ax+bにリスト(傾きと切片)の値を入れて式を書いて、そのできた式に検査数xを渡せば陽性者数yが予測できます。

因みに、近似の多項回帰式を1次でやりましたが、deg=5でやってみますと、

# フィッティング (5次近似式)
fit5 = np.polyfit(n1, n2, deg=5)
fit5

array([ 7.73878014e-23, -9.75502661e-17,  2.34671726e-11, -2.02461192e-06,
        1.03421075e-01, -6.56222138e+01])

全部で6個の値が出てきました。5次方程式の各定数5個と、最後の定数項の1個。



予測する

さて、あとはy=ax+bの式に得られたaとbを入れれば式が完成しますが、numpy.poly1d()関数の引数に先ほどのリストを渡せば、わざわざ式を書かずに済ませられます。

# 近似式に傾きと切片のリストを代入してモデル式を作る
func = np.poly1d(fit)

# できたモデル式にデータフレームの最終行の検査件数を入れると…
func(df.iloc[-1,0])

1854.5509455218305

これがデータフレームの最終行にある検査件数から予測したその時の陽性者数です。実際の陽性者数は、

# 最終行の陽性者数
df.iloc[-1,1]

1864

ほぼ同じ。
予測精度の高いモデルができた。

そしたらば、検査件数を100万にすると、どれぐらい陽性者数が出るか見てみましょ。

# 検査件数100万のときの陽性者数予測値
func(1000000)

41241.006590238554

約4万1千人の陽性者数が出ると予測された。検査100万で陽性4万ですので陽性率は約4%です。

近似式の傾き(リストの最初の値)が陽性率を表しているということになります。



たとえばイギリスでは、１日に100万を超える検査をして3万を超える陽性者数が出ている。かたや日本では、上のデータフレームを見て分かるように、１日に3万～6万ほどの検査数で、陽性者数は2000人ほど。
陽性者数だけ見ると3万と2千で10倍以上違いますが、率で見るとあまり違いはありません。



一次近似式のグラフを描く

バラついた散布図に近似式を加えてグラフ化してみます。

# 陽性者数/検査件数のグラフ
plt.figure(facecolor='lightblue')

# 散布図
plt.scatter(df.iloc[:,0], df.iloc[:,1])

# 近似式直線
ystart, ylast = func(0), func(200000)
plt.plot([0, 200000],[ystart, ylast], color="r")

plt.title('陽性者数/検査件数-近似式グラフ')
plt.xlabel('検査件数')
plt.ylabel('陽性者数')
plt.grid()

plt.savefig('陽性者数-検査件数-近似直線グラフ.jpg')
plt.show()

f:id:chayarokurokuro:20210719174042j:plain



numpy.polyfit()の引数deg=1にしたので近似式が一次式になっているのですけど、degの数値を2,3,4と上げて行けば曲線の多項式で近似できます。しかし、曲線の式にした所で果たして近似と言えるかは…。過学習というやつ。

以上のことをscikit-learnでやるならば、

from sklearn.linear_model import LinearRegression

model = LinearRegression()

とか使うとできます。



おまけ 月別平均のグラフ

月ごとの陽性率をグラフで描いてみます。

df.resample('m').mean()で月ごと平均値を出してからやります。

月ごと集計

# 月ごと平均値集計
_df = df.resample('m').mean()
# 陽性率の列追加
_df['陽性率'] = _df.iloc[:,1] / _df.iloc[:,0]

# 表示
_df
PCR 検査実施件数(単日) PCR 検査陽性者数(単日) 陽性率
日付
2020-02-29 62.400000 8.400000 0.134615
2020-03-31 1045.310345 61.931034 0.059247
2020-04-30 3942.517241 418.517241 0.106155
2020-05-31 3968.290323 80.258065 0.020225
2020-06-30 4033.166667 58.266667 0.014447
2020-07-31 10134.258065 560.225806 0.055280
2020-08-31 19922.645161 1032.258065 0.051813
2020-09-30 18915.000000 503.033333 0.026594
2020-10-31 17733.064516 567.193548 0.031985
2020-11-30 25787.233333 1571.066667 0.060924
2020-12-31 41094.903226 2770.677419 0.067421
2021-01-31 60757.451613 4922.096774 0.081012
2021-02-28 48725.892857 1494.607143 0.030674
2021-03-31 50038.580645 1361.322581 0.027205
2021-04-30 64456.966667 3882.533333 0.060235
2021-05-31 80550.064516 4915.709677 0.061027
2021-06-30 64418.300000 1765.766667 0.027411
2021-07-31 60743.000000 1793.000000 0.029518

月ごと平均値の陽性率時系列グラフ

# 陽性率の時系列グラフ
plt.figure(facecolor='lightpink')
plt.plot(_df.iloc[:, 2])
plt.title('月ごと平均値時系列陽性率')
plt.xlabel('日付')
plt.ylabel('陽性率')
plt.grid()
plt.gcf().autofmt_xdate() # 横軸重なり解消
plt.savefig('月ごと平均値時系列陽性率.jpg')
plt.show()

f:id:chayarokurokuro:20210719174200j:plain

2020年の最初のほうは検査数が激しく少ないので陽性率が高く出ている。2021年の1月、2月は冬なので高い。

コロナウィルスは2週間に一回ほどの割合で変異するらしいけど、3ヶ月ぐらいで周期性があるように見えます。なんだろ？

月ごと平均値陽性率散布図

# 月ごと平均値の陽性者数/検査件数の散布図
plt.figure(facecolor='lightgreen')

plt.scatter(_df.iloc[:,0], _df.iloc[:,1])

plt.title('月ごと平均値の陽性者数/検査件数グラフ')
plt.xlabel('検査件数')
plt.ylabel('陽性者数')
plt.grid()

plt.savefig('月ごと平均値陽性率散布図.jpg')
plt.show()

f:id:chayarokurokuro:20210719174232j:plain

ほぼほぼ右肩上がりの一直線ですかね。もしくは、右下の4点は夏場かな、冬場用と夏場用で直線を2本引くか。いずれにしても検査数を上げれば陽性者数も上がる。



おわりに

独立変数は今回は検査数の１つでしたので単回帰分析です。複数の独立変数から回帰分析するのは重回帰分析で、numpyだけでやるにはnumpy.linalg.lstsqを使うとできる(参照 : Numpyだけで重回帰分析 | 分析ノート )

重回帰分析はscikit-learnなら内蔵するデータを使ってボストンの住宅価格を予測するお決まりのアレです。



さて、コロナの陽性者数は検査数に依存します。陽性者数だけ見て一喜一憂しても仕方ないです。
検査数が毎日一定なら分かりますけどね。

100万人ぐらい人を集めて2,3日の期間限定で一気にPCR検査をやれば、陽性者数は爆上がりするはず。そして一気に検査数を元の水準に減らせば、10人増えただの100人増えただの言ってるマスコミの茶番にさすがに気づくんじゃないかな。茶番と言っても経済悪化させて人殺してるわけですからね。冗談じゃないよ。

目立ちたがりのYouTuberが面白がってやってくれないかな。 日本はそんな検査能力ないか。



以上。