俺のOneNote

俺のOneNote

データ分析が仕事な人のOneNote愛とか、分析小話とか。

pandasだけでWEBスクレイピングする

一般的にスクレイピングを使用とすると、requestsでHTMLなりjsonなりを拾ってくるのが一般的かと思われます。

Python, Requestsの使い方 | note.nkmk.me

ただ、tableタグで構造化されている場合、pandasだけでデータを拾ってこれるので大変お手軽です。

例えばwikiでオリンピックの各国メダル数をみると、HTMLのテーブル形式に整理されています。

f:id:kopaprin:20200430141023p:plain


f:id:kopaprin:20200430141141p:plain

以下、URLを指定してpandas.read_htmlするだけです。

import pandas as pd
df_list = pd.read_html("https://ja.wikipedia.org/wiki/%E8%BF%91%E4%BB%A3%E3%82%AA%E3%83%AA%E3%83%B3%E3%83%94%E3%83%83%E3%82%AF%E3%81%A7%E3%81%AE%E5%9B%BD%E3%83%BB%E5%9C%B0%E5%9F%9F%E5%88%A5%E3%83%A1%E3%83%80%E3%83%AB%E7%B7%8F%E7%8D%B2%E5%BE%97%E6%95%B0%E4%B8%80%E8%A6%A7")

戻り値はlistです。

print(type(df_list))
>>> <class 'list'>

自分が欲しい要素を確認して完了です。

df_list[0].head()
>>>
国・地域    夏季参加回数  金 銀 銅 合計  冬季参加回数  金.1.1.1  合計.1   参加回数    金.2.2.2  総数
0  アフガニスタン (AFG) 14 0  0  2  2  0  0  0  0  0  14 0  0  2  2
1  アルジェリア (ALG)    13 5  4  8  17 3  0  0  0  0  16 5  4  8  17
2  アルゼンチン (ARG)    24 21 25 28 74 19 0  0  0  0  43 21 25 28 74
3  アルメニア (ARM)   6  2  6  6  14 7  0  0  0  0  13 2  6  6  14
4  オーストララシア (ANZ) [ANZ]    2  3  4  5  12 0  0  0  0  0  2  3  4  5  12

BeautifulSoupセレクタ等の知識も不要で、HTMLから正規表現で抽出する必要もなく、便利です。
クローリングする場合にも、もしtableタグの内容であればpandas経由すると簡単にテーブルデータ化できることが多々あります。

Power BIとの出会いが人生を変えた話。

この記事は Microsoft Power BI Advent Calendar 2019の12/22 担当分です

f:id:kopaprin:20191222213307p:plain

お話の要点

タイトルとおり、「Power BI」との出会いが僕の人生を変えた話。
ただのポエムなので、Qiitaでなくこちらに書きました。

  • 新しい技術・知識の習得に取り組むのはいいぞ、という話
  • コミュニティはいいぞ、という話
  • 清水の舞台から飛び降りるといいぞ、という話

あんまりPower BI関係なくてすみません。

1 出会い編

当時の僕は、「オープンデータ」を可視化して共有できるプラットフォームを作りたいと思っていました。
もちろん、当時はRESASのようなWEBサイトは無く、HTML/CSS、js、PHPSQLを独学でイチから学びました。
エンジニアではなかった僕はもちろんうまくいくわけがなく、途方に暮れていました。
そんな中で会ったのがBIツール 「Power BI」でした。

データを簡単に取り込み、
WEBで共有でき、
インタラクティブに可視化できる。

理想としてたプラットフォームの姿がそのままそこにありました。
データサイエンス分野で生きていくとしたら、BIツールは絶対に習熟すべき!
そう確信したのを今でもよく覚えています。

2 コミュニティ参加編

独学で学んでいましたが、当時は充実した学習コンテンツもなくQiitaのような情報提供のブログもほとんどない状況でした。
英語も得意でない僕のほぼ唯一の情報源がyoshida taikiさんの吉田の備忘録でした。

吉田の備忘録に助けられつつ、必死にreferenceを読みながら学習していく中、 勉強会というものが開催されていることを知りました。
われらがPBI神が主催されているPowerBI勉強会コミュニティです。

何か学びのきっかけになるかもしれないと考えた僕はさっそく参加応募。 2016/11/26(土)のPowerBI勉強会が人生初の社外コミュニティ参加になりました。

初めて勉強会というものに出たら、同じ課題、ツール、技術を追う人たちが 社外で営利関係なく情報交換・情報共有している事実がとても新鮮で、ちょっとしたカルチャーショックになりました。
今まで社内の事しか考えていなかった僕は、初めて自分の社会一般の市場価値について意識しました。

3 登壇編

勉強会で魅力的な人にたくさんお会いし感銘を受けた僕。
その日以降変わったことが色々ありました。

  • 将来は自分も第一線で活躍している皆さんと肩を並べて話してみたいなと思ったこと
  • 社内だけでなく、社外も踏まえて自分の知識とスキルを磨こうと思ったこと
  • もじもじしてないで色々活動していこうと思ったこと

そんな意識でPower BI 勉強会 #3 のLTに手を挙げたのが1回目。
社外発信など初めてだったので、「清水の舞台から飛び降りてみた」。そんな感覚で申し込んだLTでした。
その後も、色々な場面で登壇など情報発信の機会をいただき、コミュニティ活動の楽しさに惹かれていくのでした。

4 MVP受賞&転職&出版編

コミュニティで登壇したり、必死に技術向上に励んだりした結果、人生の転機になったことに色々恵まれました。

MVPは家庭事情で長く続けることができませんでしたが、コミュニティ活動でアウトプットしてきた経験や、たくさんの尊敬できる仲間に恵まれたことは大きな財産になっています。
正直、技術コミュニティに参加していなかったら今の僕は無かったと思います。

転職も大きな決断でしたが、自分の技術がどこまで通用するのか計る大きなきっかけになりました。
自分ではたいしたスキル・経験ではないと思っていても、意外と通用することも多いんだな。と思えるようになりました。

共著出版の機会をいただけたのも、こういった活動を続けてこれた結果からかな、とも思います。
経験も知識も技術もまだまだ未熟ですが、機会に恵まれるようどんどんインプット・アウトプットを続けていきたいと思います。

5 今後編

Power BI と出会って人生が変わった今後。

  • 新しい技術・知識の習得に取り組むのはいいぞ、という話
    今まで、自分は頭が悪いので、難しいことは理解できない・習得できないと思っていました。
    無論、人より頭が弱いので時間はかかりますが、まじめに向き合えば何でも人並み程度には習得できると理解できたのは大きな収穫。
    今後もデータサイエンス分野の技術向上に邁進したい所存です。

  • コミュニティはいいぞ、という話
    社内に籠ると、情報収集・発信の視野が狭くなり、簡単に井の中の蛙になりそうです。
    コミュニティ活動、ついったらんど等、外部との関係を積極的にとっていきたい所存です。

  • 清水の舞台から飛び降りるといいぞ、という話
    自分の成長は、常に清水の舞台を飛び降りることから始まりました。
    逆に言えば、未知なことに挑戦しなければ一生成長することはなさそうです。
    リスク承知で、手を伸ばせそうなことには積極的に手を伸ばしていきたい所存です。

おわりに

ここまで書いてみると、Power BI に出会ったところから人生変わったなぁと感じつつ、変わったのはコミュニティで凄い人たちにたくさん会ってきたからかもしれません。
皆さんの活躍を見ていると、いつも「自分全然雑魚やん、がんばらんと・・・」と思わせてくれます。

自分の周囲をスゴイ人たちで囲んでおくのが、人生を変える一番大事な事なのではないかなぁ。  

Python機械学習タスクにおける欠損補完の精度をいろいろ検証してみた件。

社内機械学習コンペにて、欠損補完の精度およびシミュレーションが個人的課題になったので、
振り返りの自習&覚書です。

※この記事の欠損処理方法は、理論的に正しくない・妥当でない場合が多々あるかもしれません。  中の人の勉強不足なので、さらっと流していただくか、正しい認識をご教示いただけると喜びます。  処理内容間違ってたらぜひお知らせくださいm(__)m

※なお、本記事は以下Colaboratory公開ノートブックの焼き直しです。

colab.research.google.com

今回の自習内容は、sklearn や fancyimpute 等、欠損値・行列補完アルゴリズムによる欠損値補完処理により、機械学習による推計精度にどの程度影響するか、様々なパターンをシミュレーションすることで、今後の機械学習タスクの参考にしようとするものです。

まずは必要なライブラリ準備。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import fancyimpute
from sklearn.datasets import load_boston
import warnings
warnings.filterwarnings('ignore') # 実行に影響ないアラートを無視
np.random.seed(123)

統計科学では、欠損値は主に以下3分類で考えられるケースが多いです。

欠損値の種類 内容
MCAR (Missing Completely At Random) 完全にランダムに欠損
MAR (Missing At Random) 観測データに依存する欠損
MNAR (Missing Not At Random) 欠損データに依存する欠損

完全にランダムな欠損(MCAR)
データに依存せずに発生する欠損値です。

観測データに依存する欠損(MAR)
欠損有無が、他のデータと相関する欠損値です。
(例えば、男性は欠損値が多い、夜の時間は欠損値が多い etc)

欠損データに依存する欠損(MNAR)
欠損有無が、欠損データそのものと相関している欠損値です。データを見ただけでは判別が難しいのが特徴です。
(例えば、データが大きいほど欠損値の確率が増える etc)

参考:Suny side up!:欠損値があるデータの分析(清水裕士)

この記事では上記3種の欠損を意図的に発生させ、
各種欠損値処理の方法により、 target をどの程度正確に予測できるかを検証します。

1 sklearn Boston データセット

今回利用するデータセットはsklearnのサンプルデータbostonデータセットです。

boston = load_boston()
target = pd.DataFrame(boston.target)
boston_df = pd.DataFrame(boston.data, columns=boston.feature_names)
target.columns = ["target"]
boston_df = pd.concat([boston_df, target],axis=1)
print("bostonデータセット:",boston_df.shape)
display(boston_df.head(3))

f:id:kopaprin:20190808013948j:plain

1-1 欠損値データセットの作成

bostonデータセット全体の506レコード中、任意の100レコードについて、欠損値生成を行った3種のデータを用意します。

欠損を発生させるカラムは LSTATです。( 給与の低い職業に従事する人口の割合 (%))

これは比較的 target と相関が高く、データが適度に散布されているため、このカラムを選定しています。

plt.figure(figsize=(8,5),facecolor="white")
sns.distplot(boston_df["LSTAT"])
plt.title("LSTAT")
plt.show()

sns.relplot(data=boston_df, x="LSTAT", y="target",height=5,aspect=1.5)
plt.title("LSTAT * target")
plt.show()

f:id:kopaprin:20190808014216p:plainf:id:kopaprin:20190808014220p:plain

MCAR

完全無作為に欠損値を生成したDataFrameを boston_df_mcar として格納します。

np.random.seed(123)
rand_mcar = np.sort(np.random.choice(boston_df.index,100,replace=False))

def make_null(df,index,col_name): # 欠損値代入の関数です。以下、この関数を流用します。
  df1 = df.drop(index)
  df2 = df.iloc[index, :]
  df_c = df2[col_name]
  df2[col_name] = np.nan
  result_df = pd.concat([df1,df2]).sort_index()
  return result_df, df_c

boston_df_mcar , mcar_correct = make_null(boston_df, rand_mcar, "LSTAT")

MAR

その他観測地と相関のある欠損生成を行います。 少し簡易的ですが、カラム AGE の上位200レコードから80レコード、下位200レコードから20レコード分の欠損値を無作為に生成します。

top200 = boston_df.sort_values("AGE",ascending=False).head(200)
under200 = boston_df.sort_values("AGE",ascending=True).head(200)

rand_top = np.random.choice(top200.index,80,replace=False)
rand_under = np.random.choice(under200.index,20,replace=False)
rand_mar = np.hstack((rand_top, rand_under))

boston_df_mar , mar_correct = make_null(boston_df, rand_mar, "LSTAT")

MNAR

欠損データそのものに依存した欠損生成を行います。
ここでは極端ですが、 カラム LSTAT の上位100レコードを欠損値にします。

np.random.seed(789)

top_mnar100 = boston_df.sort_values("LSTAT",ascending=False).head(100)

rand_mnar = np.random.choice(top_mnar100.index,100,replace=False)

boston_df_mnar , mnar_correct = make_null(boston_df, rand_mnar, "LSTAT")

1-2 欠損データの形状確認

これで、3種の欠損パターンを生成したデータセットが完成しました。 内容を少し確認します。

print("3つのデータセットの形状確認")
print(boston_df_mcar.shape)
print(boston_df_mar.shape)
print(boston_df_mnar.shape)

print("欠損処理した「LSTAT」カラムのデータ確認")

plt.figure(figsize=(10,4),facecolor="white")
sns.distplot(boston_df["LSTAT"])
sns.distplot(boston_df_mcar.dropna()["LSTAT"])
sns.distplot(boston_df_mar.dropna()["LSTAT"])
sns.distplot(boston_df_mnar.dropna()["LSTAT"])

plt.legend(["master","MCAR","MAR","MNAR"])
plt.show()

3つのデータセットの形状確認 (506, 14) (506, 14) (506, 14) 欠損処理した「LSTAT」カラムのデータ確認 f:id:kopaprin:20190808014558p:plain

LSTATカラムは MCAR < MAR < MNAR の順にいびつな形状になっているようです。

1-3トレーニングデータと予測対象データの作成

学習させるデータと、予測するデータを分割します。
3つのデータセットがありますが、予測対象にするデータはそれぞれ同じインデックスとし、比較対象できるように設計します。

np.random.seed(999)
rand_i = np.sort(np.random.choice(boston_df.index,100,replace=False))

mcar_train = boston_df_mcar.drop(rand_i)
mcar_test = boston_df_mcar.iloc[rand_i,:]

mar_train = boston_df_mar.drop(rand_i)
mar_test = boston_df_mar.iloc[rand_i,:]

mnar_train = boston_df_mnar.drop(rand_i)
mnar_test = boston_df_mnar.iloc[rand_i,:]

作成したデータを以下のとおり整理します。

データセット 内容
target trainの予測用正解データ
answer testの正解データ(精度検証用)
mcar_train MCARによるデータセット(train)
mcar_test MCARによるデータセット(test)
mar_train MARによるデータセット(train)
mar_test MARによるデータセット(test)
mnar_train MNARによるデータセット(train)
mnar_test MNARによるデータセット(test)

2 欠損値補完

欠損値の補完方法には、いくつかの種類があります。

  1. リストワイズ
    1つでも欠損値のあるサンプル(または変数)は分析から除外 ※厳密にはサンプル削除はリストワイズとは呼ばないかもしれません。
     ここでは便宜上、両者をリストワイズ(rowとcolumnで区分)と呼ばせていただきます。

  2. ペアワイズ
    分析対象にする変数に欠損値が含まれていなければすべて分析対象に含める

  3. 単一代入
    欠損を何らかのデータで補完する

  4. 多重代入 欠損値に異なる値を推定した複数のデータを生成し、その分析結果を統合し補完する

参考:欠測データの統計解析 (統計解析スタンダード)


ペアワイズは機械学習に向いておらず、多重代入は簡易に実行するパッケージが整備されていません。
そこで、本のノートブックでは、リストワイズと単一代入をいくつかのパターンで組み合わせ、どのように精度が変化するか検証します。

2-1 trainデータの欠損処理

trainデータは、主に2種類の処理方法があります。

  1. リストワイズで欠損値を含むcolumn、またはサンプル(row)を削除する

  2. 欠損値を補完する、または欠損値の情報を埋め込む

「1」のケースは、欠損データと比較してtrainデータの量が多く、欠損値の種類がMCARであることが自明なときに行うことが多いと思われます。
一部の欠損したサンプル(row)を削除しても、学習には支障がないと判断した場合、こちらの処理のほうが適切かもしれません。
(補完する処理は、一般的にはバイアスを生み、精度としては実データより悪くなることが想定されるためです)

ただし、以下のようなケースでリストワイズによりcolumn, rowを削除するのは適切だとは言えません。

- 欠損がMAR, MNARの(可能性がある)ケースでrowを削除する
欠損しているデータに何らかの傾向があるため、安易に欠損しているレコード(row)を削除すると、データセットに偏りが生じる可能性があります。

- 欠損が存在するcolumnが、targetの精度に強く影響するにも関わらずcolumnを削除する
欠損が生じているcolumnそのものを削除することも手法としてはありますが、targetの精度に影響するような重要なデータである場合、column丸ごと捨ててしまうのは適切であるとは言えません。
このケースでは、多少バイアスが発生したとしても、補完するなどの処理をしたほうが有効である場合が多々あります。
また、そもそもcolumnの数が少ない場合などは、なるべくデータを有効活用するよう、安易に削除しないほうが良いと思われます。


「2」のケースは、元々のデータ数が少ない場合や、欠損値がMAR, MNARの可能性が捨てきれない場合の対応がほとんどです。

単一代入方法については、今回は以下2種を試すこととします。

項目 内容
代表値代入法 平均値・中央値・最頻値など、欠損があるcolumnの統計表を代入する方法。
(確率的)回帰代入法 欠損していない変数を説明変数として欠損値を予測する方法。



  • 代表値代入法
    簡易にデータ補完できる反面、データのバイアスがかなり大きくなることが懸念される手法です。

  • 回帰代入法
    統計モデル(機械学習モデル)を使用し、欠損値を予測します。他のデータを使用して予測するため、比較的バイアスを抑えられることが予測される。


また、この中には分類されませんが、Pythonパッケージの fancyimputeには欠損補完用の行列補完アルゴリズムがいくつか用意されています。

今回はこのパッケージ内のアルゴリズムもいくつか選定し、試してみることにします。

fancyimpute

2-2 testデータの欠損処理

testデータの欠損処理方法については、注意すべき点が2点あります。

  • リストワイズ(row)は不可
    モデルを学習するためのtrainデータでは、どのデータを学習に用いるかは設計者が自由に調整できます。
    この点から、欠損が発生しているサンプル(row)を削除することは問題ありません。
    しかし予測対象のサンプル(row)を任意に削除することはできません。
    「このサンプルは予測しません」ということは実務上でも競技上でも許されません。

  • 単一代入の平均値や学習データはtrainを使用
    代入するデータに使用する統計量やモデルは、trainデータから学習させなければなりません。 ※ 予測対象データは、課題によって1レコードずつ与えらるデータかもしれません。 また、予測対象のデータを学習に使用すると、予測対象が違うたびにモデルも異なった内容になってしまい、実務上は好ましくありません。
    したがって、代入に使用するデータはすべてtrainデータを基にすることが原則となっています。

