# 井出 剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。 #
# # # ## 井出本の基本は、データに対するモデルを使って負の対数尤度を求め、それを異常検出関数として使うことです。 #
# # # ## Davisの体重（weight）のヒストグラムを見ると、分布が左右対称ではないことに気づきます。 # 井出本では、ガンマ分布を使ってこの非対称な分布の異常度を検出しています。 #
## いつものように必要なパッケージを読込ます。 #
# # 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') # In[3]: # 2.2 carパッケージのDavisを使って例題を試す r('library(car)') r('library(MASS)') # ## Davisの体重のヒストグラムを再度プロットしてみます。 # 1個だけ160のところにありますが、2章ではこれが体重と身長を入れ間違えたものと結論づけています。 #
# # In[4]: # Rのデータをpandaのデータフレームに変換する davis = RDf2PandaDf("Davis") # weightの分布をみる # 1個だけ、160以上の値を示している davis.weight.hist(bins=20) plt.show() # ## 2章と同様にガンマ分布の対数尤度を求めてみます。 #
## 井出本では、ガンマ分布を以下の様に定義しています。 # $$ # \mathcal{G} (x | k, s) = \frac{1}{s \Gamma(k)}(\frac{x}{s})^{k-1} exp(- \frac{x}{s}) # $$ #
## これから対数尤度Lは、以下の様になります。 # $$ # L(k, s | \mathcal{D}) = \sum_{n=1}^N \left [ -ln\{s \Gamma(k) \} + (k-1) ln\frac{x^{(n)}}{s} - \frac{x^{(n)}}{s} \right ] # $$ #
## この式をパラメータkとsで微分し、０と等しいとします。 #
## $-ln\{s \Gamma(k) \}$が$-ln(s) - ln(\Gamma(k))$となりますから、sの偏微分は以下の様になります。 # # $$ # \begin{eqnarray} # \frac{\partial L}{\partial s} & = & \sum_{n=1}^N \left [ -\frac{1}{s} + (k-1)\frac{\partial}{\partial s} ln \frac{x^{(n)}}{s} + \frac{x^{(n)}}{s^2} \right ] \\ # & = & \sum_{n=1}^N \left [ -\frac{1}{s} - \frac{(k-1)}{s} + \frac{x^{(n)}}{s^2} \right ] \\ # & = & \sum_{n=1}^N \left [ -\frac{k}{s} + \frac{x^{(n)}}{s^2} \right ] # \end{eqnarray} # $$ # この式が0となることから、 # $$ # \hat{s} = \frac{1}{\hat{k} N} \sum_{n=1}^N x^{(n)} = \frac{\hat{\mu}}{\hat{k}} # $$ #
## 同様にkの偏微分は以下の様になります。 # $$ # \frac{\partial L}{\partial k} = \sum_{n=1}^N \left [ -\frac{1}{\Gamma(k)}\frac{\partial \Gamma(k)}{\partial k} + ln \frac{x^{(n)}}{s}\right ] # $$ # しかし、kがガンマ関数が入っているので、簡単には$\hat{k}$は求まりません。 #
# # ## 井出本では、最尤推定でk, sを求める代わりにモーメント法を使ってkとsを推定しています。 # # 1次のモーメントは以下の様に計算されます。$\Gamma(k+1) = k \Gamma(k)$を関係を使っています。(以下の式は未フォロー) # # $$ # \lt x \gt = \int_0^{\infty} dx \, x \, \mathcal{G}(x | k,s) = s\frac{\Gamma(k+1)}{\Gamma(k)} \int_0^{\infty} dx \, x \, \mathcal{G}(x | k+1,s) = ks # $$ # 同様に２次のモーメントは、以下の様になります。 # # $$ # \lt x^2 \gt = k(k+1)s^2 # $$ #
## これとデータから求められる１次の２次モーメントの関係 # # $$ # \lt x \gt = \frac{1}{N} \sum_{n=1}^N x^{(n)} # $$ # # $$ # \lt x^2 \gt = \frac{1}{N} \sum_{n=1}^N {x^{(n)}}^2 # $$ #
## この関係から # $$ # \hat{k}_{mo} = \frac{\hat{\mu}^2}{\hat{\sigma}^2}, \hat{s}_{mo} = \frac{\hat{\sigma}^2}{\hat{\mu}} # $$ #
# # ## davisのweight(体重）でモーメント法を使ってk, sを求めてみます。 #
# # # In[5]: # ガンマ分布の異常検出 N = len(davis.weight) mu = davis.mean()['weight'] si = davis.std()['weight']*(N-1)/N kmo = (mu/si)^2 smo = si^2/mu print kmo, smo # ## Rのfitdistrを使ってdavisのweight(体重）のガンマ分布を推定してみます。 # shapeがkの値、rateが1/sの値として返されます。 #
## モーメント法とfitdistrの結果は、そこそこ近い結果となっています。 #
# # # In[6]: r('ml <- fitdistr(Davis$weight, "gamma")') # In[7]: kml = r('ml$estimate["shape"]') sml = 1/r('ml$estimate["rate"]') print kml print sml # ## このノートの特徴は、Sageとpythonを使って井出本と同じ結果を求めることにあります。 #
## scipyのstatsのgamma.fitを使って、kとsを求めてみます。 # gamma.fitでは、floc=0を指定しないとfitdistrと同じ結果になりません。 #
# # In[8]: # SciPyのstatsを使ってガンマ分布の最尤推定値を求める、floc=0を指定するとfitdistrと同じ結果になる from scipy import stats fit_alpha, fit_loc, fit_beta = stats.gamma.fit(np.array(davis.weight), floc=0) print(fit_alpha, fit_loc, fit_beta) # ## Sage6.7からヒストグラムがサポートされました。 #
## 通常のヒストグラムと正規化（normed=True）したヒストグラムが表示でき、 # 他のプロットデータと重ね合わせることができます。 #
# # In[9]: # Pandasのデータフレームからweightを取り出す weight = list(davis.weight.values) # In[10]: # Sage6.7で導入されたhistogramを使用、binsのデフォルトは10 histogram(weight, figsize=4, bins=24) # In[11]: # 正規化したグラフも可能 hg = histogram(weight, figsize=4, bins=24, normed=True) hg.show() # In[12]: def _gamma(x, k, s): return 1/(s*gamma(k))*(x/s)^(k-1)*exp(-x/s) # In[13]: gp = plot(_gamma(x, 22.48548, 2.926333), [x, 0, 180], rgbcolor='red') (gp+hg).show() # ## 以下の様に正規分布関数_gaussを定義して、ヒストグラム、正規分布、ガンマ分布の推定結果を # 重ね合わせてみます。 #
# # In[14]: def _gauss(x, mu, sig2): return 1/sqrt(2*pi*sig2)*exp(-1/(2*sig2)*(x - mu)^2) # In[15]: # 正規分布として推定した結果（2章の結果を利用） sig2 = 226.7 mu = 65.8 ga = plot(_gauss(x, mu, sig2), [x, 0, 180], rgbcolor='green') (hg+gp+ga).show(figsize=4) # ## 最初、flocを指定しないでフィッティングをしたため、fitdistrと結果が異なり # 悩んでいました。しかしfloc=0とすることで、fitdistrと同じ結果が得られることが # 分かり、floc制約をしない場合、どのような推定結果がでるのか調べて見ました。 #
## 今回、可視化はmatplotを使いました。floc=0の設定をしない方がより良くフィット # していることが分かります。 #
# # In[16]: # locも含めた分布のフィッティング結果 samp = np.array(davis.weight) fit_alpha, fit_loc, fit_beta = stats.gamma.fit(samp) print fit_alpha, fit_loc, fit_beta x = np.linspace(20,180,200) plt.figure() # fitted distribution pdf_fitted = stats.gamma.pdf(x, fit_alpha, loc=fit_loc, scale=fit_beta) plt.title('Gamma distribution') plt.plot(x,pdf_fitted,'r-') plt.hist(samp,bins=24, normed=1,alpha=.3) plt.show() # In[17]: # fitdistrによるフィッティング結果 fit_alpha = 22.48548 fit_loc = 0 fit_beta = 2.926333 x = np.linspace(20,180,200) plt.figure() # fitted distribution pdf_fitted = stats.gamma.pdf(x, fit_alpha, loc=fit_loc, scale=fit_beta) plt.title('Gamma distribution') plt.plot(x,pdf_fitted,'r-') plt.hist(samp,bins=24, normed=1,alpha=.3) plt.show() # ## データに異常データが混ざっている場合の例として、以下の様に2種類の正規分布が混ざっている場合に # ついて考えてみます。 #
# # In[18]: # 訓練データに異常データがまざっている場合（2個の正規分布の混合分布を考える） x = var('x') N = 1000 mu0 = 3 mu1 = 0 sig0 = 0.5 sig1 = 3 pi0 = 0.6 pi1 = 0.4 # 分布を表示 g0 = plot(_gauss(x, mu0, sig0^2), [x, -5, 6], rgbcolor='green') g1 = plot(_gauss(x, mu1, sig1^2), [x, -5, 6], rgbcolor='blue') g01 = plot(_gauss(x, mu0, sig0^2)+ _gauss(x, mu1, sig1^2), [x, -5, 6], rgbcolor='red') (g0 + g1 + g01).show(figsize=4) # ## 分布は以下の様な混合正規分布として表されます。 # $$ # p(x) = \pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) + \pi_1 \mathcal{N}(x | \mu_1, \sigma_1^2) # $$ #
## パラメータを$ \theta = (\pi_0, \mu_0, \sigma_0^2, \pi_1, \mu_1, \sigma_1^2)$にまとめると、 # 対数尤度は、以下の様になります。 # $$ # L(\theta, \mathcal{D}) = \sum_{n=1}^N ln \left \{ \pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) + \pi_1 \mathcal{N}(x | \mu_1, \sigma_1^2) \right \} # $$ #
## これを$\mu_0$で偏微分すると、以下の様になります。 # $$ # 0 = \frac{\partial L}{\partial \mu_0} = \sum_{n=1}^N \frac{\pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) }{\pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) + \pi_1 \mathcal{N}(x | \mu_1, \sigma_1^2)} \frac{(x^{(n)} - \mu_0)}{\sigma_0^2} # $$ # 同様に$\sigma_0^2$で偏微分すると、以下の様になります。 # $$ # 0 = \frac{\partial L}{\partial \sigma_0^2} = \sum_{n=1}^N \frac{\pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) }{\pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) + \pi_1 \mathcal{N}(x | \mu_1, \sigma_1^2)} \frac{ \{-(x^{(n)} - \mu_0)^2 + \sigma_0^2 \}}{2} # $$ # ここで、データが$\pi_i$のどちらの集合に属するかその期待値を帰属度として以下の様に定義します。 # $$ # q_i^{(n)} = \frac{\pi_i \mathcal{N}(x | \mu_i, \sigma_i^2) }{\pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) + \pi_1 \mathcal{N}(x | \mu_1, \sigma_1^2)} # $$ #
## $\mu_i, \sigma_i^2$の推定値は、以下の様に求まります。 # $$ # \hat{\mu}_i = \frac{\sum_{n=1}^N q_i^{(n)} x^{(n)}}{\sum_{n'=1}^N q_i^{(n')}} # $$ # $$ # \hat{\sigma}_i^2 = \frac{\sum_{n=1}^N q_i^{(n)} (x^{(n)} - \mu_i)^2}{\sum_{n'=1}^N q_i^{(n')}} # $$ #
## また、$\pi_i$は、以下の様になります。 # $$ # \hat{\pi}_i = \frac{1}{N} \sum_{n=1}^N q_i^{(n)} # $$ #
# # ## 混合正規分布をEMアルゴリズムを使って求めます。 #
## EMアルゴリズムでは、以下の手順でパラメータの値を求めます。 #
# 以下の2つの正規分布をもつテストデータ1000個を生成します。 # $$ # (\pi_0, \pi_1) = (0.6, 0.4), (\mu_0, \sigma_0) = (3, 0.5), (\mu_1, \sigma_1) = (0, 3) # $$ #
# # In[19]: # テストデータの生成 T0 = RealDistribution('gaussian', sig0, seed=1) d0 = [T0.get_random_element()+mu0 for i in range(pi0*N)] h0 = histogram(d0, bins=20, normed=True, color='green') T1 = RealDistribution('gaussian', sig1, seed=1) d1 = [T1.get_random_element()+mu1 for i in range(pi1*N)] h1 = histogram(d1, bins=20, normed=True, color='blue') # プロットして確認 (h0+h1+g0+g1+g01).show(figsize=4) # ## 各パラメータの初期値を以下の様に設定します。 # $$ # (\pi_0, \pi_1) = (0.5, 0.5), (\mu_0, \sigma_0) = (5, 1.0), (\mu_1, \sigma_1) = (-5, 5.0) # $$ #
# # In[20]: # データを結合 d = d0 + d1 # pi, sig, muの初期値 pi0 = 0.5 pi1 = 0.5 mu0 = 5.0 mu1 = -5.0 sig0 = 1.0 sig1 = 5.0 # ## EMアルゴリズムを10回繰り返します。 #
## 高々10回の繰り返しで、とても精度良くパラメータが推定されています。 #
# # In[21]: # EMアルゴリズムを10回繰り返す x = np.array(d) for i in range(10): piN0 = pi0*stats.norm.pdf(x, mu0, sig0); piN1 = pi1*stats.norm.pdf(x, mu1, sig1) qn0 = piN0/(piN0 + piN1); qn1 = piN1/(piN0 + piN1) pi0 = sum(qn0)/N; pi1 = sum(qn1)/N mu0 = sum(qn0*x)/(N*pi0); mu1 = sum(qn1*x)/(N*pi1) sig0 = sqrt(sum(qn0*(x - mu0)*(x - mu0))/(N*pi0)) sig1 = sqrt(sum(qn1*(x - mu1)*(x - mu1))/(N*pi1)) # In[22]: # 結果出力 print pi0, pi1 print mu0, mu1 print sig0, sig1 # In[23]: x = var('x') g01 = plot((_gauss(x, mu0, sig0^2)+ _gauss(x, mu1, sig1^2))/2, [x, -9, 9], rgbcolor='red') hd = histogram(d, bins=50, normed=True) (hd + g01).show(figsize=4) # In[24]: # 帰属度qn0, qn1の分布 q0 = list_plot(list(qn0), rgbcolor='red') q1 = list_plot(list(qn1)) (q0 + q1).show(figsize=4) # ## 最後にScikit-learnの混合正規分布パッケージsklearn.mixtureを使って、 # davisの体重・身長のデータに対して、混合正規分布を求めてみます。 #
## GMMは、異常なデータに影響を受けるため、12番目のデータを除いた集合で、 # 混合正規分布を求めます。 #
## GMMの引数のn_componentsが求める分布の数です。RのMclustを使うとベストな分布数も計算できますが、 # ここでは2として計算します。 #
# # In[25]: # 混合正規分布 # SageでのGMMは、http://www15191ue.sakura.ne.jp:8000/home/pub/57/ を参照してください。 from sklearn import mixture # 分類器の生成 classifier = mixture.GMM(n_components=2, covariance_type='full') # In[26]: # GMMを求める # 12番目のデータを除く（体重と身長を間違えたデータ） davis_e11 = davis[davis.index != 11] data = davis_e11[['weight', 'height']] # 混合正規分布を求める classifier.fit(data) # ## 大切なのは、GMMによる分類ではなく、その結果を可視化して特徴を把握することです。 #
## GMMによってデータがどのグループに属するか分類した結果は、predict関数で取得し、 # davis_e11のpredにセットします。 #
## 次にseabornでpredで色を変えて分布をプロットします。 # たったこれだけで、GMMの結果を確認することができます。 #
# # In[28]: # 分類結果をpredにセット # davis_e11['pred']= classifier.predict(data)で、Try using .loc[row_indexer,col_indexer] = value insteadのエラー # 以下のように変更してもエラーがでます。 davis_e11.loc[:, 'pred'] = classifier.predict(data) # In[29]: # 分類結果の表示 # 身長（height）と体重（weight）の関係をみる sns.lmplot('weight', 'height', data=davis_e11, hue='pred', fit_reg=False ) plt.show() # ## GMMの結果は、$\mu$がmeans_、$\sigma$がcovars_、$\pi$がweights_変数にセットされています。 #
## Sageのcontour_plot関数を使って正規分布をプロットしてみます。 #
# # In[30]: # クラス分けされたデータ cls0 = davis_e11[davis_e11.pred == 0][['weight', 'height']].values cls1 = davis_e11[davis_e11.pred == 1][['weight', 'height']].values # In[31]: # 正規分布の推定値 mu0 = vector(classifier.means_[0]) mu1 = vector(classifier.means_[1]) sig0 = matrix(classifier._get_covars()[0]) sig1 = matrix(classifier._get_covars()[1]) show(mu0); show(mu1) show(sig0); show(sig1) # In[32]: # ベクトルのガウス分布を定義 def _gauss(v, mu, sigma): d = len(v); sigma_inv = sigma.inverse(); sigma_abs_sqrt = sigma.det().sqrt(); val = -(v - mu) * sigma_inv * (v - mu).column()/2; a = (2*pi)**(-d/2) * sigma_abs_sqrt**-0.5 return a * e**val[0]; # In[33]: x, y = var('x y') # データのプロット pt_plt = Graphics() for pt in cls1: pt_plt += point(pt, rgbcolor='red', zorder=2) for pt in cls0: pt_plt += point(pt, rgbcolor='blue', zorder=2) # 正規分布を表示 r_cnt_plt = contour_plot(lambda x, y : _gauss(vector([x, y]), mu1, sig1), [x, 20, 130], [y, 140, 200], contours = 1, cmap=['red'], fill=False) b_cnt_plt = contour_plot(lambda x, y : _gauss(vector([x, y]), mu0, sig0), [x, 20, 130], [y, 140, 200], contours = 1, cmap=['blue'], fill=False) (pt_plt + r_cnt_plt + b_cnt_plt).show(figsize=4) # ## 混合正規分布の異常度は、以下の様になります。 # $$ # a(x') = -ln \left \{ \sum_{k=1}^K \hat{\pi}_k \mathcal{N}(x' | \hat{\mu}_k, \hat{\Sigma}_k) \right \} # $$ #
# # In[34]: #異常度の計算 # pi_nを求める pi_vec = vector(classifier.weights_); show(pi_vec) # davis全データに対する事後分布を求める a = [- log(pi_vec[0]*_gauss(vector(x), mu0, sig0)+pi_vec[0]*_gauss(vector(x), mu1, sig1)).n() for x in davis[['weight', 'height']].values ] # In[35]: list_plot(a, figsize=4) # In[ ]: