「Pythonではじめる異常検知入門」を寄り道写経 ～ 第３章「異常度と評価指数」
第３章「異常度と評価指数」
書籍の著者 笛田薫 先生、江崎剛史 先生、李鍾賛 先生
この記事は、テキスト「Pythonではじめる異常検知入門」の第３章「異常度と評価指数」の通称「寄り道写経」を取り扱います。
今回はデータの可視化に寄り道しました。
ではテキストを開いて異常検知の旅に出発です🚀
このシリーズは書籍「Pythonではじめる異常検知入門」（科学情報出版、「テキスト」と呼びます）の異常検知の理論・数式とPythonプログラムを参考にしながら、テキストにはプログラムの紹介が無いけれども気になったテーマ、または、テキストのプログラム以外の方法を試したいテーマを「実験的」にPythonコード化する寄り道写経ドキュメンタリーです。
はじめに
テキスト「Pythonではじめる異常検知入門」のご紹介
テキストは、2023年4月に発売された異常検知の入門書です。
数式展開あり、Python実装ありのテキストなのです。
Jupyter Notebook 形式のソースコードと csv 形式のデータは、書籍購入者限定特典として書籍掲載のURLからダウンロードできます。
引用表記
この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「Pythonではじめる異常検知入門－基礎から実践まで－」初版、著者 笛田薫／江崎剛史／李鍾賛、オーム社
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
第３章 異常度と評価指数
Jupyter Notebook 形式（拡張子 .ipynb）でPythonコードを書きます。
テキストのこの章は主に「異常度$${a(x^{\prime})}$$の概念」「混同行列」「$${F}$$値」「ROC曲線とAUC」を掲載しています。
この記事はテキストのご近所さんに当たる以下の寄り道写経に取り組みます。
① 図3-1 の散布図の改造
numpyの乱数生成器とseabornで改造します。
② 図3-4 のヒストグラムの改造
「線は平均…、標準偏差…の正規分布」を追加描画します。
インポート
### インポート
# 数値・確率計算
import pandas as pd
import numpy as np
import scipy.stats as stats
# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'Google Colab をご利用の場合は「描画」箇所を以下のように差し替えてください。
!pip install japanize_matplotlib
import matplotlib.pyplot as plt
import japanize_matplotlib
import seaborn as sns
3-2-1 異常度算出の例１：データ間の距離を参考に正常と異常を考える
テキストのデータ作成および散布図描画コードを引用させていただき、次の内容を変更いたしました。
コードの★印が変更箇所です。
・乱数生成に numpy の「Random Generator」を使います。
・標準正規分布乱数ではなく、正規分布乱数を生成します。
・散布図描画に seaborn を使います。
### 図3-1 異常の例１：車のガソリンの使用量と走行距離の関係を確認した図
## データの作成
# 初期値設定
n_data = 20 # データ数
rng = np.random.default_rng(seed=1) # 乱数生成器★
# データを生成するための関数
def generate_data(n_data):
# 平均20、標準偏差15の正規分布に従う乱数★
x = rng.normal(loc=20, scale=15, size=n_data)
# 6xに平均30、標準偏差30の正規分布に従う乱数を加算★
y = 6*x + rng.normal(loc=30, scale=30, size=n_data)
return x, y
### データを生成
x, y = generate_data(n_data) # データの生成
y[1] = y[1] - 200 # 外れ値の作成
### データを可視化★
fig = plt.figure(figsize=(5, 3))
sns.scatterplot(x=x, y=y, s=80, alpha=0.7)
plt.xlabel('ガソリンの使用量[L]', fontsize=13)
plt.ylabel('走行距離[km]', fontsize=13)
plt.grid(lw=0.5)
plt.show()【実行結果】
乱数生成方法を変更したのでテキストと異なる描画になっています。
テキストは散布図の可視化後、この走行距離データから離れて、正常データ・異常データ・異常度の分布の可視化に進みます。
テキストの可視化コードの寄り道写経はしないことにします。
3-2-2 異常度算出の例２：正規分布を仮定して正常と異常を考える
テキストの図3-4「血圧測定の分布」につながるデータ作成等の一連のコードを寄り道写経いたします。
テキストのシナリオは、ある人の20日間の血圧が正規分布に従っており、95%信頼区間を目安にして異常値を把握してみよう、というものです。
データの登録
テキストのデータを引用して、pandas のデータフレームに格納します。
### データの登録 テキストのデータを引用
values = [132, 133, 133, 130, 136,
134, 133, 132, 134, 131,
134, 135, 134, 132, 132,
133, 135, 136, 141, 133]
data1 = pd.DataFrame({'day': range(1, 21),
'value': values})
print('data1.shape: ', data1.shape)
display(data1.head())【実行結果】
20日間の血圧の時系列データです。
時系列プロット
折れ線グラフで時系列の変化を見てみましょう。
### 時系列プロット
plt.figure(figsize=(8, 3))
sns.lineplot(data=data1, x='day', y='value', lw=0.9)
plt.xticks(data1['day'])
plt.xlabel('日にち')
plt.ylabel('最高血圧[mmHg]')
plt.grid(lw=0.5);【実行結果】
19日目の血圧が突出している感じがします。
基本統計量の算出
pandas の関数を用いて、平均・分散・標準偏差を算出します。
分散・標準偏差は標本サイズ$${n}$$で割り算しています（不偏統計量ではありません）。
### 統計量の算出
x_mean = data1['value'].mean()
x_var = data1['value'].var(ddof=0)
x_std = data1['value'].std(ddof=0)
print(f'平均 : {x_mean:8.3f}')
print(f'分散 : {x_var:8.3f}')
print(f'標準偏差: {x_std:8.3f}')【実行結果】
ちなみに pandas の describe 関数の標準偏差 std は不偏標準偏差です。
print('describe() : ', data1['value'].describe().loc['std'])
print('不偏標準偏差: ', data1['value'].std(ddof=1))【実行結果】
95%信頼区間の算出
scipy.stats を用いており、テキストとほぼ同じコードです。
### テキスト52ページの95%信頼区間
conf_int95 = stats.norm.interval(confidence=0.95, loc=x_mean, scale=x_std)
print(f'95%信頼区間 [{conf_int95[0]:.3f}, {conf_int95[1]:.3f}]')【実行結果】
図3-4 のヒストグラムの描画
今までのコードで算出した平均、標準偏差、95%信頼区間を用いて、テキストのヒストグラムに ①正規分布の確率密度関数曲線と、② 95%信頼区間の区切り線を追加描画します。
### 図3-4 Aさんの20日間の血圧測定の結果の分布
## 描画用データの作成
# ビンの設定
bins = np.linspace(128, 142, 15)
# valueの平均と標準偏差をパラメータにする正規分布の確率密度関数の算出
xval = np.linspace(bins.min(), bins.max(), 1001)
yval = stats.norm.pdf(x=xval, loc=x_mean, scale=x_std)
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 3))
twinx1 = ax.twinx() # 正規分布の確率密度関数の軸
# ヒストグラムの描画
sns.histplot(data=data1, x=data1['value'], bins=bins, ec='white', ax=ax)
ax.grid(lw=0.5)
# 95%信頼区間の両端の描画
ax.axvline(conf_int95[0], color='tab:red', ls='--')
ax.axvline(conf_int95[1], color='tab:red', ls='--')
ax.set(title=f'95%信頼区間 [{conf_int95[0]:.1f}, {conf_int95[1]:.1f}]',
xlabel='最高血圧[mmHg]', ylabel='回数')
# 正規分布の確率密度関数の描画
sns.lineplot(x=xval, y=yval, lw=0.8, ax=twinx1)
twinx1.set(yticks=[]);【実行結果】
青線が正規分布の確率密度関数です。
ヒストグラムの分布は正規分布に近似してそうですね！
また95%信頼区間から外れるのは、横軸141～142の19日目だけのようです。
第３章の寄り道写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
この記事が参加している募集
この記事が気に入ったらサポートをしてみませんか？