![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
2015/08/16からのアクセス回数 6947
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www15191ue.sakura.ne.jp:8000/home/pub/58/
また、私の公開しているSageのサーバ（http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。
井出 剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。
井出本では、Rを使って例題を実行していますが、ここではSage, Pythonをベースに例題を試します。
井出本の基本は、データに対するモデルを使って負の対数尤度を求め、それを異常検出関数として使うことです。
Davisの体重（weight）のヒストグラムを見ると、分布が左右対称ではないことに気づきます。 井出本では、ガンマ分布を使ってこの非対称な分布の異常度を検出しています。
いつものように必要なパッケージを読込ます。
sageへの入力:
# RとPandasで必要なユーティリティ
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# RUtilにRとPandasのデータフレームを相互に変換する関数を追加+GgSaveを追加（2015/07/12）
load(DATA + 'RUtil.py')
sageへの入力:
# 2.2 carパッケージのDavisを使って例題を試す
r('library(car)')
r('library(MASS)')
[1] "MASS" "car" "jsonlite" "ggplot2" "stats" "graphics" "grDevices" "utils" [9] "datasets" "methods" "base"
Davisの体重のヒストグラムを再度プロットしてみます。 1個だけ160のところにありますが、2章ではこれが体重と身長を入れ間違えたものと結論づけています。
sageへの入力:
# Rのデータをpandaのデータフレームに変換する
davis = RDf2PandaDf("Davis")
# weightの分布をみる
# 1個だけ、160以上の値を示している
ggplot(davis, aes(x='weight')) + geom_histogram(binwidth=5.0,color="grey")
<ggplot: (8738330208357)>
sageへの入力:
GgSave('fig_2.1.png', fac=0.8)
Saving 11.0 x 8.0 in image.
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を求めてみます。
sageへの入力:
# ガンマ分布の異常検出 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
19.1928236809 3.42836474164
Rのfitdistrを使ってdavisのweight(体重）のガンマ分布を推定してみます。 shapeがkの値、rateが1/sの値として返されます。
モーメント法とfitdistrの結果は、そこそこ近い結果となっています。
sageへの入力:
r('ml <- fitdistr(Davis$weight, "gamma")')
shape rate 22.4854793 0.3417247 ( 2.2317648) ( 0.0342978)
sageへの入力:
kml = r('ml$estimate["shape"]')
sml = 1/r('ml$estimate["rate"]')
print kml
print sml
shape
22.48548
rate
2.926333
このノートの特徴は、Sageとpythonを使って井出本と同じ結果を求めることにあります。
scipyのstatsのgamma.fitを使って、kとsを求めてみます。 gamma.fitでは、floc=0を指定しないとfitdistrと同じ結果になりません。
sageへの入力:
# 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)
(22.485303592007426, 0, 2.9263558630975797)
今回バージョンアップしたSage6.7からヒストグラムがサポートされました。
通常のヒストグラムと正規化（normed=True）したヒストグラムが表示でき、 他のプロットデータと重ね合わせることができます。
sageへの入力:
# Pandasのデータフレームからweightを取り出す weight = list(davis.weight.values)
sageへの入力:
# Sage6.7で導入されたhistogramを使用、binsのデフォルトは10 histogram(weight, figsize=5, bins=24)
sageへの入力:
# 正規化したグラフも可能 hg = histogram(weight, figsize=5, bins=24, normed=True) hg.show()
sageへの入力:
def _gamma(x, k, s):
return 1/(s*gamma(k))*(x/s)^(k-1)*exp(-x/s)
sageへの入力:
gp = plot(_gamma(x, 22.48548, 2.926333), [x, 0, 180], rgbcolor='red') (gp+hg).show()
以下の様に正規分布関数_gaussを定義して、ヒストグラム、正規分布、ガンマ分布の推定結果を 重ね合わせてみます。
sageへの入力:
def _gauss(x, mu, sig2):
return 1/sqrt(2*pi*sig2)*exp(-1/(2*sig2)*(x - mu)^2)
sageへの入力:
# 正規分布として推定した結果（2章の結果を利用） sig2 = 226.7 mu = 65.8 ga = plot(_gauss(x, mu, sig2), [x, 0, 180], rgbcolor='green') (hg+gp+ga).show(figsize=5)
最初、flocを指定しないでフィッティングをしたため、fitdistrと結果が異なり 悩んでいました。しかしfloc=0とすることで、fitdistrと同じ結果が得られることが 分かり、floc制約をしない場合、どのような推定結果がでるのか調べて見ました。
今回、可視化はmatplotを使いました。floc=0の設定をしない方がより良くフィット していることが分かります。
sageへの入力:
# 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)
GgSave("5.png", 0.7)
4.37074130051 36.4314864684 6.71934493135 Saving 8.0 x 6.0 in image.
sageへの入力:
# 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)
GgSave("6.png", 0.7)
Saving 8.0 x 6.0 in image.
データに異常データが混ざっている場合の例として、以下の様に2種類の正規分布が混ざっている場合に ついて考えてみます。
sageへの入力:
# 訓練データに異常データがまざっている場合（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=5)
分布は以下の様な混合正規分布として表されます。 $$ 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) $$
sageへの入力:
# テストデータの生成
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=5)
各パラメータの初期値を以下の様に設定します。 $$ (\pi_0, \pi_1) = (0.5, 0.5), (\mu_0, \sigma_0) = (5, 1.0), (\mu_1, \sigma_1) = (-5, 5.0) $$
sageへの入力:
# データを結合 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回の繰り返しで、とても精度良くパラメータが推定されています。
sageへの入力:
# 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))
sageへの入力:
# 結果出力 print pi0, pi1 print mu0, mu1 print sig0, sig1
0.602760455466 0.397239544534 2.98730105891 -0.292164057442 0.512590277559 2.85810740491
sageへの入力:
# g0とg1を合わせて正規化したグラフ
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=5)
sageへの入力:
# 帰属度qn0, qn1の分布 q0 = list_plot(list(qn0), rgbcolor='red') q1 = list_plot(list(qn1)) (q0 + q1).show(figsize=5)
最後にScikit-learnの混合正規分布パッケージsklearn.mixtureを使って、 davisの体重・身長のデータに対して、混合正規分布を求めてみます。
GMMは、異常なデータに影響を受けるため、12番目のデータを除いた集合で、 混合正規分布を求めます。
GMMの引数のn_componentsが求める分布の数です。RのMclustを使うとベストな分布数も計算できますが、 ここでは2として計算します。
sageへの入力:
# 混合正規分布 # SageでのGMMは、http://www15191ue.sakura.ne.jp:8000/home/pub/57/ を参照してください。 from sklearn import mixture # 分類器の生成 classifier = mixture.GMM(n_components=2, covariance_type='full')
sageへの入力:
# GMMを求める # 12番目のデータを除く（体重と身長を間違えたデータ） davis_e11 = davis[davis.index != 11] data = davis_e11[['weight', 'height']].values # 混合正規分布を求める classifier.fit(data)
GMM(covariance_type='full', init_params='wmc', min_covar=0.001, n_components=2, n_init=1, n_iter=100, params='wmc', random_state=None, thresh=None, tol=0.001)
大切なのは、GMMによる分類ではなく、その結果を可視化して特徴を把握することです。
GMMによってデータがどのグループに属するか分類した結果は、predict関数で取得し、 davis_e11のpredにセットします。
次にggplotでpredで色を変えて分布をプロットします。 たったこれだけで、GMMの結果を確認することができます。
sageへの入力:
# 分類結果をpredにセット davis_e11['pred'] = classifier.predict(data) # 分類結果の表示 # 身長（height）と体重（weight）の関係をみる ggplot(davis_e11, aes(x='weight', y='height', colour='pred')) + geom_point()
<ggplot: (8738328197473)>
sageへの入力:
GgSave('fig_3.8.png', fac=0.8)
Saving 11.0 x 8.0 in image.
GMMの結果は、\(\mu\)がmeans_、\(\sigma\)がcovars_、\(\pi\)がweights_変数にセットされています。
Sageのcontour_plot関数を使って正規分布をプロットしてみます。
sageへの入力:
# クラス分けされたデータ cls0 = davis_e11[davis_e11.pred == 0][['weight', 'height']].values cls1 = davis_e11[davis_e11.pred == 1][['weight', 'height']].values
sageへの入力:
# 正規分布の推定値 mu0 = vector(classifier.means_[0]) mu1 = vector(classifier.means_[1]) sig0 = matrix(classifier._get_covars()[0]) sig1 = matrix(classifier._get_covars()[1]) view(mu0); view(mu1) view(sig0); view(sig1)
163.092600076 & 55.5095275006 \\ 55.5095275006 & 55.8057145758 48.9419981283 & 32.6024590071 \\ 32.6024590071 & 45.2848005417
sageへの入力:
# ベクトルのガウス分布を定義
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];
sageへの入力:
x, y = var('x y')
# データのプロット
pt_plt = Graphics()
for pt in cls1:
pt_plt += point(pt, rgbcolor='red')
for pt in cls0:
pt_plt += point(pt, rgbcolor='blue')
# 正規分布を表示
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=5)
混合正規分布の異常度は、以下の様になります。 $$ a(x') = -ln \left \{ \sum_{k=1}^K \hat{\pi}_k \mathcal{N}(x' | \hat{\mu}_k, \hat{\Sigma}_k) \right \} $$
sageへの入力:
#異常度の計算 # pi_nを求める pi_vec = vector(classifier.weights_); view(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 ]
sageへの入力:
list_plot(a, figsize=5)
皆様のご意見、ご希望をお待ちしております。