2-2 欠損補完データの生成

上記踏まえ、3種のデータフレームにを基にした欠損処理後データフレームを作成します。

削除方法 内容 適用データ
A リストワイズ(row) 欠損が発生しているレコード(row)の削除 trainのみ
B 平均値代入 trainデータの平均を代入 train, test
C 回帰代入(LinearRegression) 線形モデルによる推計値を代入 train, test
D 回帰代入(RandomForest) ランダムフォレストによる推計値を代入 train, test
E 回帰代入(LightGB) LightGBMによる推計値を代入 train, test
F 行列補完(KNN) fancyimpute.KNN(最近傍法)を利用した行列補完 train, test
G 行列補完(softimpute) fancyimpute.softimpute(SVD分解の反復処理)を利用した行列補完 train, test


まず、train,testデータセットのリストを作成し、同リストを引数として渡すと、代入したリストを返す関数を作成していきます。

mcar_list = [mcar_train, mcar_test]
mar_list = [mar_train, mar_test]
mnar_list = [mnar_train, mnar_test]

A リストワイズ(row)

欠損が発生しているレコード(row)を削除します。

def listwise_row(df_list):
  for i,df in enumerate(df_list):
    if i == 0 :
      result_df = df.dropna()
    else:
      pass
  return result_df


mcar_A_train = listwise_row(mcar_list)
mar_A_train = listwise_row(mar_list)
mnar_A_train = listwise_row(mnar_list)

B 平均値代入

欠損値に平均値を代入します。

def sub_mean(df_list):
  df1 = df_list[0].fillna(df_list[0].mean())
  df2 = df_list[1].fillna(df_list[0].mean())
  return df1, df2

mcar_B_train, mcar_B_test = sub_mean(mcar_list)
mar_B_train, mar_B_test = sub_mean(mar_list)
mnar_B_train, mnar_B_test = sub_mean(mnar_list)

C 回帰代入(LinearRegression)

scikit-learnのLinearRegressionを用いた回帰代入を行います。

from sklearn.linear_model import LinearRegression

def sub_lr(df_list):
  train_master = df_list[0].dropna()
  train_miss = df_list[0][df_list[0].isnull().any(axis=1)]
  test_master = df_list[1].dropna()
  test_miss = df_list[1][df_list[1].isnull().any(axis=1)]
  model = LinearRegression().fit(train_master.drop(["LSTAT","target"],axis=1), train_master["LSTAT"])
  train_miss["LSTAT"] = model.predict(train_miss.drop(["LSTAT","target"],axis=1))
  test_miss["LSTAT"] = model.predict(test_miss.drop(["LSTAT","target"],axis=1))
  
  train = pd.concat([train_master,train_miss]).sort_index()
  test  = pd.concat([test_master,test_miss]).sort_index()
  return train, test

mcar_C_train, mcar_C_test = sub_lr(mcar_list)
mar_C_train, mar_C_test = sub_lr(mar_list)
mnar_C_train, mnar_C_test = sub_lr(mnar_list)

D 回帰代入(RandomForest)

scikit-learnのRandomForestを用いた回帰代入を行います。

from sklearn.ensemble import RandomForestRegressor

def sub_rfr(df_list):
  train_master = df_list[0].dropna()
  train_miss = df_list[0][df_list[0].isnull().any(axis=1)]
  test_master = df_list[1].dropna()
  test_miss = df_list[1][df_list[1].isnull().any(axis=1)]
  model = RandomForestRegressor(random_state=1).fit(train_master.drop(["LSTAT","target"],axis=1), train_master["LSTAT"])
  train_miss["LSTAT"] = model.predict(train_miss.drop(["LSTAT","target"],axis=1))
  test_miss["LSTAT"] = model.predict(test_miss.drop(["LSTAT","target"],axis=1))
  
  train = pd.concat([train_master,train_miss]).sort_index()
  test  = pd.concat([test_master,test_miss]).sort_index()
  return train, test

mcar_D_train, mcar_D_test = sub_rfr(mcar_list)
mar_D_train, mar_D_test = sub_rfr(mar_list)
mnar_D_train, mnar_D_test = sub_rfr(mnar_list)

E 回帰代入(LightGBM)

lightgbmを用いた回帰代入を行います。

from lightgbm import LGBMRegressor
from sklearn.model_selection import train_test_split

def sub_lgbm(df_list):
  train_master = df_list[0].dropna()
  train_miss = df_list[0][df_list[0].isnull().any(axis=1)]
  test_master = df_list[1].dropna()
  test_miss = df_list[1][df_list[1].isnull().any(axis=1)]
  
  X_train, X_val, y_train, y_val = train_test_split(train_master.drop(["LSTAT","target"],axis=1),
                                                   train_master["LSTAT"],
                                                   test_size=0.2,
                                                   random_state=1)
  
    
  model = LGBMRegressor(random_state=1).fit(X_train, y_train,
                                            eval_set=[(X_val, y_val)],
                                            eval_metric = "mean_squared_error",
                                            early_stopping_rounds=20,
                                            verbose=20)
  train_miss["LSTAT"] = model.predict(train_miss.drop(["LSTAT","target"],axis=1))
  test_miss["LSTAT"] = model.predict(test_miss.drop(["LSTAT","target"],axis=1))
  
  train = pd.concat([train_master,train_miss]).sort_index()
  test  = pd.concat([test_master,test_miss]).sort_index()
  return train, test

mcar_E_train, mcar_E_test = sub_lgbm(mcar_list)
mar_E_train, mar_E_test = sub_lgbm(mar_list)
mnar_E_train, mnar_E_test = sub_lgbm(mnar_list)

F 行列補完(KNN)

fancyimpute KNN による最近傍法を利用した行列補完

from fancyimpute import KNN

def sub_knn(df_list):
  model = KNN(k=3)
  result_list=[]
  for i,df in enumerate(df_list):
    df_c = df.copy()
    pred = model.fit_transform(df_c.drop("target",axis=1))
    df_c["LSTAT"] = pred[:,12]
    result_list.append(df_c)
  return result_list[0], result_list[1]

