PyLearn SIG 04: 機械学習によるクラス分類モデルの構築¶
前回まで、時系列信号を以下2種類の特徴量に変換する術を学びました。
- 時間領域特徴量
- 平均、分散、歪度、尖度、最小値、最大値、中央値、RMS、動的特徴量
- 周波数領域特徴量
- パワーバンド、周波数領域エントロピー、エネルギー
これらは一つ一つ、人間が解釈可能であり、信号の物理的な意味を有する特徴量となります。時系列信号から何らかの状態を判別したい場合、
- 入力: 特徴量
- 出力: 状態
となります。このように何らかの状態を分類する問題を、機械学習の分野ではクラス分類問題と呼びます。ここでは、特徴量からクラス分類問題を解くモデルを構築する方法について解説します。
0. 観測信号のロード¶
前回までは、1回の観測データのみを使用していました。今回は、150回の観測データをロードします。datディレクトリの中に、1.csv〜150.csvがあります。これらはそれぞれ、以下の状況下で測定したセンサデータです。
- クラスA2: モータを設置させたネジを一箇所緩めた場合（50例）1.csv 〜 50.csv
- クラスA1: グリスを塗らずにモータを回した場合（50例）51.csv 〜 100.csv
- クラスN: 適切な状態でモータを回した場合（50例）101.csv 〜 150.csv
- いずれもサンプリング周波数100Hz
機械学習の分野では、分類したい対象をクラスと呼びます。NはNormal（正常）、A1/A2はAnomaly（異常）から取ってきたものです。
import numpy as np
import random
# 正常クラスの信号読み込み
Sig_N=[]
for i in range(101, 151):
Sig_N.append(np.loadtxt("dat/" + str(i) + ".csv",
delimiter=",", skiprows=1, usecols=range(2,8)))
# 異常1クラスの信号読み込み
Sig_A1=[]
for i in range(51, 101):
Sig_A1.append(np.loadtxt("dat/" + str(i) + ".csv",
delimiter=",", skiprows=1, usecols=range(2,8)))
# 異常2クラスの信号読み込み
Sig_A2=[]
for i in range(1, 51):
Sig_A2.append(np.loadtxt("dat/" + str(i) + ".csv",
delimiter=",", skiprows=1, usecols=range(2,8)))
これでSig_N 〜 Sig_A2のリストに対し、6軸のセンサデータが取得されました。
1. 7チャンネルへの拡張（合成加速度の算出）¶
続いて、合成加速度を求めて7チャンネル慣性センサデータへと変換します。関数化しておきましょう。
# 合成加速度を算出する関数:
# Input: 6チャンネルの慣性センサデータ（3軸加速度＋3軸角速度）
# Output: 最後の列に合成加速度が追加された7チャンネルの慣性センサデータ
def MakingMag(signal):
M = signal[:,0] * signal[:,0] + signal[:,1] * signal[:,1] + (signal[:,2]) * (signal[:,2])
M = np.sqrt(M)
M_new = np.reshape(M, (len(M), 1))
# 7チャンネル慣性センサデータの取得
sig7 = np.concatenate([signal, M_new], 1)
return sig7
この関数を使って、150回の観測データ（6チャンネル）を、合成加速度が追加された7チャンネルデータに変換します。
# 正常データ
Sig7_N = []
for i in range(0, len(Sig_N)):
Sig7_N.append(MakingMag(Sig_N[i]))
# 異常1データ
Sig7_A1 = []
for i in range(0, len(Sig_A1)):
Sig7_A1.append(MakingMag(Sig_A1[i]))
# 異常2データ
Sig7_A2 = []
for i in range(0, len(Sig_A2)):
Sig7_A2.append(MakingMag(Sig_A2[i]))
2. 特徴量変換¶
続いて、150回の観測信号を特徴量に変換してみます。これを格納するリストとして、
- FV_N: クラスNの特徴量を保存するリスト
- FV_A1: クラスA1の特徴量を保存するリスト
- FV_A2: クラスA2の特徴量を保存するリスト
を定義します（専門用語としてこれを特徴量ベクトルと呼びます。英語でFeature vector: FV）。変数の構造は、
- FV_N[ i ]: i 回目の観測信号
- FV_N[ i ][ j ]: i 回目の観測信号の軸 j
- FV_N[ i ][ j ][ k ]: i 回目の観測信号の軸 j の特徴量 k
とします（FV_A1, A2も同じ）。
軸 j とは、それぞれ以下の意味を持ちます。
- j = 0: X軸加速度
- j = 1: Y軸加速度
- j = 2: Z軸加速度
- j = 3: X軸角速度
- j = 4: Y軸角速度
- j = 5: Z軸角速度
- j = 6: 合成加速度
k 番目の特徴量は、それぞれ以下のように定義することとします。
- k = 0: 平均、
- k = 1: 分散、
- k = 2: 歪度、
- k = 3: 尖度、
- k = 4: 最小値、
- k = 5: 最大値、
- k = 6: 中央値、
- k = 7: RMS、
- k = 8: 動的特徴量、
- k = 9: パワーバンド（0 〜 10 Hz）
- k = 10: パワーバンド（10 〜 20 Hz）
- k = 11: パワーバンド（20 〜 30 Hz）
- k = 12: パワーバンド（30 〜 40 Hz）
- k = 13: パワーバンド（40 〜 50 Hz）# ここでナイキスト周波数
- k = 14: 周波数領域エントロピー
- k = 15: エネルギー
すなわち、
- FV_N[5][3][2]: クラスNの5回目の観測信号のX軸角速度の歪度
- FV_A1[0][2][4]: クラスA1の0回目の観測信号のZ軸角速度の最小値
となります。なお、これらの特徴量のうち、pythonで用意されていないものもたくさんあるので、まずはそれを求める関数を定義しておきます。
import numpy as np
import scipy.stats as st
from sklearn.linear_model import LinearRegression
from scipy.fftpack import fft
# *****************************************************************
# RMSを求める関数
# Input: 1チャンネルの信号
def Calc_RMS(signal):
a = signal * signal # 二乗
sum_a = np.sum(a) # 総和
sqrt_a = np.sqrt(sum_a) # 平方根
RMS = np.mean(sqrt_a) # 平均値
return RMS
# *****************************************************************
# 動的特徴量を求める関数
# Input: 1チャンネルの信号
def Calc_LinearCoef(signal):
# 入力変数（信号値に対応する連番の数字）
X=np.arange(len(signal))
X=np.reshape(X, (len(X), 1))
# 出力変数
Y=signal
# 線形回帰モデルのパラメータフィッティング
ModelLR = LinearRegression() # 線形回帰モデルのセット
ModelLR.fit(X, Y) # パラメータ獲得（Xは時間、Yは信号値）
# 動的特徴量を返り値に
return ModelLR.coef_[0]
# *****************************************************************
# パワーバンドを求める関数
# signal: 時系列信号（1軸分のみ）
# sf: サンプリング周波数
# min_f: パワーバンドを算出する周波数の最小値
# max_f: パワーバンドを算出する周波数の最大値
def Calc_PowerBand(signal, sf, min_f, max_f):
# パワースペクトル算出
L = len(signal) # 信号長
freq = np.linspace(0, sf, L) # 周波数
yf = np.abs(fft(signal)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
# 指定された周波数の探索
n1, n2 = [], []
for i in range(1, int(L/2)): # <- 直流成分は探索範囲から除外
n1.append(np.abs(freq[i] - min_f))
n2.append(np.abs(freq[i] - max_f))
min_id = np.argmin(n1) # min_fと最も近い周波数が格納された配列番号取得
max_id = np.argmin(n2)# max_fと最も近い周波数が格納された配列番号取得
# 指定された周波数のパワーバンド算出
PB = np.sum(yf[min_id:max_id])
return PB
# *****************************************************************
# 周波数領域エントロピーを求める関数
# signal: 時系列信号（1軸分のみ）
# sf: サンプリング周波数
def Calc_FreqEnt(signal, sf):
# パワースペクトル算出
L = len(signal) # 信号長
freq = np.linspace(0, sf, L) # 周波数
yf = np.abs(fft(signal)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
# 手続き1
a = 0
for i in range(1, int(L/2)): # <- 直流成分除去 & ナイキスト周波数まで
a = a + yf[i]
# 手続き2 & 3
E = 0
for i in range(1, int(L/2)): # <- 直流成分除去 & ナイキスト周波数まで
E = E + (yf[i]/a) * np.log2(yf[i]/a)
E = -E
return E
# *****************************************************************
# エネルギーを算出する関数
# signal: 時系列信号（1軸分のみ）
# sf: サンプリング周波数
def Calc_Energy(signal, sf):
# パワースペクトル算出
L = len(signal) # 信号長
freq = np.linspace(0, sf, L) # 周波数
yf = np.abs(fft(signal)) # 振幅スペクトル
yf = yf * yf # パワースペクトル
# エネルギー算出
En = np.mean(yf[1:int(L/2)]) # <- 直流成分除去 & ナイキスト周波数まで
return En
個別の特徴量を算出する関数ができました。続いて、7チャンネル慣性センサデータを全ての特徴量に変換する関数を定義します。先ほど定義した FV[ i ] [ j ] [ k ]構成の出力を行います。
def FeatureVector(signal, NumClass, NumChannel, NumFeature, sf):
FV = []
for i in range(0, NumClass):
TempList = []
for j in range(0, NumChannel):
TempAry = np.zeros([1, NumFeature])
# 観測信号 i の軸 j の...
sig = signal[i][:, j]
TempAry[0, 0] = np.mean(sig) # j=0 平均
TempAry[0, 1] = np.var(sig) # j=1 分散
TempAry[0, 2] = st.skew(sig) # j=2 歪度
TempAry[0, 3] = st.kurtosis(sig) # j=3 尖度
TempAry[0, 4] = np.min(sig) # j=4 最小値
TempAry[0, 5] = np.max(sig) # j=5 最大値
TempAry[0, 6] = np.median(sig) # j=6 中央値
TempAry[0, 7] = Calc_RMS(sig) # j=7 RMS
TempAry[0, 8] = Calc_LinearCoef(sig) # j=8 動的特徴量
TempAry[0, 9] = Calc_PowerBand(sig, sf, 0, 10) # j=9 PB1
TempAry[0, 10] = Calc_PowerBand(sig, sf, 10, 20) # j=10 PB2
TempAry[0, 11] = Calc_PowerBand(sig, sf, 20, 30) # j=11 PB3
TempAry[0, 12] = Calc_PowerBand(sig, sf, 30, 40) # j=12 PB4
TempAry[0, 13] = Calc_PowerBand(sig, sf, 40, 50) # j=13 PB5
TempAry[0, 14] = Calc_FreqEnt(sig, sf) # j=14 周波数領域エントロピー
TempAry[0, 15] = Calc_Energy(sig, sf) # j=15 エネルギー
TempList.append(TempAry[0])
FV.append(TempList)
return FV
特徴量変換の関数ができたところで、実際にクラスN, A1, A2について、それぞれ50観測分のデータを特徴量に変換します。
NumChannel = len(Sig7_N[0][0, :]) # 軸の最大数
NumFeature = 16 # 扱う特徴量の数
NumClassN = len(Sig7_N) # クラスNのデータ数
NumClassA1 = len(Sig7_A1) # クラスNのデータ数
NumClassA2 = len(Sig7_A2) # クラスNのデータ数
sf = 100 # サンプリング周波数
# クラスNの特徴量変換
FV_N = FeatureVector(Sig7_N, NumClassN, NumChannel, NumFeature, sf)
# クラスA1の特徴量変換
FV_A1 = FeatureVector(Sig7_A1, NumClassA1, NumChannel, NumFeature, sf)
# クラスA2の特徴量変換
FV_A2 = FeatureVector(Sig7_A2, NumClassA2, NumChannel, NumFeature, sf)
これで、7チャンネルのデータを16種の特徴量に変換しました。すなわち、$7 \times 16 = 112$個の特徴量に変換されました。これは、1回の観測を112次元のベクトルに変換したことと等価です。1回の観測を112個の特徴で表現しました、ということです。これらの特徴量を入力として、様々な教師あり学習モデル（決定木、k近傍法、ニューラルネット、サポートベクターマシンなど）で学習させることで、クラスN, A1, A2を分類するモデルが出来上がります。
3. 特徴量の規格化¶
ここらでクラスを分類するモデルを作りたいところですが、その前にもう一つやることがあります。それは、特徴量の規格化です。この必要性を認識するために、
- X軸加速度の平均値（クラスN、0回目の観測データ）
- X軸加速度の分散（クラスN、0回目の観測データ）
- X軸加速度のエネルギー（クラスN、0回目の観測データ）
をそれぞれ表示してみます。
print("平均値: ", round(FV_N[0][0][0]))
print("分散値: ", round(FV_N[0][0][1]))
print("エネルギー: ", round(FV_N[0][0][15]))
見るとわかるように、値のスケールがかなり違います。 これではちょっとよくないことが起こります。 例えば、以下の2次元の平面を考えてください。
- 縦軸: 平均値
- 横軸: 分散値
この空間上で、「上に1進む」と「右に1進む」という二つの意味を考えると、それが有する物理的意味がまるで異なります。より具体的にいうと、
- 上に1進む → 平均値を1上げる → 平均は値のスケールがあまり大きくない → 信号のバイアスを上昇させる。
- 右に1進む → 分散値を1上げる → 分散は値のスケールが極めて大きい → 信号のばらつきをほんの僅かに上昇させる。
となります。このように、同じ「距離が1」でも、上方向か、右方向かにより、移動前と移動後の信号の違いがかなり違います。決定木、k近傍法、サポートベクターマシンのような多くの教師あり学習は、空間上の距離を重要な要素とみなしますし、ニューラルネットワークはそもそも値のスケールが大きすぎるとうまく学習を行うことができます。
このことから、特徴量は規格化して使用します。規格化には2つの方法があります。
- Min-max normalization
- 特徴量をすべて最小0、最大1とする規格化です。$f'_l$は規格化された値、$f_l$はある特徴量の$l$番目の値で、$f_{\rm \{ min, max \}}$は該当する特徴量の最小/最大値とすると、定義は以下のようになります。
- $f'_l = \frac{f_l - f_{\rm min}}{f_{\max} - f_{\rm min}} $
- Z-score normalization
- 特徴量をすべて、平均0、標準偏差1とする規格化です。$f'_l$は規格化された値、$f_l$はある特徴量の$l$番目の値で、$f_{\rm \{ mean, std \}}$は該当する特徴量の平均/標準偏差とすると、定義は以下のようになります。
- $f'_l = \frac{f_l - f_{\rm mean}}{f_{\rm std}} $
なお、$f_{\max}, f_{\rm min}, f_{\rm mean}, f_{\rm std}$は、個別のクラスごとではなく、クラスN, A1, A2すべてまとめて算出してください。というのは、目的はクラスを分類することであり、実際に何らかのシステムを運用するときは、データを入力する前にはどのクラスか明らかではありません。とすると、クラスごとに$f_{\max}, f_{\rm min}, f_{\rm mean}, f_{\rm std}$を求めても、入力データがどのクラスかわからないので、どれを使っていいのかわからなくなってしまいます。
今回は、Z-score normalizationを試してみます。このため、はじめに、各特徴量の平均と標準偏差を求めます。今回はリストの構造として、
- Zscore[ j ][ k ][ 0 ]: 軸 j の特徴量 k の平均
- Zscore[ j ][ k ][ 1 ]: 軸 j の特徴量 k の標準偏差
とすることにします。
# 観測データ数の獲得
NumDatN = len(FV_N)
NumDatA1 = len(FV_A1)
NumDatA2 = len(FV_A2)
Zscore = []
for j in range(0, NumChannel):
temp2 = []
for k in range(0, NumFeature):
temp1 = []
temp_zscore = np.zeros([1, 2])
# 3クラスの50観測、合計150観測に対して
# 軸 j , 特徴量 k を1列に並べる
for i in range(0, NumDatN):
temp1.append(FV_N[i][j][k])
for i in range(0, NumDatA1):
temp1.append(FV_A1[i][j][k])
for i in range(0, NumDatA2):
temp1.append(FV_A1[i][j][k])
# 軸 j , 特徴量 k の 平均・標準偏差を取得
Val_mean = np.mean(temp1)
Val_std = np.std(temp1)
temp_zscore[0, 0] = Val_mean
temp_zscore[0, 1] = Val_std
temp2.append(temp_zscore[0])
Zscore.append(temp2)
規格化を行うためのパラメータ（平均と標準偏差）を求めることができました。それでは、ここで得られた値を用いて規格化を行います。
# クラスNの特徴量の規格化
for i in range(0, NumClassN):
for j in range(0, NumChannel):
for k in range(0, NumFeature):
FV_N[i][j][k] = (FV_N[i][j][k] - Zscore[j][k][0]) / Zscore[j][k][1]
# クラスA1の特徴量の規格化
for i in range(0, NumClassA1):
for j in range(0, NumChannel):
for k in range(0, NumFeature):
FV_A1[i][j][k] = (FV_A1[i][j][k] - Zscore[j][k][0]) / Zscore[j][k][1]
# クラスA2の特徴量の規格化
for i in range(0, NumClassA2):
for j in range(0, NumChannel):
for k in range(0, NumFeature):
FV_A2[i][j][k] = (FV_A2[i][j][k] - Zscore[j][k][0]) / Zscore[j][k][1]
これで規格化が行われました。どのような値になったか、X軸加速度の分散値を3例ずつみてみます。
for i in range(0, 3):
print("クラスN / Xacc / 分散: ", round(FV_N[i][0][1], 2))
print("クラスA1 / Xacc / 分散: ", round(FV_A1[i][0][1], 2))
print("クラスA2 / Xacc / 分散: ", round(FV_A2[i][0][1], 2))
16000近くあったX軸加速度の分散のスケールがおとなしくなりました。続いてエネルギーを見てみます。
for i in range(0, 3):
print("クラスN / Xacc / エネルギー: ", round(FV_N[i][0][15], 2))
print("クラスA1 / Xacc / エネルギー: ", round(FV_A1[i][0][15], 2))
print("クラスA2 / Xacc / エネルギー: ", round(FV_A2[i][0][15], 2))
エネルギーは億を超えた値が出ていました。こちらもいい感じに規格化できています。
4. 特徴量ベクトルへの加工¶
最後に、特徴量ベクトルへと加工します。現在、軸 j と特徴量 kを分離するような変数の管理をしていました。これは人間にとってはわかりやすい構造ですが、教師あり学習のモデルを構築する上で、コンピュータ側にとって、このようなわかりやすさは不要です。pythonで教師あり学習のモデルを構築する際は、1回の観測に対する特徴量を、1次元の長いベクトルに変化させる必要があります。ここではその加工を行います。今回、7チャンネル、16種の特徴量を扱っているので、合計112次元の特徴量ベクトルが生み出されます。
FV_N_ML = np.zeros([NumClassN, 112]) # <- 50回の観測、112次元の特徴量（7チャンネル * 16種の特徴量）
FV_A1_ML = np.zeros([NumClassA1, 112]) # <- 50回の観測、112次元の特徴量（7チャンネル * 16種の特徴量）
FV_A2_ML = np.zeros([NumClassA2, 112]) # <- 50回の観測、112次元の特徴量（7チャンネル * 16種の特徴量）
# クラスN
for i in range(0, NumClassN):
FV_N_ML[i, :] = np.concatenate([
FV_N[i][0], # X軸加速度の特徴量 16種
FV_N[i][1], # Y軸加速度の特徴量 16種
FV_N[i][2], # Z軸加速度の特徴量 16種
FV_N[i][3], # X軸角速度の特徴量 16種
FV_N[i][4], # Y軸角速度の特徴量 16種
FV_N[i][5], # Z軸角速度の特徴量 16種
FV_N[i][6], # 合成加速度の特徴量 16種
], 0)
# クラスA1
for i in range(0, NumClassA1):
FV_A1_ML[i, :] = np.concatenate([
FV_A1[i][0], # X軸加速度の特徴量 16種
FV_A1[i][1], # Y軸加速度の特徴量 16種
FV_A1[i][2], # Z軸加速度の特徴量 16種
FV_A1[i][3], # X軸角速度の特徴量 16種
FV_A1[i][4], # Y軸角速度の特徴量 16種
FV_A1[i][5], # Z軸角速度の特徴量 16種
FV_A1[i][6], # 合成加速度の特徴量 16種
], 0)
# クラスA2
for i in range(0, NumClassA2):
FV_A2_ML[i, :] = np.concatenate([
FV_A2[i][0], # X軸加速度の特徴量 16種
FV_A2[i][1], # Y軸加速度の特徴量 16種
FV_A2[i][2], # Z軸加速度の特徴量 16種
FV_A2[i][3], # X軸角速度の特徴量 16種
FV_A2[i][4], # Y軸角速度の特徴量 16種
FV_A2[i][5], # Z軸角速度の特徴量 16種
FV_A2[i][6], # 合成加速度の特徴量 16種
], 0)
# 特徴量IDと種別を表示
FeatureName = ["平均", "分散", "歪度", "尖度", "最小", "最大", "中央値", "RMS", "動的",
"PB1", "PB2", "PB3", "PB4", "PB5", "エントロピー", "エネルギー"]
AxisName = ["Xacc", "Yacc", "Zacc", "Xang", "Yang", "Zang", "Mag"]
i=0
for j in range(0, NumChannel):
for k in range(0, NumFeature):
print("特徴量ID ", i, ": ", AxisName[j], FeatureName[k])
i = i + 1
これで特徴量ベクトルに加工できました。FV_N_MLには、50行112列のの数字が集まっています。行が50回分の個別の観測結果、列が112種の特徴量の値となります（A1, A2も同じ）。どんな値が得られているのか、imshow関数と呼ばれるヒートマップを作成する関数で可視化してみます。
c=(-2, 2) # 色の最小・最大範囲
from matplotlib import pyplot as plt
plt.figure(figsize=(50, 30)) # 描画領域の横幅, 縦幅
# クラスN
plt.subplot(3,1,1)
plt.imshow(FV_N_ML, clim=c)
plt.colorbar()
plt.title('Feature vector: Class N', size=30)
plt.xlabel('Feature ID', size=30)
plt.ylabel('Observation', size=30)
plt.xticks(np.arange(0, 112, 5), size = 15) # メモリ間隔
plt.yticks(np.arange(0, 50, 5), size = 15)
# クラスA1
plt.subplot(3,1,2)
plt.imshow(FV_A1_ML, clim=c)
plt.colorbar()
plt.title('Feature vector: Class A1', size=30)
plt.xlabel('Feature ID', size=30)
plt.ylabel('Observation', size=30)
plt.xticks(np.arange(0, 112, 5), size = 15) # メモリ間隔
plt.yticks(np.arange(0, 50, 5), size = 15)
# クラスA2
plt.subplot(3,1,3)
plt.imshow(FV_A2_ML, clim=c)
plt.colorbar()
plt.title('Feature vector: Class A2', size=30)
plt.xlabel('Feature ID', size=30)
plt.ylabel('Observation', size=30)
plt.xticks(np.arange(0, 112, 5), size = 15) # メモリ間隔
plt.yticks(np.arange(0, 50, 5), size = 15)
plt.show()
上から順に、クラスN, A1, A2です。縦軸が50回分の観測、横軸が112個の特徴量となります。これを見るだけで、正常なモータ、グリスを塗らない場合、ネジをゆるめた場合で、特徴量にどのような違いが出るか、かなりたくさんの考察ができそうです。
機械学習において重要な特徴量は、以下2つを満たすものです。
1点目: 同じクラスの観測では、毎回、似たような値を取るような特徴量: これは、上の3枚のグラフのある1枚に着目して、縦軸0〜50回の観測において、いずれも似たような色を持つ特徴量が有効である、という意味です。なぜかというと、同じクラス（正常1、正常2、正常3、…）を観測しているのに、毎回違う値が出る特徴量だと、不安定となるためです。
2点目: 異なるクラス間で、異なる値を持つ特徴量: これは、当然と言えば当然です。特徴量の値に応じてクラスを分類するわけですから、異なるクラス間で同じ値を取っている特徴量は、分類に使用することができません。
具体的には、
- だめ特徴量: 46と62（同クラス内でバラついており、異クラス間で似ている特徴量）
- よい特徴量: 36と53（同クラス内でバラついておらず、異クラス間で似ていない特徴量）
実際にこれを、散布図で確認してみます。
# N, A1, A2の順に散布図を作成
plt.figure(figsize=(14, 4)) # 描画領域の横幅, 縦幅
# ダメ特徴量
plt.subplot(1,2,1)
ax1 = 46
ax2 = 62
plt.scatter(FV_N_ML[:,ax1], FV_N_ML[:,ax2], s=50, c='blue', alpha=0.8)
plt.scatter(FV_A1_ML[:,ax1], FV_A2_ML[:,ax2], s=50, c='red', alpha=0.8)
plt.scatter(FV_A2_ML[:,ax1], FV_A2_ML[:,ax2], s=50, c='green', alpha=0.8)
plt.title("Bad: Feature space of F" + str(ax1) + " and F" + str(ax2))
plt.xlabel("Feature " + str(ax1))
plt.ylabel("Feature " + str(ax2))
plt.grid(True) #グリッド線（True:引く、False:引かない）
plt.legend(["N", "A1", "A2"]) # 凡例
# よい特徴量
plt.subplot(1,2,2)
ax1 = 36
ax2 = 53
plt.scatter(FV_N_ML[:,ax1], FV_N_ML[:,ax2], s=50, c='blue', alpha=0.8)
plt.scatter(FV_A1_ML[:,ax1], FV_A2_ML[:,ax2], s=50, c='red', alpha=0.8)
plt.scatter(FV_A2_ML[:,ax1], FV_A2_ML[:,ax2], s=50, c='green', alpha=0.8)
plt.title("Good: Feature space of F" + str(ax1) + " and F" + str(ax2))
plt.xlabel("Feature " + str(ax1))
plt.ylabel("Feature " + str(ax2))
plt.grid(True) #グリッド線（True:引く、False:引かない）
plt.legend(["N", "A1", "A2"]) # 凡例
ダメな特徴量では、クラスをうまく分類することができません。一方、よい特徴量を採用した場合は、いい感じにクラスを分類できそうなことがわかります。このような散布図を、2次元の特徴量空間と呼びます。よい特徴量をきちんと採用できた場合、かなり高精度で正常と異常を切り分けられそうです。...という感じで、かなりたくさんの考察ができそうです。きちんと研究や開発する場合は、上の3枚のマップや、散布図による特徴量空間をきちんと確認して色々と考察した方がいいと思いますが、ここではプログラムについて説明するので、先に進みます。
補足: 機械学習の分野では、次元の呪いと呼ばれる問題があります。これは、採用する特徴量が多すぎると、莫大なデータセットがなくてはうまく人工知能が学習してくれず、精度がかなり低くなってしまう、という問題です。そのため、クラス分類において意味のない特徴量を積極的に削減しようとする手続きを行います。この手続きを特徴量選択（Feature Engineering; FE）と呼びます。クラス分類問題を解く機械学習はたくさん種類がありますが、その多くは、特徴選択を実施してくれるアルゴリズムが導入されていません。例えば、ニューラルネットワークやサポートベクターマシンは、与えた特徴量を削減したりせず、すべて使用してクラス分類問題をときます。すなわち、次元の呪いを回避して性能の良いモデルを構築したい場合、分析者が特徴選択を行う必要があります。一方、決定木と呼ばれるクラス分類問題の解法は、特徴選択もアルゴリズムに組み込まれているので、多くの場合、分析者が特徴選択を行う必要がありません。
特徴量選択の手法としては、クラス内分散・クラス間分散比、Out-of-Bag Error法、Minimum References Set、ReliefF、統計的仮説検定などがありますが、機械学習にあまり慣れてないうちは、勉強することが増えてしまうので、全ての特徴量をとりあえず使う、で良いと思います。
5. ラベリング¶
機械学習である判定モデルを作るためには、入力に対する答えを用意してあげなくては、知能が形成されません。これまで入力を必死に求めていきました。次はそれに対する答えを用意していきます。
# クラスNのラベル
Output_N = []
for i in range(0, NumClassN):
Output_N.append("N")
# クラスA1のラベル
Output_A1 = []
for i in range(0, NumClassA1):
Output_A1.append("A1")
# クラスA2のラベル
Output_A2 = []
for i in range(0, NumClassA2):
Output_A2.append("A2")
# 結合
Y=np.concatenate([Output_N, Output_A1, Output_A2])
print(Y)
これで、入力に対する解答 Y ができました。同じように、問題も X としてまとめてみましょう。
X=np.concatenate([FV_N_ML, FV_A1_ML, FV_A2_ML])
このように、
- 問題: X （150問、112次元の特徴量）
- 解答: Y （150問、ラベルが解答）
という構造のデータセットができました。
6. 教師・テストデータセット分割¶
機械学習では、与えられた問題を正しく正答しようとすると知識が得られます。そのため、学習させれば全問正解してしまうこともほとんどです。しかしそれでは、未知のデータに正しく正解できるかわかりません。このことから、学習に使うデータと精度を評価するデータを分割させることが普通です。これで得られるデータセットを教師データ、テストデータと言います。
教師データ: 機械学習で知能を形成させるためのデータ。学習にはこちらを使用する。このデータの誤差を訓練誤差と呼ぶ。
テストデータ: 学習により得られた知能の精度を評価するためのデータ。こちらは学習には使用しない。このデータの誤差を汎化誤差と呼ぶ。
機械学習で何らかの自動判定を実現しようとする際、目的となるのは汎化誤差の最小化です。訓練誤差の最小化はあんまり意味がないので注意してください。訓練誤差と汎化誤差の関係をより深く知りたい人は、PyLearnMLCR 05: 精度評価指標とモデル選択の基礎を参照してください。
以下に、教師データとテストデータを分割するコードを示します。train_test_split関数を使用します。test_sizeはテストデータの割合、random_stateは乱数のシード（変更すると結果が変わる、同じ値を使用すると毎回同じ結果が出る）です。
from sklearn.model_selection import train_test_split
(X_train, X_test, Y_train, Y_test) = train_test_split(X, Y, test_size=0.3, random_state=1)
print('教師データサイズ: ', np.shape(Y_train))
print('テストデータサイズ: ', np.shape(Y_test))
print('教師データの中身: ', Y_train)
print('テストデータの中身: ', Y_test)
7. 機械学習によるモデル構築¶
7.1.1 決定木の構築¶
決定木によりクラスN, A1, A2を判定するモデルを構築してみます。 決定木の詳細はPyLearnMLCR 01: 決定木を参照してください。以下は決定木を構築するコードです。複雑な木構造が欲しい場合には、min_samples_leafに小さい整数を、単純な木構造が欲しい場合には大きい整数を与えてください。これは、木の一番下のサンプル数をその値以上にする、という制約をかけるパラメータです。
from sklearn import tree
np.random.seed(1) # 擬似乱数シード： 毎回同じ乱数を出す
TreeModel = tree.DecisionTreeClassifier(criterion='gini', min_samples_leaf=1) # モデル構造の定義
TreeModel = TreeModel.fit(X_train, Y_train) # 学習開始
どのような木構造が得られたのか、可視化してみます。可視化は以下のコードで実現できます。
import pydotplus
dot_data = 'Tree.dot'
tree.export_graphviz(TreeModel, out_file=dot_data) # モデル入力
graph = pydotplus.graphviz.graph_from_dot_file(dot_data)
# 保存
graph.write_png('Tree.png') # pngで保存
graph.write_pdf('Tree.pdf') # pdfで保存
# jupyterで表示
from IPython.display import Image
Image(graph.create_png())
木構造を見ると、X17とX47があることがわかります。これは、17番目の特徴量と47番目の特徴量が使用されていることを意味します。
- X17: Y軸加速度の分散
- X47: Z軸加速度のエネルギー
value = [a, b, c]とは、左から順に、A1, A2, Nのクラスの観測数です。ソートされてしまうのでこういう順番になります。木構造から解釈すると、以下の傾向が見て取れます。
正常 クラスNの特徴
- Y軸加速度の分散が小さい。
異常1 クラスA1の特徴
- Y軸加速度の分散が大きく、かつ、Z軸加速度のエネルギーが大きい。
異常2 クラスA2の特徴
- Y軸加速度の分散が大きく、かつ、Z軸加速度のエネルギーが小さい。
さて、本当にこの構造なのか散布図で確認してみます。
# N, A1, A2の順に散布図を作成
ax1 = 17
ax2 = 47
plt.scatter(FV_N_ML[:,ax1], FV_N_ML[:,ax2], s=50, c='blue', alpha=0.8)
plt.scatter(FV_A1_ML[:,ax1], FV_A1_ML[:,ax2], s=50, c='red', alpha=0.8)
plt.scatter(FV_A2_ML[:,ax1], FV_A2_ML[:,ax2], s=50, c='green', alpha=0.8)
plt.title("Good: Feature space of F" + str(ax1) + " and F" + str(ax2))
plt.xlabel("Feature " + str(ax1))
plt.ylabel("Feature " + str(ax2))
plt.grid(True) #グリッド線（True:引く、False:引かない）
plt.legend(["N", "A1", "A2"]) # 凡例
今回は教師データだけでなく、テストデータも含めて可視化している点に注意してください。この結果を見ると、先ほど整理した分岐条件と一致していることがわかります。ただ、ちょっと微妙な特徴量が採用されています。クラスNとクラスA2の距離が近いためです。無論、綺麗に分離はできますが、未知のデータを入力したとき、クラスA2がクラスNの領域に入ってしまったり、あるいはその逆のことが容易に起きそうです。このようなことを避けたい場合は、異なるクラス間の距離が近い特徴量を除去した上で、決定木の学習をかける、などの対応が必要です。
ところで、112個もの特徴量を与えたのに、実際に使用した特徴量はたった二つになってしまいました。決定木は、「必要最小限の情報のみできれいにクラスを分割すること」を目標としたアルゴリズムが導入されているためです。そのため、今回の分類問題においては、この2つで十分きれいに分割できる、という判断がされたわけです。単に問題が単純なことが原因で、決定木が簡単な問題しか解けないわけではありません。
勝手に特徴量を捨てるな！と思う場合には、ニューラルネットワークやサポートベクターマシンの適用を検討しましょう。しかしその場合は、人間に対するわかりやすさは消失するので、注意してください。
7.1.2 決定木の汎化能力の測定¶
先ほどの例では、リーフノード（一番下の四角）にクラスが混ざっていることはありませんでした。したがって、訓練誤差は0、教師データをすべて正しく分類できた、となります。それでは次に、学習に使用しなかったテストデータを活用して汎化能力（未知のデータに対する対応能力）を測定してみます。
テストデータに対する推定結果を得るには、以下のように記述します。
Y_test_est = TreeModel.predict(X_test)
print(Y_test_est)
続いて、推定されたクラス（Y_test_est）と実際のクラス（Y_test）の結果を照合してみます。
print(Y_test_est == Y_test)
最後の1つだけ間違えていることがわかります。どの問題をどう間違えたのか、確認してみます。
print("実測: ", Y_test[-1])
print("推定: ", Y_test_est[-1])
リストで-1を指定すると、最後の要素を取り出すことができます。今回は、A1をA2だと誤答したことがわかりました。ところで、もうちょっと詳しく精度を調べたい場合もあります。このとき、混合行列と呼ばれる結果を算出します。
- 行成分: 実測（上から順にクラスA1, A2, N）
- 列成分: 推定（左から順にクラスA1, A2, N）
混合行列はこのような並びで作られます。対角成分に数字が大きい場合、推定クラスをそれだけ正しく判定できたことを意味します。また、各セルの値を見ることで、どのクラスをどのクラスと間違えたのか、構築したモデルの正誤に対する個性を詳細に知ることができます。
from sklearn.metrics import confusion_matrix
print(confusion_matrix(Y_test, Y_test_est))
1例を除き、対角成分に集まっていることを確認できました。これらはすべて正答です。対角成分以外にある1は、上から1つ目なので実測がクラスA1、左から2つ目なので推定がクラスA2、つまりクラスA1をクラスA2と誤答したという意味があります。
なお、適合率、再現率、F値と呼ばれる精度評価方法もあります。適合率はそのクラスを正しく正答できた割合、再現率はそのクラスを取り漏らさなかった割合、F値は適合率と再現率の相乗平均として定義されます。これは、以下のコードで算出できます。
from sklearn.metrics import precision_recall_fscore_support
precision_recall_fscore_support(Y_test, Y_test_est)
左から順に、クラスA1, A2, Nです。上から順に、適合率、再現率、F値となります。0ほど悪く、1ほど良いとなります。もし、ある状態の検出に対する正しさを重要視したいのであれば適合率を、ある状態の取り漏らしを無くすることを重要視したいのであれば、再現率を重視します。適合率を重視するとちょっと怪しいものは検出しないモデルになり、再現率を重視するとちょっと怪しいものでもどんどん検出するモデルになります。すごい高価な宝石だけを拾うか、石かもしれないものをどんどん拾うか、そのスタンスの違いが表現された指標です。どっちでも良い場合はF値重視でやればokです。
7.1.3 意図的にダメな特徴量を使用した決定木（深い木）¶
先ほど、62と46番目の特徴量がダメであることを散布図で調べました。この特徴量を使って無理やり決定木を作るとどうなるのか、確認してみます。
# ダメな特徴量だけを取り出し（62と46）
X_train_new = np.concatenate([X_train[:, 62:63], X_train[:, 46:47]], 1)
np.random.seed(1) # 擬似乱数シード： 毎回同じ乱数を出す
TreeModel2 = tree.DecisionTreeClassifier(criterion='gini', min_samples_leaf=1) # モデル構造の定義
TreeModel2 = TreeModel2.fit(X_train_new, Y_train) # ダメな特徴量で学習
import pydotplus
dot_data2 = 'Tree2.dot'
tree.export_graphviz(TreeModel2, out_file=dot_data2) # モデル入力
graph2 = pydotplus.graphviz.graph_from_dot_file(dot_data2)
# 保存
graph2.write_png('Tree2.png') # pngで保存
graph2.write_pdf('Tree2.pdf') # pdfで保存
# jupyterで表示
from IPython.display import Image
Image(graph2.create_png())
# テストデータもダメな特徴量だけ取り出し
X_test_new = np.concatenate([X_test[:, 62:63], X_test[:, 46:47]], 1)
# テストデータでの推定
Y_test_est2 = TreeModel2.predict(X_test_new)
# 混合行列
print(confusion_matrix(Y_test, Y_test_est2))
# 適合率、再現率、F値算出
precision_recall_fscore_support(Y_test, Y_test_est2)
精度はどれも6割くらいです。あまり良くないですね。
7.1.4 意図的にダメな特徴量を使用した決定木（浅い木）¶
ダメな特徴量を使っているのだから、精度が低いのは当然です。しかし、ときとして良い特徴量がほとんどない中で、まあまあの精度を狙わなくてはいけないときがあります。
実は、先ほどのモデルの精度が悪かった理由は、2つあります。1つ目はダメな特徴量を使っていること、2つ目はモデルが複雑すぎることです。先ほどのモデルは、min_samples_leafに1が設定されています。つまり、一番下のノードのサンプルサイズに対して、1を許容するという意味です。こうすると、教師データを無理やりすべて正答しようとしてしまいます。こうすると、教師データには正答するが、テストデータには正答しない、いわゆる過学習と呼ばれる現象が起きてしまいます。教師データの精度を上げるすぎることはよくないのです。この問題を解消するために、min_samples_leafに25くらいの数字を入れてみましょう。
np.random.seed(1) # 擬似乱数シード： 毎回同じ乱数を出す
TreeModel3 = tree.DecisionTreeClassifier(criterion='gini', min_samples_leaf=25) # モデル構造の定義
TreeModel3 = TreeModel3.fit(X_train_new, Y_train) # ダメな特徴量で学習
import pydotplus
dot_data3 = 'Tree3.dot'
tree.export_graphviz(TreeModel3, out_file=dot_data3) # モデル入力
graph3 = pydotplus.graphviz.graph_from_dot_file(dot_data3)
# 保存
graph3.write_png('Tree3.png') # pngで保存
graph3.write_pdf('Tree3.pdf') # pdfで保存
# jupyterで表示
from IPython.display import Image
Image(graph3.create_png())
リーフノードに25サンプル以上ないといけないという制約が与えらたので、木が単純になりました。リーフノードを見ると、教師データを完全に分離しきれていません。そのため、教師データを無理やり正答させることで、未知のデータには対応できなくなるという、過学習が解消されることが期待されます。実際にこの精度を見てみます。
# テストデータでの推定
Y_test_est3 = TreeModel3.predict(X_test_new)
# 混合行列
print(confusion_matrix(Y_test, Y_test_est3))
# 適合率、再現率、F値算出
precision_recall_fscore_support(Y_test, Y_test_est3)
精度はまだあまり良いとは言えませんが、深い木よりは明らかに精度が良くなっています。このように、ダメな特徴量を使っていても、学習を工夫することで、精度を高めることが可能です。もしあなたが、精度が悪いという場面に直面したとき、以下2点、検討してください。
- ダメな特徴量を使っていないか？ もっといい特徴量はないか？
- モデルが複雑すぎないか？
そうすると、まあまあ良い感じのモデルを狙えるようになります。
7.2 他の機械学習の適用¶
今回扱っている課題はクラス分類問題ですので、決定木以外でも様々な方法でモデルを構築することができます。興味がある場合は、以下を参照してください。
8. モデルのセーブとロード¶
実運用環境下で使用するには、構築したモデルをセーブすることが必要です。まずはモデルの保存から。
# モデルの保存
import pickle
filename = 'TreeModel3.sav'
pickle.dump(TreeModel3, open(filename, 'wb'))
これでモデルを保存できました。3番目のモデルです。自分のディレクトリにあるか確認してみてください。次に、モデルのロードです。
# モデルのロード
filename = 'TreeModel3.sav' # ファイル名を指定
Load_model = pickle.load(open(filename, 'rb'))
これでモデルがロードできました。実際に使えるか確認してみます。TreeModel3は2つの入力（特徴量）を入れると動くので、この2つの数字を入れてみます。今回は適当に、1つ目を1、2つ目を2にしてみます。
Res = Load_model.predict([[1, 2]]) # <- ロードしたモデルに、適当な特徴量を2つ入力
print(Res)
クラスA1と判定されました。62番目の特徴量と46番目の特徴量が1と2のときは、異常1と推定するようです。
9. まとめ¶
本日は、時間・周波数領域特徴量から、クラス分類問題を決定木により実現する方法を述べました。色々とややこしい手続きが多いですが、慣れるとすぐにできるようになります。