FrontPage

2015/08/16からのアクセス回数 6227

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://www15191ue.sakura.ne.jp:8000/home/pub/58/

また、私の公開しているSageのサーバ（http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。

井出　剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。

井出本では、Rを使って例題を実行していますが、ここではSage, Pythonをベースに例題を試します。

3章　非正規データからの異常検知 †

井出本の基本は、データに対するモデルを使って負の対数尤度を求め、それを異常検出関数として使うことです。

↑

分布が左右対称でない場合 †

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.

th_fig_2.1.jpg

↑

ガンマ分布の対数尤度 †

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を求める代わりにモーメント法を使って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の体重を使って計算してみる †

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を使ってガンマ分布を推定 †

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
↑

scipyのstatsを使ってガンマ分布を推定 †

このノートの特徴は、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)

sage0.png

sageへの入力:

# 正規化したグラフも可能
hg = histogram(weight, figsize=5, bins=24, normed=True)
hg.show()

sage1.png

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()

sage2.png

↑

2章の正規分布との重ね合わせ †

以下の様に正規分布関数_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)

sage3.png

↑

stats.gamma.fitの精度 †

最初、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.

th_fig-5.jpg

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.

th_fig-6.jpg

↑

訓練データに異常データがまざっている場合 †

データに異常データが混ざっている場合の例として、以下の様に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)

sage4.png

分布は以下の様な混合正規分布として表されます。 $$ 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アルゴリズムを使って求めます。

EMアルゴリズムでは、以下の手順でパラメータの値を求めます。

  1. パラメータ\(\theta\)の初期値をセット
  2. 潜在変数（ここではどの正規分布に属するか）の期待値\(q_i^{(n)}\)が最大になるようなパラメータを求める
  3. 求まった潜在変数からパラメータの値を再設定する
  4. 値が収束するまで２）に戻る
↑

テストデータの生成 †

以下の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)

sage5.png

↑

初期値のセット †

各パラメータの初期値を以下の様に設定します。 $$ (\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アルゴリズムの実行 †

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)

sage6.png

sageへの入力:

# 帰属度qn0, qn1の分布
q0 = list_plot(list(qn0), rgbcolor='red')
q1 = list_plot(list(qn1))
(q0 + q1).show(figsize=5)

sage7.png

↑

パッケージを使って混合正規分布（GMM）を求める †

最後に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.

th_fig_3.8.jpg

↑

正規分布の形もプロット †

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)

sage8.png

↑

混合正規分布データの異常度 †

混合正規分布の異常度は、以下の様になります。 $$ 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)

sage9.png

↑

コメント †

選択肢 投票
おもしろかった 14  
そうでもない 0  
わかりずらい 0  

皆様のご意見、ご希望をお待ちしております。


(Input image string)

添付ファイル: fileth_fig-6.jpg 1085件 [詳細] fileth_fig-5.jpg 1108件 [詳細] fileth_fig_3.8.jpg 1009件 [詳細] fileth_fig_2.1.jpg 1064件 [詳細] filesage9.png 1034件 [詳細] filesage8.png 1028件 [詳細] filesage7.png 1199件 [詳細] filesage6.png 1084件 [詳細] filesage5.png 1086件 [詳細] filesage4.png 1119件 [詳細] filesage3.png 1055件 [詳細] filesage2.png 1222件 [詳細] filesage1.png 1027件 [詳細] filesage0.png 1014件 [詳細]
トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2023-11-08 (水) 08:31:40 (170d)
SmartDoc