mcar_F_train, mcar_F_test = sub_knn(mcar_list)
mar_F_train, mar_F_test = sub_knn(mar_list)
mnar_F_train, mnar_F_test = sub_knn(mnar_list)

G 行列補完 (softimpute)

fancyimpute softimputeによる行列補完

from fancyimpute import BiScaler, SoftImpute


def sub_soft(df_list):
  norm = BiScaler()
  model = SoftImpute()
  result_list = []
  for i, df in enumerate(df_list):
    df_c = df.copy()
    pred = model.fit_transform(df_c.drop("target",axis=1).values)
    df_c["LSTAT"] = pred[:,12]
    result_list.append(df_c)
  return result_list[0], result_list[1]


mcar_G_train, mcar_G_test = sub_soft(mcar_list)
mar_G_train, mar_G_test = sub_soft(mar_list)
mnar_G_train, mnar_G_test = sub_soft(mnar_list)

2-3 train, test全パターンデータの整理

# train_list
mcar_train_list = [mcar_A_train, mcar_B_train, mcar_C_train, mcar_D_train,
                   mcar_E_train, mcar_F_train,mcar_G_train]
mar_train_list = [mar_A_train, mar_B_train, mar_C_train, mar_D_train,
                  mar_E_train, mar_F_train, mar_G_train]
mnar_train_list = [mnar_A_train, mnar_B_train, mnar_C_train, mnar_D_train,
                   mnar_E_train, mnar_F_train, mnar_G_train]

# test_list
mcar_test_list = [mcar_B_test, mcar_C_test, mcar_D_test, mcar_E_test,
                  mcar_F_test, mcar_G_test]
mar_test_list = [mar_B_test, mar_C_test, mar_D_test, mar_E_test,
                  mar_F_test, mar_G_test]
mnar_test_list = [mnar_B_test, mnar_C_test, mnar_D_test, mnar_E_test,
                  mnar_F_test, mnar_G_test]

3 学習・予測・精度検証

この時点で、MCAR, MAR, MNARで作成した欠損を補完したデータフレームが、それぞれ補完方法の違いによりA~Hの種類で準備しました。

ここから、MCAR, MAR, MNARごとに学習・予測・精度検証を行い、どのような欠損にはどのような補完方法が適当なのかシミュレーションしていきます。

なお、精度検証で利用する機械学習アルゴリズムは線形モデル(LinearRegression)とツリーモデル(RandomForest)とします。

評価指標はMSE(平均二乗誤差)とし、MSEが少ないほうが精度が良いとされます。

3-1 精度検証用の予測モデルの作成と精度抽出

まず、以下のとおり使用ライブラリをインポートします。

from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error

線形モデル(LinearRegression)により、全train, testデータの組み合わせによるMSE精度をリスト化します。

# 線形モデルによる精度リストの抽出

lr = LinearRegression()

def predict(train_list, test_list, model):
  mse_list = []
  for i,train in enumerate(train_list):
    train_X = train.drop("target",axis=1)
    target = train["target"]
    model.fit(train_X,target)
    for j,test in enumerate(test_list):
      test_X = test.drop("target",axis=1)
      true = test["target"]
      pred = model.predict(test_X)
      metric = mean_squared_error(true, pred)
      mse_list.append(metric)
  return mse_list

MCAR_LR = predict(mcar_train_list, mcar_test_list, lr)
MAR_LR = predict(mar_train_list, mar_test_list, lr)
MNAR_LR = predict(mnar_train_list, mnar_test_list, lr)

以下で、MCAR, MAR, MNARそれぞれのデータフレームに整えます。

MCAR_LR_df = pd.DataFrame(np.array(MCAR_LR).reshape(7,6))
MCAR_LR_df.columns = ["test_mean","test_LR","test_RFR","test_LGBM","test_KNN","test_soft"]
MCAR_LR_df.index = ["train_listwise","train_mean","train_LR","train_RFR","train_LGBM","train_KNN","train_soft"]

MAR_LR_df = pd.DataFrame(np.array(MAR_LR).reshape(7,6))
MAR_LR_df.columns = ["test_mean","test_LR","test_RFR","test_LGBM","test_KNN","test_soft"]
MAR_LR_df.index = ["train_listwise","train_mean","train_LR","train_RFR","train_LGBM","train_KNN","train_soft"]

MNAR_LR_df = pd.DataFrame(np.array(MNAR_LR).reshape(7,6))
MNAR_LR_df.columns = ["test_mean","test_LR","test_RFR","test_LGBM","test_KNN","test_soft"]
MNAR_LR_df.index = ["train_listwise","train_mean","train_LR","train_RFR","train_LGBM","train_KNN","train_soft"]

続いて、ツリーモデル(ランダムフォレスト)を使用して同様の手順を繰り返します。

rfr = RandomForestRegressor(random_state=1)

def predict(train_list, test_list, model):
  mse_list = []
  for i,train in enumerate(train_list):
    train_X = train.drop("target",axis=1)
    target = train["target"]
    model.fit(train_X,target)
    for j,test in enumerate(test_list):
      test_X = test.drop("target",axis=1)
      true = test["target"]
      pred = model.predict(test_X)
      metric = mean_squared_error(true, pred)
      mse_list.append(metric)
  return mse_list

MCAR_RFR = predict(mcar_train_list, mcar_test_list, rfr)
MAR_RFR = predict(mar_train_list, mar_test_list, rfr)
MNAR_RFR = predict(mnar_train_list, mnar_test_list, rfr)
MCAR_RFR_df = pd.DataFrame(np.array(MCAR_RFR).reshape(7,6))
MCAR_RFR_df.columns = ["test_mean","test_LR","test_RFR","test_LGBM","test_KNN","test_soft"]
MCAR_RFR_df.index = ["train_listwise","train_mean","train_LR","train_RFR","train_LGBM","train_KNN","train_soft"]

MAR_RFR_df = pd.DataFrame(np.array(MAR_RFR).reshape(7,6))
MAR_RFR_df.columns = ["test_mean","test_LR","test_RFR","test_LGBM","test_KNN","test_soft"]
MAR_RFR_df.index = ["train_listwise","train_mean","train_LR","train_RFR","train_LGBM","train_KNN","train_soft"]

