2020-06-14

Issue Creation Meetup を活用して課題設定できるデータ分析者を目指す！

先日、参加したIssue Creation Meetupが有意義だったので内容を投稿します。
次回、参加を検討される方の参考になれば幸いです😃

はじめに

データサイエンティストに求められる能力の中に「課題設定力」があります。
プロのデータサイエンティストによるとこの能力はデータ分析プロジェクトの半分を占めるくらい重要です。

※以下の動画3:35あたりを参照。（問題設計≒課題設定と置き換えています）

www.youtube.com

Issue Creation Meetupはまさにこの課題設定力を鍛えるのための勉強会です。
データ分析系の勉強会として新しい試みではないでしょうか？

connpass.com

勉強会の内容

以下のお題について、チームで知見を出し、発表するというものでした。

お題：「あなたはゲームアプリ運営会社の社員です。
　　　　　あるサービスの離脱率を下げる施策を考えよ。」

参加者：5人/チーム×2チーム
役割：チーム内でファシリテーター、発表者、書記を決める
制限時間：1時間

与えられたデータは以下です。

ファイル	カラム	行数
users_sample.csv	インストール日、ユーザーID、性別、年齢、デバイスタイプ、アドネットワーク	4872
dau_sample.csv	ログ（日付け）、ユーザーID	20066
action_sample.csv	ログ（日付け）、ログ（時刻）、ユーザーID、デバイスタイプ、アクション名、アクションID	88865

主宰者が実務で使ったものを、企業秘密に触れない範囲で加工したものです。
なので、ビジネスで起こる現象をある程度反映しているデータになっています。

主宰者より、以下の基本方針が示されました。

■現状の確認
　・KPIの基礎集計 - インストール数／DAU／継続率など
■問題設計
　・仮説出し、検証
　・目的変数、説明変数の定義
　・結果を踏まえての施策案

与えられた情報は以上です。

1時間後に何かしらの施策を発表しなければいけません。
kaggleなどのデータ分析コンペでは、きっちり設定された問題に対して、モデリングや評価方法のスキルを競います。それに対して、Issue Creation Meetupでは、ビジネスでありがちなふわっとした課題に対して課題の設定と施策の決定が求められます。

CRISP-DMにおける棲み分けをざっくり書くと以下のようになるかと。

f:id:banquet-kuma:20200611175804p:plain — 求められる能力の違い

どちらの能力がより重要という訳ではなく、バランスよく鍛えていくことが大切だと思います。

赤字で書いた内容を1時間で行うのはかなり大変です😅

チームでの活動内容

自己紹介もそこそこに、どういう方針を取るか協議しました。

チーム内にゲーム業界のドメイン知識がある人がいなかったので、
データを片っ端から可視化してて、それを基に皆で仮説を出し合うことにしました。
本来、データを可視化する前に仮説を立てて、その仮説を立証するための可視化を行うのがベターと思いますが、とにかく手を動かそうと言うことで意見が一致しました。

以下の時系列データを視化しました。

・離脱率
・アクティブユーザー数
・アクション別のアクション数
・デバイス別のインストール数
・アドネットワーク別の流入数

さて、お気付きでしょうか？
離脱率ってどう定義するのでしょうか？結構難しいです。
協議していると時間が経ってしまうことと、与えられたデータから簡単に算出できなさそうだったので、1日のアクティブユーザー数の推移で代替することにしました。

つまり、「1日のアクティブユーザー数が減っている　≒　離脱率が上昇している」と見なしました。アクティブユーザー数の推移を可視化すると以下のようになりました。

f:id:banquet-kuma:20200612093032p:plain — アクティブユーザー数の推移

8月に盛り上がり（夏休み？）を見せるものの、8月中旬から下がっています。
このことから、昨今の離脱率の上昇は本当であると結論付けました。

なお、お題の文言から、離脱率を下げる施策を考えていましたが、
同じようなヒジネス価値を生むのであれば、違う手段を取ってもOKでした。
（後で気付いた）

例えば、「1日のアクティブユーザー数を増やすにはどうすれば良いか？」とういうふうに思考を転換しても良いです。
この点については、主宰者のyokkunsさん（DATUM STUDIO株式会社の取締役CAO）より懇親会で解説がありました。

次に我々のチームで着目したのが、アクション数の推移です。
アクションには次の3つがあります。

・message
　ユーザー同士がメッセージを送る行為（チャット機能みたいなもの）

・asking_help
　ユーザーが他のユーザーに救援を要請する行為（協力して敵に立ち向かう）

・battle
　ユーザーが対戦ゲームを行う行為（このゲームの本筋な楽しみ方）

それぞれのアクション数を可視化すると以下のようになりました。

f:id:banquet-kuma:20200612214504p:plain — 各アクション数の推移

グラフを見て分かることは、battleがmessageとasking_helpに比べて少ないということです。battleだけ抜き出すと以下のようになります。

f:id:banquet-kuma:20200612214722p:plain — battle数の推移

8月中旬のピーク以降減っています。
このことから、入会したユーザーの（ゲーム本来の楽しみ方である）battleへの移行率が低いことが離脱率上昇の要因と考えました。

それでは、なぜbattleへの移行率が低いのでしょうか？

もう一度、各アクション数のグラフ見て意見を出し合いました。

・messageとbattleは8月中旬以降減少
・asking_helpは8月中旬以降も減っていない
・messageとbattleは同じような位置にピークが立つ
・messageとasking_helpのピークに遅れてbattleのピークが立つように見える
　（①、②の位置関係）

f:id:banquet-kuma:20200612220724p:plain — message , asking_help , battleのピーク位置の関係,

これらの事実から、ユーザーの行動を以下のように推測しました。

f:id:banquet-kuma:20200612224814p:plain — ユーザーの平均的な行動

ただし、message & asking_helpとbattleの間に結構なタイムラグがあります。
このことから、battleについては、「〇回以上行ったユーザーのログが残る」などの仮定が必要になりそうです。（この仮定に一番無理があるかも😅　）

何故、battleが継続しないか？
battleのピーク（①、②）の前後でasking_helpの数は減っておらず、
救援要請が続いていることから、ゲームの難易度が難しすぎるのでは無いかと考えました。つまり、、

f:id:banquet-kuma:20200613204512p:plain — インストール → 離脱までの流れ

この仮説を検証する施策として以下を提案しました。

「ゲームの難易度別にA/Bテストを行い、離脱率の変化を見る」

意見をここまで纏めたところでタイムアップとなりました😂.

懇親会での裏話

発表の後はRemo上での懇親会に突入しました。
一線級のデータサイエンティストに質問し放題というおいしい懇親会です😆

■発表に対してのフィードバック
　・全体を見てから、各論に議論を移すのは良かった。
　　ピークの日にどんな出来事があったなど、詳細にデータを見ていくのは次のステップ
　　でOK。
　・実は離脱の原因を特定するのは難しい。何の前触れもなく、突然消失するユーザーが
　　一定数いる。
　　だから発想を変えて、新規入会者を増やすという手を考えるのもあり。
　　（ビジネス価値は変わらない）

■質疑応答
Q1.何故この勉強会を企画したかのか？
A1.データ分析できる様々なツールやkaggleのおかげで、"How to" は浸透した。
　しかし、課題設定できる分析者が圧倒的に足りていない。
　課題設定は文脈を読まないとできないタスクであり、今のところ機械化できない。
　逆に言うと、それができないと将来的に淘汰される。

　この話を私のイメージでまとめると以下のようになります。

f:id:banquet-kuma:20200613214654p:plain — どこを目指すか

Q2.こんな未経験からの参入者は嫌だ
A2.手法に拘る人
　データをよく見る前に、機械学習の〜手法にデータを突っ込もうとする人。
　こういう人は矯正するのに時間がかかる。

Q3.今回、時間制限が厳しかったが、1時程度で知見を求められることはあるのか？
A3.クライアントからというより、上司から「あのデータどうなった？何か言えそう？」
　みたいな無茶振りが飛んでくることはある（笑）

所感

短い時間で初対面の人と知見を出さないといけないため、参加する前はかなりハードルが
高いと感じていました。しかし、一つの目標に向かって一致団結すると意外となんとかなるものです。

時間制限の観点で大事と思ったことは以下です。

・とにかく思ったことを発言する（言うか言わないか迷っている時間がもったいない）。
・一方で、人が発言している間は静かに聞く。
・良いと思った意見はとにかく肯定する。
　「良いね！やってみよう！👍」みたいな感じに。
・短時間で可視化するためにExcelを使い倒せた方が良い。
・時間配分を意識して、余裕を持って意見をまとめる。
　「残り〇分になりました。ここまでの意見を要約すると〜」みたいな促しをする。

課題設定力はすぐに身に付くものではないので、こういった機会を活用して継続的に鍛錬していきたいです。
課題設定力を伸ばしたい方にはオススメの勉強会です！
参加者の皆様、特に主宰いただいたTakashiMinodaさん、yokkunsさんありがとうございました！😃

2020-05-19

殺人事件はランダムに起こっていると言えるか？

はじめに

統計を勉強するにあたり、
教科書の演習問題を解くだけではモチベーションが上がらないので、実際のデータに当てはめて考察してみました。
テーマは「指数分布、中心極限定理、ガンマ分布」です😀

データセットの入手

指数分布とガンマ分布は、いずれもランダムに起きる事象を扱う確率分布です。
なので、ランダムに起きていそうな現象を探しました。

Google Cloud 一般公開データセットを見ていたところ、
シカゴの犯罪に関するデータがあったので、それを使ってみることにしました。
犯罪のカテゴリーとしては、ランダムに起こりそうなものを選びたかったので、「HOMICIDE（殺人）」を選びました。

以下のリンクからGoogle Couldにログインします。

cloud.google.com

（データセットのソース　↓

City of Chicago | Data Portal | City of Chicago | Data Portal）

「chicago」で検索すると、以下のデータセットが表示されます。
「Chicago Crime Data」をクリックします（2020/5/21現在の表示）

f:id:banquet-kuma:20200521191917p:plain — Chicago Crime Data

次に「データセットを表示」をクリックします。

f:id:banquet-kuma:20200521192350p:plain — Chicago Crime Data

画面上部の空白にSQLのクエリを書いて「実行」を押します。
すると、所望のデータだけを抽出することができます。

f:id:banquet-kuma:20200521192827p:plain — 「HOMICIDE」を抜き出すクエリ

今回は、以下のクエリにより
2001年1月1日～最終更新（2020年5月9日15時過ぎ）までに起きた殺人事件の発生時刻を抜き出しました。トータルで10223件です。

SELECT 
    date
FROM 
    `bigquery-public-data.chicago_crime.crime` 
WHERE 
    primary_type = "HOMICIDE"
ORDER BY
    date

抽出し終わったら、「クエリを保存」を押し、csvに保存します。

f:id:banquet-kuma:20200521200045p:plain — csvで保存

ちなみに、上記の操作にはGoogle Cloud Platformというシステムを使っています。
有料ですが、12 か月間 $300 分の無料トライアルが用意されているので、この程度の操作であれば無料で使うことができます😃

