#!/usr/bin/env python
# coding: utf-8
#
#
# 井出 剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。 #
# # # ## 井出本の基本は、データに対するモデルを使って負の対数尤度を求め、それを異常検出関数として使うことです。 #
## M次元、N個の観測データ$\mathcal{D}$が以下の様な多次元正規分布に従っていると、仮定します。 # $$ # \mathcal{N}(x|\mu, \Sigma) = \frac{|\Sigma|^{-1/2}}{(2 \pi )^{M/2}} exp \left \{ -\frac{1}{2}(x - \mu)^T \Sigma^{-1} (x - \mu) \right \} # $$ #
## 対数尤度は、以下の様に計算されます。 # $$ # L(\mu, \Sigma | \mathcal{D}) = ln \prod_{n=1}^{N} \mathcal{N} (x^{(n)} | \mu, \Sigma) = \sum_{n=1}^N ln \mathcal{N} (x^{(n)} | \mu, \Sigma) # $$ #
## これに多次元正規分布$\mathcal{N}$の定義を代入すると、対数尤度Lは以下の様になります。 # $$ # L(\mu, \Sigma | \mathcal{D}) = -\frac{M N}{2}ln(2 \pi) - \frac{N}{2} ln | \Sigma | # -\frac{1}{2} \sum_{n=1}^N (x^{(n)} - \mu)^T \Sigma^{-1} (x^{(n)} - \mu) # $$ # #
## これを$\mu, \Sigma$で偏微分して、尤度Lを最大にする$\mu, \Sigma$を求めます。 #
## 最初にLを$\mu$で偏微分すると、以下の様になります。 # $$ # \frac{\partial L(\mu, \Sigma | \mathcal{D})}{\partial \mu} = \sum_{n=1}^N \Sigma^{-1} (x^{(n)} - \mu) # $$ #
## 上記の式が0になる$\hat{\mu}$は以下の様に求まります。 # $$ # \hat{\mu} = \frac{1}{N} \sum_{n=1}^N x^{(n)} # $$ #
## 次に、$\Sigma$での偏微分です。これには、以下の関係式を使います。 # $$ # -ln | \Sigma | = ln | \Sigma^{-1} | # $$ # と、 # $$ # (x^{(n)} - \mu)^T \Sigma^{-1} (x^{(n)} - \mu) = Tr \left \{ \Sigma^{-1} (x^{(n)} - \mu) (x^{(n)} - \mu)^T \right \} # $$ # と、以下の公式を使います。 # $$ # \frac{\partial}{\partial A} Tr(A B) = \frac{\partial}{\partial A} Tr(B A) = B^T # $$ # $$ # \frac{\partial}{\partial A} ln |A| = (A^{-1})^T # $$ #
## 偏微分の結果は、以下の様になります。共分散行列$\Sigma$が対象行列であるので、$\Sigma = \Sigma^T$を使っています。 # $$ # \begin{eqnarray} # \frac{\partial L(\mu, \Sigma | \mathcal{D})}{\partial (\Sigma^{-1})} & = & \frac{N}{2}ln | \Sigma^{-1} | # - \frac{1}{2} \Sigma^{-1} (x^{(n)} - \mu) (x^{(n)} - \mu)^T \\ # & = & \frac{N}{2} \Sigma^T - \frac{1}{2}\sum_{n=1}^N (x^{(n)} - \mu) (x^{(n)} - \mu)^T \\ # & = & \frac{N}{2} \Sigma - \frac{1}{2}\sum_{n=1}^N (x^{(n)} - \mu) (x^{(n)} - \mu)^T # \end{eqnarray} # $$ #
## この値が0になるような$\hat{\Sigma}$は、以下の様に求まります。 # $$ # \hat{\Sigma} = \frac{1}{N} \sum_{n=1}^N (x^{(n)} - \mu) (x^{(n)} - \mu)^T # $$ #
## 負の対数尤度を基に、異常度$a(x')$を以下の様に定義します。これは、マハラノビス距離と呼ばれ、「ばらつきが大きい方向の変動は大目にみる」効果があるそうです。 # $$ # a(x') = ({x'} - \mu)^T \Sigma^{-1} ({x'} - \mu) # $$ #
## しかし、実際のRでのプログラムをみると、ちょっと違っています。実行例2.5(Davisに対する異常度の計算部分)は、以下の計算をしています。 # $$ # a(x') = \sum_{m=1}^M ({x'} - \mu)^2 \Sigma^{-1} # $$ #
## 値としては、Wikipediaにある共分散行列が対角行列であるときの以下の式に近いと思われます。 # $$ # d(\vec{x}, \vec{y}) = \sqrt{\sum_{i=1}^P \frac{(x_i - y_i)^2}{\sigma_i^2}} # $$ #
# # ## 最初に必要なライブラリーやパッケージをロードしておきます。jupyter用に新しくなったRUtil.pyも使います。 #
# # In[1]: get_ipython().run_cell_magic('HTML', '', '\n') # In[2]: # RとPandasのデータフレームを相互に変換する関数を読み込む # Rの必要なライブラリ r('library(ggplot2)') r('library(jsonlite)') # python用のパッケージ import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns get_ipython().run_line_magic('matplotlib', 'inline') # jupyter用のdisplayメソッド from IPython.display import display, Latex, HTML, Math, JSON # sageユーティリティ load('script/sage_util.py') # Rユーティリティ load('script/RUtil.py') # ## Rのcarパッケージに入っているデータDavisを使って、多次元データの異常度を求めてみましょう。 #
# # In[3]: # 2.2 carパッケージのDavisを使って例題を試す # carパッケージのインストールには、時間が掛かりますので注意してください。 r('library(car)') # ## Rのデータをpandaのデータフレームに変換して、データの性格を調べます。 # 最初にinfoとdescribeで個数や欠損値、平均、分散などの統計情報を求めます。 #
# # In[5]: # Rのデータをpandaのデータフレームに変換する davis = RDf2PandaDf("Davis") # データの個数と欠損値の有無をみる davis.info() # In[7]: # データの情報を取り出す davis.describe() # ## 次にDavisの体重のヒストグラムと体重(weight)と身長(height)の散布図をプロットし、おおざっぱな性質を調べます。 #
# # # In[10]: # weightの分布をみる davis.weight.hist(bins=20) plt.show() # In[12]: # 身長(height)と体重(weight)の関係をみる sns.lmplot('weight', 'height', data=davis, fit_reg=False ) plt.show() # ## Sageとpandas(numpy)の機能を使って共分散を求めみましょう。 #
## 平均mxを求めてデータの平均からの差Xcを求めます。pandasではこれらをmeanとマイナス記号で処理することができます。 #
## 数式では共分散行列$\hat{\Sigma}$は、以下の様に定義されていますが、 # $$ # \hat{\mu} = \frac{1}{N}\sum_{n=1}^N x^{(n)}, # \hat{\Sigma} = \frac{1}{N} \sum_{n=1}^N (x^{(n)} - \hat{\mu})(x^{(n)} - \hat{\mu})^T # $$ # Xcが列ベクトルではなく、行ベクトルなので、以下の様に計算しています。 # $$ # \hat{\Sigma} = \frac{1}{N} Xc^T Xc # $$ #
# # In[13]: # pandaデータフレームで、weightとheightの平均からのずれXcを計算する X = davis[['weight', 'height']] mx = X.mean() Xc = X - mx # ## numpyの固有の理由だと思いますが、マトリックスの積に.dot()関数を使います。この辺はちょっと違和感があります。 #
## Pnadasの機能を使うと共分散行列は一発で計算できるのですが欠損値の影響なのか値が一致しません。 #
# # In[15]: # Pandaで標本共分散行列を求めるには print Xc.T.dot(Xc) /len(Xc) # もっと簡単なのはcov関数を使う(でも値が微妙に異なるのは、欠損値があるからか?) print Xc.cov() # ## Sageのmatrixを使った異常度の計算方法を以下に示します。 #
## 理論をベースとしているSageと実験値を処理するためのPandas, # numpyをどうやって融合するのかがSage使いの腕の見せ所です。 #
## と言ってもやっていることは簡単です。マトリックスの要素単位の処理をするには、 # Sageのmatrixにnumpy()メソッドで、numpyのマトリックスに変換するだけです。 #
## Rやnumpyでは多くのデータを一括して処理することが多く、要素単位の処理に四則演算子が # 使われますが、Sageではマトリックスやベクトルの定義に沿って計算が行われます。 #
# # In[16]: # Sageのmatrixを使った計算 # 中心化したデータをSageのマトリックスにする Xc =matrix((X - mx).values) # In[17]: # 標本共分散行列を求める Sx = Xc.T*Xc/ Xc.nrows(); Sx # In[19]: # (XcとSxの逆行列の積)の各要素とXcの各要素を掛け合わせる # ここでのトリックは、numpyの*が要素毎の積であることを利用 # numpy()関数を使ってマトリックスをnumpyのマトリックスに変換する a_prev = (Xc*Sx.inverse()).numpy() * Xc.numpy() # 行単位の和を求める a = [sum(x) for x in a_prev] # リストプロット list_plot(a, figsize=4) # ## pandaとnumpyの機能を使って異常値を計算します。ここでは、$\Sigma^{-1} Xc$がsolveの解であることを # 使っています。また、次数の和には、pythonのsumを使っています。 #
## 出力結果では、1つだけ突出した異常データが見つかります。 #
## 全体的な処理から以降ではPandaを使った分散行列の方式を使います。 #
# # In[21]: # Pandaのみを使った計算 # 中心化したデータの値をXcにセット Xc =(X - mx).values # 標本共分散行列を求める Sx = (X - mx).cov().values # 異常値aを求めるΣ^-1 * Xcは、solveの解を利用 a_prev = np.linalg.solve(Sx, Xc.T).T * Xc # 行単位の和を求める a = [sum(x) for x in a_prev] # リストプロット list_plot(a, figsize=4) # ## 異常値の検出では、全体に対する個々の観測値の異常を検出するのに使用しましたが、 # 個々の観測値の変数値の異常を検出するには、マハラノビス=タグチ法を使用します。 #
## M変数のなかからいくつかの変数を選び、1変数当たりの異常度を以下の様に検出します。 # # $$ # SN_q = -10 log_{10} \left \{ \frac{1}{N'} \sum_{n=1}^{N'} \frac{1}{a_q(x'^{(n)})/M_q} \right \} # $$ #
## マハラノビス=タグチ法の実行例2.6をSageを使って試してみましょう。 #
## MASSパッケージのroadデータには、アメリカ26州の交通事故死亡者数deaths、運転者数drivers、 # 人口密度popden、郊外地区の道路長rural、1月における1日の最高気温の平均値temp、1年度との燃料消費量 # fuelの6変数を記録しています。 #
## データの前段階でroadをdriversで割って対数化する必要があるのですが、pandaで処理する方法が見つからず、 # Rで前処理することにしました(残念)。 #
## 取り込んだデータは、infoで個数と欠損値の有無を調べ、内容を調べます。 # $rowという変な名前が存在するので、stateに置き替えます。 #
# # In[23]: # マハラノビス=タグチ法 # MASSパッケージのroadデータを使って計算 r('library(MASS)') # roadをPandaのデータフレームに変換 # road / road$driversとlogを使った対数化がPandaで処理する方法がみつからない r('X <- road / road$drivers') r('X <-log(X[, -2] + 1)') road = RDf2PandaDf("X") road.info() # In[24]: # _rowという変な名前が存在するので、headで内容をチェック、_rowが州の名前だと判明 road.head() # In[25]: # 州のカラムにstateをセットする road.columns = ['state', 'deaths', 'fuel', 'popden', 'rural', 'temp'] # Calif州のデータをピックアップする road[road.state == 'Calif'] # ## 異常度の計算と同様に異常度を求め、次元数で割って1変数当たりの異常度を求めます。 #
# # In[26]: # pandaデータフレームで、平均からのずれXcを計算する X = road[['deaths', 'fuel', 'popden', 'rural', 'temp']] mx = X.mean() Xc = X - mx Xc.shape[1] # In[27]: # Pandaのみを使った計算 # 中心化したデータをnumpyのマトリックスにする Xc =(X - mx).values # 標本共分散行列を求める Sx = (X - mx).cov().values # 異常値aを求めるΣ^-1 * Xcは、solveの解を利用 a_prev = np.linalg.solve(Sx, Xc.T).T * Xc # 行単位の和を求め、1変数当たりの異常度を計算する a = [sum(x)/Xc.shape[1] for x in a_prev] # リストプロット list_plot(a, figsize=4) # ## 最も異常度の大きい、CalifのデータをXcから取り出し、SN比を計算し、棒グラフで表示します。 #
## 表示順はdeaths fuel popden rural tempで、図2.7(b)と並びが異なりますが、 # 値は同じ結果になっています。 #
## この結果から、異常に寄与している変数は、通事故死亡者数deaths、1年度との燃料消費量 # fuel、人口密度popdenであることが分かります。 # この結果は、車が多く使用される過密な大都市で交通事故死亡者が増加するという一般的な # 常識と一致します。 #
# # In[28]: # XcからCalifのデータを取得 Xc_prime = Xc[4] # 要素毎の計算をしたいので、numpyの形式に変換 SN1 = 10*np.log10(Xc_prime^2/Sx.diagonal()) print SN1 # In[29]: # 左からdeaths fuel popden rural tempで、図2.7(b)と並びが異なる bar_chart(SN1, figsize=4) # In[ ]: