overflow33の日記

python 機械学習 系の記事を書いて行きたい所存

引き続きCOVIDの状況分析

結論

・日本は、他の国と比べるとまだマシ。
・ヨーロッパでは、感染者数の増加ペースは飽和したが、依然油断ができない。
アメリカは地獄。トランプの頭はおかしい。
・感染が拡大してしまった国では、収束が見通せない。
・今後、ブラジル、ロシアが心配になってきた。

グラフの見方

見やすさのため縦軸は、対数。
横軸は、2020/01/22を Day 0 とした。
(Day80が2020/04/11、Day88が2020/04/19)

日本

f:id:overflow33:20200419212612p:plain
日本
ここ10日ほど、回復者数に対して死者数の増え方が急になっている。
どうやら、TVで報道されているように、医療崩壊が始まっているというのは事実だろう。
医療崩壊は、回復者数と死者数との差および、位相差(タイミングのずれ)によく現れる。
宣言の効果で、来週の増加ペースがどれだけ落ちるかが見どころだが、このグラフを見ると期待が持てない。

他の国と比べると死者数の増加ペースは遅く、政府、国民はよくやっている方ではないだろうか。
それでもなお、ウイルスの感染力が強力で太刀打ちできていない。
本格的な収束には、いよいよ強権の発動が必要かもしれない。

今更だが、緊急事態宣言があと10日早ければ、韓国と同じコースを辿れたのではないかと…。
やはり、拡大が始まってからでは遅すぎたのかもしれないが、結論を出すには時期尚早である。

アメリ

f:id:overflow33:20200419213002p:plain
アメリ
まさに世紀末。
この世の終わりかというくらいの蔓延期である。
このタイミングで、経済活動の再開を提示できるトランプ大統領サイコパスかと…。

中南米

f:id:overflow33:20200419212701p:plain
メキシコ
f:id:overflow33:20200419213250p:plain
ブラジル
ペースは鈍いが、拡大中。
暖かくなると感染力が弱くなるのは本当かも。
今後、ノーガードを決め込んだブラジルがどうなるか心配。

ヨーロッパ

f:id:overflow33:20200419212526p:plain
イタリア
f:id:overflow33:20200419212815p:plain
スペイン
f:id:overflow33:20200419212353p:plain
フランス
f:id:overflow33:20200419212400p:plain
ドイツ
ペースは鈍ってきたが、順当に感染拡大中。
収束は遠い。

f:id:overflow33:20200419212947p:plain
イギリス
f:id:overflow33:20200419213123p:plain
オランダ
医療崩壊とはこのこと。
現時点で、一日当たりで見ると、回復者<<死者数。
運び込まれたときには、瀕死か、もう無くなっているケースが多数あると思われる。

f:id:overflow33:20200419212742p:plain
ロシア
感染拡大中。
プーチンが檄を飛ばしたらしいが、効果の程は不明。
今後、アメリカ同様の地獄になるのでは…?。

イラン

f:id:overflow33:20200419212459p:plain
イラン
言われていたほど、悲惨なことにはなっていない。
アメリカよりはマシ。
拡大も阻止はできているが、なかなか感染者数が減らない。

インド

f:id:overflow33:20200419212432p:plain
インド
こちらも、着実に感染拡大中。
今後の爆発が心配される。

韓国

f:id:overflow33:20200419212623p:plain
韓国
理想的。

中国

f:id:overflow33:20200419212312p:plain
中国
死者数の追加報告が合った分、グラフの最後が跳ねている。
まあ、まだ隠されてるのかもしれないが…。


詳しい分析の手法は、前回のブログを参照。
前回と異なる点は以下。
前回:7日間移動合計
今回:7日間移動平均
値を一日あたりに直したかどうかです。
overflow33.hatenablog.com

データのソース