指数分布

単位時間にλ回起きるランダムな現象の発生間隔を表します。

$f(x) =λe^{-λx}　(x：時間)・・・(1)$

殺人事件がランダムに起こっているのであれば、その発生間隔が指数分布に従うはずです。先ほどのcsvを使って、殺人事件の発生間隔を算出してプロットします。

f:id:banquet-kuma:20200522145223p:plain

発生間隔のプロット

Pythonを使って行います。

import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.stats import expon
from scipy.stats import gamma
import japanize_matplotlib
import seaborn as sns
sns.set(font='MS Gothic')

# データの読み込み
data = pd.read_csv("データセットから抽出したcsvのパス")

# date列から"UTC"を削除する
UTC_del = lambda x: x.replace(" UTC","")
data["date"] = data["date"].map(UTC_del)

#date列をdatetime型に変換する
data["date"] = pd.to_datetime(data["date"])

# date列でソートする（不要かも）
data = data.sort_values('date')

# date列の1階差分を取る
data_diff = data.copy()
data_diff["delta time"] = data_diff["date"].diff()

# 0行目は欠損するので除く
data_diff = data_diff.iloc[1:,:]

# delta time列を分に変換する
trans_min = lambda x: x.total_seconds()/60
data_diff["delta time (min)"] = data_diff["delta time"].map(trans_min)

# 発生間隔の最短、最長、平均を求める
print("殺人事件が起こらなかった最短時間：{}分".format(round(data_diff["delta time (min)"].min(),1)))
print("殺人事件が起こらなかった最長時間：{}分".format(round(data_diff["delta time (min)"].max(),1)))
print("殺人事件が起こらなかった平均時間：{}分".format(round(data_diff["delta time (min)"].mean(),1)))

殺人事件が起こらなかった最短時間：0.0分
殺人事件が起こらなかった最長時間：12631.0分
殺人事件が起こらなかった平均時間：995.7分

平均すると16.6時間に1件の殺人が発生しているようです。

# 殺人件数が発生する平均的な時間間隔
inv_lambda = data_diff["delta time (min)"].mean()

# 指数分布を発生させるための時間軸
x =  np.arange(0, 6000, 1)

# ヒストグラムを指数分布の理論曲線を描く
plt.figure(figsize=(20, 10))

# 実データのヒストグラム
plt.hist(data_diff["delta time (min)"], 
         bins=int(round(np.log2(len(data_diff))+1)), 
         density=True,label="殺人事件の発生間隔の頻度")

# 指数分布の理論曲線（λ = 1分当たりの殺人事件発生回数の頻度）
plt.plot(x,expon.pdf(x,scale=inv_lambda), 
         color="red",label="λ＝[1分あたりの発生回数]とした指数分布")

plt.xlabel("発生間隔(分)",fontsize=18)
plt.ylabel("確率密度",fontsize=18)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.xlim(data_diff["delta time (min)"].min(),6000)
plt.legend(fontsize=18,loc="best")
plt.grid(b=True, which='major', axis='both')

f:id:banquet-kuma:20200522164637p:plain — 発生間隔の頻度と指数分布の比較

指数分布ぽくなっていますね。

ヒストグラムを描く際の階級数には注意が必要です。
階級数の取り方で、ヒストグラムの見え方が如何様にも変わるからです💦
今回は以下のスタージェスの公式を使いました。

$k≒1+log_{2}N　(k：階級数、N：観測数)$

さて、上のグラフの理論曲線（赤線）を見て何か疑問に思われないでしょう？
「なぜ、発生間隔0（分）で確率密度が最大になるのか？」と。
これは、ある殺人事件が起こった次の瞬間に殺人事件が起こる確率が最も高いという
ことです。

感覚的には、殺人事件が1件が起こったら、しばらく平和な時間が続きそうじゃないですか？😅（というか、続いて欲しい）
同じような疑問を感じられた方がいましたら、是非以下の動画をご覧ください。
まさにその点が解説されています。

www.youtube.com

しっくりくる説明はたくみ先生にお任せするとして、
数式の観点からこの疑問を解決します。

ある殺人事件が起こってから、ｓの間殺人事件が発生しなかったという前提を置きます。その時にさらにｔ経過しても殺人事件が発生しない確率を考えます。
この確率は以下のように書けます。

$\displaystyle{P(X≧s+t|X≧s)}$

$\displaystyle{=\frac{P(X≧s+t,X≧s)}{P(X≧s)}}$

$\displaystyle={\frac{P(X≧s+t)}{P(X≧s)}}$

$\displaystyle{=\frac{e^ {-λ(s+t)}}{e^ {-λs}}}$

$=e^ {-λt}$

$=P(X≧t)$

ある殺人事件が発生してから、tの間殺人事件が発生しない確率と等しくなりました。
つまり、ランダムに起こる事象に対しては「sの間事象が発生しなかったから、次のtの間も発生しにくいとか、発生しやすいということは全く言えない」ということです。
これを「無記憶性」と言います。この性質があるので、ある殺人事件が起こった次の瞬間に別の殺人事件が起こる可能性が最も高いわけです。

f:id:banquet-kuma:20200522223854p:plain — 無記憶性の説明

統計学的には、
「天災は~~忘れたころに~~忘れる前にやってくる」と言うのが正しそうです。

分布関数を求めてみます。

plt.plot(x, expon.cdf(x,scale=inv_lambda)*100, color="red"
                ,alpha=0.5, label="殺人事件の発生確率".format(round(inv_lambda,1)))

plt.xlabel("発生間隔(分)",fontsize=24)
plt.ylabel("分布関数×100",fontsize=24)
plt.xticks(fontsize=22)
plt.yticks(fontsize=22)
plt.xlim(data_diff["delta time (min)"].min(),6000)
plt.legend(fontsize=28,loc="best")
plt.grid(b=True, which='major', axis='both')

f:id:banquet-kuma:20200523171402p:plain — 分布関数

print("1時間以内に殺人が起きる確率：{} %".format(round(expon.cdf(60,scale=inv_lambda)*100,1)))
print("24時間以内に殺人が起きる確率：{} %".format(round(expon.cdf(60*24,scale=inv_lambda)*100,1)))

1時間以内に殺人が起きる確率：5.8 %
24時間以内に殺人が起きる確率：76.5 %

中心極限定理

せっかくなので、中心極限定理が成り立っているか確認します。

中心極限定理とは、

母集団がどんな分布でも、サンプルサイズが十分大きければ、
標本平均の確率分布が正規分布に従うと見なしてよいという定理です。

式で書くと以下のようになります。

$\displaystyle{標本平均　\frac{(X_{1}+X_{2}+・・・X_{n})}{n}　～　N(μ , \frac{σ^2}{n})}・・・(2)$

$X：殺人事件の発生間隔、n：サンプルサイズ、μ：Xの母平均、σ^2：Xの母分散$

殺人事件の発生間隔の母平均、母分散は不明なので、
不偏推定量である標本平均と不偏分散を使います。

標本数を100に固定して、
サンプルサイズn = 1、10、100、500についてヒストグラムと(2)式から求めた正規分布を
書くと以下のようになります。

# 殺人事件の間隔をランダムサンプリングして標本平均を求める関数
# sample_size：サンプルサイズ、number_of_samples：標本数

def cal_mean(sample_size , number_of_samples):
    # 標本平均を入れるリストを作成
    mu_list = []
    for i in range(number_of_samples):
        # indexをランダムに生成
        random_index = random.sample(list(data.index), k=sample_size)
        # 平均を算出
        mu = data_diff["delta time (min)"][random_index].mean()
        # リストに格納する
        mu_list.append(mu)
    
    return mu_list

# 標本数は100に固定する
number_of_samples = 100

# サンプルサイズ1の時の標本平均値
n1 = cal_mean(1 , number_of_samples)

# サンプルサイズ10の標本平均値
n10 = cal_mean(10 , number_of_samples)

# サンプルサイズ100の標本平均値
n100 = cal_mean(100 , number_of_samples)

# サンプルサイズ500の標本平均値
n500 = cal_mean(500 , number_of_samples)

# 時間軸
X = np.arange(0,3000,0.1)

#標本平均（＝母平均の不偏推定量）
mu = data_diff["delta time (min)"].mean()

# 不偏標準偏差（＝母標準偏差の不偏推定量）
std = data_diff["delta time (min)"].std()

# サンプルサイズ1の時の描画
Y1 = norm.pdf(X,mu,std/np.sqrt(1))
plt.plot(X,Y1,color='#123456',linestyle="-",label="N(μ,$σ^{2}$/1)")

plt.hist(n1, 
         bins=int(round(np.log2(len(n1))+1)), 
         color='#123456',
         density=True,label="サンプルサイズ = 1",
         alpha=0.3)

# サンプルサイズ10の時の描画
Y10 = norm.pdf(X,mu,std/np.sqrt(10))
plt.plot(X,Y10,color='#4169e1',linestyle="-",label="N(μ,$σ^{2}$/10)")

plt.hist(n10, 
         bins=int(round(np.log2(len(n10))+1)), 
         color='#4169e1',
         density=True,label="サンプルサイズ = 10",
         alpha=0.3)

# サンプルサイズ100の時の描画
Y100 = norm.pdf(X,mu,std/np.sqrt(100))
plt.plot(X,Y100,color='#1300e6',linestyle="-",label="N(μ,$σ^{2}$/100)")

plt.hist(n100, 
         bins=int(round(np.log2(len(n10))+1)), 
         color='#1300e6',
         density=True,label="サンプルサイズ = 100",
         alpha=0.3)

# サンプルサイズ500の時の描画
Y500 = norm.pdf(X,mu,std/np.sqrt(500))
plt.plot(X,Y500,color='#e6bf00',linestyle="-",label="N(μ,$σ^{2}$/500)")

plt.hist(n500, 
         bins=int(round(np.log2(len(n10))+1)), 
         color='#e6bf00',
         density=True,label="サンプルサイズ = 500",
         alpha=0.3)


plt.xlabel("発生間隔の平均値(分)",fontsize=24)
plt.ylabel("確率密度",fontsize=24)
plt.xticks(fontsize=22)
plt.yticks(fontsize=22)
plt.xlim(data_diff["delta time (min)"].min(),3000)
plt.legend(fontsize=28,loc="best")
plt.grid(b=True, which='major', axis='both')

f:id:banquet-kuma:20200523185825p:plain — 中心極限定理

確かに、サンプルサイズが大きくなるにつれてヒストグラムが(2)式の正規分布に近づいていることが分かります。サンプルサイズ10ですでに峰が見え始めています。そして、サンプルサイズが増加するに従って、全データを使って算出した標本平均995.7分の辺りにデータが集中してくることが分かります。

最後に10件の殺人事件が起こるのに要する時間がガンマ分布にあてはまるかどうか見てみます。

ガンマ分布

ガンマ分布は平均的に期間 $β$ ごとに1回起きるランダムな事象がα回起こるまでに要する
時間の分布を表します。

$\displaystyle{Ga(α , β)=}$