MNAR_RFR_df = pd.DataFrame(np.array(MNAR_RFR).reshape(7,6))
MNAR_RFR_df.columns = ["test_mean","test_LR","test_RFR","test_LGBM","test_KNN","test_soft"]
MNAR_RFR_df.index = ["train_listwise","train_mean","train_LR","train_RFR","train_LGBM","train_KNN","train_soft"]

上記までの処理で、精度検証のデータがすべて整いました。

3-2 精度の可視化と評価

最後に、精度を可視化・分析するためのコードを記載していきます。

def bar_data(df1,df2):
  df_list=[df1, df2]
  result = pd.DataFrame()
  for i,df in enumerate(df_list):
    s_df = df.stack().reset_index()
    s_df.columns = ["b","c","mse"]
    if i == 0:
      s_df["a"] = "Model_LR : "
    else :
      s_df["a"] = "Model_RFR : "
    s_df["predict_model"] = s_df["a"]+":"+s_df["b"]+":"+s_df["c"]
    s_df = s_df[["predict_model","mse"]]
    result = pd.concat([result, s_df]).sort_values(by="mse")
  return result

MCAR_result = bar_data(MCAR_LR_df,MCAR_RFR_df)
MAR_result = bar_data(MAR_LR_df,MAR_RFR_df)
MNAR_result = bar_data(MNAR_LR_df,MNAR_RFR_df)
def barplot(df,title):
  plt.figure(figsize=(5,30),facecolor="white")
  plt.title(title)
  sns.barplot(data = df, x="mse", y="predict_model")
  plt.show();


def heatmap(df,title):
  plt.figure(figsize=(10,7),facecolor="white")
  plt.title(title)
  sns.heatmap(df,annot=True,fmt="1.1f",cmap="coolwarm",)
  plt.show();

MCAR、MAR、MNARごとの精度ランキング

MCARの精度は、trainデータをsoftimputeで補完し、LRやLGBMなど、機械学習系のモデルでtestデータを補完したデータセットが精度上位にきています。
testデータをmeanで補完したモデルは総じて成績が悪く、安易に使用しないほうがよさそうです。

また、targetの予測モデルはLRよりもRFRが優秀です。

barplot(MCAR_result,"MCAR_result")

f:id:kopaprin:20190808020036p:plain

MARのデータセットではtrainデータをリストワイズ(row)で除去したモデルが優秀です。 欠損が別column(ここでは AGE )と相関しているため、無理に補完するよりも、その行そのものがtrainデータに混ざらない方が汎化性能が高まるようです。

こちらも testデータをmeanで置き換えるのは成績が悪く、最終のtarget予測モデルはRFRが高精度です。

barplot(MAR_result,"MAR_result")

f:id:kopaprin:20190808020332p:plain

MNARは特徴的で、trainデータをmeanで置き換えるモデルが比較的上位にきています。 (統計学の世界では禁則に近いtrain, testともにmean置換が精度で3位!)

平均値置換でなぜ精度が上がるのか考察が深まっていませんが、様々はパターンで検証の余地がありそうです。

barplot(MNAR_result,"MNAR_result")

f:id:kopaprin:20190808020435p:plain

LR, RFRのヒートマップと精度差分

ここでは、精度のヒートマップを用いて傾向を可視化しています。 また、LRとRFRの精度を比較し、相対的にどちらのモデルのほうが精度が高いか可視化しています。
上段がLRで予測したモデルの精度、中断がRFRで予測したモデルの精度、下段が両モデルの差分となっています。

MCARでは、3段目がほとんど赤いヒートマップとなっているため、RFRのほうが各パターンにおける精度が高いようです。
LRではtestデータをLRで予測したモデルの精度が高い傾向です。
RFRではtrainデータをsoftimputeで予測したモデルの精度が高い傾向です。

heatmap(MCAR_LR_df,"MCAR_LR: heatmap")
heatmap(MCAR_RFR_df,"MCAR_RFR: heatmap")
heatmap(MCAR_LR_df - MCAR_RFR_df,"MCAR : LR - RFR : heatmap")

f:id:kopaprin:20190808020600p:plain f:id:kopaprin:20190808020631p:plain f:id:kopaprin:20190808020657p:plain

MARでは、RFRで予測したモデルのうち、trainをリストワイズ除去したモデルの精度がやはり高いことがわかります。
LRの予測モデルでは、train_meanやtrain_softが優秀です。

heatmap(MAR_LR_df,"MAR_LR: heatmap")
heatmap(MAR_RFR_df,"MAR_RFR: heatmap")
heatmap(MAR_LR_df - MAR_RFR_df,"MAR : LR - RFR : heatmap")

f:id:kopaprin:20190808020751p:plain f:id:kopaprin:20190808020818p:plain f:id:kopaprin:20190808020844p:plain

MNARは、LRモデルの予測値はどれも安定している一方、RFRは精度に大きなばらつきがあります。
バイアスにより、一部モデルでは過学習が起こっている可能性が推察できます。

ただし、testデータの補完方法によっては、かなり高い精度で推計できていることもわかります。

heatmap(MNAR_LR_df,"MNAR_LR: heatmap")
heatmap(MNAR_RFR_df,"MNAR_RFR: heatmap")
heatmap(MNAR_LR_df - MNAR_RFR_df,"MNAR : LR - RFR : heatmap")

f:id:kopaprin:20190808020945p:plain f:id:kopaprin:20190808021010p:plain f:id:kopaprin:20190808021033p:plain

4 考察

今回のシミュレーションにより、様々な欠損補完方法と、その精度検証を行うことができました。


総括としては以下のとおりです。

  • 欠損の発生種別により、補完方法、targetの推計方法で大きなばらつきがあり、「これをやれば正解」というモデルは無さそうである
  • 禁則とされていた平均値置換も、場合により高い精度を得る可能性がある

また、課題として以下のようなことがあげられます。

  • 今回は総500レコードの小規模なデータセットであり、大規模データでも同様の事象になるとは限らない
  • 欠損値補完に使用したモデルや、target予測に使用したモデル、いずれもパラメーター調整を行っていない
    パラメーター調整を行えば、精度は大きく変化する可能性がある

以上、結論としては、機械学習タスクによる欠損値補完処理は、いくつかのパターンを試行・検証することが最適なのではないかと考えられます。

Google Colaboratoryのスクラッチコードセルが大変便利だったので是非皆さんに使っていただきたい件

つい、

print(hogehoge)

とか

df.head()

とか