2019 Novel Coronavirus COVID-19 (2019-nCoV) Data Repository by Johns Hopkins CSSE
ArcGIS Dashboards
GitHub - CSSEGISandData/COVID-19: Novel Coronavirus (COVID-19) Cases, provided by JHU CSSE

COVID-19 感染から回復、死亡までの日数比較

概要と結論

日本では、COVID-19に罹ってから、回復、死亡するまでにどのくらい時間がかかっているのか、考えてみた。
各国の1週あたりの感染者数、回復者数、死亡者数の推移から、医療現場の状況を推察した。
分析では、感染者一人ひとりのデータを追うのは大変なので、日別の感染者数から、大体の目安を付けることにする。

・日本では、感染確認から回復確認までに約14日程度(10-20日)、感染確認から死亡確認まで26日程度(20-30日)
医療崩壊が起こると、1週あたりの回復者と死者の割合が 1:1 に近づき、回復者<死者となる場合がある

詳細

データ処理の内容

各日の確認数から 前4日、後3日の合計7日(=1週)の合計値を計算し、プロットした。
1日当たりだと、雑音が多く土日の影響も出てしまうため、移動合計を利用した。
縦軸の数値は、その日を基準にしたときの、1週あたりの人数を表している。

見やすさのため縦軸は、対数。
横軸は、2020/01/22を Day 0 とした。

日本の場合

f:id:overflow33:20200405205945p:plain

日本の場合、第一波と第二波が見られるのが特徴である。
第一波に注目する。
感染者数のピークが、約Day 10、収束が 約Day 17 である。
回復者数のピークは、約Day 25、収束が 約Day 34 である。
死者数のピークは、約Day 37、収束が 約Day 42 である。

ピークと収束の間隔から
感染確認から回復確認までの日数は、だいたい14-15日。
感染確認から死亡確認までの日数は、だいたい 25-27日。
となる。

第二波に関しては、まだピークも収束も迎えておらず、推定が難しい。

イタリア、スペイン、アメリカの場合

f:id:overflow33:20200405210022p:plainf:id:overflow33:20200405210042p:plainf:id:overflow33:20200405210057p:plain

控えめに言ってヤバい。
なにがヤバいって、1weekあたりの死者数と回復者数がほぼ同じ。
回復者と死者が同相で増えていく。
さらに、アメリカでは、回復者よりも死亡者の位相が早い(最近になって、追いついては来たが)。
病院に来て、医者が診たときには、重症もしくは瀕死となった方が多数いることが容易に想像できる。
完全に医療崩壊である。

オランダの場合

f:id:overflow33:20200405210457p:plain
これはもっとヤバい。
少し前まで感染確認後、平均10日程度で、かなりの割合で人が亡くなっている!
しかも、回復者は死亡者より少ない状況が続いている。
感染者数の増加ペースが鈍っていることだけが救いである。
(もはや検査すらできていない可能性もあるが、そうでないことを願う…。)

韓国、ドイツの場合

f:id:overflow33:20200405210012p:plainf:id:overflow33:20200405210113p:plain
感染者数>回復者数>死者数
となっており、医療が機能してる様子が伺える。
ドイツの感染者数は非常に多いが、落ち始めており、この山を越えれば大丈夫だろう。
あと1-2週が山に見える。
韓国は、非常に上手く制御されている。
感染者数も少なく、死者の増加も見られない。
各国が韓国を手本にする理由がよく分かる。

中国の場合

f:id:overflow33:20200405210203p:plain
死者が隠されていないか、心配になるグラフである。
感染拡大後、死者のピークが先に来ているが、回復者数も一気に増えていく。
感染者が重症化する前に医療が機能して、命を助けた結果であると信じたい。
この数字が本当なら、巨大な隔離施設を作った甲斐があったといえる。

データのソース

2019 Novel Coronavirus COVID-19 (2019-nCoV) Data Repository by Johns Hopkins CSSE
ArcGIS Dashboards
GitHub - CSSEGISandData/COVID-19: Novel Coronavirus (COVID-19) Cases, provided by JHU CSSE

解析コード

main.py

from os import path

import matplotlib.pyplot as plt

import loadData
import myLoadModule as mLM

from matplotlib import rcParams
rcParams['font.family'] = 'IPAexGothic'
rcParams['font.size'] = 20

def main():
    folder = "data"
    deathfileName = "200404_death.csv"
    confirmedfileName = "200404_confirmed.csv"
    recoveredfileName = "200404_recovered.csv"

    saveFolder = "fig_200404"

    deathfilePath = path.join(folder, deathfileName)
    confirmedfilePath = path.join(folder, confirmedfileName)
    recoveredfilePath = path.join(folder, recoveredfileName)

    deathData = loadData.data()
    deathData.load(deathfilePath)

    confirmedData = loadData.data()
    confirmedData.load(confirmedfilePath)

    recoveredData = loadData.data()
    recoveredData.load(recoveredfilePath)

    countryList = confirmedData.countryList
    DataList = [deathData, confirmedData, recoveredData]
    labelList = ['死亡者数', '確認数', '回復数']
    my_plot(countryList, DataList, labelList, saveFolder)


def my_plot(contryList, DataList, labelList, saveFolder=None):

    dpi = 96
    width = 1920*0.8
    height = 1080*0.8

    for contry in contryList:
        for Data in DataList:
            fig = plt.figure(1, figsize=(width / 100, height / 100), dpi=dpi)
            ax = fig.add_subplot(1, 1, 1)

            plot_y = Data.getCountryDataPerDay7(contry)
            ax.semilogy(plot_y)
            # ax.plot(plot_y)
            ax.grid(True)
            ax.set_title(contry)
            ax.set_ylabel('[7日間移動合計 人数]')
            ax.set_xlabel('[日]')
        ax.legend(labelList)

        plt.tight_layout()

        if saveFolder is None:
            plt.show()
        else:
            mLM.checkFolder(saveFolder)
            saveFilePath = path.join(saveFolder, contry.replace('*', '') + '.png')
            fig.savefig(saveFilePath)

        plt.close(fig)

if __name__ == "__main__":
    main()

loadData.py

from os import path

import pandas as pd

class data:
    def __init__(self):
        self.filePath = ""
        self.df = pd.DataFrame([])
        self.countryDf = pd.DataFrame([])
        self.countryDfDay = pd.DataFrame([])
        self.countryDfDay7 = pd.DataFrame([])
        self.countryList = []

    def load(self, filePath):
        self.filePath = filePath
        self.df = pd.read_csv(filePath)

        # print(self.df)
        self.make_countryData()
        self.makePerDayData()

    # 国別のデータにする
    def make_countryData(self):
        col = "Country/Region"
        self.countryList = sorted(list(set(self.df[col].values)))

        for country in self.countryList:
            idx = (self.df[col] == country)
            oneCountry = self.df[idx].sum()
            self.countryDf = pd.concat([self.countryDf, oneCountry], axis=1)
        # 緯度経度や場所の情報を取り除く
        self.countryDf.drop(index=['Lat', 'Long', 'Province/State', 'Country/Region'], inplace=True)
        # 列名を国名に付け替える
        self.countryDf.columns = self.countryList
        # print(self.countryDf)

    def makePerDayData(self):
        self.countryDfDay = self.countryDf.diff()
        self.countryDfDay7 = self.countryDfDay.rolling(7, center=True).sum()
        print(self.countryDfDay7)

    def getCountryDataPerDay(self, country):
        if country in self.countryList:
            return self.countryDfDay[country].values
        else:
            return None

    def getCountryDataPerDay7(self, country):
        if country in self.countryList:
            return self.countryDfDay7[country].values
        else:
            return None

    def getCountryData(self, country):
        if country in self.countryList:
            return self.countryDf[country].values
        else:
            return None

myLoadModule.py

from os import path
from os import mkdir

def checkFolder(dir):
    if path.exists(dir):
        pass
    else:
        mkdir(dir)

線形判別(LDA: Linear Discriminant Analysis) とは

線形判別(LDA: Linear Discriminant Analysis) とは

特徴量空間で、各クラスをガウス分布へ当てはめ、当てはめた分布を直線で分けるとき、どの方向に直線を引くのが良いか求めるものである。

フィッシャーの線形判別と検索すれば解説が出てくる。
例えば、以下。
フィッシャーの線形判別分析法 - Qiita

この説明では、引数(priors: 事前確率)の意味が理解できないため、以下では、sci-kit learn の解説ページをベースに説明する。**

sci-kit learn 公式の解説

1.2. Linear and Quadratic Discriminant Analysis — scikit-learn 0.22.2 documentation

ナイーブベイズ推定において、分布をガウス分布で近似し、境界を特徴量の線形結合の形(線形)に限定したものと解釈する。

問題設定

ある特徴量ベクトルXが与えられたとき、Xがクラスkに属する事後確率は、ベイズの定理から、次のように求められる。
P(y=k | X) = \frac{P(X | y=k) P(y=k)}{P(X)} = \frac{P(X | y=k) P(y = k)}{ \sum_{l} P(X | y=l) \cdot P(y=l)}
ここで、X: 特徴量ベクトル, y: Xの属するクラス, k: クラス である。

ここで、P(y= k | X)が最大になるようなkを求めれば良い。
(ここまでは、ナイーブベイズ推定と同じ)

各クラスの分布をガウス分布と仮定

各クラスの分布P(X | y)ガウス分布であると仮定する。
P(X | y=k) = \frac{1}{(2\pi)^{d/2} |\Sigma_k|^{1/2}}\exp\left(-\frac{1}{2} (X-\mu_k)^t \Sigma_k^{-1} (X-\mu_k)\right)
ここで、\mu_k: クラスkの平均ベクトル,  \Sigma_k: クラスkの共分散行列

各クラスの共分散行列が一致すると仮定

線形判別では、さらに、すべてのクラスにおいて、クラス内共分散行列(分散共分散行列)が一致する(\Sigma_k = \Sigma)と仮定すると、特徴量空間で線形の判別となる。

クラスの境界が線形となることの証明

線形となることは、2つのクラスk,lに対して、確率が同じになる境界が、特徴量の線形結合で表せることを示せばよい。

2つのクラスk,lに対して、確率が同じになるのは、次の場合である。
\frac{P(y=k | X)}{P(y=l | X)}=1

両辺に対数を取って、ベイズの定理、各クラスの分布P(X | y)ガウス分布であること、各クラスの共分散行列が同じであることを用いると、以下の式が得られる。
\log\left(\frac{P(y=k|X)}{P(y=l|X)}\right)=\log\left(\frac{P(X|y=k)P(y=k)}{P(X|y=l)P(y=l)}\right)=0 \Leftrightarrow
(\mu_k-\mu_l)^t\Sigma^{-1} X =\frac{1}{2} (\mu_k^t \Sigma^{-1} \mu_k - \mu_l^t \Sigma^{-1} \mu_l)- \log\frac{P(y=k)}{P(y=l)}

右辺(\frac{1}{2} (\mu_k^t \Sigma^{-1} \mu_k - \mu_l^t \Sigma^{-1} \mu_l)- \log\frac{P(y=k)}{P(y=l)})、および、Xの係数(\mu_k-\mu_l)^t\Sigma^{-1}は、定数であり、クラスk,lの境界は特徴量空間で線形である。
*P(y=k), P(y=l)は、各クラスがどの程度発生するかを表す事前確率である。

計算するもの

以上から線形判別では、各クラスに共通で当てはめができる共分散行列\Sigmaを求めれば良い。

ソルバー ’lsqr’ では、各クラスの共分散行列と平均の積\Sigma_l \mu_lと、共通の共分散行列と平均の積\Sigma \mu_lとの二乗誤差の総和が最も小さくなるものを利用する。

ソルバー 'eigen' では、基本的にクラス内分散をクラス間分散の比が大きくなるものを計算している。
これは、フィッシャーの線形判別として解説されているものである。

参考:
http://www.svcl.ucsd.edu/courses/ece271B-F09/handouts/Dimensionality3.pdf

’svd’では、'eigen'での計算の効率化のため、固有値を計算せずに、特異値分解(SVD: Singular Value Decomposition)を利用して、直接解を求めている(ようである 未確認)。


**
なお、事前確率(priors)を入力しない場合、学習データに含まれるデータの割合がそのまま事前確率となる。
学習データとテストデータでデータに含まれるラベルの割合が大きく異なることが予想できる場合、テストデータに合わせて事前確率を入力するべきである。

sci-kit Learn のiris(アヤメ)データを分類付きでCSVへ保存する

import pandas as pd
import numpy as np
from sklearn.datasets import load_iris

# データセットの読み込み
iris = load_iris()
# データセットに正解ラベルを追加する
# iris.data や iris.target は <numpy.ndarray> 型
# np.append は次元を揃える -> reshape(-1, 1) -> -1 は要素数の自動設定
data = np.append(iris.data, iris.target.reshape(-1, 1), axis=1)
# 特徴量の名前を取得(列ラベルにする) **
colLabels = iris.feature_names.copy()
# data に targetを追加したので、列ラベルにもtargetを追加する
colLabels.append("target")
# pandas DataFrameクラスへ格納
df = pd.DataFrame(data, columns=colLabels)

# ラベルの名前も DataFrame へ追加する
target_names = iris.target_names
# 追加する列を確保
df['target_name'] = ""
# 種類毎にtargetのインデックス(idx)を探し、DataFrameのidxを指定して追加する
for i, target_name in enumerate(target_names):
    idx = df['target'] == i
    df.loc[idx, 'target_name'] = str(target_name)

# データセットの確認
print(df.head(5))
# # CSV形式で保存
saveFilePath = 'iris.csv'
df.to_csv(saveFilePath, encoding="shift-jis")

iris(アヤメ)のデータセットを読み出し、ラベルと種類の名前を追加し、csv形式でファイルへ保存するスクリプトです。
機械学習もまずは、データを読み出して眺めるところから始めましょう。

読み込んだデータの整理には、pandas が非常に便利です。
今回は、csv形式で保存するために正解となるラベル(target)やラベル名(target_name)を列に加えて保存しました。
学習時には、間違って正解ラベルを学習する特徴量に加えないよう、注意してください。

深いコピー(Deep Copy)と浅いコピー(Shallow Copy)

**で示した部分の解説
colLabel = iris.feature_names.copy()
.copy() で深いコピーする。

.copyなしの浅いコピーでは、appendしたとき iris.feature_names の中身まで変わってしまう。
浅いコピーは、C言語ならポインタ渡しを行っているのと同じで、メモリのアドレスのみがコピーされ、メモリ領域を共有する。
深いコピー(.copy())は、実態の複製であり、ここでは iris.feature_names と同じ内容を持った新しいメモリ領域が colLabel に割り当てられる。

C Matlab python での 2次元配列 比較

2次元配列の書き方

C言語

int x[3][3] = {{0}, {0}};
int i, j;
for (i=0; i<3; i++){
    for(j=0; j<3; j++){
        x[i][j] = i + j;
    }
}
for (i=0; i<3; i++){
    for(j=0; j<3; j++){
        printf("%d ", x[i][j]);
    }
    printf("\n");
}

output

0 1 2
1 2 3
2 3 4

Matlab

x = zeros(3, 3)
for i = 1: 3
    for j = 1: 3
        x(i, j) = i + j
    end