$\displaystyle{\frac{1}{Γ(α)}\frac{1}{β}(\frac{x}{β})^{α-1}exp(-x/β) ,x>0}$

$\displaystyle{Γ(α)}$ は以下で定義される。

$\displaystyle{Γ(α)=\int_{0}^{∞} y^{\alpha-1}\exp(-y)dy}$

ここで、 $α$ =10、 $β$ =995.7分（発生間隔の平均）です。

# dateの10階差分を取る
data_diff10 = data.copy()
data_diff10["delta time"] = data_diff10["date"].diff(periods=10)

# 10行目までは欠損するので削除する
data_diff10 = data_diff10.iloc[10:,:]

# delta time列をトータル分に変換する
trans_min = lambda x: x.total_seconds()/60
data_diff10["delta time (min)"] = data_diff10["delta time"].map(trans_min)

# 時間軸
x =  np.arange(0, 40000, 1)

plt.figure(figsize=(20, 10))

# 殺人事件が10件発生するのに要する時間
plt.hist(data_diff10["delta time (min)"], 
         bins=int(round(np.log2(len(data_diff10))+1)), 
         density=True,label="殺人事件10件の発生間隔の頻度")

# ガンマ分布の理論曲線
plt.plot(x, gamma.pdf(x,a=10,scale=inv_lambda),
         color="red",label="Ga(α,β)")

plt.xlabel("殺人事件10件の発生間隔(分)",fontsize=24)
plt.ylabel("確率密度",fontsize=24)
plt.xticks(fontsize=22)
plt.yticks(fontsize=22)
plt.xlim(data_diff10["delta time (min)"].min(),30000)
plt.legend(fontsize=28,loc="best")
plt.grid(b=True, which='major', axis='both')

f:id:banquet-kuma:20200523232803p:plain — ガンマ分布とヒストグラムの比較

うーん、結構ガンマ分布から外れています。
（ということは、指数分布からもそれなりに外れていたということか。）

おそらく、 $β$ = 一定としてしまっていることに無理があるのではないかと思います🤔

$β$ は2001年1月1日～ 2020年5月9日の間の平均的な殺人事件の発生間隔であり、
ランダムにこの頻度で殺人事件が発生していれば、ヒストグラムと $\displaystyle{Ga(α , β)}$ は一致するはずです。

しかし、実際には上記の期間には様々なことが起こっています。
9.11テロ、リーマンショック、大統領の交代　etc.

殺人事件の発生に影響を及ぼす経済状況の変化や法律の改正があったかもしれません。
したがって、完全にランダムに殺人事件が発生していたとは言えないように思います。
つまり、殺人事件が起きやすい期間と起きにくい期間があったのではないかと。

数式で表すと、

$\displaystyle{Ga(α , β(t))=}$

$\displaystyle{\frac{1}{Γ(α)}\frac{1}{β(t)}(\frac{x}{β(t)})^{α-1}exp(-x/β(t)) ,x>0}$

$β$ が刻々と変わる時間の関数なのではないかと。

この $β(t)$ の出し方は勉強が進んだら検討してみたいと思います。
（そもそも数学的に関数を導出できるのか否かすら分かりません😅）

そもそも、ランダムであることって数学的に厳密に定義するとどうなるんだろう🙄？

最後に

統計を勉強をするにあたっては、世の中の現象がどんな確率分布に従っているか見てみる
と面白いです。教科書どおりの数式にピッタリあてはまることは少ないかもしれせん。でも、逆にそれが勉強意欲をそそってくれますｗ
オープンデータが豊富にあり、PythonやRでサクッと可視化できる時代になったのでこれらを勉強に活用しない手はないですね🤗

参考図書

現代数理統計学の基礎 (共立講座数学の魅力)

現代数理統計学の基礎 (共立講座数学の魅力) | 達也, 久保川 |本 | 通販 | Amazon

2020-05-16

エリア51への行き方　詳細版

ミリタリー好きな人のために、エリア51への行き方を記載します。
メジャーな観光地を訪れるだけでは満足できない人におススメのスポットです。
ラスベガスに行かれた際は是非お立ち寄りを😎

エリア51

正式名称グルーム・レイク空軍基地。
アメリカ空軍がステレス機など最新鋭の兵器開発を行っている基地です。2013年まで公式に存在を認めていなかったトップシークレットな基地です。ラスベガスから北西に200kmほどの砂漠地帯にあり、グルーム湖という塩湖の上にあります。

F-117 ナイトホークやSR-71 ブラックバードが開発されたことで有名な基地です🤩

f:id:banquet-kuma:20200510211533p:plain — F-117ナイトホーク

f:id:banquet-kuma:20200510210505p:plain — SR-71 ブラックバード

ja.wikipedia.org

また、ロズウェル事件で墜落したとされるUFOから発見された宇宙人が解剖されたと噂される基地です。かなり怪しい話ですが（笑）詳しく知りたい方は以下をお読みください。

f:id:banquet-kuma:20200510204930p:plain — エリア51 世界でもっとも有名な秘密基地の真実

www.amazon.co.jp

エリア51への行き方

まずラスベガスでレンタカーを借りましょう！（というか、一般人はその方法しかありません。マッカラン国際空港から通勤用の便が出ているらしいですが、当然関係者以外は乗れません）。

f:id:banquet-kuma:20200510205935p:plain

エリア51をメインディッシュとして、ラスベガスから右回りに、95号線→6号線→375号線経由で行きます。このルートを一周するとだいたい800km！東京～広島くらいの距離です😀

f:id:banquet-kuma:20200513010147p:plain — 右回りルート（Las Vegas →95号→6号→375号→93号→Las Vegas）

途中、デス・バレー国立公園に寄りましたが、まさに「死の谷」という感じで日本で見らない景色が見られて非常におススメです。「猿の惑星」や「スター・ウォーズ」のロケ地として有名で、世界最高気温56.7°Cが観測された場所でもあります。脱水症状に注意しましょう💦

tg.tripadvisor.jp

95号線沿いにはゴールドラッシュ時代に建設された町で、ゴーストタウンになっているスポットがあり、見て周るのも面白いです。

f:id:banquet-kuma:20200516211331p:plain — ゴーストタウン

さて、話は6号線から375号線へ分岐する地点へ一気に飛びます（笑）。
ここからが本題。

f:id:banquet-kuma:20200514224341p:plain — 6号線から375号線への分岐点

ネバダ州道375号線は全米一交通量の少ない州道であり、Extraterrestrial Highway（地球外高速道路）と呼ばれています。よくUFO（未確認飛行物体）が目撃されるんでそう呼ばれるみたいです🛸

375号線をエリア51の最寄り町Rachelを目指してひた走ります。
気持ち良いくらいの1本道、そして全く車とすれ違わない。。
200km/h出しても捕まりません。
ガソリンスタドが100km先までないような地域がざらにあるので、ガソリンの残量に注意しましょう！こんなところでガス欠したらマジ洒落にならない🚗

f:id:banquet-kuma:20200514224724p:plain — 道路で寝れます

途中、LOW FLYING AIR CRAFT（低空飛行の飛行機に注意）の標識が！！
注意しようがないですがｗ
軍用機のテスト飛行でもやっているんでしょうか！？
たぶん、通行人が目撃しているUFOて最新鋭のステレス機か何かなんじゃないかな🤔

f:id:banquet-kuma:20200514225139p:plain — 何が飛んでくる？

シールが貼られまくったExtraterrestrial Highwayの標識が見えてきたらRachelはすぐそこです。私も貼ろうとしたんですが、全く届かなかったです😅

f:id:banquet-kuma:20200515233624p:plain — Extraterrestrial Highwayの標識

Rachelには「Little Alien/リトル・エイリアン」というカフェがあり、世界中からUFOマニアが集まってきます。食事が取れ、宿泊できます。

f:id:banquet-kuma:20200515234107p:plain — Little Alien

f:id:banquet-kuma:20200516213539p:plain — お出迎え

コーラでも飲んで一服しながら、エリア51までの道を確認しましょう。
店でエリア51までの詳細な地図をもらえます😆（確か無料だったはず）。

f:id:banquet-kuma:20200515233408p:plain — エリア51までの詳細な地図

赤い点線の左側は立ち入り禁止区域なので、絶対に超えないようにしましょう👽
明日は「エリア51に最も近づける地点（point21）」まで行きます。

Little Alienの周囲には何もなく、星空がめちゃくちゃ綺麗なので、宿泊をおススメします✨

f:id:banquet-kuma:20200515234907p:plain — Rachelの夕暮れ

翌朝は今まで聞いたことがないくらいの爆音で目が覚めました🤣
実はRachelの近くには、エリア51以外にネリス空軍演習場とネバダ核実験場があります。
おそらく、戦闘機から出るソニックブームか爆撃訓練の音だと思います✈✈

f:id:banquet-kuma:20200515002104p:plain — Rachel周辺の軍事サイト

Little Alienの脇には放射線量計があり、現在は人体に問題ないレベルですが、
核実験が盛んだった冷戦時代は放射線量が高かったようです💦

いよいよ、エリア51に接近します！
昨日もらった地図を見ると2とおりの行き方があることが分かります。

Rachelから375号線を南下して、

(a) Black Mailbox（point7)を右折する
(b) グルームレイクロードへの分岐（point10）を右折する

(b)の道をおススメします。私は(a)の道を選んで迷いました💦
Black Mailboxは誰が何のために設置したか全く不明なMailboxです。
投函すると宇宙人と文通できるとかできないとか✉

f:id:banquet-kuma:20200516000302p:plain — 白いBlack Mailbox

f:id:banquet-kuma:20200516001206p:plain — point10の分岐

グルームレイクロードへの分岐を右折したら、エリア51に最も近づける地点（point21）
まで延々直進します。これが結構長い、、んですが、スピードを出し過ぎてうっかり制限区域に入ってしまうとヤバいので、安全運転で行きましょう！！

f:id:banquet-kuma:20200516001645p:plain — エリア51に続く道

見てのとおり、グルームレイクロードは未舗装なので、オフロードタイプのレンタカーを借りることをおススメします🚜

運が良ければ、周辺の町とエリア51を行き来する通勤バス（White Bus）に遭遇します。
窓が黒塗りなのが不気味ですね。。通勤者の顔を晒さないためでしょうか？
町の中にあるバス停の場所も極秘みたいです。こういう場所で働く人って、周囲に仕事の話を一切できないんでしょうね😅

f:id:banquet-kuma:20200516110006p:plain — White Bus

15kmほど走るとpoint21に到着します。ここがエリア51に最も近づける地点です。
左右にWARNINGの標識があるのが目印です。ゲートがないので、うっかり超えてしまいそうですが、絶対に超えてはいけません。私は他の基地で誤って基地内をさまよい、憲兵に拘束されたことがあります（前科はついていません）🤣。ここは、他のどの基地よりもヤバい事態になること間違いなしです。

f:id:banquet-kuma:20200516104628p:plain

ここを超えるな！写真を撮るな！
　・$1000の罰金 or 禁錮6か月　、もしくは両方
　・殺傷能力のある武器を使用する可能性があるよ

と書かれています😱

