sage/入門機械学習による異常検出勉強ノート2章
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2015/08/16からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/56/
また、私の公開しているSageのサーバ（http://www15191ue.sak...
アップロードし、実行したり、変更していろいろ動きを試すこ...
井出 剛著の「入門機械学習による異常検出」(以降、井出本と...
井出本では、Rを使って例題を実行していますが、ここではSage...
** 2章 正規分布に従うデータからの異常検出 [#y0b5dad5]
井出本の基本は、データに対するモデルを使って負の対数尤度...
*** 多変数正規分布に基づく異常検知 [#f276fc2b]
M次元、N個の観測データ\(\mathcal{D}\)が以下の様な多次元正...
$$
\mathcal{N}(x|\mu, \Sigma) = \frac{|\Sigma|^{-1/2}}{(2 \p...
$$
対数尤度は、以下の様に計算されます。
$$
L(\mu, \Sigma | \mathcal{D}) = ln \prod_{n=1}^{N} \mathca...
$$
これに多次元正規分布\(\mathcal{N}\)の定義を代入すると、対...
$$
L(\mu, \Sigma | \mathcal{D}) = -\frac{M N}{2}ln(2 \pi) - ...
$$
これを\(\mu, \Sigma\)で偏微分して、尤度Lを最大にする\(\mu...
最初にLを\(\mu\)で偏微分すると、以下の様になります。
$$
\frac{\partial L(\mu, \Sigma | \mathcal{D})}{\partial \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 ...
$$
と、以下の公式を使います。
$$
\frac{\partial}{\partial A} Tr(A B) = \frac{\partial}{\pa...
$$
$$
\frac{\partial}{\partial A} ln |A| = (A^{-1})^T
$$
偏微分の結果は、以下の様になります。共分散行列\(\Sigma\)...
$$
\begin{eqnarray}
\frac{\partial L(\mu, \Sigma | \mathcal{D})}{\partial (\S...
& = & \frac{N}{2} \Sigma^T - \frac{1}{2}\sum_{n=1}^N (x^{...
& = & \frac{N}{2} \Sigma - \frac{1}{2}\sum_{n=1}^N (x^{(n...
\end{eqnarray}
$$
この値が０になるような\(\hat{\Sigma}\)は、以下の様に求ま...
$$
\hat{\Sigma} = \frac{1}{N} \sum_{n=1}^N (x^{(n)} - \mu) ...
$$
負の対数尤度を基に、異常度\(a(x')\)を以下の様に定義します...
$$
a(x') = ({x'} - \mu)^T \Sigma^{-1} ({x'} - \mu)
$$
しかし、実際のRでのプログラムをみると、ちょっと違っていま...
$$
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...
$$
** データを使って異常度を計算する [#l6f7f976]
Rのcarパッケージに入っているデータDavisを使って、多次元デ...
*** 準備 [#xa9fde68]
いつものように必要なパッケージを読込ます。
sageへの入力:
#pre{{
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
sageへの入力:
#pre{{
# 2.2 carパッケージのDavisを使って例題を試す
# carパッケージのインストールには、時間が掛かりますので注...
r('library(car)')
}}
#pre{{
[1] "car" "jsonlite" "ggplot2" "stats" "gra...
[9] "methods" "base"
}}
*** データの性格を知る [#q1d144e8]
Rのデータをpandaのデータフレームに変換して、データの性格...
最初にinfoとdescribeで個数や欠損値、平均、分散などの統計...
sageへの入力:
#pre{{
# Rのデータをpandaのデータフレームに変換する
davis = RDf2PandaDf("Davis")
# データの個数と欠損値の有無をみる
davis.info()
}}
#pre{{
<class 'pandas.core.frame.DataFrame'>
Int64Index: 200 entries, 0 to 199
Data columns (total 5 columns):
height 200 non-null int64
repht 183 non-null float64
repwt 183 non-null float64
sex 200 non-null object
weight 200 non-null int64
dtypes: float64(2), int64(2), object(1)
memory usage: 9.4+ KB
}}
sageへの入力:
#pre{{
# データの情報を取り出す
davis.describe()
}}
#pre{{
height repht repwt weight
count 200.000000 183.000000 183.000000 200.000000
mean 170.020000 168.497268 65.622951 65.800000
std 12.007937 9.467048 13.776669 15.095009
min 57.000000 148.000000 41.000000 39.000000
25% 164.000000 160.500000 55.000000 55.000000
50% 169.500000 168.000000 63.000000 63.000000
75% 177.250000 175.000000 73.500000 74.000000
max 197.000000 200.000000 124.000000 166.000000
}}
次にDavisの体重のヒストグラムと体重（weight）と身長（heig...
sageへの入力:
#pre{{
# weightの分布をみる
# 1個だけ、160以上の値を示している
ggplot(davis, aes(x='weight')) + geom_histogram(binwidth=...
}}
#pre{{
<ggplot: (8730188077121)>
}}
sageへの入力:
#pre{{
GgSave('fig_2.1.png', fac=0.8)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(th_fig_2.1.jpg);
sageへの入力:
#pre{{
# 身長（height）と体重（weight）の関係をみる
ggplot(davis, aes(x='weight', y='height')) + geom_point()
}}
#pre{{
<ggplot: (8730187375813)>
}}
sageへの入力:
#pre{{
GgSave('fig_2.5.png', fac=0.8)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(th_fig_2.5.jpg);
*** Sageとpandasの機能を使って共分散行列を求める [#i02f72...
Sageとpandas(numpy)の機能を使って共分散を求めみましょう。
平均mxを求めてデータの平均からの差Xcを求めます。pandasで...
数式では共分散行列\(\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{\...
$$
Xcが列ベクトルではなく、行ベクトルなので、以下の様に計算...
$$
\hat{\Sigma} = \frac{1}{N} Xc^T Xc
$$
sageへの入力:
#pre{{
# pandaデータフレームで、weightとheightの平均からのずれXc...
X = davis[['weight', 'height']]
mx = X.mean()
Xc = X - mx
}}
numpyの固有の理由だと思いますが、マトリックスの積に.dot()...
Pnadasの機能を使うと共分散行列は一発で計算できるのですが...
sageへの入力:
#pre{{
# Pandaで標本共分散行列を求めるには
print Xc.T.dot(Xc) /len(Xc)
# もっと簡単なのはcov関数を使う（でも値が微妙に異なるのは...
print Xc.cov()
}}
#pre{{
weight height
weight 226.720 34.2040
height 34.204 143.4696
weight height
weight 227.859296 34.375879
height 34.375879 144.190553
}}
*** Sageのマトリックスを使った異常度の計算 [#ae7fa507]
Sageのmatrixを使った異常度の計算方法を以下に示します。
理論をベースとしているSageと実験値を処理するためのPandas,
numpyをどうやって融合するのかがSage使いの腕の見せ所です。
と言ってもやっていることは簡単です。マトリックスの要素単...
Sageのmatrixにnumpy()メソッドで、numpyのマトリックスに変...
Rやnumpyでは多くのデータを一括して処理することが多く、要...
使われますが、Sageではマトリックスやベクトルの定義に沿っ...
sageへの入力:
#pre{{
# Sageのmatrixを使った計算
# 中心化したデータをSageのマトリックスにする
Xc =matrix((X - mx).values)
}}
sageへの入力:
#pre{{
# 標本共分散行列を求める
Sx = Xc.T*Xc/ Xc.nrows(); Sx
}}
#pre{{
[ 226.7199999999998 34.203999999999986]
[34.203999999999986 143.4696000000001]
}}
sageへの入力:
#pre{{
# (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=5)
}}
&ref(sage0.png);
*** Pandaを使った異常度の計算 [#cfbefb8d]
pandaとnumpyの機能を使って異常値を計算します。ここでは、\...
使っています。また、次数の和には、pythonのsumを使っていま...
出力結果では、１つだけ突出した異常データが見つかります。
全体的な処理から以降ではPandaを使った分散行列の方式を使い...
sageへの入力:
#pre{{
# 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=5)
}}
&ref(sage1.png);
** マハラノビス=タグチ法 [#b1694e66]
異常値の検出では、全体に対する個々の観測値の異常を検出す...
個々の観測値の変数値の異常を検出するには、マハラノビス=タ...
M変数のなかからいくつかの変数を選び、1変数当たりの異常度...
$$
SN_q = -10 log_{10} \left \{ \frac{1}{N'} \sum_{n=1}^{N'}...
$$
*** MASSパッケージのroadデータで実験 [#q3b28eee]
マハラノビス=タグチ法の実行例2.6をSageを使って試してみま...
MASSパッケージのroadデータには、アメリカ２６州の交通事故...
人口密度popden、郊外地区の道路長rural、1月における1日の最...
fuelの６変数を記録しています。
データの前段階でroadをdriversで割って対数化する必要がある...
Rで前処理することにしました（残念）。
取り込んだデータは、infoで個数と欠損値の有無を調べ、内容...
_rowという変な名前が存在するので、stateに置き替えます。
sageへの入力:
#pre{{
# マハラノビス=タグチ法
# 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()
}}
#pre{{
<class 'pandas.core.frame.DataFrame'>
Int64Index: 26 entries, 0 to 25
Data columns (total 6 columns):
_row 26 non-null object
deaths 26 non-null float64
fuel 26 non-null float64
popden 26 non-null float64
rural 26 non-null float64
temp 26 non-null float64
dtypes: float64(5), object(1)
memory usage: 1.4+ KB
}}
sageへの入力:
#pre{{
# _rowという変な名前が存在するので、headで内容をチェック
road.head()
}}
#pre{{
_row deaths fuel popden rural temp
0 Alabama 1.9638 0.5614 0.3401 0.3491 0.3310
1 Alaska 1.5911 0.4470 0.0357 0.4294 1.3157
2 Arizona 2.0098 0.5390 0.1239 0.3094 0.5326
3 Arkanas 2.0740 0.5902 0.3145 0.5842 0.4411
4 Calif 1.7888 0.1046 0.0999 0.1168 0.0660
}}
sageへの入力:
#pre{{
# 州のカラムにstateをセットする
road.columns = ['state', 'deaths', 'fuel', 'popden', '...
# Calif州のデータをピックアップする
road[road.state == 'Calif']
}}
#pre{{
state deaths fuel popden rural temp
4 Calif 1.7888 0.1046 0.0999 0.1168 0.066
}}
*** 1変数当たりの異常度とCalifのSN比を求める [#i9423026]
異常度の計算と同様に異常度を求め、次元数で割って1変数当た...
sageへの入力:
#pre{{
# pandaデータフレームで、平均からのずれXcを計算する
X = road[['deaths', 'fuel', 'popden', 'rural', 'temp']]
mx = X.mean()
Xc = X - mx
Xc.shape[1]
}}
#pre{{
5
}}
sageへの入力:
#pre{{
# 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
# 行単位の和を求め、１変数当たりの異常度を計算する
a = [sum(x)/Xc.shape[1] for x in a_prev]
# リストプロット
list_plot(a, figsize=5)
}}
&ref(sage2.png);
最も異常度の大きい、CalifのデータをXcから取り出し、SN比を...
表示順はdeaths fuel popden rural tempで、図2.7(b)と並...
値は同じ結果になっています。
この結果から、異常に寄与している変数は、通事故死亡者数dea...
fuel、人口密度popdenであることが分かります。
この結果は、車が多く使用される過密な大都市で交通事故死亡...
常識と一致します。
sageへの入力:
#pre{{
# XcからCalifのデータを取得
Xc_prime = Xc[4]
# 要素毎の計算をしたいので、numpyの形式に変換
SN1 = 10*np.log10(Xc_prime^2/Sx.diagonal()); SN1
}}
#pre{{
array([-14.80328085, 12.82597259, -6.26565335, -0.2057...
}}
sageへの入力:
#pre{{
# 左からdeaths fuel popden rural tempで、図2.7(b)と並...
bar_chart(SN1, figsize=5)
}}
&ref(sage3.png);
** コメント [#t17dfc02]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- このページを載せていただきありがとうございます。興味を...
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2015/08/16からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/56/
また、私の公開しているSageのサーバ（http://www15191ue.sak...
アップロードし、実行したり、変更していろいろ動きを試すこ...
井出 剛著の「入門機械学習による異常検出」(以降、井出本と...
井出本では、Rを使って例題を実行していますが、ここではSage...
** 2章 正規分布に従うデータからの異常検出 [#y0b5dad5]
井出本の基本は、データに対するモデルを使って負の対数尤度...
*** 多変数正規分布に基づく異常検知 [#f276fc2b]
M次元、N個の観測データ\(\mathcal{D}\)が以下の様な多次元正...
$$
\mathcal{N}(x|\mu, \Sigma) = \frac{|\Sigma|^{-1/2}}{(2 \p...
$$
対数尤度は、以下の様に計算されます。
$$
L(\mu, \Sigma | \mathcal{D}) = ln \prod_{n=1}^{N} \mathca...
$$
これに多次元正規分布\(\mathcal{N}\)の定義を代入すると、対...
$$
L(\mu, \Sigma | \mathcal{D}) = -\frac{M N}{2}ln(2 \pi) - ...
$$
これを\(\mu, \Sigma\)で偏微分して、尤度Lを最大にする\(\mu...
最初にLを\(\mu\)で偏微分すると、以下の様になります。
$$
\frac{\partial L(\mu, \Sigma | \mathcal{D})}{\partial \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 ...
$$
と、以下の公式を使います。
$$
\frac{\partial}{\partial A} Tr(A B) = \frac{\partial}{\pa...
$$
$$
\frac{\partial}{\partial A} ln |A| = (A^{-1})^T
$$
偏微分の結果は、以下の様になります。共分散行列\(\Sigma\)...
$$
\begin{eqnarray}
\frac{\partial L(\mu, \Sigma | \mathcal{D})}{\partial (\S...
& = & \frac{N}{2} \Sigma^T - \frac{1}{2}\sum_{n=1}^N (x^{...
& = & \frac{N}{2} \Sigma - \frac{1}{2}\sum_{n=1}^N (x^{(n...
\end{eqnarray}
$$
この値が０になるような\(\hat{\Sigma}\)は、以下の様に求ま...
$$
\hat{\Sigma} = \frac{1}{N} \sum_{n=1}^N (x^{(n)} - \mu) ...
$$
負の対数尤度を基に、異常度\(a(x')\)を以下の様に定義します...
$$
a(x') = ({x'} - \mu)^T \Sigma^{-1} ({x'} - \mu)
$$
しかし、実際のRでのプログラムをみると、ちょっと違っていま...
$$
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...
$$
** データを使って異常度を計算する [#l6f7f976]
Rのcarパッケージに入っているデータDavisを使って、多次元デ...
*** 準備 [#xa9fde68]
いつものように必要なパッケージを読込ます。
sageへの入力:
#pre{{
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
sageへの入力:
#pre{{
# 2.2 carパッケージのDavisを使って例題を試す
# carパッケージのインストールには、時間が掛かりますので注...
r('library(car)')
}}
#pre{{
[1] "car" "jsonlite" "ggplot2" "stats" "gra...
[9] "methods" "base"
}}
*** データの性格を知る [#q1d144e8]
Rのデータをpandaのデータフレームに変換して、データの性格...
最初にinfoとdescribeで個数や欠損値、平均、分散などの統計...
sageへの入力:
#pre{{
# Rのデータをpandaのデータフレームに変換する
davis = RDf2PandaDf("Davis")
# データの個数と欠損値の有無をみる
davis.info()
}}
#pre{{
<class 'pandas.core.frame.DataFrame'>
Int64Index: 200 entries, 0 to 199
Data columns (total 5 columns):
height 200 non-null int64
repht 183 non-null float64
repwt 183 non-null float64
sex 200 non-null object
weight 200 non-null int64
dtypes: float64(2), int64(2), object(1)
memory usage: 9.4+ KB
}}
sageへの入力:
#pre{{
# データの情報を取り出す
davis.describe()
}}
#pre{{
height repht repwt weight
count 200.000000 183.000000 183.000000 200.000000
mean 170.020000 168.497268 65.622951 65.800000
std 12.007937 9.467048 13.776669 15.095009
min 57.000000 148.000000 41.000000 39.000000
25% 164.000000 160.500000 55.000000 55.000000
50% 169.500000 168.000000 63.000000 63.000000
75% 177.250000 175.000000 73.500000 74.000000
max 197.000000 200.000000 124.000000 166.000000
}}
次にDavisの体重のヒストグラムと体重（weight）と身長（heig...
sageへの入力:
#pre{{
# weightの分布をみる
# 1個だけ、160以上の値を示している
ggplot(davis, aes(x='weight')) + geom_histogram(binwidth=...
}}
#pre{{
<ggplot: (8730188077121)>
}}
sageへの入力:
#pre{{
GgSave('fig_2.1.png', fac=0.8)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(th_fig_2.1.jpg);
sageへの入力:
#pre{{
# 身長（height）と体重（weight）の関係をみる
ggplot(davis, aes(x='weight', y='height')) + geom_point()
}}
#pre{{
<ggplot: (8730187375813)>
}}
sageへの入力:
#pre{{
GgSave('fig_2.5.png', fac=0.8)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(th_fig_2.5.jpg);
*** Sageとpandasの機能を使って共分散行列を求める [#i02f72...
Sageとpandas(numpy)の機能を使って共分散を求めみましょう。
平均mxを求めてデータの平均からの差Xcを求めます。pandasで...
数式では共分散行列\(\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{\...
$$
Xcが列ベクトルではなく、行ベクトルなので、以下の様に計算...
$$
\hat{\Sigma} = \frac{1}{N} Xc^T Xc
$$
sageへの入力:
#pre{{
# pandaデータフレームで、weightとheightの平均からのずれXc...
X = davis[['weight', 'height']]
mx = X.mean()
Xc = X - mx
}}
numpyの固有の理由だと思いますが、マトリックスの積に.dot()...
Pnadasの機能を使うと共分散行列は一発で計算できるのですが...
sageへの入力:
#pre{{
# Pandaで標本共分散行列を求めるには
print Xc.T.dot(Xc) /len(Xc)
# もっと簡単なのはcov関数を使う（でも値が微妙に異なるのは...
print Xc.cov()
}}
#pre{{
weight height
weight 226.720 34.2040
height 34.204 143.4696
weight height
weight 227.859296 34.375879
height 34.375879 144.190553
}}
*** Sageのマトリックスを使った異常度の計算 [#ae7fa507]
Sageのmatrixを使った異常度の計算方法を以下に示します。
理論をベースとしているSageと実験値を処理するためのPandas,
numpyをどうやって融合するのかがSage使いの腕の見せ所です。
と言ってもやっていることは簡単です。マトリックスの要素単...
Sageのmatrixにnumpy()メソッドで、numpyのマトリックスに変...
Rやnumpyでは多くのデータを一括して処理することが多く、要...
使われますが、Sageではマトリックスやベクトルの定義に沿っ...
sageへの入力:
#pre{{
# Sageのmatrixを使った計算
# 中心化したデータをSageのマトリックスにする
Xc =matrix((X - mx).values)
}}
sageへの入力:
#pre{{
# 標本共分散行列を求める
Sx = Xc.T*Xc/ Xc.nrows(); Sx
}}
#pre{{
[ 226.7199999999998 34.203999999999986]
[34.203999999999986 143.4696000000001]
}}
sageへの入力:
#pre{{
# (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=5)
}}
&ref(sage0.png);
*** Pandaを使った異常度の計算 [#cfbefb8d]
pandaとnumpyの機能を使って異常値を計算します。ここでは、\...
使っています。また、次数の和には、pythonのsumを使っていま...
出力結果では、１つだけ突出した異常データが見つかります。
全体的な処理から以降ではPandaを使った分散行列の方式を使い...
sageへの入力:
#pre{{
# 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=5)
}}
&ref(sage1.png);
** マハラノビス=タグチ法 [#b1694e66]
異常値の検出では、全体に対する個々の観測値の異常を検出す...
個々の観測値の変数値の異常を検出するには、マハラノビス=タ...
M変数のなかからいくつかの変数を選び、1変数当たりの異常度...
$$
SN_q = -10 log_{10} \left \{ \frac{1}{N'} \sum_{n=1}^{N'}...
$$
*** MASSパッケージのroadデータで実験 [#q3b28eee]
マハラノビス=タグチ法の実行例2.6をSageを使って試してみま...
MASSパッケージのroadデータには、アメリカ２６州の交通事故...
人口密度popden、郊外地区の道路長rural、1月における1日の最...
fuelの６変数を記録しています。
データの前段階でroadをdriversで割って対数化する必要がある...
Rで前処理することにしました（残念）。
取り込んだデータは、infoで個数と欠損値の有無を調べ、内容...
_rowという変な名前が存在するので、stateに置き替えます。
sageへの入力:
#pre{{
# マハラノビス=タグチ法
# 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()
}}
#pre{{
<class 'pandas.core.frame.DataFrame'>
Int64Index: 26 entries, 0 to 25
Data columns (total 6 columns):
_row 26 non-null object
deaths 26 non-null float64
fuel 26 non-null float64
popden 26 non-null float64
rural 26 non-null float64
temp 26 non-null float64
dtypes: float64(5), object(1)
memory usage: 1.4+ KB
}}
sageへの入力:
#pre{{
# _rowという変な名前が存在するので、headで内容をチェック
road.head()
}}
#pre{{
_row deaths fuel popden rural temp
0 Alabama 1.9638 0.5614 0.3401 0.3491 0.3310
1 Alaska 1.5911 0.4470 0.0357 0.4294 1.3157
2 Arizona 2.0098 0.5390 0.1239 0.3094 0.5326
3 Arkanas 2.0740 0.5902 0.3145 0.5842 0.4411
4 Calif 1.7888 0.1046 0.0999 0.1168 0.0660
}}
sageへの入力:
#pre{{
# 州のカラムにstateをセットする
road.columns = ['state', 'deaths', 'fuel', 'popden', '...
# Calif州のデータをピックアップする
road[road.state == 'Calif']
}}
#pre{{
state deaths fuel popden rural temp
4 Calif 1.7888 0.1046 0.0999 0.1168 0.066
}}
*** 1変数当たりの異常度とCalifのSN比を求める [#i9423026]
異常度の計算と同様に異常度を求め、次元数で割って1変数当た...
sageへの入力:
#pre{{
# pandaデータフレームで、平均からのずれXcを計算する
X = road[['deaths', 'fuel', 'popden', 'rural', 'temp']]
mx = X.mean()
Xc = X - mx
Xc.shape[1]
}}
#pre{{
5
}}
sageへの入力:
#pre{{
# 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
# 行単位の和を求め、１変数当たりの異常度を計算する
a = [sum(x)/Xc.shape[1] for x in a_prev]
# リストプロット
list_plot(a, figsize=5)
}}
&ref(sage2.png);
最も異常度の大きい、CalifのデータをXcから取り出し、SN比を...
表示順はdeaths fuel popden rural tempで、図2.7(b)と並...
値は同じ結果になっています。
この結果から、異常に寄与している変数は、通事故死亡者数dea...
fuel、人口密度popdenであることが分かります。
この結果は、車が多く使用される過密な大都市で交通事故死亡...
常識と一致します。
sageへの入力:
#pre{{
# XcからCalifのデータを取得
Xc_prime = Xc[4]
# 要素毎の計算をしたいので、numpyの形式に変換
SN1 = 10*np.log10(Xc_prime^2/Sx.diagonal()); SN1
}}
#pre{{
array([-14.80328085, 12.82597259, -6.26565335, -0.2057...
}}
sageへの入力:
#pre{{
# 左からdeaths fuel popden rural tempで、図2.7(b)と並...
bar_chart(SN1, figsize=5)
}}
&ref(sage3.png);
** コメント [#t17dfc02]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- このページを載せていただきありがとうございます。興味を...
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