df.shape

とかで汚らしくしてしまう迷える子羊たちは、 ぜひGoogle Colaboratoryのスクラッチコードセルを使っていただきたい。 というお話でした。

※ winショートカットはCtrl + alt + N

colab.gif

Android版OneNoteでWindows10の付箋が同期できるぞ

こちらTwitterでも書き込んだのですが、AndroidOneNoteで「付箋」なるものがリリースされてました。

さてさて、実際にPCで同期されるか試してみます。

Windows10の付箋ってどこにあんの?

ここにありました↓

f:id:kopaprin:20190117004922p:plain

Sticky Notesがアプリ名ですね。
「付箋」で検索してもヒットしないので注意。

起動していろいろ作ってみるとこんな感じ。 f:id:kopaprin:20190117005220p:plain

AndroidOneNoteでもきちんと共有されます。
f:id:kopaprin:20190117005626p:plain

付箋ユーザーには便利なのかな?

2019年の抱負をば。

皆様、新年あけましておめでとうございます。

f:id:kopaprin:20190101235143j:plain

昨年は、このブログを開設したり、MVP受賞したり、転職したり、引越ししたり、激動の1年となりました。

昨年をもろもろ振り返り、現在のスペックを記録し、今年1年の抱負でも書き残しておこうと思います。

1年後の自分、がんばれよ。

【2018年の成果と挫折、スペック】

  • Microsoft MVP受賞
    よくやった俺氏。しかし家庭内事情で休日活動がだいぶ制限されているぞ。今後どうやって活動しよう?

  • kaggle expart未達成
    そもそもsubmit全然してないという。ここ数か月はいい調子で取り組めているので、あれこれ考えずまずは取り組むべし。

  • TOEIC 700点
    そもそも受験も勉強もしてねぇ。やってるのはPodcastリスニングだけ。読み・書き絶望的。

  • Tableau
    本業でお世話になることになったTableau。とりあえずexpart名乗って恥じないレベルにしないとMSMVP for DataPlatform受賞歴の名が泣く。
    ってもBIは面白いのでたぶん大丈夫。

  • OneNote
    あれ、このブログってOneNote情報メインじゃなかったっけ?いつから中の人のポエムが中心になったのか。

【2019年の抱負】

今年のことば守破離

基礎力を固めることを中心にしなければならないことが多そうな1年なので、守破離が私の基本指針になりそうです。

  • kaggle expart の取得
    ここは昨年以上に精力的に取り組んでいきたいところ。社内kaggle部も作ったし、下地はできつつある。

  • Tableau pro
    なにか明確な目標があるわけでもないが、Tableauを使いこなす技術の習得。

  • Power BI pro
    MVPの継続というほうが分かりやすいが、諸事情で活動しづらくなってしまったのと、そもそもMVPが目的ではなくて、BIを楽しむことが目的なので、引き続きPower BI 技術を磨いて発信。

  • 書籍執筆
    これは今進行中プロジェクトをしっかりやれという自戒。

  • 統計検定1級取得
    TOEICは一旦捨て置き、こちらに焦点を当てる1年にします。解析学線形代数学、数理統計学の沼に飛び込む。

おわりに

33歳になりました。

基本をおろそかにしておくと、口だけ薄っぺらアナリストの出来上がりになりそうで怖いです。

そのための守破離の「守」が今年の姿勢です。きっちり基礎を叩き上げたいところ。

特に数学面、プログラミング技術面あたりの知識・技術を確立させたい所存。(本当は英語も叩き込みたいところであるが・・・)

2020年の抱負を書く時は、多少の自信と実績を引っ提げて読み返せると良いですなぁ。

使わないと時代遅れ?BIが実現するデータ分析の高度化・民主化・効率化の話

BI(ビジネスインテリジェンス)ツールということば、
かなり広まってきたような気がします。

一方で、
「BIツールなんて難しそう」
「普段の仕事に高度なツールなんて不要」
「そもそもメリットが分からない」
という方々も少なくないのではないかと。

本記事は、Microsoft Power BI Advent Calendar 2018として、
Microsoftが提供しているBIツール「Power BI」を例にとり、
そのメリットや課題について、全くBIを使ったことがない人にもわかるようにお伝えしていこうと思います。


データ分析の高度化・民主化・効率化

BIを活用する理由はこの3点に集約されるというのが私の考えです。
これらの言葉を聞いて少しでも関心を持ったら少しでもBIの利用・運用について試してみて損はないと思います。

以下、具体的にどのようなことがあげられるか見ていきましょう。

データ分析の高度化

グラフの可視化、ピボットテーブルなど、初歩的なデータ分析においてExcelが最もビジネスにおいて利用されてきた(されている)のではないでしょうか。
これについて、Excelが良い・悪いという話ではありません。
Excelの小回りの良さはビジネスのあらゆる場面で重宝しますし、多様な関数、VBAにより、柔軟に分析ニーズに応えることも可能です。

また、Excel標準機能では難しい統計解析・可視化についても、以下のような統計解析・可視化ツールも無償・有償で取り入れることができます。 以下ツールはいずれも大変優秀な機能を有しています。

appsource.microsoft.com

norimune.net

bellcurve.jp

このような応用範囲の広さと深さ、そして手軽さから、Excelは使い続けられてきました。

そのほか、データ分析の高度化という意味だけで言えば、統計解析を得意とするR言語や、機械学習ブームで人気急上昇中のPythonといったプログラミング言語を利用する方針もあります。

R: The R Project for Statistical Computing

Welcome to Python.org

このような様々な選択肢がある中、BIツールにおいても単なる可視化にとどまることなく、統計解析、機械学習といった解析的な技術が利用できるようになってきています。
Power BIでは以下のような機能や流れがあります。


Python、Rの統合

Power BI では、PythonやRのスクリプトを実行できるようになっています。
Pythonプレビュー機能

powerbi.microsoft.com

docs.microsoft.com

恥ずかしながら、Power BI × Python や Power BI × R などについて、以下QiitaやYoutubeにて解説していますので気になる方はぜひご覧ください。

qiita.com

www.youtube.com

Python・Rの統合により、Python・Rで提供されている高度な統計解析・機械学習・可視化ライブラリの恩恵を受けることができます。
さらに、「Rなんて書けないよ!」という方についても、Rの統計解析・可視化スクリプトを自動で実行してくれるカスタムビジュアルも複数提供されています。