一見、無防備な境界なのですが、写真向かって右側の丘の上からがっつり監視されています。

f:id:banquet-kuma:20200516112952p:plain — 監視

よく見ると物体検知器のようなものもあります。

f:id:banquet-kuma:20200516113330p:plain — 物体検知？

この場所には世界中から軍事マニアが集まってきて、
なかば観光地化（？）しているので、ここからの写真撮影は黙認されているようです。

この先で世界最高の頭脳が、
人類がまだ見ぬ技術の開発を行っていると思うとワクワクします🤩
が、凡人は思いを馳せるだけですw

最後に

point21までいくツアーがあるみたいです笑
こういう場所は自分で行ったほうが楽しいと思うんだけどな～

www.veltra.com

2020-04-28

DLG Crossでデータマネジメントの学び方のヒントを得る！

はじめに

先日、「製造現場におけるデータマネジメントを考える」というタイトルでブログを書きましたが、データマネジメントの学び方のヒントが得られそうな勉強会があったので参加してみました。

参加したのは、データラーニングギルド*1が主催するデータ分析人材の交流イベントです。今回のテーマは今話題のデータマネジメントに関してでした。

最近話題になった以下のツイートにある職業だとデータエンジニアがデータマネジメントに直結するのですが、AIスペシャリストであっても、データサイエンティストであっても、データマネジメントされている前提で活躍できると思う訳です。

今後、データ分析を活用する企業が増えていくと同時にデータマネジメントできる人材の価値も更に上がっていくと思われます。

LinkedInの発表している「需要の高まる職業トップ15」の2020年版はこんな感じのようですが、要するにまあ「フロントエンドからバックエンドまで色々できてクラウドとDevOps周りもやれてデータドリブンでビジネスに貢献できるエンジニアは安泰ですよ」ということですかねｗ😁https://t.co/Jfnzn91d3a pic.twitter.com/hnS09Q2Gfy
— 勝又健太｜雑食系エンジニア｜参加者数ランキング日本第5位のオンラインサロン主催 (@poly_soft) 2020年4月8日

data-learning-guild.connpass.com

勉強会で学んだこと

今回は5件の発表がありました。
なんと嬉しいことに資料に加えて当日の録画が公開されていますので、内容をまだ確認されていない方は必見です！非常に示唆に富む内容になっています。それぞれの発表についれ感じたことを記載します。

「データマネジメントを学ぶ上での課題と環境の提供方法」

（By 村上智之　DLGギルド長）

今後、データマネジメントできる人材の価値が上がっていくことは間違いないので、なんとしてもデータマネジメントを学び、実践しなければいけません（笑）

そこで壁になるのが、「卵が先が鶏が先か？」問題です。実務経験が無いと案件に携われないし、案件に携われないと実務経験が積めないという未経験者に良くあるジレンマですね🤣

更に、実務でないまでも個人でデータマネジメントを学ぶこと自体が難しいんです。
模擬的に課題を捏造するにしてもkaggleのような分析コンペで用いられるデータは既にマネジメントされた後のデータですもんね。

村上氏の発表では、未経験者がいかにデータマネジメントの実践経験を積めば良いか述べられています。気になる方は見てみてください。

www.youtube.com

データマネジメントを個人が気軽に学べる機会を - Speaker Deck

「はじめてのデータパイプライン」

（By kotaro1）

こちらの発表は全くの未経験から約1年でデータ分析基盤の構築ができるようになった方のお話です。
上達していく過程を文学少女が野球少女になってチームで活躍できるまでの過程に例えており、非常に分かり易いです。まさに、村上氏が発表された内容を体現された方ですね！

よく感じることなのですが、学ぶ意欲があっても何から着手して良いか分からない！てことありますよね。特に未経験からキャリアチェンジした者にとって、0 → 1にするステップが難しい。。データ基盤構築についてもまさにこれが言えます。kotaro1さんが如何にして0を1,2,3..にしたか知りたい方は動画をご確認ください（またこのパターン）！

www.youtube.com

はじめてのデータパイプライン - Speaker Deck

「GUIプレパレーションツールの光と闇」

（By voovovo）

データ分析基盤構築のコンサルをバリバリやられている方のお話です。
発表の中でAlteryxというプレパレーションツールの紹介がありました。プレパレーションツールとは、ETL処理、機械学習、空間情報処理がGUIベースでできるツールのことです。プレパレーションツールを導入する目的をを仰々しく言うと、市民データサイエンティストの創造です。皆がPythonやR、SQLの知識がなくてもデータ分析基盤を構築できるようになります。

プレパレーションツールにはAlteryx以外に色々なベンダーがあるようで、主要ベンダーが一同に介して議論を繰り広げた様子がこちらにレポートされていますので、是非ご覧ください。「競合ベンダー6社が集結！第5回関西Tableauユーザー会（ETL祭り）に参加しました」ベンダー間のバチバチした感じと、ユーザーからの容赦ない質問が面白いです。Alteryx担当者のコメントに余裕を感じるのは私だけでしょうか（笑）

この発表の特筆すべき点は、プレパレーションツール導入の闇についても述べられている点です。この手の話はベンダーからなかなか引き出せないのではないでしょうか！？そして、深い闇に陥らないための対策についても述べられています（素晴らしい！）。組織論にも発展する内容となっており、サラリーマンは皆思い当たるようなシチュエーションの話題も出てきます。

「いなかったら、強くなれ、、」深い。語り尽くせないので、動画を見てください。

www.youtube.com

GUIプリパレーションツールに潜む光と闇(voovovo) - Speaker Deck

データの民主化とデータパイプラインマネジメント

（By 増田貴志）

フリーランスでデータアーキテクトをやられている方の「データの民主化」についての発表です。（個人的には「データ分析の民主化」と言っていない点が肝かと思いました）

データパイプライン完成後に発生する問題が「なる早問題」。。
なまじパイプラインができたものだから、各方面からデータ加工やら集計の依頼が舞い込む訳ですね。。
これに対処するためにデータエンジニアがやらなければいけない仕事が「データの民主化」です。

簡単に言うと、

●　誰でもデータを確認できる
●　誰でもデータをETLできる

という状態を作ることです。上述した「GUIプレパレーションツールの光と闇」の発表にも通ずるものがありますね。

「データの民主化」のメリットとしては、

●　データサイエンスチームに変な依頼が来ないようになる（SQLクエリの微修正とか）
　　→　本来の分析業務に集中できる
●　転職者がいても、データが属人化していないので引き継ぎがスムーズ
　　→　人材流出による損失を最小限にできる（データエンジニアは昨今引く手あまたです！）

データの民主化については、以下のレポートも参考になります。
「データ分析者を社内で増やすには？メルカリから学ぶ「ゆるふわBI」の取り組み」

さて、「データの民主化」のための一つの取り組みとして増田さんが提唱されているのが、「データ分析のためのSQLコード規約」です。
私はSQL初心者なのですが、SQLてPythonと違って、大文字/小文字の区別なく動いちゃったりと自由度が高いんですね。なので、スパゲッティになり易いです。規約に従うことで、可読性が上がり、「データの民主化」を加速することができます。

データサイエンティストやデータエンジニアを目指してSQLの学習を開始される方は意識してみると後々お得なことがありそうです！

qiita.com

最後に、
「自分達のためでなく、人のためにデータウェアハウスを作ることこそが、民主化に当たってのデータエンジニアの仕事」という言葉が心に残りました。

www.youtube.com

20200419_データの民主化とデータパイプライン - Google スライド

データマネジメントを伴わない経営は、破綻する

〜2つのデータ分析プロジェクトからまなぶ残酷な真実〜

（Bt yuzutas0）

さて、こちらの過激なタイトルの発表、「データマネジメントが30分でわかる本」の著者の発表です。

私は書籍を読んでからこの発表を聞いたのですが、データマネジメントに見識が浅い方は動画を見てから書籍を読むことをおススメします！動画の分かり易さが半端ないです😆

yuzutas0さんの発表は独特のスタイルで、

「A. 創業5年のベンチャー企業、
　B. 創業50年の大企業、
　一方はデータ分析をビジネス価値につなげられ、片方はデータ分析を始めることすらできなかった・・さて成功したのはどちらの企業でしょうか？」

と言う挑戦的な質問から始まります。みなさんは、どちらが成功したと思いますか？
私は直感でA.と答えました。実は、この問いに答えること自体がデータマネジメントに通ずるものがあるんです。

f:id:banquet-kuma:20200428160017p:plain:w600 — データマネジメントがないとこうなります

成功した企業と失敗した企業を分けた要因は何だったのか・・？
気になる方は是非以下をcheck !

www.youtube.com

データマネジメントなき経営は、破綻する。 #dl_guild / 20200419 - Speaker Deck

最後に

全ての動画を見ると、3時間くらいかかりますが、データ分析基盤構築やデータマネジメントに関心がある方は見て損はないです。こんなコンテンツが無料で見られる時代て、、初学者にとってはありがたいです😂

*1:データ分析人材のためのオンラインサロン。 https://data-learning.com/guild

2020-04-12

製造現場におけるデータマネジメントを考える

はじめに

データ分析業に転職してから、機械学習や統計学の手法を中心に学んで来たのですが、それってデータありきの話なんですよね。
データ分析をビジネスに活用するためには、まずデータそのものの準備や適切な管理が必要になります。Garbage in Garbage outという言葉もありますもんね。
在宅勤務で時間に余裕ができた今、データマネジメントについて学んでみることにしました。

勉強材料

データマネジメントを学ぶ書籍としては、
「データマネジメント知識体系ガイド（DMBOK）」というバイブルがありますが、600P超えで如何せん初心者にはハードルが高いです💦。
と思っていたところ、この本を実務者の経験を踏まえて要約した書籍がありました。
その名も「データマネジメントが30分でわかる本」。この本の優れている点はタイトルに書かれているとおり30分でデータマネジメントの概要が分かる点です。

章ごと所要時間別に要約が記載されています。これは非常に便利で、全体をザックリ把握してから詳細を読み進めることができます。私はまず、各章の30秒コースを全て読み、その後3分コース→30分コース→3時間コースと読み進めました。初学者の方はDMBOKを読む前に一読することをおすすめします。

●【30秒コース】一言で！
●【3分コース】なにそれ
●【30分コース】どうして
●【3時間コース】ケーススタディ

f:id:banquet-kuma:20200412012436p:plain:w600 — データマネジメントが30分でわかる本

書籍の購入はこちらから↓
Amazon CAPTCHA

製造現場におけるデータマネジメントを考えてみる

私は元半導体の生産技術者なので、良くも悪くも事例を製造現場に置き変えて考えてしまいます。
上記書籍はWeb系の事業会社を題材にしたものですが、過去の経験を踏まえて、メーカーの製造現場でのデータマネジメントについて考えたいと思います。ネタバレしない範囲で、各章ごとに学んだことを述べたいと思います。

第1章　データアーキテクチャ：

データ分析の目的はあくまで、ビジネスで価値を出すことです。データとビジネスを紐つける青写真（データアーキテクチャ）が必要になります。