end

python (list)

# x = [[0]*3]*3  # NG
x = [[0]*3 for i in range(3)]
for i in range(3):
    for j in range(3):
        x[i][j] = i + j
print(x)

output

[[0, 1, 2], [1, 2, 3], [2, 3, 4]]

python (numpy)

import numpy as np
x = np.zeros([3, 3])
for i in range(3):
    for j in range(3):
        x[i, j] = i + j
print(x)

output

[[0. 1. 2.]
 [1. 2. 3.]
 [2. 3. 4.]]

Cとmatlabは1次元配列と同じ要領で、2次元配列を作ることができる。
python では、リスト内包表記で初期化する必要がある。
pythonの(2次元)リストは、文字列なども格納できるため便利だが、数値のみを格納する場合は、素直にnumpyを使った方が簡単である。

複雑なテーブルを作成する場合

pythonでは、 pandas の DataFrame 利用することを推奨する。
Matlabでは、table が利用できる。( 2013b 以降 )


解説:
python の2次元 list を作ろうとする場合、初期化で

x = [[0]*3]*3  # NG

としたくなるが、この場合、外側の[~]*3で中の3つの要素(アドレス)がコピーされる。
(この場合、x = [[x1, x2, x3], [x1, x2, x3], [x1, x2, x3]] のような状態)
ここに、例えば"x1=5"と数値を入れると、3つのx1がすべて5になる。
(x = [[5, x2, x3], [5, x2, x3], [5, x2, x3]] )
上の例では、次の結果となる。

python (list)

x = [[0]*3]*3  # NG
for i in range(3):
    for j in range(3):
        x[i][j] = i + j
print(x)

output

[[2, 3, 4], [2, 3, 4], [2, 3, 4]]

参考:
Pythonのリスト(配列)を任意の値・要素数で初期化 | note.nkmk.me

C Matlab python での 配列 比較

1次元配列の書き方

C言語

int i;
int x[11]={0};
for(i=0; i<=10; i++;){
    x[i] = i;
}

Matlab

x = zeros(10)
for i = 1: 11
    x(i) = i - 1
end

python (list)

x = [0]*11
for i in range(11):
    x[i] = i

python (numpy)

import numpy as np
x = np.zeros(11)
for i in range(11):
    x[i] = i

xは < class 'list' > 型ではなく、 < class 'numpy.ndarray' > 型になる。


いずれも、xに 0 から 10 までの整数を配列として代入するプログラムである。
動作を高速化するため、先に配列分のメモリを確保する書き方になっている。


Matlabpythonでは、動作が遅くなるが、順次メモリを確保させる書き方ができる。
xの要素が一つ増えるたびにメモリ確保が行われ、動作が遅くなるので注意。

Matlab

for i = 1: 11
    x(i) = i - 1
end

python (list)

x = []
for i in range(11):
    x.append(i)

python (numpy)

import numpy as np
x = np.array([])
for i in range(11):
    x = np.append(x, i)

C言語でメモリを動的に確保するには、malloc.hを利用する必要がある。
pythonmatlabは楽でいいですね。

C Matlab python での if文 比較

if文(条件分岐)の書き方

C言語

int x,y;
x = 0;
if (x == 0){
    y = x + 1;
}else if(x == 1){
    y = x - 1;
}else{
    y = x;
}

Matlab

x = 0;
if x == 0
    y = x + 1;
else if x == 1
    y = x - 1;
else
    y = x;
end

python

x = 0
if x == 0:
    y = x + 1
elif x == 1:
    y = x - 1
else:
    y = x

C言語

条件式の括弧 ( ) と、分岐時に実行される処理に大括弧 { } が必要。

Matlab

条件分岐が終了したことを示すために end が必要。
if文の後ろの セミコロン(;) は不要。

python

if や elif, else の条件式の最後に コロン(:) が必要
else if を elif と書く。