docs.microsoft.com

BI上での統計解析・機械学習は日進月歩で進化しています。
そして、その流れをさらに加速させる動きがあります。
Power BIとMicrosftが提供するAI関連サービスとの統合です。


Azure Cognitive Servicesや Azure Machine Learningとの統合

2018年11月14日、Power BI Blogにて以下の案内がありました。

powerbi.microsoft.com

おそらく、多くの人が待ち望んでいるであろう、MicrosoftのBIツールとAIツールの統合を進めている趣旨の記事です。
すでに学習されたデータを手元のデータに適用できるAzure Cognitive Services、
手元のデータを学習させ、自身で様々な機械学習モデルが構築できるAzure Machine Learning、
これらがPower BIと統合すれば、BI環境はさらに高度化していくことは言うまでもありません。
それがほぼノンコーディングで実現できるのであれば、BIだけでなくAIの民主化にも寄与するものと思われます。

次は、データ分析の民主化の話です。


データ分析の民主化

データ分析の民主化。 煩わしい表現をしていますが、要はだれでも高度なデータ分析ができるようになることです。

インタラクティブなレポートの威力

BIは、データ分析結果の元ファイルを配布することを前提としていません。 分析結果は、レポートをWEB(クラウド)にアップロードしてステークホルダー間で共有することが基本です。
さらに、そのレポートはインタラクティブに動作し、共有者が自分の見たい視点で操作することができます。
これは以下サンプルから実際に体験してみましょう。

上記元データは、UCI Machine Learning RepositoryのAdultデータを利用しています。

このように、分析者が作ったレポートについて、分析・可視化などをしたことがないような人に対しても、 データドリブンな判断を促すことが可能になります。


WEB上のレポート、インタラクティブな操作が可能にするデータドリブンな判断

WEB上でインタラクティブなレポートを共有することについて、一番効果を発揮するのが、打合せや会議におけるブレーンストーミングなどではないでしょうか?

BIが日常の会議で活用されていない場合、以下のようになっていませんか?

f:id:kopaprin:20181201015830j:plain

データについて出てきた用意していない質問は、持ち帰りで分析、再度のレポート作成。
または、会議で言いっぱなしで特にその後検証することもなく終わる・・・。

これではせっかくの会議にかけた時間がもったいないですね。

BIを運用できるようになるとこう変わるかもしれません。

f:id:kopaprin:20181201020027j:plain

おそらく会議の成果やレポート作成にかかる手間の差は歴然ではないでしょうか? BI導入によるデータ分析環境の民主化は、データドリブンな判断を促すだけでなく、効率化にもつながります。

次は効率化の視点です。


データ分析の効率化

先ほどみたように、レポートがインタラクティブに動作することにより、レポートを何度も作成したり、Excelでグラフをつくり、PowerPointにまとめ・・・といったような作業の必要がほとんどなくなる可能性があります。

そのほかにも、BIにより効率化しそうなことは多々あります。

多様なデータソースから取得・自動更新

Power BIでのデータ分析の基本は、別に格納されているデータソースからデータを取得し、データモデルとして加工・可視化することです。

データソースは、Desktop上のcsvファイルからオンプレミスデータベース、クラウド、各種WEBサービスなど、多岐にわたります。

docs.microsoft.com

様々なデータソースから取得したデータは、一度宛先を指定すれば、自動的に更新が可能です。 ※更新の詳細は、以下リンク先を確認してください。

docs.microsoft.com

これにより、例えば月次、週次でまとめていたレポートなどを一切の手間なく、更新し続けることが可能になります。

ダッシュボードによる様々なレポートモニタリング

Power BI Serviceには、ダッシュボード機能があります。 これは、作成、WEBで共有したレポートの中から、自身に関係のあるビジュアル等を取捨選択し、1か所でまとめて管理できる機能です。 f:id:kopaprin:20181201022128p:plain

これにより、自身がチェックしたいデータやレポートだけを一元管理できる仕組みとなっています。 様々なデータの場所を覚え、いちいちチェックしに行く必要はありません。


ビッグデータの整形

複雑なビッグデータを可視化しやすく整形し、モデルとして取り込むのは非常に大きなリソースが必要なケースがありました。 それを簡素化し、クリック操作で分析しやすくするPower BI Dataflowsがリリースされました。(プレビュー)

docs.microsoft.com

これをうまく活用すれば、専門の技術者が何日も要していたデータモデルへの整形を、分析者が簡易に実行できるようになる可能性があります。


BIの運用に向けて

Power BIは無償にて始めることができますが、組織内共有など、ビジネスで本格的に扱うとなれば有料のサブスクリプションが必要になります。
他のBIツールについても、金額の大小はあれど、費用がかかるものがほとんどなはずです。
費用面は導入に向けて大きなハードルになると思いますが、そのほかにも乗り越えるべきハードルがあります。

BIは単なる可視化ツールではなく、インサイトを生むためのツール

BIツールは魅力的なツールであり、様々な可視化メニューが存在します。
一歩間違えると、「可視化」が目的になり、「インサイトを得る」ことが欠損してしまう可能性があります。
魅力的なビジュアルは重要ですが、5W1Hを踏まえたレポート設計が必要です。

  • Who(だれが見る・誰のためのレポートか)
  • When(いつ利用するレポートか)
  • Where(どこで利用されるレポートか)
  • What(なにを表現したいレポートか)
  • Why(なぜこのレポートを作るのか)
  • How(どのように使ってもらうレポートか)

これらを明確に応えられるレポートは、おそらくインサイトを生むことができるBIレポートです。

BI設計者の育成

設計者というと大げさですが、BIを導入し、活用していくためには既存情報システム部門だけでの取り組みでは不十分かもしれません。
BIが使われるかどうかは、さまざまなIT基板以上に経営・マーケティングなど、現場レベルで利用価値のある内容である必要があります。
現場の声をレポートに落とし込める設計者を育成することが欠かせません。

powerbi.microsoft.com


おわりに

昨今感じているBIを運用するうえで知っておいてほしいこと、ぜひ伝えたいことは記載できたつもりです。
しかし、Power BI をはじめ様々なBIツールの機能はこれだけでは語りつくせないことは承知しています。
この記事を読んで、少しでもBIに興味をもっていただければ本望です。