例えば、製造現場に於いては、以下のことがイメージされます。
いずれも「生産性向上」のための施策になります。

●　工程ごとの進捗率が自動的に可視化される
●　工程ごとの不良率が自動的に可視化、保存される
●　装置のログが自動的に可視化、保存される
　（不良率とログは同じダッシュボードで見られるようにする。相関があるかもしれないので）

過去を振り返ると、データを現場で探すのに時間を要していたので、まずその時間を無くすことを目指したいと思います。可視化は従来Excelのグラフで行っていましたが、BIツールを使って色々な角度からデータを確認できるようにしたいです。

f:id:banquet-kuma:20200409172335p:plain:w600 — 製造工程のデータアーキテクチャ

第2章　データストレージとオペレーション：

データから価値を生みたいのであれば、データの管理を属人化させるのではなく、管理するシステムをしっかり作りましょうね！ということと理解しました。私がかつていた現場と目指す姿との比較は以下になります。

	現状	目指す姿
データの保管場所	現場のPCと工程担当者のデスクPC	データベース
データの運搬	USB	データアーキテクチャ

第3章　データ統合と相互運用性

ゴールに到達するためには、抽出(Extract)、変換(Transform)、取り込み(Load)（ETLと呼ばれる）を業務の要件を満たすように行う必要があります。そのためには、まず要件を定義しないといけません。
例えば、不良率の可視化の例だと、不良が多発した場合緊急的にラインを停止する必要があります。そのため、不良率は製品ごとにデータの更新、可視化を行う必要があります。1日の生産が終わった後に実はめちゃくちゃ不良率高かったことが分かっても後の祭りです。

第4章　データモデリングとデザイン

事業が大規模になってくると、データ量が増え、データどうしの関連が複雑になってきます。データの関連性がこんがらがるとデータ活用に於いてリスクでしかありません（不要なデータを混入させたり、必要なデータを消したり）。データどうしの関連を整理する必要があります。DMBOKではデータの粒度を次の3段階に分けて整理しています。

● 概念モデル
基本的で重要な概念のみが表現されたもの（Ex. そのラインで作られる製品群）

● 論理モデル
詳細な用件が記載されたもの。ただし、実装上の課題は考慮しない
（例）
　・その製品の工程一覧
　・工程ごとの担当者一覧
　・工程ごとの製造装置一覧

● 物理モデル
具体的な技術的ソリューションを指し、
論理モデルをデータベース上のテーブルとして表現した関連図（ER図）
例えばある製品データについてのER図を作ってみました。

f:id:banquet-kuma:20200409182330p:plain — 製品データのER図（イメージ）

この章で重要と感じた言葉に「モデルストーミング」という言葉があります。
これは、皆で紙やホワイトボードにデータモデリングの案を書き出してブラッシュアップする作業です。
色々な人の意見を聞いてコミュニケーションを上手く取りながら進めて行くのが重要そうです（どんな仕事でもそうですが）。また、大企業になると工場毎が別会社みたいなものなので、モデルストーミングするメンバーを以下のように分けたほうが良さそうです。別々の工場で同じような工程がある場合、データを共有できれば不良原因に早く辿り着けるかもしれません（他の工場で過去に同じような問題に直面し、既に解決済である可能性がある）。

f:id:banquet-kuma:20200409185215p:plain:w600 — ある製品のデータモデル

第5章　マスターデータ管理

データは一貫した信頼できる状態（マスターデータ）で管理することを学びました。
Excelファイルなどでよくあるのですが、複数似たような名前のファイルがあり、どのファイルが最新か分からず利用者が混乱する場面です（笑）マスターデータで管理することで無駄なコストを無くすことができます。書籍にはマスターデータの管理を行うのに便利なツールについて記載さています。

第6章　ドキュメントとコンテンツ管理

さてここまでは、構造化データに注目してきましたが、データには非構造化データもあります。非構造化データは管理の目が届きにくく、忘れさられがちなので注意が必要です。非構造化データには以下のようなものがあります。

● 製品の設計図
● 製品の組みたて図
● 作業要領書
● 作業手順書

非構造化データはWordに書かれたドキュメントをイメージしました。
このようなデータにはメーカーのノウハウが濃縮されているため、構造化データよりも重要かもしれません。過去にサイバー攻撃で狙われたのもこの手のデータが多いように思います😱。

第7章　データセキュリティ

ある意味この章が一番重要かもしれません。情報流出は顧客の信頼を失いますし、市場競争力も落ちます。
この章には私が前から気になっていたことが書かれていました。「クラウドでのデータ保護はオンプレと同じレベルで可能か？」についてです。作者の見解は書籍で確認いただきたいですが、私はエッジAIのクラウドAIに対する利点はセキリュティが強固なことと思っていました。しかし、そう決めるけるのは早計と思い始めました。この点については追加調査するつもりです。

第8章　データ品質

データ品質は利用者のニーズを満たしているか否か（＝ビジネスに活用できるか否か）で決まることを理解しました。例えば、検査データを例に挙げると、製品規格がμ±3σだとして、検査装置の測定ばらつきの標準偏差が3σだと話になりません。また、測定精度が要求を満たしていても欠損が多ければ問題になり得ます。極端な例ですが、この場合データ品質が低いということになりますね。

f:id:banquet-kuma:20200410213801p:plain:w600 — データ品質が低い場合に起こること

書籍にはデータ品質が低い場合の様々なリスクが載っています。
また、なるほどと思ったのですが、はじめに各データの品質を定義してやる必要があります。ゴールを明確にしておき、過剰品質を避けます。前述の検査データの例だと、測定ばらつきを0にすることはできません。なので、測定ばらつきのせいで、規格外れを規格内、規格内を規格外れと判定する確率を0にはできません。しかし、ほとんど起こり得ないことを想定してデータ品質を追求することはコストの観点から得策でないです。どこまで追求するのか初めに線引きをしておかなければいけません。

第9章　データウェアハウジングとビジネスインテリジェンス

意思決定を支えるデータを、適切なタイミングで意思決定者に提示するシステムを構築することと理解しました。単にシステムを構築するのではなく、現場のビジネスロジックに沿った可視化やレポーティングが行われなければいけません。筆者は利用者のニーズを正しく捉えるために5W1Hの観点でデータの活用を整理しており、非常に参考になりました。私が担当していた業務に置き換えてみると以下のようになります。

f:id:banquet-kuma:20200411184811p:plain:w400 — 歩留り改善のためのBI構築に向けての5W1H

データフローの設計も重要になります。従来は生産技術者が装置や現場のPCからデータを抜き取って指標を集計、Excelで可視化していましたが、以下のようなデータフローを構築すればその作業を無くすことができ、創造的な業務に時間を割くことができます。システムを構築する際の費用対効果を的確に見積もれ、周囲を牽引できる人が社内で重宝されると思います。

f:id:banquet-kuma:20200411205649p:plain:w700 — 歩留り改善に向けてのデータフロー

第10章　メタデータ管理

メタデータ管理とは「データを説明するデータを管理」することを指します。製造工程から出てくるデータを読み解くには専門知識が要るので生データは工程担当者しか見ません。また、万人が理解できる必要もないと考えています。なので、工程担当者がデータの意味を理解していれば良いとも考えられるのですが、、
市場不具合などが発生すると製造データを顧客に開示する必要が発生し、メタデータの管理が杜撰だと追及される可能性があります（だから、整備しておくわけではないですが）。ISO 9001の観点からも重要と言えそうです。

第11章　データガバナンス

これまでの章の内容を実践するために、社内のルールを定めて組織的な意思決定を行うことです。
ROIを最大化するために、以下が重要となります。

●　課題の特定
闇雲にデータ基盤を構築してもビジネスに活用できる訳ではありません。初めに現場や意思決定者と共に課題を明文化して「あるべき姿」を決めましょう。

●　案件の分割
半導体製造には多数の工程があるため、データ基盤の構築も複雑になると思われます。工場内の各ラインについて、上記のデータフローを構築していきますが、初めから完璧なものを目指すと挫折します（笑）。初めは、取り組み易い現場、製品、工程から着手しましょう。

●　段階的にリリースする
前述の内容と被るのですが、手戻りを無くすにはスモールスタートが重要ということです。基盤構築の各段階に細かくホールドポイントを設けて関係者に使い勝手などの意見を聞きながら進めましょう。

●　定性情報の活用
これはWeb業界の事例と、私が志向している製造現場の事例で話が違ってくると思いますが、データに疑問があれば、数字だけから読み解こうとするのではなく、製造員や装置メーカーにヒアリングしようということ捉えました。その事象が起こっているのが現場なので、実地でのヒアリングから有益な情報が得られることが多いです。

最後に

書籍を読了して、
データ活用でビジネス価値を創出したいなら、データマネジメントから！
ということが理解できました。

車を使って移動時間を節約したかったら、良いエンジンも大切ですが、整備された道と精製された燃料が必要ですよね？車がポンコツでもこの二つをしっかり作れば生産性が上がる事例もあると思います（例えが下手ですみません）。生産技術の場合、データからインサイトを得るのはまだ機械化できない要素が多いです。なので、データの収集、抽出、変換を自動化し、その時間を節約することで、インサイトを得るまでの時間を短縮する（ひいては意思決定までの時間を短縮する）ことに価値があると感じます。

言うは易く行うは難しで、結局データマネジメントには以下の能力が求められのではないかと思いました。やはりこの壁を超えないとですね。。

「あなたもたったの3ヶ月でデータサイエンティストに！！……なれるなんて本気で思ってるの？　～データラーニングスクール設立に込めた3つの想い～」　から引用

・ビジネス上の課題を見つけて、良き問いをしっかりと立てられること
・その問いの対して適切な解決策を提案できる幅広い知識と経験、バランス感覚
・自分の持てる武器を使ってしっかりと問題の解を導き出すこと
・そこから導き出された答を現場に浸透させてビジネス価値に転換すること

note.com

2020-04-08

COVID-19の感染者数をPlotly Expressを使って動くヒートマップで表示する

はじめに

先日ブログの題材に取り上げたCOVID-19の感染者数をPlotly Expressを使ってサックと可視化してみました。

やりたいこと

感染者数のヒートマップを地図上にプロットして、日付け毎に動かしたい。
以下のようなグラフの濃淡を日毎に変えたいわけです。濃淡は感染者数を表ます。

f:id:banquet-kuma:20200406211544p:plain — 地図上に表示したヒートマップ

Plotly Express

Pythonの可視化ライブラリPlotlyのwrapperで、簡単にjupyter notebook上でインタラクティブに可視化できます（アニメーションもOK！）。時系列データをグラフ上で簡単に動かせるの特長です。感覚としては、matplotlibやseabornと同じくらい簡単に書けて、しかもインタラクティブなグラフを作れます。
pipでインストールしておきます。

pip install plotly-express

plotly.com

Plotly ExpressをJupyterLab上で使う

JupyterLabの拡張機能にplotly-extensionなるものがあり、これをインストールしておく必要があります。npmが必要になるので事前にインストールしておきます。

jupyter labextension install @jupyterlab/plotly-extension

www.npmjs.com

github.com

アニメーションの作り方

上記のPlotly Expressのサイトにおあつらえむきなグラフがあるので丸パクします（笑）以下の4行だけでアニメーションが作れます。

import plotly.express as px

df = px.data.gapminder()
fig = px.choropleth(df, locations="iso_alpha", color="lifeExp",\
　　 hover_name="country", animation_frame="year", range_color=[20,80])
fig.show()

streamable.com

上記のコードのうち、
以下の項目をCOVID-19の感染者数のデータに差し替えます。

● df : 各国の日毎の感染者数のデータフレーム
● locations：ISO 3166-1 alpha-3に基づく国名コード
● color：濃淡を表す変数（今回の場合、感染者数）
● hover_name：地図上にカーソルを当てた時に表示する変数
● animation_frame：アニメーションを動かす軸（今回の場合、日付け）
● range_color：濃淡バーの表示範囲

感染者数を表すデータフレームの中に国名コードの列を作ってやる必要があります。感染者数のデータは前回のブログ同様以下から入手しました。
data.humdata.org

コードを書く

#  必要なライブラリのインポート
import plotly.express as px
import pandas as pd
import numpy as np

# データの読み込み
Confirmed_df = pd.read_csv("time_series-ncov-Confirmed0122-0321.csv",header=1)

# 列名を変える
Confirmed_df.rename(columns={'#adm1+name': 'State', '#country+name': 'Country',"#affected+infected+value+num":"Infected person(per day)"},inplace=True)

#  "#date"列をdatetime64[ns]型 に変更
Confirmed_df["#date"] = pd.to_datetime(Confirmed_df["#date"], format='%Y-%m-%d')

print(Confirmed_df.head())

	State	Country	#geo+lat	#geo+lon	#date	Infected person(per day)
0	nan	Afghanistan	33	65	2020/3/21	24
1	nan	Afghanistan	33	65	2020/3/20	24
2	nan	Afghanistan	33	65	2020/3/19	22
3	nan	Afghanistan	33	65	2020/3/18	22
4	nan	Afghanistan	33	65	2020/3/17	22

"Infected person(per day)"の列が1日の感染者数を表ます。

このデータフレームに国名コードを紐付けるために国名コードを取得する必要があります。以下から落とすことができます。

https://gist.github.com/tadast/8827699#file-countries_codes_and_coordinates-csv

# 国名コードを読み込む
ISO_alpha_3_df = pd.read_csv("countries_codes_and_coordinates.csv",header=0)

# コードが"AFG"という風に" "で囲まれているので、削除する。空白も削除しておく。
f_replace = lambda x:x.replace('"','').replace(' ','')
ISO_alpha_3_df = ISO_alpha_3_df.applymap(f_replace)

print(ISO_alpha_3_df.head())

	Country	Alpha-2 code	Alpha-3 code	Numeric code	Latitude (average)	Longitude (average)
0	Afghanistan	AF	AFG	4	33	65
1	Albania	AL	ALB	8	41	20
2	Algeria	DZ	DZA	12	28	3
3	AmericanSamoa	AS	ASM	16	-14.3333	-170
4	Andorra	AD	AND	20	42.5	1.6

# 国名と国名コードを紐つけた辞書を作る
dict_country = dict(zip(ISO_alpha_3_df["Country"], ISO_alpha_3_df["Alpha-3 code"]))

print(dict_country)

{'Afghanistan': 'AFG',
'Albania': 'ALB',
'Algeria': 'DZA',
'AmericanSamoa': 'ASM',
'Andorra': 'AND',
'Angola': 'AGO',
'Anguilla': 'AIA',
'Antarctica': 'ATA',
'AntiguaandBarbuda': 'ATG',
'Argentina': 'ARG',
'Armenia': 'ARM',
'Aruba': 'ABW',
・・・}

この時点で厄介な問題に直面しました。
国名コードのソースであるISO_alpha_3_dfでの国名（"Country"列）と感染者数のソースであるConfirmed_dfの国名（"Country"列）に表記揺れがあるため、dict_countryから国名コードを取得できない国が発生します。どの国が漏れるか調べます。

# 国名コードが取れない国一覧
print(set(Confirmed_df.Country)-set(dict_country.keys()))

{'Antigua and Barbuda',
'Bahamas, The',
'Bosnia and Herzegovina',
'Burkina Faso',
'Cabo Verde',
'Cape Verde',
'Central African Republic',
'Congo (Brazzaville)',
'Congo (Kinshasa)',
'Costa Rica',
"Cote d'Ivoire",
'Cruise Ship',
'Czechia',
'Dominican Republic',
'East Timor',
'El Salvador',
'Equatorial Guinea',
'Eswatini',
'Gambia, The',
'Holy See',
'Iran',
'Korea, South',
'Kosovo',
'Moldova',
'New Zealand',
'North Macedonia',
'Papua New Guinea',
'Saint Lucia',
'Saint Vincent and the Grenadines',
'San Marino',
'Saudi Arabia',
'South Africa',
'Sri Lanka',
'Taiwan*',
'Tanzania',
'Trinidad and Tobago',
'US',
'United Arab Emirates',
'United Kingdom'}

感染者数のソース(Confirmed_df)にある166ヵ国中、36ヵ国が漏れることが分かりました。注目すべき国は手動で辞書に国名コードを登録しました。それ以外の国は諦めました。

# 表記揺れで辞書に登録されない国は手動で登録する
dict_country["US"]="USA"
dict_country["United Kingdom"]="GBR"
dict_country["Iran"]="IRN"
dict_country["Taiwan*"]="TWN"
dict_country["Korea, South"]="KOR"

感染者数のソース(Confirmed_df)に国名コード列を作ります

# Confirmed_dfにiso_alpha列を作る
# Confirmed_dfの"Country"列の国名がdict_countryのキーにあれば、
# dict_countryの値をConfirmed_dfの"iso_alpha"列に加える

def func(x,dict_country=dict_country):
    if x in list(dict_country.keys()):
        return dict_country[x]
    else:
        return None
        
Confirmed_df["iso_alpha"] = Confirmed_df["Country"].apply(func)

print(Confirmed_df.head())

	State	Country	#geo+lat	#geo+lon	#date	Infected person(per day)	iso_alpha
0	nan	Afghanistan	33	65	2020-03-21 00:00:00	24	AFG
1	nan	Afghanistan	33	65	2020-03-20 00:00:00	24	AFG
2	nan	Afghanistan	33	65	2020-03-19 00:00:00	22	AFG
3	nan	Afghanistan	33	65	2020-03-18 00:00:00	22	AFG
4	nan	Afghanistan	33	65	2020-03-17 00:00:00	22	AFG

これで日時/感染者数/国名コードが1つのデータフレームに紐付きました。

欠損値を確認します。

Confirmed_df.isnull().sum()

State 9420
Country 0
#geo+lat 0
#geo+lon 0
#date 0
Infected person(per day) 0
iso_alpha 2040
dtype: int64

"state"（州）と"iso_alpha"（国コード）列に多くの欠損があります。
"iso_alpha" 列の欠損は前述したとおり、感染者数のソースと国コードのソースに表記揺れが生じるために発生しています。"iso_alpha"列に欠損があるとエラーが発生するため欠損行を削除しておきます。

# 'iso_alpha'列にNaNがある行を削除する
Confirmed_df.dropna(subset=['iso_alpha'], axis=0, inplace=True)

国別にデータフレームを作成し、縦につなげます。

# 国コードが取得できた国名の一覧をリストにする
country_list = list(set(Confirmed_df["Country"]))

# 各国別のデータフレームを格納するリストを作る
dataframe_list = []

# インデックス（'pandas.core.indexes.datetimes.DatetimeIndex'をstr型に変更する)
# これは、Plotly Expressのanimation_frame指定する型はstr or int 型である必要があるため
f_strftime = lambda x:x.strftime('%Y-%m-%d')

# 国別の日毎の感染者数のデータフレームを作ってdataframe_list に格納する
for i in range(len(country_list)):
    df2 = Confirmed_df[Confirmed_df["Country"] == country_list[i]].sort_values(by="#date")
    df2 = pd.DataFrame(df2.groupby(['#date'])['Infected person(per day)'].sum())
    df2['Infected person(per day)'] = df2['Infected person(per day)'].diff()
    df2.dropna(subset=['Infected person(per day)'],inplace = True)
    df2["Date"] = df2.index.map(f_strftime)
    df2["Country"] = country_list[i]
    df2["iso_alpha"] = dict_country[country_list[i]]
    df2.reset_index(inplace=True, drop=True)
    dataframe_list.append(df2)

#各国のデータフレームを縦に繋ぐ
df_all = pd.concat(dataframe_list, axis=0)

感染者数そのもので濃淡を付けると、中国以外の国の変化が分かりにくいので、
感染者数の対数をとって変化を見やすくします。エラーがでないようにするために、感染者数に微小な数（1.0e-10）を足しておきます。

df_all["log(Infected person(per day))"]=np.log10(df_all["Infected person(per day)"]+1.0e-10)

print(df_all.head().to_markdown())

	Infected person(per day)	Date	Country	iso_alpha	log(Infected person(per day))
0	0	2020-01-23	Togo	TGO	-10
1	0	2020-01-24	Togo	TGO	-10
2	0	2020-01-25	Togo	TGO	-10
3	0	2020-01-26	Togo	TGO	-10
4	0	2020-01-27	Togo	TGO	-10

ようやく可視化できるデータフレームが完成しまた。感染者が0の日は、"log(Infected person(per day)"の列が-10になっていますが、濃淡バーを0以上に設定するので影響ないと思われます。最後にPlotly Expressで可視化します。

# グラフの作成
fig = px.choropleth(df_all, locations="iso_alpha", color="log(Infected person(per day))", \
         hover_name="Country", animation_frame="Date", range_color=[0,4],width=1000, height=800)
fig.show()

以下のような動くヒートマップが表示されます。

streamable.com

1/23時点では色がついている（＝感染者が発生している）のは中国だけですが、2月下旬になると米国、イタリア、イランあたりに色が付き始めますね。地理的な要因（中国に近い国から感染が広がる）は無さそうです。灰色は、国名コードを取得できていな国です。Plotly Expressの優れている点は地図上にカーソルを持っていくと"hover_name"で指定した列の値が確認できる点です。

最後に

Plotly Expressは可視化自体は簡単なのですが、グラフの種類によっては前処理が面倒なものもあります。簡単なグラフはmatplotlibやseabornと同じ操作感で書けて、カーソルをあてることで値を確認できるので、便利です。
以下はkaggleのデータでEDAした時の様子です。折れ線グラフにカーソルをあてると値を確認できます

streamable.com

TableauやPower BIと棲み分けて、効率的な可視化を行っていきたいです。

2020-04-06

COVID-19のデータから感染者数のグレンジャー因果性を調べる

はじめに

世界を危機に陥れてるCOVID-19ですが、
昨今勉強を進めている時系列分析で何か知見が得られないか調べてみることにしました。
しかしながら、初学者であるため間違った解釈をしている可能性があります。
どしどし、ご指摘お願いいたしますm(_ _)m

やりたいこと

VARモデルを使って中国の感染者数から韓国、日本への感染者数に対してグレンジャー因果性があるかを検定する。検定には、Pythonのstatsmodelsを使いました。
www.statsmodels.org

VARモデル

ARモデル（自己回帰モデル）を多変量に拡張したものであり、以下の式で表されます。

$\begin{align}\begin{aligned}Y_t = \nu + A_1 Y_{t-1} + \ldots + A_p Y_{t-p} + u_t\\u_t \sim {\sf Normal}(0, \Sigma)\end{aligned}\end{align}$ 　　・・・1

$Y_t$ ：t時点の目的変数ベクトル

$Y_{t-i}$ ：t-i時点の目的変数ベクトル

$\nu$ ：n×1定数ベクトル

$A_p$ ：n×n定数ベクトル

$u_t$ ：ベクトルホワイトノイズ

Varモデルの優れている点は、自身の過去のデータだけでなく関係ありそうな他変数の過去のデータも回帰に使えることです。例えば、仮に中国、韓国、日本の感染者数にそれぞれ相関があると仮定します。
その場合、1式は以下のように書き下せます。

3変量のVAR(p)モデル

t時点の中国の感染者数=

$\nu_A + A_{11}×(t-1時点の中国の感染者数)+A_{12}×(t-1時点の韓国の感染者$
$数)+A_{13}×(t-1時点の日本の感染者数)+$
$A_{21}×(t-2時点の中国の感染者数)+A_{22}×(t-2時点の韓国の感染者+$
$数)+A_{23}×(t-2時点の日本の感染者数)+・・・+$
$A_{p1}×(t-p時点の中国の感染者数)+A_{p2}×(t-p時点の韓国の感染者数)+$
$A_{p3}×(t-p時点の日本の感染者数)+$
$u_{ta}$
・・・2-A

t時点の韓国の感染者数=

$\nu_B + B_{11}×(t-1時点の中国の感染者数)+B_{12}×(t-1時点の韓国の感染者$
$数)+B_{13}×(t-1時点の日本の感染者数)+$
$B_{21}×(t-2時点の中国の感染者数)+B_{22}×(t-2時点の韓国の感染者+$
$数)+B_{23}×(t-2時点の日本の感染者数)+・・・+$
$B_{p1}×(t-p時点の中国の感染者数)+B_{p2}×(t-p時点の韓国の感染者数)+$
$B_{p3}×(t-p時点の日本の感染者数)+$
$u_{tb}$
・・・2-B

t時点の日本の感染者数=

$\nu_C +C_{11}×(t-1時点の中国の感染者数)+C_{12}×(t-1時点の韓国の感染者$
$数)+C_{13}×(t-1時点の日本の感染者数)+$
$C_{21}×(t-2時点の中国の感染者数)+C_{22}×(t-2時点の韓国の感染者+$
$数)+C_{23}×(t-2時点の日本の感染者数)+・・・+$
$C_{p1}×(t-p時点の中国の感染者数)+C_{p2}×(t-p時点の韓国の感染者数)+$
$C_{p3}×(t-p時点の日本の感染者数)+$
$u_{tc}$
・・・2-C

沖本本pp.74によると

「モデルに含める変数は分析者の興味によって決めるものであり、目的に応じて様々です。例えば、国際株式市場の関係を分析したい場合は、日本、アメリカ、イギリスの株式収益率などからなるベクトルを考えます。（中略）重要なことは、分析者が、動学的関係に興味のある変数を全てモデルに含めることであるが、（中略）」

と書かれており、どんな変数をベクトルに加えるかは分析者の主観によるところもあるようです。

また、Varモデルで気を付けないといけないことは、非定常なデータには適用できないということです。例えば、Pythonで時系列分析を扱う際によく用いられるstatesmodelsのVarモデルのページには以下の注意書きがあります。非定常過程やトレンドを持つ過程は定常過程に変換してからVarモデルを適用しないと、不適切な結果を招くと書かれています。

f:id:banquet-kuma:20200405015818p:plain — Varモデルの注意書き

www.statsmodels.org

Varモデルの解説はこちらのサイトが分かり易かったです。
tjo.hatenablog.com

グレンジャー因果性

背後に明確な理論（経済・ファイナンス理論など）を仮定せずとも時系列データだけから因果性の有無を判断できる指標です。

沖本本pp.80には以下のように定義されています。

定義4.1(グレンジャー因果性)
現在と過去のxの値だけに基づいた将来のxの予測と、現在と過去のxとyの値に基づいた将来のxの予測をして、後者のMSE(Mean Squared Error ;平均2乗誤差)の方が小さくなる場合、 $y_t$ から $x_t$ へのグレンジャー因果性(Granger causality)が存在すると言われる。

定義4.2（一般的なグレンジャー因果性）
$x_t$ と $y_t$ をベクトル過程とする。また、xとyの現在と過去の値を含む、時点tにおいて利用可能な情報の集合を $Ω_t$ とし、 $Ω_t$ から現在と過去のyを取り除いたものを $Ω̃ _t$ とする。このとき、 $Ω̃ _t$ に基づいた将来のxの予測と、 $Ω_t$ に基づいた将来のxの予測を比較して、後者のMSEのほうが小さくなる場合、 $y_t$ から $x_t$ へのグレンジャー因果性が存在するといわれる。ここで、MSEの大小は行列の意味での大小であることに注意されたい。

例えば、日本の感染者数から中国の感染者数へのグレンジャー因果性が存在しないということは、2-A式に於いて $A_{13}=A_{23}=・・・=A_{p3}=0$ と同値になります。つまり、中国の感染者数を表す式に於いて日本の感染者数に関する係数が全て0になるということです。これは、以下の手順で判断できます。

沖本本pp.81

n変量VAR(p)モデルにおけるグレンジャー因果性検定の手順
(1)VARモデルにおける $Y_{kt}$ のモデルをOLSで推定し、その残渣平方和を $SSR_1$ とする。
(2)VARモデルにおける制約を課したモデルをOLSで推定し、その残渣平方和を $SSR_0$ とする。
(3)F統計量を $F=\frac{(SSR_0-SSR_1)/r}{SSR_1/(T-np-1)}$ で計算します。rはグレンジャー因果検定に必要な制約の数である。
(4) $rF$ を $χ^2(r)$ の95%点と比較し、 $rF$ に方が大きくければ、ある変数(群)から $Y_{kt}$ へのグレンジャー因果性は存在し、小さければグレンジャー因果性は存在しないと結論する。

グレンジャー因果性検定に於いて、
帰無仮説( $H_{0}$ )は「各変数群がその他の変数へ対してグレンジャー因果性がない」
になります。つまり、帰無仮説が棄却されればグレンジャー因果性があると判断されるわけです。

(2)の「制約を課した」という意味は、
2-A式に於いて $A_{12}=A_{22}=・・・=A_{p2}=0$ や $A_{13}=A_{23}=・・・=A_{p3}=0$ を課すという意味です。つまり、t時点の中国の感染者数をt-1〜t-p時点の中国の感染者数だけで表すと言うことです。日本、韓国の過去の感染者数を考慮したモデルでの残渣 $SSR_1$ に対して、制約を課したモデル（中国の過去の感染者数だけを考慮したモデル）での残渣 $SSR_0$ が十分に大きければ中国の感染者数と日本、韓国の感染者数に因果性があると判断するわけです。

グレンジャー因果性には長短・短所があり、まとめると以下のようになります。

長所
●定義が非常に明確である
●時系列データを用いて容易に検定できる

短所
●因果性が存在する必要条件であるが、十分条件でない
●真の因果性の方向とグレンジャー因果性の方向が一致するとは限らない
●定性的概念であり、関係の強さが測れない

因果関係の方向がはっきりしている場合を除いては、ある変数が他の変数の予測に有用かどうかと言う観点で解釈するのが良いみたいです。

データを確認する

各国の感染者数のデータは以下のサイトからcsvで入手しました。

data.humdata.org

1/22〜3/21に於ける、中国、韓国、日本、米国、イタリア、ドイツ、スペイン累積感染者数を可視化します。

# 必要なライブラリのインポート
import statsmodels as sm
import statsmodels.api as sap
from statsmodels.graphics import tsaplots
from statsmodels.tsa import stattools
from statsmodels.tsa.api import VAR, DynamicVAR
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
sns.set()

# データの読み込み
Confirmed_df = pd.read_csv("time_series-ncov-Confirmed0122-0321.csv",\
header=1,index_col="#date", parse_dates=True)

##affected+infected+value+num(感染者数)をint型にキャストしておく
Confirmed_df["#affected+infected+value+num"] = Confirmed_df["#affected+infected+value+num"].astype(int)

# 国ごとにデータを取得する
# 日付け（昇順）でソートしておく
China_pre_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "China"].sort_index()
US_pre_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "US"].sort_index()
Italy_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "Italy"].sort_index()
Spain_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "Spain"].sort_index()
Germany_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "Germany"].sort_index()
Korea_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "Korea, South"].sort_index()
Japan_cumsum = Confirmed_df[Confirmed_df["#country+name"] == "Japan"].sort_index()

# ChinaとUSは複数の州のデータがあるので国全体の合計にする
China_cumsum = China_pre_cumsum.groupby(['#date'])['#affected+infected+value+num'].sum()
US_cumsum = US_pre_cumsum.groupby(['#date'])['#affected+infected+value+num'].sum()

# 原系列の表示(感染者数の累積和)
plt.figure(figsize = (20,10))
ax = plt.subplot()

ax.plot(China_cumsum,label = "Total Confirmed in China")
ax.plot(US_cumsum,label = "Total Confirmed in US")
ax.plot(Italy_cumsum["#affected+infected+value+num"],label = "Total Confirmed in Italy")
ax.plot(Spain_cumsum["#affected+infected+value+num"],label = "Total Confirmed in Spain")
ax.plot(Germany_cumsum["#affected+infected+value+num"],label = "Total Confirmed in Germany")
ax.plot(Korea_cumsum["#affected+infected+value+num"],label = "Total Confirmed in Korea")
ax.plot(Japan_cumsum["#affected+infected+value+num"],label = "Total Confirmed in Japan")


plt.title('Coronavirus COVID-19',fontsize = 20)
plt.xlabel('date',fontsize = 20)
plt.ylabel('Total number of confirmed people',fontsize = 20)
plt.tick_params(labelsize=20)
plt.setp(ax.get_xticklabels(), rotation=90)
plt.legend(loc = "upper left",fontsize = 20)

f:id:banquet-kuma:20200405113852p:plain — 1/22 〜 3/21の各国の累積感染者数

今回使用したデータでは、中国と米国は州別にデータが入力されていたため、国全体の合計に変換しました。3/22時点では、中国の累積感染者数が最も多いですね。
日本、韓国だけを抜き出したのが以下のグラフです。
感染者数は増加傾向にあるものの、米国、イタリア、ドイツ、スペインに比べると明らかに傾きが緩いです。

f:id:banquet-kuma:20200405114325p:plain — 累積感染者数から日本、韓国のみ抜粋

次に、上記のデータに対して1階差分を取り日毎の感染者数を算出します。
1/22は前日のデータがないため、データが消失します。

# 日毎の感染者数を求める(1階差分を取る)
# 日毎の感染者数が負になっている日はデータの誤りとして除いている
China_diff = China_cumsum.diff().dropna()
China_diff = China_diff[China_diff>=0]
US_diff = US_cumsum.diff().dropna()
US_diff = US_diff[US_diff>=0]
Italy_diff = Italy_cumsum["#affected+infected+value+num"].diff().dropna()
Italy_diff = Italy_diff[Italy_diff>=0]
Spain_diff = Spain_cumsum["#affected+infected+value+num"].diff().dropna()
Spain_diff = Spain_diff[Spain_diff>=0]
Germany_diff = Germany_cumsum["#affected+infected+value+num"].diff().dropna()
Germany_diff = Germany_diff[Germany_diff>=0]
Korea_diff = Korea_cumsum["#affected+infected+value+num"].diff().dropna()
Korea_diff = Korea_diff[Korea_diff >=0]
Japan_diff = Japan_cumsum["#affected+infected+value+num"].diff().dropna()
Japan_diff  = Japan_diff[Japan_diff>=0]

#グラフの表示は累積感染者数と同じなので割愛

f:id:banquet-kuma:20200405122507p:plain — 1/22 〜 3/21の日毎の感染者数

中国のデータに於いて、2/13に巨大なピークがあります。これは湖北省が感染者の認定基準を変えたことによるものと思われます。こういう場合データをどう扱うか悩みますが（汗）、一旦このまま分析を進めることにします。

f:id:banquet-kuma:20200405124342p:plain — 中国湖北省に於ける認定基準の変更

ja.m.wikipedia.org

日毎の感染者数で韓国、日本だけを抜き出したのが以下のグラフ。
韓国はピークを過ぎているが、日本は増加トレンドにあるように見えます。

f:id:banquet-kuma:20200405181311p:plain — 日毎の感染者数から日本、韓国のみ抜粋

単位根検定を行う

前述したとおり、Varモデルは定常過程を前提としているため、過程が定常か調べておく必要があります。
これは単位根検定で調べることができます。単位根検定には主に以下の3つがあるようです。原理は理解できていないので割愛しますが、対立仮説と帰無仮説の対応は以下のとおりです。

名称	帰無仮説( $H_{0}$ )	対立仮説( $H_{1}$ )
Dickey-Fuller(DF)検定	単位根AR(1)過程である	定常AR(1)過程である
拡張DF(ADF)検定	単位根AR(p)過程である	定常AR(p)過程である
Phillips-Perron(PP)検定*	単位根AR(p)過程である	定常AR(p)過程である

（*ホワイトノイズの分散不均一性まで考慮した検定手法）

今回は日毎の中国、韓国、日本のデータに対してADF検定を行います。
結果は帰無仮説と対立仮説のモデルの置き方によって違ってくるため、3つのモデルを考えます。

-	帰無仮説( $H_{0}$ )	対立仮説( $H_{1}$ )	説明	regression
場合1	$𝑦_𝑡=𝑦_{𝑡−1}+𝑢_𝑡$	$𝑦_𝑡=𝜌𝑦_{𝑡−1}+𝑢_𝑡$ , \|𝜌\|<1	トレンド項なし、定数項なし	nc
場合2	$𝑦_𝑡=𝑦_{𝑡−1}+𝑢_𝑡$	$𝑦_𝑡=c+𝜌𝑦_{𝑡−1}+𝑢_𝑡$ , \|𝜌\|<1	トレンド項なし、定数項あり	c
場合3	$𝑦_𝑡=c+𝑦_{𝑡−1}+𝑢_𝑡$	$𝑦_𝑡=c+𝜌𝑦_{𝑡−1}+𝛿𝑡+𝑢_𝑡$ , \|𝜌\|<1	1次のトレンド項あり、定数項あり	ct

表中のregressionは以下で行うstattools.adfullerのオプションのregressionに対応します。
𝛿𝑡は1次のトレンドを表ます。

# 中国の日毎の感染者数に対してADF検定を行う
# 出力値は先頭からt値、p値、使用されたラグの数、計算に使用されたデータ数、各棄却域に対するt値、AICの最大値
# トレンド項あり（2次まで）、定数項あり
ctt = stattools.adfuller(China_diff, regression="ctt")
# トレンド項あり（1次）、定数項あり
ct = stattools.adfuller(China_diff, regression="ct")
# トレンド項なし、定数項あり
c = stattools.adfuller(China_diff, regression="c")
# トレンド項なし、定数項なし
nc = stattools.adfuller(China_diff, regression="nc")

# 韓国、日本も同様に行う

日毎の感染者数に対する単位根検定のp値をまとめると以下のようになりました。

国	場合1	場合2	場合3	結論
中国	0.051	1.0e-4	4.4e-05	場合2,3を仮定すると定常過程と言える
韓国	0.26	0.55	0.70	いずれの場合も単位根過程でないとは言えない
日本	0.74	0.80	2.7e-7	場合3を仮定すると定常過程と言える

韓国のデータは定常過程と言えないので除外し、中国の感染者数と日本の感染者数のグレンジャー因果性だけを検定することにします。

Varモデルの構築

中国と日本の日毎の感染者数のデータでからVarモデルを構築します。

# 中国と日本のデータで共通する日付けのデータだけを抜き出す
diff_CJ = pd.merge(China_diff, Japan_diff, on='#date', how='inner')
# 列名を変更する
diff_CJ.columns = ["China_diff","Japan_diff"]

#Varモデルを作成する
model_diff = VAR(diff_CJ)

フィッティングにあたってモデルのラグ（1式のp）を決めてやる必要があります。
現時点から何時点過去に遡って回帰式を構成するかを決める作業です。
これは、情報量基準(AIC、BIC）を指標に決めます。

# 最適なラグを探索する。ラグ10まで
model_diff.select_order(10).summary()

f:id:banquet-kuma:20200405210124p:plain:w300 — 各ラグに於ける情報量基準の計算

ラグ1まで考慮した場合に情報量基準が最小になりました。
実際は以下のようにmaxlags=10とic="aic"と指定してやることでラグ10までの範囲でAICが最小となるラグを自動的選んでフィッティングしてくれます。

# AIC基準で最適なラグを選択したモデルへのあてはめ
result_diff = model_diff.fit(maxlags=10, ic="aic")

#あてはめの結果を表示
result_diff.summary()

結果のサマリーを表示します。

f:id:banquet-kuma:20200406012134p:plain — あてはめの結果

サマリーの中に記載されている係数を用いて、
中国の感染者数と日本の感染者数は以下のように書けます。

t時点の中国の感染者数=

$const+ L1.China\_diff×(t-1時点の中国の感染者数)+$
$L1.Japan\_diff×(t-1時点の日本の感染者数)$

t時点の日本の感染者数=

$const+ L1.China\_diff×(t-1時点の中国の感染者数)+$
$L1.Japan\_diff×(t-1時点の日本の感染者数)$

それぞれの係数 $L1.China\_diff$ 、 $L1.Japan\_diff$ は、
中国の感染者数(China_diff)については、
「Results for equation China_diff」に記載されている値を
日本の感染者数(Japan_diff)については、
「Results for equation Japan_diff」に記載されている値を使います。

黒枠内に示した回帰係数のp値を見ると、
日本の感染者数に於いては、 $L1.China\_diff$ が有意になっておらず
（つまり、帰無仮説 $L1.China\_diff$ =0を棄却できない）、日本の感染者数に及ぼす
中国の感染者数の影響が有意でないことが示唆されました。

ここまでの結果で、中国の感染者数から日本の感染者数に対してグレンジャー因果性がないことが示唆されますが、一応グレンジャー因果性検定をやってみます。

result_diff.test_causality("Japan_diff", "China_diff", kind = "f").summary()

f:id:banquet-kuma:20200405215817p:plain:w300 — グレンジャー因果性検定の結果

帰無仮説：
「中国の感染者数から日本の感染者数に対してグレンジャー因果性が存在しない」を棄却できませんでした。したがって、昨日の中国の感染者数と今日の日本の感染者数にグレンジャー因果性がないと結論付けました。

この結果は以下を考えると妥当という気がします（分析手法が正しければですが、、）。
●1日という短期間で中国 → 日本間でウイルスの媒介が波及しそうにないこと
●感染が発覚した人のみのデータを使っているので相関が上手く拾えていない可能性があること

最後に

COVID-19の実データを使ってグレンジャー因果性検定をやってみましたが、かなり難しさを感じました。テキストに付随するデータだとすんなり結論を導けるのですが、実データとなると「そもそもこのモデルを使って良いデータなのか？前提条件を満たしているのか？」と言った疑念が浮かんできます（笑）。実務やコンペなどで経験を積んで正しく分析できるようになりたいです。

参考文献

手法やコードは以下の書籍を参考にしました。

説明がコンパクトに纏まっていて非常に分かり易いです。
www.amazon.co.jp

Pythonによる時系列分析のコードが豊富です。
www.amazon.co.jp

はじめに

勉強会の内容

チームでの活動内容

懇親会での裏話

所感

はじめに

データセットの入手

指数分布

発生間隔のプロット

ガンマ分布

最後に

参考図書

エリア51への行き方

最後に

はじめに

勉強会で学んだこと

「データマネジメントを学ぶ上での課題と環境の提供方法」

（By 村上智之 DLGギルド長）

「はじめてのデータパイプライン」

（By kotaro1）

「GUIプレパレーションツールの光と闇」

（By voovovo）

データの民主化とデータパイプラインマネジメント

（By 増田貴志）

データマネジメントを伴わない経営は、破綻する

〜2つのデータ分析プロジェクトからまなぶ残酷な真実〜

（Bt yuzutas0）

最後に

はじめに

勉強材料

製造現場におけるデータマネジメントを考えてみる

第1章 データアーキテクチャ：

第2章 データストレージとオペレーション：

第3章 データ統合と相互運用性

第4章 データモデリングとデザイン

第5章 マスターデータ管理

第6章 ドキュメントとコンテンツ管理

第7章 データセキュリティ

第8章 データ品質

第9章 データウェアハウジングとビジネスインテリジェンス

第10章 メタデータ管理

第11章 データガバナンス

最後に

はじめに

やりたいこと

Plotly Express

Plotly ExpressをJupyterLab上で使う

アニメーションの作り方

コードを書く

最後に

はじめに

やりたいこと

VARモデル

グレンジャー因果性

データを確認する

単位根検定を行う

Varモデルの構築

最後に

参考文献

（By 村上智之　DLGギルド長）

第1章　データアーキテクチャ：

第2章　データストレージとオペレーション：

第3章　データ統合と相互運用性

第4章　データモデリングとデザイン

第5章　マスターデータ管理

第6章　ドキュメントとコンテンツ管理

第7章　データセキュリティ

第8章　データ品質

第9章　データウェアハウジングとビジネスインテリジェンス

第10章　メタデータ管理

第11章　データガバナンス