南西フォーラムAIミニセミナー

AIを活用した故障検知器の作り方

〜時系列信号を対象としたAIの開発・構築手法の事例〜

2019.02.25 於 さがみはら産業創造センター 〒252-0131 神奈川県相模原市緑区西橋本５丁目４−２１

講師: 大前佑斗
東京高専 電気工学科
助教 博士（工学）
連絡先: yuto.omae@gmail.com

0. はじめに

本日はお忙しいところ、お集まりいただきありがとうございました。本日は以下の点について解説します。

目次

    1. 人工知能とは何か
    1. 古典的方法による人工知能の構成方法
      2.1 加速度・角速度データ
      2.2 データの読み込み
      2.3 特徴量空間と人工知能
      2.4 特徴量変換
      2.5 教師データとテストデータ
      2.6 次元の呪い
      2.7 有効な特徴量の探索（自分で探そう）
      2.8 有効な特徴量の探索（アルゴリズムで探そう）
      2.9 人工知能の構築
      　2.9.1 決定木
      　2.9.2 サポートベクターマシン
      　2.9.3 ニューラルネットワーク
      2.10 精度評価: 適合率（Precision）　VS　再現率（Recall）
    1. 深層学習による人工知能の構築
      3.1 深層学習による特徴量学習
      3.2 深層学習の特徴量と古典的な特徴量はどちらが有効なのか？
      3.3 深層学習による特徴量を利用した人工知能の構築
    1. おわりに
  • Appendix. 環境構築

プログラムリスト

  • Programming 01. センサデータのロード
  • Programming 02. 波形の表示
  • Programming 03. 特徴量変換（平均と分散）
  • Programming 04. 教師・テストデータセット作成
  • Programming 05. 目で見て有効な特徴量を探す
  • Programming 06. Out-of-Bag Error法による特徴量重要度評価
  • Programming 07. Out-of-Bag Error法において最適となる特徴量空間の可視化
  • Programming 08. 最適な2つの特徴量のみを抽出
  • Programming 09. 決定木の構築
  • Programming 10. 決定木による識別境界の可視化
  • Programming 11. ニューラルネットワークの学習
  • Programming 12. ニューラルネットワークによる識別境界の可視化
  • Programming 13. ニューラルネットワークによる識別境界の可視化（動画）
  • Programming 14. 線形・非線形サポートベクターマシンの学習
  • Programming 15. サポートベクターマシンによる識別境界の可視化
  • Programming 16. 過学習を起こさせた非線形サポートベクターマシンによる識別境界の可視化
  • Programming 17. テストデータによる精度評価（適合率、再現率、F値）
  • Programming 18. 深層学習による特徴量学習
  • Programming 19. 深層学習により抽出された特徴量の可視化（動画）
  • Programming 20. Out-of-Bag Error法による特徴量重要度評価（深層学習の特徴量 VS 古典的な特徴量）
  • Programming 21. 深層学習とニューラルネットワークの連携（一般的な深層学習）（動画）
  • Programming 22. 深層学習とサポートベクターマシンの連携（動画）
  • Programming 23. 深層学習と決定木の連携（動画）

これらを学ぶことで、以下のような具体的課題に適用可能です。

  1. 心拍から心臓の状態を判断する。
  2. 加速度・角速度データから機械の状態を判断する。
  3. 加速度・角速度データから歩き方がおかしい人を判定する。
  4. マイクの音声データから、吃音症かどうかを判定する。
  5. ...など

1. 人工知能とは何か

ご存知の方も多いと思いますが、人工知能とはただの最適化数学の応用でしかありません。以下の例を考えてみてください。

  • 今、何時間勉強したら、テストで何点を取れるのか知りたいとします。
  • 勉強時間を$x$、テストの点数を$y$とします。
  • 勉強した方がテストの点数が良さそうなので、$y=ax+b$というモデルで予測を試みます。（★）
  • $a, b$の決め方によって、予測モデルの性能が左右されます。
  • モデルの性能の良さを定義します。(▲)
  • 実際にデータをたくさん集めて、当てはまりの良い$a, b$を探します。（■）
  • 最終的に、$a, b$が確定し、勉強時間$x$からテストの点数$y$を知ることができます。
  • 例えば、$a, b = 5, 10$が最適なのであれば、$y=5x+10$が最終的に出来上がったモデルです。

最近あちこちで人工知能、人工知能と呼ばれていますが、上のの手続きをやっているにすぎません。 人工知能の分野では、(★)、(▲)、(■)に、以下のような専門用語がついています。


記号 呼び方 意味 具体的な方法
★ モデル 知能の形状みたいなもの ニューラルネット、SVM、決定木、深層学習など
▲ 評価関数 頭の良さの定義づけ 平均二乗誤差、適合率、再現率、F値など
■ 学習 　頭が良くなるように$a, b$を変更 　最急勾配降下法、モーメンタム、遺伝的アルゴリズムなど


人工知能に興味があって自分で勉強している人などは、具体的な方法に示したことを知っていると思います。 これを点で理解していると、どれがなんなのかわかりませんが、点と点を結んだ線として理解すると、★、▲、■の関係になります。 人工知能そのものの研究をしている人たちは、★、▲、■に新しい方法を提案しています （今回はあまり関係ありませんが、大前は(▲)を専門としています）。

人工知能を作ったことがある人もない人も、作る場合には、

  • 自分はどんなモデルを採用していて、
  • 頭の良さ（評価関数）をどのように定義していて、
  • どのような学習を行わせているのか

という点は、最低限意識する必要はあります。

2. 古典的方法による人工知能の構成方法

今回は具体例として、慣性センサから機械の異常状態を検出するという課題を扱ってみます。

2.1 加速度・角速度データ

一般的な慣性センサは、3軸加速度と3軸角速度を取得することができます。 センサデータを扱うときは、サンプリング周波数をきちんと記録することが必要です。 サンプリング周波数とは、1秒間に何個データを測定するかという意味です。 例えば、サンプリング周波数が100Hzの場合は、1秒間に100個データを取る、となります。 データ点数が500個あれば、5秒間の測定データだとわかるわけです。 しかし、サンプリング周波数を記録していなければ、何秒のデータなのかわからなくなってしまうので、データの取り直しになってしまいます。

2.2 データの読み込み

今回、事例としてあるモータの横に、慣性センサをつけて、サンプリング周波数100Hzとして、150例のデータを取得しました。 取得したデータの内訳は以下の通りです。

  • クラスN: 適切な状態でモータを回した場合（50例）
  • クラスA1: グリスを塗らずにモータを回した場合（50例）
  • クラスA2: モータを設置させたネジを一箇所緩めた場合（50例）

人工知能の分野において、分類したい状態のことをクラスと呼びます。今回で言えば、N, A1, A2がクラスに相当するわけです。 そして、入力したセンサデータからどのクラスに属するのか分類する問題を、クラス分類問題と呼びます。

  • Signal_N: クラスNの波形（正常）

    • 50本のデータ、3軸加速度・3軸加速度、サンプリング100Hzで1秒分の信号、$50\times6\times100サイズ$

  • Signal_A1: クラスA1の波形（異常1）

    • 50本のデータ、3軸加速度・3軸加速度、サンプリング100Hzで1秒分の信号、$50\times6\times100サイズ$

  • Signal_A2: クラスA2の波形（異常2）

    • 50本のデータ、3軸加速度・3軸加速度、サンプリング100Hzで1秒分の信号、$50\times6\times100サイズ$

CSVデータの中には、以下のようなセンサデータが入っているイメージです。

Seq Xacc Yacc Zacc Xang Yang Zang
1 1.3 2.4 3.2 3.6 2.4 5.1
2 3.4 2.5 1.3 3.2 2.2 5.7
... ... ... ... ... ... ...
100 0.1 0.3 1.4 1.2 4.2 5.7

Programming 01. センサデータのロード

In [1]:
import numpy as np
import random

# 正常クラスの信号読み込み
Signal_Class_N_temp=[]
for i in range(101, 151):
    Signal_Class_N_temp.append(np.loadtxt("sensdat/dat/" + str(i) + ".csv", delimiter=",", skiprows=1, usecols=range(1,8)))

# 異常1クラスの信号読み込み
Signal_Class_A1_temp=[]
for i in range(51, 101):
    Signal_Class_A1_temp.append(np.loadtxt("sensdat/dat/" + str(i) + ".csv", delimiter=",", skiprows=1, usecols=range(1,8)))

# 異常2クラスの信号読み込み
Signal_Class_A2_temp=[]
for i in range(1, 51):
    Signal_Class_A2_temp.append(np.loadtxt("sensdat/dat/" + str(i) + ".csv", delimiter=",", skiprows=1, usecols=range(1,8)))

# 信号長を1秒に固定
SigLength = 100
Signal_Class_N = []
Signal_Class_A1 = []
Signal_Class_A2 = []
for i in range(0,50):
    Signal_Class_N.append(Signal_Class_N_temp[i][0:SigLength,:])
    Signal_Class_A1.append(Signal_Class_A1_temp[i][0:SigLength,:])
    Signal_Class_A2.append(Signal_Class_A2_temp[i][0:SigLength,:])

Programming 02. 波形の表示

ロードしたセンサデータ波形を表示してみます。1クラス1信号のみです。

  • 左側: クラスN（正常）　　中央: クラスA1（異常1）　　右側: クラスA2（異常2）
In [3]:
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(20, 4))
ViewID=0

# クラスN 正常
plt.subplot(1,3,1)
plt.plot(Signal_Class_N[ViewID][:,1:4])
plt.title('Class N')
plt.xlabel('Sequence')
plt.ylabel('Amplitude')
plt.legend(['Xacc', 'Yacc', 'Zacc', 'Xang', 'Yang', 'Zang'])
plt.ylim(-2000, 13000)

# クラスA1 異常
plt.subplot(1,3,2)
plt.plot(Signal_Class_A1[ViewID][:,1:4])
plt.title('Class A1')
plt.xlabel('Sequence')
plt.ylabel('Amplitude')
plt.legend(['Xacc', 'Yacc', 'Zacc', 'Xang', 'Yang', 'Zang'])
plt.ylim(-2000, 13000)

# クラスA2 異常
plt.subplot(1,3,3)
plt.plot(Signal_Class_A2[ViewID][:,1:4])
plt.title('Class A2')
plt.xlabel('Sequence')
plt.ylabel('Amplitude')
plt.legend(['Xacc', 'Yacc', 'Zacc', 'Xang', 'Yang', 'Zang'])
plt.ylim(-2000, 13000)
Out[3]:
(-2000, 13000)

2.3 特徴量空間と人工知能

この後の処理は、深層学習か、古典的な機械学習かで、微妙に処理が異なってきますが、いずれにせよ、人工知能の中身にあるイメージを理解することは大変重要です。
人工知能をブラックボックスで使わないための入り口は、以下の言葉を理解することだと考えています。

  • 人工知能の正体は$n$次元の空間に生成された、$n-1次元$の部分空間であること。

ややこしいと思うので、$n=2, 3$で考えてみます。

  • 人工知能の正体は
    • $2$次元の空間に生成された、$1次元$の部分空間であること。
    • $3$次元の空間に生成された、$2次元$の部分空間であること。

もうちょっと噛み砕くと、以下のようになります。

  • 人工知能の正体は
    • 面に生成された、線であること。
    • 空間に生成された、面であること。

やはりわかりにくいと思うので、イラストで考えてみます。

上の例では、2次元の空間に1次元の識別境界を生成しています。この1次元こそがAIであることはイメージがついたでしょうか。 これを抽象化させると、

  • 人工知能の正体は$n$次元の空間に生成された、$n-1次元$の部分空間であること。

となります。

2.4 特徴量変換

それでは、$6\times100$の信号を特徴量に変換してみます。上で述べたとおり、特徴量には、

  • 時間領域特徴量（平均、分散など）
  • 周波数領域特徴量（特定周波数帯のスペクトル密度など）
  • 深層学習由来の特徴量

の3種類があることがわかりました。周波数成分まで踏み込むとややこしくなるので、ひとまず時間領域特徴量の特に平均と分散だけで考えてみます（実際にやる場合は、高周波成分のスペクトル密度なども、特徴量として採用したほうがいいです）。

  • 平均値（Mean）

    • 信号の平均的な値を表します。

  • 分散値（Variance）

    • 信号のばらつき具合を表します。

この2つを採用するので、具体的に算出する特徴量は、以下の12個となります。

  • Mean(Xacc): X軸加速度の平均値
  • Mean(Yacc): Y軸加速度の平均値
  • Mean(Zacc): Z軸加速度の平均値
  • Mean(Xang): X軸角速度の平均値
  • Mean(Yang): Y軸角速度の平均値
  • Mean(Zang): Z軸角速度の平均値
  • Var(Xacc): X軸加速度の分散値
  • Var(Yacc): Y軸加速度の分散値
  • Var(Zacc): Z軸加速度の分散値
  • Var(Xang): X軸角速度の分散値
  • Var(Yang): Y軸角速度の分散値
  • Var(Zang): Z軸角速度の分散値

これにより、600個の数値で表現される一回の観測結果が、12個の数字で表現することができます。

Programming 03. 特徴量変換（平均と分散）

以下のコードは、信号を平均・分散の特徴量に変換するコードです。

In [4]:
NumSignal = 50 # 信号の数（1クラス50信号）
NumAx = 6 # 軸の数（3軸加速度・角速度）

# クラスN（正常データ）の平均値と分散値を算出
Features_N=[[]]
for i in range(0, NumSignal):
    for j in range(1, NumAx+1):
        Features_N[i].append(np.round(np.mean(Signal_Class_N[i][:,j]),3))
    for j in range(1, NumAx+1):
        Features_N[i].append(np.round(np.var(Signal_Class_N[i][:,j]),3))
    if i != NumSignal-1:
        Features_N.append([])
        
# クラスA1（異常1データ）の平均値と分散値を算出
Features_A1=[[]]
for i in range(0, NumSignal):
    for j in range(1, NumAx+1):
        Features_A1[i].append(np.round(np.mean(Signal_Class_A1[i][:,j]),3))
    for j in range(1, NumAx+1):
        Features_A1[i].append(np.round(np.var(Signal_Class_A1[i][:,j]),3))
    if i != NumSignal-1:
        Features_A1.append([])
        
# クラスA2（異常２データ）の平均値と分散値を算出
Features_A2=[[]]
for i in range(0, NumSignal):
    for j in range(1, NumAx+1):
        Features_A2[i].append(np.round(np.mean(Signal_Class_A2[i][:,j]),3))
    for j in range(1, NumAx+1):
        Features_A2[i].append(np.round(np.var(Signal_Class_A2[i][:,j]),3))
    if i != NumSignal-1:
        Features_A2.append([])


# 各クラス、一本の信号だけ特徴量を表示
Factor = ['1 :Xacc平均', '2 :Yacc平均', '3 :Zacc平均', # 3軸加速度の平均
          '4 :Xang平均', '5 :Yang平均', '6 :Zang平均', # 3軸角速度の平均
          '7 :Xacc分散', '8 :Yacc分散', '9 :Zacc分散', # 3軸加速度の分散
          '10:Xang分散', '11:Yang分散', '12:Zang分散'] # ３軸角速度の分散

print(' *** 各クラスの0番目の信号の特徴量（12次元） ***')    
print('特徴量名   | クラスN   | クラスA1   | クラスA2')
for i in range(0,12):
    print(Factor[i], '  | ', Features_N[0][i], ' | ', Features_A1[0][i], ' | ', Features_A2[0][i])

 *** 各クラスの0番目の信号の特徴量（12次元） ***
特徴量名   | クラスN   | クラスA1   | クラスA2
1 :Xacc平均   |  -94.94  |  -51.5  |  -59.55
2 :Yacc平均   |  -73.1  |  -28.76  |  -79.43
3 :Zacc平均   |  9853.74  |  9705.2  |  9907.4
4 :Xang平均   |  -16.63  |  -15.42  |  -23.03
5 :Yang平均   |  -1.17  |  -15.58  |  -26.54
6 :Zang平均   |  5.09  |  3.7  |  6.51
7 :Xacc分散   |  11833.816  |  519735.37  |  415307.428
8 :Yacc分散   |  11627.19  |  615169.502  |  173216.065
9 :Zacc分散   |  165086.432  |  436522.46  |  1506606.56
10:Xang分散   |  7346.353  |  19951.044  |  107218.409
11:Yang分散   |  7348.821  |  60853.064  |  112902.468
12:Zang分散   |  135.622  |  2105.51  |  3939.17

各クラス（正常、異常1、異常2）の0番目の信号の特徴量だけを表示してみました。

  • 特徴量に変換する前: $100\times6$次元の行列（サンプリング周波数100Hzにおける1秒の信号、6軸）
  • 特徴量に変換した後: $12$次元のベクトル

というように、600次元の信号が12次元の特徴量に変化したことはわかるでしょうか。このように、人工知能を構築するには、低次元化を行う必要があります。

2.5 教師データとテストデータ

人工知能の分野では、データセットを2つ、教師データ、テストデータに分割します。 これはそれぞれ以下の意味を持ちます。

  • 教師データ：
    • 人工知能が学習を行うためのデータセット
  • テストデータ：
    • 人工知能の性能を評価するためのデータセット

人工知能を構築する様々な手法は、通常、教師データを正しく正答させるような学習を行います。ただし、教師データを良く当てるようなモデルを作ったとしても、それが直ちに良いモデルだということができません。人間の学力調査において、教材の問題はよく解けるが、期末テストではボロボロであることを通常、頭が良いとは言いません。教師データの精度が良いとはすなわち、教材の問題をよく解けることと等価で、期末テストでどのくらい解けるかを意味するわけではありません。人工知能の分野でも、期末テストを用意して、本当の頭の良さを測ってあげる必要があります。このためのデータセットのことを、テストデータと言います。

今回は150信号（1クラス50信号、3クラス）なので、以下のようにデータセットを分割します。

  • 教師データ： 90信号（1クラス30信号、3クラス）
  • テストデータ： 60信号（1クラス20信号、3クラス）

Programming 04. 教師・テストデータセット作成

以下のコードは、教師データ、テストデータを取得するコードです。それぞれ以下の意味を持ちます。

  • X_train: 教師データの特徴量（12次元、90信号分）
  • Y_train: 教師データの正答ラベル（クラスラベルN/A1/A2、それぞれ30ずつ、合計90信号分）
  • X_test: テストデータの特徴量（12次元、60信号分）
  • Y_test: テストデータの正答ラベル（クラスラベルN/A1/A2、それぞれ20ずつ、合計60信号分）
In [5]:
# 正答ラベルの付与
ClassLabel_N=[]
for i in range(0, NumSignal):
    ClassLabel_N.append("N")

ClassLabel_A1=[]
for i in range(0, NumSignal):
    ClassLabel_A1.append("A1")

ClassLabel_A2=[]
for i in range(0, NumSignal):
    ClassLabel_A2.append("A2")

# 教師データセット（1クラス0〜29番目の30信号、3クラスのため、合計90信号）
X_train = np.concatenate([Features_N[0:30], Features_A1[0:30], Features_A2[0:30]])
Y_train = np.concatenate([ClassLabel_N[0:30], ClassLabel_A1[0:30], ClassLabel_A2[0:30]])

# テストデータセット（1クラス30〜49番目の20信号、3クラスのため、合計60信号）
X_test = np.concatenate([Features_N[30:49], Features_A1[30:49], Features_A2[30:49]])
Y_test = np.concatenate([ClassLabel_N[30:49], ClassLabel_A1[30:49], ClassLabel_A2[30:49]])

2.6 次元の呪い

さて、モータの異常状態を検出したいとき、どんな情報が必要でしょうか。横揺れのばらつき、縦揺れのばらつき、高周波成分の大きさ……、色々と必要そうです。 当然、情報は多ければ多いほど、有利です。 人間が何らかの知的な判断を行うときも、情報は少ないよりも多い方が有利です。 そう考えると、たくさんの情報、すなわち、たくさんの特徴量を採用して人工知能を構築することは、一見良さそうに思えます。

実はこれ、一概に正しいとは言えないのです。 より具体的に言うと、たくさんの特徴量を採用することには、メリットとデメリットが混在しているのです。 これをイメージした図を以下に示します。

このように、人工知能に与える情報を増やしすぎたとき、教師データが莫大にければ良い精度が出せなくなる現象を、次元の呪いと言います。 すなわち、今回用意した12次元の特徴量をすべて使用するのではなく、何らかの形で12次元の特徴量の中から、不要なもの、有効なものを分別させることが必要になります。

大量の特徴量の中から人工知能の精度向上に寄与する優れた特徴量を探索すると言う問題は、Feature Engineering と呼ばれる人工知能の学術分野となります。これは、数ある人工知能の研究領域の中で、極めて重要な部分として知られており、世界中の研究者がいろいろな議論を行なっています。代表的な方法としては、

  • 目で見て探す（ → 2.7節）
  • Out-of-Bag Error法（ → 2.8節）
  • クラス内分散・クラス間分散比
  • ベイズ誤差率
  • ReliefF
  • t検定

などがあります。すべて説明すると大変なので、上の2つだけ、以下で解説します。

2.7 有効な特徴量の探索（自分で探そう）

まず、目で見て良い特徴量を探す、ということをやってみます。 ただ単に、散布図を書くだけです。

Programming 05. 目で見て有効な特徴量を探す

In [6]:
fig = plt.figure(figsize=(12, 4))

plt.subplot(1,2,1)
plt.scatter(X_train[ 0:30,0], X_train[ 0:30,1], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[30:60,0], X_train[30:60,1], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[60:90,0], X_train[60:90,1], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
plt.title("Mean(Xacc)-Mean(Yacc)")
plt.xlabel("Mean(Xacc)")
plt.ylabel("Mean(Yacc)")
plt.legend(["Class N", "Class A1", "Class A2"], loc="upper right")
plt.grid(True)

plt.subplot(1,2,2)
plt.scatter(X_train[ 0:30,4], X_train[ 0:30,5], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[30:60,4], X_train[30:60,5], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[60:90,4], X_train[60:90,5], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
plt.title("Mean(Yang)-Mean(Zang)")
plt.xlabel("Mean(Yang)")
plt.ylabel("Mean(Zang)")
plt.legend(["Class N", "Class A1", "Class A2"], loc="upper right")
plt.grid(True)

試しに以下の2つの特徴量を用いて、3クラスの散布図を描いてみました。

  • 左側の散布図（X軸加速度の平均値 - Y軸加速度の平均値）
  • 右側の散布図（Y軸角速度の平均値 - Z軸角速度の平均値）

どちらの空間も、Class A1（異常1）の検出はできそうな感じはしますね。 このような作業をいろいろな特徴量の組み合わせで行うと、Class N（正常）、A1（異常1）、A2（異常2）の3クラスをうまく分離できる特徴量が何となく見えてきます。

とはいえ、目で見て探すのは大変ですし、今回は2次元の空間（あるいはせいぜい3次元）しか可視化することができないので、限界もあります。例えば今回は、平均、分散というわずか2つの特徴量しか採用していないため、12次元（2種類$\times$6軸）となりますが、歪度や周波数領域の特徴量を採用したりすると、100次元、1000次元と、どんどん特徴量の数が増えてしまいます。このような大量の特徴量から、良い特徴量を目で見て探すことは、ほぼ不可能になってしまいます。

2.8 有効な特徴量の探索（アルゴリズムで探そう）

ここでは、Out-of-Bag Error法と呼ばれる、有効な特徴量の探索方法について述べます。 これは、以下の手順を踏みます。

  • 12次元すべての特徴量を用いて、人工知能を構築する。
  • 重要度を調べたいある1つの特徴量の値をランダムに乱す。
  • 誤差の増加度を算出する。
  • これを12個の特徴量すべてについて実施する。
  • 誤差の増加度が高い特徴量を重要な特徴量と解釈する。

注意: かなりざっくり書いています。詳細を知りたい人はブートストラップ法などを勉強すると良いです。

Programming 06. Out-of-Bag Error法による特徴量重要度評価

pythonでは、scikit-learn上でこの手法がパッケージ化されているので、簡単に実装できます。

In [7]:
from sklearn.ensemble import RandomForestClassifier
RF = RandomForestClassifier(n_estimators=300, random_state=1)
RF.fit(X_train, Y_train)
FeatureImportance = RF.feature_importances_

label = ["Mean(Xacc)", "Mean(Yacc)", "Mean(Zacc)", "Mean(Xang)", "Mean(Yang)", "Mean(Zang)", 
         "Var(Xacc)", "Var(Yacc)", "Var(Zacc)", "Var(Xang)", "Var(Yang)", "Var(Zang)"]
left = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
plt.barh(left, width=FeatureImportance, tick_label=label, align="center")
plt.xlabel("Value")
plt.ylabel("Feature's name")
plt.title("Feature importances")
plt.show()

値が大きいほど、その特徴量は有効であることを示しています。 このように、目で見て探す方法とは異なり、一瞬で、有効な特徴量は何か、定量的に知ることができます。

今回は12次元の特徴量なので、あまり恩恵を感じにくいのですが、100次元、1000次元といった莫大な特徴量の中から有効な特徴量を探索しなければならないときに、大変大きな恩恵を感じることができます。

一般的には、重要度が高い特徴量上位6〜7本を使おう、というように、ある程度の次元を保持させます。 ただし今回はわかりやすさ重視のため、2個の特徴量のみを使用することにします。すなわち、

  • 1番目に重要な特徴量: Var(Yacc) Y軸加速度の分散
  • 2番目に重要な特徴量: Var(Xacc) X軸加速度の分散

だけを使ってみます。さて、実際にこの2つの特徴量で正常、異常1、異常2の状態が分類できそうか、散布図により確認してみます。

Programming 07. Out-of-Bag Error法において最適となる特徴量空間の可視化

In [8]:
import matplotlib.ticker as ptick

fig = plt.figure(figsize=(12, 4))

ax=plt.subplot(1,2,1)
plt.scatter(X_train[ 0:30,6], X_train[ 0:30,7], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[30:60,6], X_train[30:60,7], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[60:90,6], X_train[60:90,7], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
plt.title("Var(Xacc)-Var(Yacc)")
plt.xlabel("Var(Xacc)")
plt.ylabel("Var(Yacc)")
plt.legend(["Class N", "Class A1", "Class A2"], loc="upper right")
plt.grid(True)
ax.yaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.xaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.ticklabel_format(style="sci",  axis="y",scilimits=(0,0))
ax.ticklabel_format(style="sci",  axis="x",scilimits=(0,0))

# 左下を拡大した散布図
ax=plt.subplot(1,2,2)
plt.scatter(X_train[ 0:30,6], X_train[ 0:30,7], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[30:60,6], X_train[30:60,7], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[60:90,6], X_train[60:90,7], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
plt.title("Var(Xacc)-Var(Yacc)")
plt.xlabel("Var(Xacc)")
plt.ylabel("Var(Yacc)")
plt.legend(["Class N", "Class A1", "Class A2"], loc="upper left")
plt.grid(True)
plt.xlim(0, 900000)
plt.ylim(0, 900000)
ax.yaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.xaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.ticklabel_format(style="sci",  axis="y",scilimits=(0,0))
ax.ticklabel_format(style="sci",  axis="x",scilimits=(0,0))
  • 左側: Out-of-Bag Error法における最適な特徴量空間における散布図
  • 右側: 左下が見にくいので、拡大したもの

見てみると、Class N（正常）、A1（異常1）、A2（異常2）がよく分離されていることが確認できます。Out-of-Bag Error法により、有効な特徴量を発見できたことが確認されました。そのため、有効な特徴量を取り出します。

注意: 今回は簡単な例なので、わざわざこんなことをしなくても、分散に差がありそう、ということは信号を見ればわかる人はわかるかもしれませんね。

Programming 08. 最適な2つの特徴量のみを抽出

In [9]:
X_train_optimal = X_train[:,6:8]
X_test_optimal = X_test[:,6:8]
# pythonの仕様につき、6:8とすると、６列目以上８列目未満、すなわち、6列目と7列目の2列分が取り出されます。
# 6列目: X軸加速度の分散
# 7列目: Y軸加速度の分散

2.9 人工知能の構築

分離に有効な特徴量を2つ発見できたので、試しにこの2つの特徴量を用いて、人工知能を構築してみます。いろんな種類がありますが、今回は以下の手法を用いてみます。

  • 決定木

    • If-Then節のルールベースにより、クラス分離基準を知ることができます。メリットとして、なぜそのように判断されるのか、人間が解釈しやすいという特徴があります。一方、空間内に軸に対して垂直・平行な識別境界しか生成できないので、後述する方法よりもやや性能が落ちる傾向にあります。また、有効な特徴量を自動探索してくれるアルゴリズムも備わっています。

  • ニューラルネットワーク

    • 人間の脳を模倣したモデルとして、有名です。大変良い性能を持つ可能性がある一方で、識別境界の形状を規定するパラメータが存在しないので、明らかにおかしい識別境界を生成してしまう場合があります。教師データが少ない場合に良い性能を出そうとする仕組みが導入されていないので、教師データが足りない場合に使うことはオススメしません。なお、有効な特徴量を自動探索するアルゴリズムは備わっていません。

  • サポートベクターマシン

    • 特徴量空間上で分離したいクラスに対して、あるクラスと別のクラスのちょうど真ん中に識別境界を生成しようとします。歴史的に、ニューラルネットワークの不安定さを解消するモデルとして、サポートベクターマシンが登場し、脚光を浴びました。ニューラルネットワークとは対照的に、少ない教師データでも良い識別境界を獲得しようとする仕組みが導入されているので、教師データが足りない場合に使うと良いです。なお、有効な特徴量を自動探索するアルゴリズムは備わっていません。

これらはすべて、n次元の空間に対し、$n-1$次元の識別境界を生成する手法となります。数理的な解釈としては、それ以外のなにものでもありません。今回は特徴量を2つ使用するので、2次元の空間に対し、1次元の線を生成する手法となります。

2.9.1 決定木

それでは、決定木を構築するプログラミングについて説明します。

Programming 09. 決定木の構築

In [10]:
# 決定木の学習
from sklearn import tree
Model_Tree_01 = tree.DecisionTreeClassifier(min_samples_leaf=1, random_state=15)
Model_Tree_01 = Model_Tree_01.fit(X_train_optimal, Y_train)

# 木構造の可視化
import pydotplus
dot_data = 'Tree.dot'
tree.export_graphviz(Model_Tree_01, out_file=dot_data)
graph = pydotplus.graphviz.graph_from_dot_file(dot_data)
from IPython.display import Image
Image(graph.create_png())
Out[10]:

valueというのは、ソートされたクラスラベル、左から順に（異常1, 異常2, 正常）の数となります。 一番上をルートノードと言います。今回は教師データが30信号$\times$3クラス=90信号分なので、value=[30, 30, 30]となっています。 これを、特徴量とその閾値による条件がTrue（正しい）かFalse（正しくない）の分岐により、90信号がどんどん分岐されていきます。

  • X[0]: X軸加速度の分散値
  • X[1]: Y軸加速度の分散値

を意味します。ちょっとわかりにくいと思うので、構築した決定木の分岐基準を言語化してみます。

  • X軸加速度の分散が5.5万以下　ならば クラスN（正常）
  • X軸加速度の分散が5.5万より大きい　かつ　Y軸加速度の分散が39万以下　ならば　クラスA2（異常2）
  • X軸加速度の分散が5.5万より大きい　かつ　Y軸加速度の分散が39万より大きい　ならば　クラスA1（異常1）

このように、人間がわかりやすい判断基準を得ることができました。 これでもわかりやすいのですが、特徴量空間に実際にどのような識別境界が得られたのか可視化してみます。

Programming 10. 決定木による識別境界の可視化

In [11]:
Xmin, Ymin, Xmax, Ymax = 0, 0, 900000, 900000 # 空間の最小最大値
resolution = 10000
x_mesh, y_mesh = np.meshgrid(np.arange(Xmin, Xmax, resolution),
                             np.arange(Ymin, Ymax, resolution))
MeshDat = np.array([x_mesh.ravel(), y_mesh.ravel()]).T
z = Model_Tree_01.predict(MeshDat) # 決定木に判別させる
z = np.reshape(z, (len(x_mesh), len(y_mesh))) # 空間のためにデータ整形

# 決定木の判定結果
ax = plt.subplot(1,1,1)
plt.scatter(x_mesh[z=='N'], y_mesh[z=='N'], s=5, alpha=0.3, c='blue')
plt.scatter(x_mesh[z=='A1'], y_mesh[z=='A1'], s=5, alpha=0.3, c='orange')
plt.scatter(x_mesh[z=='A2'], y_mesh[z=='A2'], s=5, alpha=0.3, c='red')

# 教師データの散布図
plt.scatter(X_train[ 0:30,6], X_train[ 0:30,7], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[30:60,6], X_train[30:60,7], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[60:90,6], X_train[60:90,7], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
plt.title("Var(Xacc)-Var(Yacc)")
plt.xlabel("Var(Xacc)")
plt.ylabel("Var(Yacc)")
plt.legend(["Class N", "Class A1", "Class A2"], loc="upper left")
plt.grid(True)
plt.xlim(0, 900000)
plt.ylim(0, 900000)
ax.yaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.xaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.ticklabel_format(style="sci",  axis="y",scilimits=(0,0))
ax.ticklabel_format(style="sci",  axis="x",scilimits=(0,0))

一つ一つのプロットが、教師データです。カラフルな領域が、教師データを使用して生成した識別境界となります。 2次元の空間上（面）に1次元の識別境界（線）ができていますね。 それはともかく、決定木により生成した識別境界は、教師データをすべて正しく分類できています。 とはいえ、X軸方向の揺れのばらつきが小さければ、Y軸方向の揺れのばらつきがいくら大きくても正常と判定されるようです。 ちょっと微妙な感じもしますね。

2.9.2 ニューラルネットワーク

続いてニューラルネットワークです。こちらも設定すべきハイパーパラメータがたくさんありますが、今回は学習回数のみをセットします。 そして、以下の設定でモデルを学習させてみます。

  • 入力層（第1層）：　2次元（X軸加速度の分散、Y軸加速度の分散）
  • 隠れ層（第2層）：　2次元
  • 出力層（第3層）：　3次元（正常、異常1、異常2、それぞれの確信度みたいな3つの数字）
  • 学習回数：　5000回
  • 活性化関数：　シグモイド関数

簡単にいうと、入力された2つの特徴量が、別の2次元に変換され、正常、異常1、異常2のどれらしいかを3次元ベクトルとして出力する2段階の関数となります。変換自体は高校3年生の数学で学ぶ程度のことで、ただの足し算、掛け算、割り算、指数関数の組み合わせなので、かなり簡単です。パラメータの学習は大学1年生くらいの数学で学ぶ内容がベースです（偏微分を使います）。

直にニューラルネットワークをプログラムで書かせる場合、200行程度必要ですが、pythonでは簡単にプログラムを書くことができます。

Programming 11. ニューラルネットワークの学習

In [12]:
from sklearn.neural_network import MLPClassifier
import warnings
warnings.filterwarnings('ignore')
Model_NN_01 = MLPClassifier(hidden_layer_sizes=(2, ), solver="sgd", activation="logistic", learning_rate_init=0.01, 
                             random_state=1, max_iter=5000, alpha=0.1, tol = pow(10, -10), verbose=False)
Model_NN_01 = Model_NN_01.fit(X_train_optimal/100000, Y_train)

これで学習が終わりです。決定木のときと同様に、2次元の空間に対し、どんな1次元が生成されたのか、確認してみます。

Programming 12. ニューラルネットワークによる識別境界の可視化

In [13]:
Xmin, Ymin, Xmax, Ymax = 0, 0, 900000, 900000 # 空間の最小最大値
resolution = 10000
x_mesh, y_mesh = np.meshgrid(np.arange(Xmin, Xmax, resolution),
                             np.arange(Ymin, Ymax, resolution))
MeshDat = np.array([x_mesh.ravel(), y_mesh.ravel()]).T
MeshDat = MeshDat / 100000
z = Model_NN_01.predict(MeshDat)
z = np.reshape(z, (len(x_mesh), len(y_mesh)))

# ニューラルネットの判定結果
ax = plt.subplot(1,1,1)
plt.scatter(x_mesh[z=='N'], y_mesh[z=='N'], s=5, alpha=0.3, c='blue')
plt.scatter(x_mesh[z=='A1'], y_mesh[z=='A1'], s=5, alpha=0.3, c='orange')
plt.scatter(x_mesh[z=='A2'], y_mesh[z=='A2'], s=5, alpha=0.3, c='red')

# 教師データの散布図
plt.scatter(X_train[ 0:30,6], X_train[ 0:30,7], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[30:60,6], X_train[30:60,7], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[60:90,6], X_train[60:90,7], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
plt.title("Var(Xacc)-Var(Yacc)")
plt.xlabel("Var(Xacc)")
plt.ylabel("Var(Yacc)")
plt.legend(["Class N", "Class A1", "Class A2"], loc="upper left")
plt.grid(True)
plt.xlim(0, 900000)
plt.ylim(0, 900000)
ax.yaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.xaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.ticklabel_format(style="sci",  axis="y",scilimits=(0,0))
ax.ticklabel_format(style="sci",  axis="x",scilimits=(0,0))

今回はどうでしょうか。決定木での識別境界のデメリットが大きく緩和されていますね。 どちらの分散も小さいときのみ、正常判定がなされるという、割と妥当な識別境界を得ることができました。

ニューラルネットワークの学習は、数式内部にある頭の良さを操作する変数を、微分により、プラスマイナスどちらにずれしたらいいか？という判断を行って、実際にその方向に変化させる、という処理の反復により実現されます。すなわちちょっとずつ頭が良くなる、ということです。これを動画してみます。

Programming 13. ニューラルネットワークによる識別境界の可視化（動画）

以下のコードは、ニューラルネットワークの学習状況、すなわち識別境界の変化を動画ファイルにして、出力するコードです。印刷された紙では特に意味がありませんので、ご注意ください。

In [14]:
import matplotlib.animation as animation

NN_Frame = []
NN_Movie = []
fig = plt.figure()
ax = plt.subplot(1,1,1)
for i in range(1,50):
    NN_Movie.append(MLPClassifier(hidden_layer_sizes=(2, ), solver="sgd", activation="logistic", 
                                  learning_rate_init=0.01, random_state=1, max_iter=i*10, alpha=0.1,
                                  tol = pow(10, -10), verbose=False))
    NN_Movie[i-1].fit(X_train_optimal/100000, Y_train)
    
    z = NN_Movie[i-1].predict(MeshDat)
    z = np.reshape(z, (len(x_mesh), len(y_mesh)))

    # ニューラルネットの判定結果
    im1 = plt.scatter(x_mesh[z=='N'], y_mesh[z=='N'], s=5, alpha=0.3, c='blue')
    im2 = plt.scatter(x_mesh[z=='A1'], y_mesh[z=='A1'], s=5, alpha=0.3, c='orange')
    im3 = plt.scatter(x_mesh[z=='A2'], y_mesh[z=='A2'], s=5, alpha=0.3, c='red')

    # 教師データの散布図
    im4 = plt.scatter(X_train[ 0:30,6], X_train[ 0:30,7], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
    im5 = plt.scatter(X_train[30:60,6], X_train[30:60,7], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
    im6 = plt.scatter(X_train[60:90,6], X_train[60:90,7], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
    
    plt.title("Var(Xacc)-Var(Yacc)")
    plt.xlabel("Var(Xacc)")
    plt.ylabel("Var(Yacc)")
    plt.legend(["Class N", "Class A1", "Class A2"], loc="upper left")
    plt.grid(True)
    plt.xlim(0, 900000)
    plt.ylim(0, 900000)
    ax.yaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
    ax.xaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
    ax.ticklabel_format(style="sci",  axis="y",scilimits=(0,0))
    ax.ticklabel_format(style="sci",  axis="x",scilimits=(0,0)) 
    
    # フレームを積む
    NN_Frame.append([im1, im2, im3, im4, im5, im6]) 
                    
ani = animation.ArtistAnimation(fig, NN_Frame, interval=500)
ani.save('NN_Movie.mp4', writer="ffmpeg")

2.9.3 サポートベクターマシン

続いて、サポートベクターマシンに移ります。こちらは、クラス間の真ん中を通るように識別境界を生成するモデルです。この特性により、教師データが少なくても比較的良い識別境界を生成することができます。サポートベクターマシンには、大きく分けて2つのモデルがあります。

  • 線形サポートベクターマシン: Linear SVM

    • 空間に対し、まっすぐな識別境界を引くモデルとなります。今回は2次元空間なので、1次元の直線という形で識別境界が生成されます。識別境界を曲げることができないので精度が悪くなりがちというデメリットがある一方、教師データ内に正解してはいけない外れ値が含まれていたとき、それを無理に正解しようとしないので、安定するというメリットがあります。このような性質を、人工知能の分野では低バリアンス性といいます。

  • 非線形サポートベクターマシン Non-linear SVM

    • 空間に対し、曲がった識別境界を生成できるモデルとなります。今回は2次元空間なので、1次元の曲線という形で識別境界を生成します。複雑な識別境界を生成できるため良い精度が得られる場合があるというメリットがある一方、分析者が識別境界の形状を制御しにくかったり、過学習が発生しやすかったり、本来当ててはいけない問題も解くような境界を得てしまうというデメリットがあります。このような性質を、人工知能の分野では高バリアンス性と言います。

    どっちがいいんだということは自分が作りたいモデルによって異なるので、今回は両方作ってみます。

Programming 14. 線形・非線形サポートベクターマシンの学習

In [15]:
from sklearn import svm

# 線形サポートベクターマシン
SVMModel_Lin01 = svm.SVC(kernel='linear', C=0.1)
SVMModel_Lin01 = SVMModel_Lin01.fit(X_train_optimal, Y_train)

#非線形サポートベクターマシン
SVMModel_RBF01 = svm.SVC(kernel='rbf', C=10, gamma=0.00001)
SVMModel_RBF01 = SVMModel_RBF01.fit(X_train_optimal/1000, Y_train)

これで学習が終わりです。これまでと同じように、どのような1次元が生成されたのか、みてみます。

Programming 15. サポートベクターマシンによる識別境界の可視化

In [16]:
Xmin, Ymin, Xmax, Ymax = 0, 0, 900000, 900000 # 空間の最小最大値
resolution = 10000
x_mesh, y_mesh = np.meshgrid(np.arange(Xmin, Xmax, resolution),
                             np.arange(Ymin, Ymax, resolution))
MeshDat = np.array([x_mesh.ravel(), y_mesh.ravel()]).T

# メッシュデータの推定
z_Lin01 = SVMModel_Lin01.predict(MeshDat) # モデル1
z_RBF01 = SVMModel_RBF01.predict(MeshDat/1000) # モデル2

# データ整形
z_Lin01 = np.reshape(z_Lin01, (len(x_mesh), len(y_mesh))) 
z_RBF01 = np.reshape(z_RBF01, (len(x_mesh), len(y_mesh))) 

fig = plt.figure(figsize=(12, 4))

# 線形SVMの可視化
ax = plt.subplot(1,2,1)
plt.scatter(x_mesh[z_Lin01=='N'], y_mesh[z_Lin01=='N'], s=5, alpha=0.3, c='blue')
plt.scatter(x_mesh[z_Lin01=='A1'], y_mesh[z_Lin01=='A1'], s=5, alpha=0.3, c='orange')
plt.scatter(x_mesh[z_Lin01=='A2'], y_mesh[z_Lin01=='A2'], s=5, alpha=0.3, c='red')

plt.scatter(X_train[ 0:30,6], X_train[ 0:30,7], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[30:60,6], X_train[30:60,7], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[60:90,6], X_train[60:90,7], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
plt.title("Linear SVM: Var(Xacc)-Var(Yacc)")
plt.xlabel("Var(Xacc)")
plt.ylabel("Var(Yacc)")
plt.legend(["Class N", "Class A1", "Class A2"], loc="upper left")
plt.grid(True)
plt.xlim(0, 900000)
plt.ylim(0, 900000)
ax.yaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.xaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.ticklabel_format(style="sci",  axis="y",scilimits=(0,0))
ax.ticklabel_format(style="sci",  axis="x",scilimits=(0,0))

# 非線形SVMの可視化
ax = plt.subplot(1,2,2)
plt.scatter(x_mesh[z_RBF01=='N'], y_mesh[z_RBF01=='N'], s=5, alpha=0.3, c='blue')
plt.scatter(x_mesh[z_RBF01=='A1'], y_mesh[z_RBF01=='A1'], s=5, alpha=0.3, c='orange')
plt.scatter(x_mesh[z_RBF01=='A2'], y_mesh[z_RBF01=='A2'], s=5, alpha=0.3, c='red')

plt.scatter(X_train[ 0:30,6], X_train[ 0:30,7], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[30:60,6], X_train[30:60,7], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[60:90,6], X_train[60:90,7], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
plt.title("Non-linear SVM: Var(Xacc)-Var(Yacc)")
plt.xlabel("Var(Xacc)")
plt.ylabel("Var(Yacc)")
plt.legend(["Class N", "Class A1", "Class A2"], loc="upper left")
plt.grid(True)
plt.xlim(0, 900000)
plt.ylim(0, 900000)
ax.yaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.xaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.ticklabel_format(style="sci",  axis="y",scilimits=(0,0))
ax.ticklabel_format(style="sci",  axis="x",scilimits=(0,0))
  • 左側: 線形サポートベクターマシン
  • 右側: 非線形サポートベクターマシン

どうでしょうか。意見が割れそうです。なんとなく、非線形サポートベクターマシンの方が良さそうですね（右側）。あと、ニューラルネットワークよりも良さそうな識別境界が得られていますね。先ほどサポートベクターマシンの説明で行なった、あるクラスと別のクラスの、ちょうど真ん中に識別境界が生成されている、ということがわかるでしょうか。一般に、サポートベクターマシンの方がニューラルネットワークよりも良い識別境界が得られると言われる所以が、この考え方となります。

ところで、先ほど非線形サポートベクターマシンは複雑な識別境界を引けるので、デメリットとして過学習が生じやすいと言いましたね。実際に、わざと過学習を起こさせたモデルを作ってみます。

Programming 16. 過学習を起こさせた非線形サポートベクターマシンによる識別境界の可視化

In [17]:
# 非線形サポートベクターマシンの学習（過学習）
SVMModel_RBF02_OverLearning = svm.SVC(kernel='rbf', C=10, gamma=0.001)
SVMModel_RBF02_OverLearning = SVMModel_RBF02_OverLearning.fit(X_train_optimal/1000, Y_train)

Xmin, Ymin, Xmax, Ymax = 0, 0, 900000, 900000 # 空間の最小最大値
resolution = 10000
x_mesh, y_mesh = np.meshgrid(np.arange(Xmin, Xmax, resolution),
                             np.arange(Ymin, Ymax, resolution))
MeshDat = np.array([x_mesh.ravel(), y_mesh.ravel()]).T
MeshDat = MeshDat / 1000
z = SVMModel_RBF02_OverLearning.predict(MeshDat)
z = np.reshape(z, (len(x_mesh), len(y_mesh)))

# ニューラルネットの判定結果
ax = plt.subplot(1,1,1)
plt.scatter(x_mesh[z=='N'], y_mesh[z=='N'], s=5, alpha=0.3, c='blue')
plt.scatter(x_mesh[z=='A1'], y_mesh[z=='A1'], s=5, alpha=0.3, c='orange')
plt.scatter(x_mesh[z=='A2'], y_mesh[z=='A2'], s=5, alpha=0.3, c='red')

# 教師データの散布図
plt.scatter(X_train[ 0:30,6], X_train[ 0:30,7], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[30:60,6], X_train[30:60,7], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
plt.scatter(X_train[60:90,6], X_train[60:90,7], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
plt.title("Var(Xacc)-Var(Yacc)")
plt.xlabel("Var(Xacc)")
plt.ylabel("Var(Yacc)")
plt.legend(["Class N", "Class A1", "Class A2"], loc="upper left")
plt.grid(True)
plt.xlim(0, 900000)
plt.ylim(0, 900000)
ax.yaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.xaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.ticklabel_format(style="sci",  axis="y",scilimits=(0,0))
ax.ticklabel_format(style="sci",  axis="x",scilimits=(0,0))

どうでしょうか。何か様子がおかしいことがわかると思います。

  • 異常1（Class A1, 黄色）の領域が広い
  • 異常2（Class A2, 赤色）の領域が狭い

果たして、実際にこのモデルを使用したとして、異常2をうまく検出できるでしょうか？
おそらくですが、できない場合が多いと思います。 ですがこのモデル、学習に使用した教師データをすべて正しく回答しています。これを見ればわかるように、人工知能の分野においては、学習に使用した教師データに正答することは、あまり意味のあることではないのです。
ここではわざと、サポートベクターマシンで過学習を発生させました。しかし、決定木でもニューラルネットワークでも過学習は発生させられます。なので、サポートベクターマシンは他のモデルより劣っているとは決して思わないようにしてください。

2.10 精度評価: 適合率（Precision）　VS　再現率（Recall）

分類問題における精度評価の指標として、一般に、適合率、再現率、F値を用います。

  • 適合率 Precision（最悪:0, 最良:1）
    • ターゲットクラスにおいて、予測されたクラスが正しい確率のこと。この値が高いほど、出力された値が信用してもいいことを意味する。
  • 再現率 Recall（最悪:0, 最良:1）
    • ターゲットクラスにおいて、実際のクラスを正しく予測した確率のこと。この値が高いほど、ターゲットクラスを取り逃がさないことを意味する。
  • F値（最悪:0, 最良:1）
    • ターゲットクラスにおける、適合率と再現率の相乗平均のこと。

適合率と再現率、一見似ていますがまったく別物の指標です。少し、頭の中で以下の事例を考えてみましょう。 糖尿病か否かを判定するモデルがあります。8人の被験者に対し、このモデルの実測と予測を比較したところ、以下の結果となりました。0が糖尿病ではない、1が糖尿病である、ということを意味します。「適合率・再現率」は「高い・低い」、どうなるでしょうか。

  • 実測: 0 0 0 0 1 1 1 1
  • 予測: 0 0 0 0 0 0 1 0

ターゲットクラスを「糖尿病である」すなわち「1」と考えてみます。こうした場合、予測モデルが糖尿病だと出力した値は、100%信用できますね。したがって、「このモデルの適合率は高い」と判断できます。一方、実際の糖尿病患者は4名おり、3名取り漏らしています。ゆえに、「このモデルの再現率は低い」と判断できます。では、以下の事例ではどうでしょうか。

  • 実測: 0 0 0 0 1 1 1 1
  • 予測: 1 1 1 1 1 1 1 1

超極端な例です。この予測モデルの糖尿病だという回答は、全然信用できません。ゆえに「このモデルの適合率は低い」と判断されます。一方、実際の糖尿病患者を取り漏らしているでしょうか。一切取り漏らしていません。ゆえに「このモデルの再現率は高い」と判断できます。

では、適合率と再現率、どちらが高いモデルがいいのでしょうか。これは扱う問題の性質によって異なります。具体例として、火災報知器を例に考えてみましょう。あれは適合率と再現率、どちらが高く、どちらが低いでしょうか。火災報知器がブーブー鳴っても、実際に火事であることはあまりない気がします。しかし火災報知器は、実際の火事を取り漏らすことはほとんどありません。つまり、適合率が低く、再現率が高い、というのが正解です。詳しいことはわかりませんが、ああいう類の問題は、再現率を極めて高くしよう、という思想が備わっています。このように、適合率と再現率、重視すべきはいずれかという問題は、状況によって異なることを覚えておくといいでしょう。

ところで、察しのいい人は気がつくかもしれませんが、再現率を高めようとすると、適合率が下がります。逆に、適合率を高めようとすると、再現率が下がります。実は、適合率と再現率にはトレードオフの関係があります。科学というのはなかなか辛辣で、うまい話はないのです。 ...とは言ったものの、火災報知器のような特殊な例を除いて、適合率も再現率も、どちらも重要です。したがって、特にどちらを大事にしたいという考えがない場合には、適合率と再現率の相乗平均であるF値をモデルの精度評価指標として採用することが普通です。論文などでどのモデルがいいのか議論したい場合には、F値を使うといいでしょう。

とはいえ、何でもかんでもF値を見るのではなく、適合率・再現率もしっかりと観察するようにしましょう。知能が有する個性を明瞭に知る一つの手がかりとなります。

ところで、先ほど教師データの精度には重要な意味はないと言いました。通常、人工知能の精度評価には学習には使用しなかったテストデータを使用します。それでは、テストデータを用いて、今まで構築してきたモデルの適合率、再現率、F値を算出してみます。

Programming 17. テストデータによる精度評価（適合率、再現率、F値）

In [18]:
from sklearn.metrics import precision_recall_fscore_support

#　決定木の精度評価
Y_Tr = Model_Tree_01.predict(X_test_optimal) # <- 今回はテストデータを使用
Res_Tr = precision_recall_fscore_support(Y_Tr, Y_test, average='weighted')
Res_Trw = []
Res_Trw.append(np.round(Res_Tr[0], 3))
Res_Trw.append(np.round(Res_Tr[1], 3))
Res_Trw.append(np.round(Res_Tr[2], 3))
print("決定木:           適合率P =",Res_Trw[0], "/再現率R =",Res_Trw[1], "/F値 =",Res_Trw[2])

# ニューラルネット 
Y_NN = Model_NN_01.predict(X_test_optimal/100000) # <- 今回はテストデータを使用
Res_NN = precision_recall_fscore_support(Y_NN, Y_test, average='weighted')
Res_NNw = []
Res_NNw.append(np.round(Res_NN[0], 3))
Res_NNw.append(np.round(Res_NN[1], 3))
Res_NNw.append(np.round(Res_NN[2], 3))
print("ニューラルネット: 適合率P =",Res_NNw[0], "/再現率R =",Res_NNw[1], "/F値 =",Res_NNw[2])

# 線形SVM
Y_SVML = SVMModel_Lin01.predict(X_test_optimal) # <- 今回はテストデータを使用
Res_SVML = precision_recall_fscore_support(Y_SVML, Y_test, average='weighted')
Res_SVMLw = []
Res_SVMLw.append(np.round(Res_SVML[0], 3))
Res_SVMLw.append(np.round(Res_SVML[1], 3))
Res_SVMLw.append(np.round(Res_SVML[2], 3))
print("線形SVM:          適合率P =",Res_SVMLw[0], "/再現率R =",Res_SVMLw[1], "/F値 =",Res_SVMLw[2])

# 非線形SVM
Y_SVMR = SVMModel_RBF01.predict(X_test_optimal/1000) # <- 今回はテストデータを使用
Res_SVMR = precision_recall_fscore_support(Y_SVMR, Y_test, average='weighted')
Res_SVMRw = []
Res_SVMRw.append(np.round(Res_SVMR[0], 3))
Res_SVMRw.append(np.round(Res_SVMR[1], 3))
Res_SVMRw.append(np.round(Res_SVMR[2], 3))
print("非線形SVM:        適合率P =",Res_SVMRw[0], "/再現率R =",Res_SVMRw[1], "/F値 =",Res_SVMRw[2])

# 非線形SVM　過学習
Y_SVMROL = SVMModel_RBF02_OverLearning.predict(X_test_optimal/1000) # <- 今回はテストデータを使用
Res_SVMROL = precision_recall_fscore_support(Y_SVMROL, Y_test, average='weighted')
Res_SVMROLw = []
Res_SVMROLw.append(np.round(Res_SVMROL[0], 3))
Res_SVMROLw.append(np.round(Res_SVMROL[1], 3))
Res_SVMROLw.append(np.round(Res_SVMROL[2], 3))
print("非線形SVM(過学習):適合率P =",Res_SVMROLw[0], "/再現率R =",Res_SVMROLw[1], "/F値 =",Res_SVMROLw[2])

決定木:           適合率P = 0.945 /再現率R = 0.93 /F値 = 0.931
ニューラルネット: 適合率P = 0.945 /再現率R = 0.93 /F値 = 0.931
線形SVM:          適合率P = 0.945 /再現率R = 0.93 /F値 = 0.931
非線形SVM:        適合率P = 0.945 /再現率R = 0.93 /F値 = 0.931
非線形SVM(過学習):適合率P = 0.919 /再現率R = 0.86 /F値 = 0.866

各指標とも、最良が1、最悪が0です。過学習させたモデルだけ、精度が少し低い状態にあります。他のモデルはすべて同い精度で、値も高いですね。本当は差をつけたかったのですが、今回は異常値検出としての問題が簡単すぎたので、すべて高い精度となってしまいました。セミナーとしてはあまり説明しやすい結果ではありませんが、自然にやるとこんな感じである、ということでしょうか。

3. 深層学習による人工知能の構築

これまで、いわゆる古典的方法による人工知能の構築方法を説明しました。次に、最近流行りの深層学習について述べます。深層学習と古典的な人工知能の間には、以下の違いがあります。

  • 古典的な人工知能の構築方法:
    • 分析者が生の信号を特徴量に変換する。（先ほどは、3軸加速度・角速度の平均と分散を使いました）
    • 有効な特徴量を探します。（先ほどは、Out-of-Bag Error法で探しました。X, Y軸加速度の分散が重要でした）
    • 特徴量空間内に、識別境界を生成する。
  • 深層学習による人工知能の構築方法:
    • 生の波形から、人工知能が有効な特徴量を生成する。
    • その特徴量空間内に、識別境界を生成する。

違いがわかったでしょうか。 深層学習のメリット、デメリットは以下の点にあります。

  • 深層学習のメリット
    • 教師データを莫大に用意できたとき、凄まじい性能となる。
    • 特徴量探索を行う必要がない。
  • 深層学習のデメリット
    • 教師データが少ないと、古典的な人工知能よりもかなり低い性能になる。
    • 特徴量学習の段階で過学習が生じることがある。
    • 学習された特徴量に物理的意味が見出せないので、なぜその判別結果になるのか説明できない。

実際に使うようになるとわかりますが、古典的な人工知能の方が良い性能を持つことも多いです。

3.1 深層学習による特徴量学習

古典的な人工知能と深層学習の違いがわかったところで、深層学習を構築してみます。ノートパソコンなので、かなり単純なモデルにしてみます。詳しい人からは怒られそうなくらい単純です。

Programming 18. 深層学習による特徴量学習

In [19]:
# 深層学習用の教師データ整形
InputDat_N=[]
InputDat_A1=[]
InputDat_A2=[]
for i in range(0, 30):
    InputDat_N.append(np.reshape(Signal_Class_N[i][0:,1:3], [1, 2*100]).ravel())
    InputDat_A1.append(np.reshape(Signal_Class_A1[i][0:,1:3], [1, 2*100]).ravel())
    InputDat_A2.append(np.reshape(Signal_Class_A2[i][0:,1:3], [1, 2*100]).ravel())
X_DL = np.concatenate([InputDat_N, InputDat_A1, InputDat_A2])
X_DL = X_DL/100

# 深層学習モデルの構築
Model_DeepLearning = MLPClassifier(hidden_layer_sizes=(5, 2, ), solver="sgd", activation="logistic", learning_rate_init=0.05, 
                             random_state=10, max_iter=700, alpha=0.01, tol = pow(10, -10), verbose=False)
Model_DeepLearning = Model_DeepLearning.fit(X_DL, Y_train)
#plt.plot(Model_DeepLearning.loss_curve_)

これで学習が終わりました。さて、どのような特徴量が学習されたでしょうか。今回は、わかりやすさ重視なので、深層学習で2つの特徴量を学習させました（普通は300個の特徴量を学習するなど、かなり数が多いです）。2つだと、2次元空間として可視化できます。実際に見てみましょう。

Programming 19. 深層学習により抽出された特徴量の可視化（動画）

In [20]:
import matplotlib.animation as animation
DL_Frame = []
DLM = []
fig = plt.figure()

# 学習回数を25回ずつ増やしながら、深層学習構築、特徴抽出を繰り返す
for i in range(1, int(700/25)):
    
    plt.clf() # ★ 動画ファイルを作成する場合にはコメントアウト

    # 深層学習モデルの構築
    DLM.append(
        MLPClassifier(hidden_layer_sizes=(5, 2, ), solver="sgd", activation="logistic", learning_rate_init=0.05, 
                      random_state=10, max_iter=(i-1)*25+1, alpha=0.01, tol = pow(10, -10), verbose=False)
    )
    DLM[i-1].fit(X_DL, Y_train)
    
    # 学習済みモデルの特徴抽出
    DLF=[]
    for j in range(0, 90):
        temp = 1/(1+np.exp(-1 * np.dot(DLM[i-1].coefs_[0].T, X_DL[j]) + DLM[i-1].intercepts_[0]))
        DLF.append(1/(1+np.exp(-1 * np.dot(DLM[i-1].coefs_[1].T, temp) + DLM[i-1].intercepts_[1])))        
    DLF = np.array(DLF)

    # 特徴量空間の可視化
    im1 = plt.scatter(DLF[0:30, 0], DLF[0:30, 1], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
    im2 = plt.scatter(DLF[30:60, 0], DLF[30:60, 1], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
    im3 = plt.scatter(DLF[60:90, 0], DLF[60:90, 1], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
    plt.title("Deep Feature Space")
    plt.xlabel("Deep Feature 1")
    plt.ylabel("Deep Feature 2")
    plt.legend(["Class N", "Class A1", "Class A2"], loc="upper right")
    plt.grid(True)
    
    # フレームを積む
    DL_Frame.append([im1, im2, im3]) 
                    
#aniDL = animation.ArtistAnimation(fig, DL_Frame, interval=500)
#aniDL.save('DL_Movie.mp4', writer="ffmpeg")

3.2 深層学習の特徴量と古典的な特徴量はどちらが有効なのか？

さて、学習した2つの特徴量はどのくらい重要なのでしょうか。古典的な人工知能の方で説明した平均や分散よりも優れた特徴量を獲得できたのでしょうか。実はこれ、（諸般の理由によりあまり適切ではありませんが。。。）先ほどやったOut-of-Bag Error法で定量的に測定することができます。 今計算している特徴量は、

  • 古典的な人工知能で算出した12個の特徴量（3軸加速度・角速度の平均と分散）：　12次元
  • 深層学習により取得した2つの特徴量：　2次元

この14個です。これら全部で、Out-of-Bag Error法を適用します：

  • 14個の特徴量すべてを使用して、人工知能を構築する。
  • 重要度を調べたいある1つの特徴量の値をランダムに乱す。
  • 誤差の増加度を算出する。
  • これを14個の特徴量すべてについて実施する。
  • 誤差の増加度が高い特徴量を重要な特徴量と解釈する。

Programming 20. Out-of-Bag Error法による特徴量重要度評価（深層学習の特徴量 VS 古典的な特徴量）

In [23]:
# 平均・分散による12次元の空間と、深層学習の2次元空間を結合させる。
X_train_DL = np.concatenate([X_train, DLF],1)

# OOB
RFDL = RandomForestClassifier(n_estimators=300, random_state=1)
RFDL.fit(X_train_DL, Y_train)
FeatureImportanceDL = RFDL.feature_importances_

# 棒グラフ
label = ["Mean(Xacc)", "Mean(Yacc)", "Mean(Zacc)", "Mean(Xang)", "Mean(Yang)", "Mean(Zang)", 
         "Var(Xacc)", "Var(Yacc)", "Var(Zacc)", "Var(Xang)", "Var(Yang)", "Var(Zang)",
         "DL Feature 1", "DL Feature 2"]
left = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
plt.barh(left, width=FeatureImportanceDL, tick_label=label, align="center")
plt.xlabel("Value")
plt.ylabel("Feature's name")
plt.title("Feature importances")
plt.show()

上ふたつが深層学習由来の特徴量の重要度です。かなり重要だと判断されるようです。ただ、よく見るとY軸加速度の分散値の方が有効みたいです。この結果からもわかるように、深層学習は古典的な人工知能構築手法よりも優れている、ということは一概に言えないわけです。

とはいえ、深層学習の特徴量は、Y軸加速度の分散という最重要な特徴量に迫ることができました。このように、深層学習は手軽にいい特徴量を獲得することができます。これはかなりいい点でしょうね。

ところで、よく深層学習と古典的な人工知能の構築は、対立する手法だと思う人も多いかもしれません。しかし、ここ1-2年の人工知能の理論研究では、平均、分散といったこれまで使われてきた特徴量と、深層学習により獲得された特徴量を混ぜた空間を構築することで、どちらか一方のみを使うよりも良い性能が得られるということがわかってきています。このように、深層学習と古典的な手法は、別に対立したりはしないのです。例えば、上の棒グラフを見ると、 - DL Feature 1: 深層学習により得た特徴量 - DL Feature 2: 深層学習により得た特徴量 - Y軸加速度の分散値: 古典的な手法による特徴量 という3次元空間で人工知能を構築すると、すごい良い性能が出そう、というのは簡単に想像できると思います。

3.3 深層学習による特徴量を利用した人工知能の構築

さて、深層学習による特徴量学習について説明しました。 深層学習はただ単に、入力の次元を下げ、空間を作る方法だということは理解できたでしょうか。 同時に、これまで解説した決定木、サポートベクターマシン、ニューラルネットワークなどが、空間上に識別境界を生成する方法だということを理解していれば、深層学習とこれまでの手法が連携できることがわかると思います。 これを目で見て確認するため、以下のモデルを構築してみます。

  • 特徴量: 深層学習、識別境界: ニューラルネットワーク（一般的な深層学習）
  • 特徴量: 深層学習、識別境界: サポートベクターマシン
  • 特徴量: 深層学習、識別境界: 決定木

Programming 21. 深層学習とニューラルネットワークの連携（一般的な深層学習）（動画）

In [24]:
import matplotlib.animation as animation
DL_Frame = []
DLM = []
fig = plt.figure()

# 学習回数を25回ずつ増やしながら、深層学習構築、特徴抽出を繰り返す
for i in range(1, int(700/25)):
    
    plt.clf() # ★注意: 動画にしたい場合には、コメントアウトすること！！
    
    # 深層学習モデルの構築
    DLM.append(
        MLPClassifier(hidden_layer_sizes=(5, 2, ), solver="sgd", activation="logistic", learning_rate_init=0.05, 
                      random_state=10, max_iter=(i-1)*25+1, alpha=0.01, tol = pow(10, -10), verbose=False)
    )
    DLM[i-1].fit(X_DL, Y_train)
    
    # 学習済みモデルの特徴抽出
    DLF=[]
    for j in range(0, 90):
        temp = 1/(1+np.exp(-1 * np.dot(DLM[i-1].coefs_[0].T, X_DL[j]) + DLM[i-1].intercepts_[0]))
        DLF.append(1/(1+np.exp(-1 * np.dot(DLM[i-1].coefs_[1].T, temp) + DLM[i-1].intercepts_[1])))        
    DLF = np.array(DLF)
    
    # 深層学習による特徴量空間上にニューラルネットワークで識別境界を生成
    DeepLearningNN = MLPClassifier(hidden_layer_sizes=(2, ), activation="logistic", learning_rate_init=0.05,
                                   max_iter = 1000, random_state=1)
    DeepLearningNN = DeepLearningNN.fit(DLF, Y_train)
    
    # ニューラルネットによる識別境界の可視化
    Xmin, Ymin, Xmax, Ymax = 0, 0, 1, 1 # 空間の最小最大値
    resolution = 0.01
    x_mesh, y_mesh = np.meshgrid(np.arange(Xmin, Xmax, resolution),
                             np.arange(Ymin, Ymax, resolution))
    MeshDat = np.array([x_mesh.ravel(), y_mesh.ravel()]).T
    z = DeepLearningNN.predict(MeshDat)
    z = np.reshape(z, (len(x_mesh), len(y_mesh)))
    im1 = plt.scatter(x_mesh[z=='N'], y_mesh[z=='N'], s=5, alpha=0.3, c='blue')
    im2 = plt.scatter(x_mesh[z=='A1'], y_mesh[z=='A1'], s=5, alpha=0.3, c='orange')
    im3 = plt.scatter(x_mesh[z=='A2'], y_mesh[z=='A2'], s=5, alpha=0.3, c='red')

    # 深層学習による特徴量空間プロットの可視化
    im4 = plt.scatter(DLF[0:30, 0], DLF[0:30, 1], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
    im5 = plt.scatter(DLF[30:60, 0], DLF[30:60, 1], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
    im6 = plt.scatter(DLF[60:90, 0], DLF[60:90, 1], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
    
    plt.title("Deep Learning + NN")
    plt.xlabel("Deep Feature 1")
    plt.ylabel("Deep Feature 2")
    plt.legend(["Class N", "Class A1", "Class A2"], loc="upper right")
    plt.grid(True)
    
    # フレームを積む
    DL_Frame.append([im1, im2, im3, im4, im5, im6]) 

plt.show()
#aniDLNN = animation.ArtistAnimation(fig, DL_Frame, interval=500)
#aniDLNN.save('DLNN_Movie.mp4', writer="ffmpeg")

Programming 22. 深層学習とサポートベクターマシンの連携（動画）

In [25]:
import matplotlib.animation as animation
DL_Frame = []
DLM = []
fig = plt.figure()

# 学習回数を25回ずつ増やしながら、深層学習構築、特徴抽出を繰り返す
for i in range(1, int(700/25)):
    
    plt.clf() # ★注意: 動画にしたい場合には、コメントアウトすること！！
    
    # 深層学習モデルの構築
    DLM.append(
        MLPClassifier(hidden_layer_sizes=(5, 2, ), solver="sgd", activation="logistic", learning_rate_init=0.05, 
                      random_state=10, max_iter=(i-1)*25+1, alpha=0.01, tol = pow(10, -10), verbose=False)
    )
    DLM[i-1].fit(X_DL, Y_train)
    
    # 学習済みモデルの特徴抽出
    DLF=[]
    for j in range(0, 90):
        temp = 1/(1+np.exp(-1 * np.dot(DLM[i-1].coefs_[0].T, X_DL[j]) + DLM[i-1].intercepts_[0]))
        DLF.append(1/(1+np.exp(-1 * np.dot(DLM[i-1].coefs_[1].T, temp) + DLM[i-1].intercepts_[1])))        
    DLF = np.array(DLF)
    
    # 深層学習による特徴量空間上にサポートベクターマシンで識別境界を生成
    DeepLearningSVM = svm.SVC(kernel='linear', C=100000)
    DeepLearningSVM = DeepLearningSVM.fit(DLF, Y_train)
    
    # サポートベクターマシンによる識別境界の可視化
    Xmin, Ymin, Xmax, Ymax = 0, 0, 1, 1 # 空間の最小最大値
    resolution = 0.01
    x_mesh, y_mesh = np.meshgrid(np.arange(Xmin, Xmax, resolution),
                             np.arange(Ymin, Ymax, resolution))
    MeshDat = np.array([x_mesh.ravel(), y_mesh.ravel()]).T
    z = DeepLearningSVM.predict(MeshDat)
    z = np.reshape(z, (len(x_mesh), len(y_mesh)))
    im1 = plt.scatter(x_mesh[z=='N'], y_mesh[z=='N'], s=5, alpha=0.3, c='blue')
    im2 = plt.scatter(x_mesh[z=='A1'], y_mesh[z=='A1'], s=5, alpha=0.3, c='orange')
    im3 = plt.scatter(x_mesh[z=='A2'], y_mesh[z=='A2'], s=5, alpha=0.3, c='red')

    # 深層学習による特徴量空間プロットの可視化
    im4 = plt.scatter(DLF[0:30, 0], DLF[0:30, 1], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
    im5 = plt.scatter(DLF[30:60, 0], DLF[30:60, 1], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
    im6 = plt.scatter(DLF[60:90, 0], DLF[60:90, 1], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
    
    plt.title("Deep Learning + SVM")
    plt.xlabel("Deep Feature 1")
    plt.ylabel("Deep Feature 2")
    plt.legend(["Class N", "Class A1", "Class A2"], loc="upper right")
    plt.grid(True)
    
    # フレームを積む
    DL_Frame.append([im1, im2, im3, im4, im5, im6]) 

plt.show()
#aniDLSVM = animation.ArtistAnimation(fig, DL_Frame, interval=500)
#aniDLSVM.save('DLSVM_Movie.mp4', writer="ffmpeg")

Programming 23. 深層学習と決定木の連携（動画）

In [26]:
import matplotlib.animation as animation
DL_Frame = []
DLM = []
fig = plt.figure()

# 学習回数を25回ずつ増やしながら、深層学習構築、特徴抽出を繰り返す
for i in range(1, int(700/25)):
    
    plt.clf() # ★注意: 動画にしたい場合には、コメントアウトすること！！
    
    # 深層学習モデルの構築
    DLM.append(
        MLPClassifier(hidden_layer_sizes=(5, 2, ), solver="sgd", activation="logistic", learning_rate_init=0.05, 
                      random_state=10, max_iter=(i-1)*25+1, alpha=0.01, tol = pow(10, -10), verbose=False)
    )
    DLM[i-1].fit(X_DL, Y_train)
    
    # 学習済みモデルの特徴抽出
    DLF=[]
    for j in range(0, 90):
        temp = 1/(1+np.exp(-1 * np.dot(DLM[i-1].coefs_[0].T, X_DL[j]) + DLM[i-1].intercepts_[0]))
        DLF.append(1/(1+np.exp(-1 * np.dot(DLM[i-1].coefs_[1].T, temp) + DLM[i-1].intercepts_[1])))        
    DLF = np.array(DLF)
    
    # 深層学習による特徴量空間上に決定木で識別境界を生成
    DeepLearningTr = tree.DecisionTreeClassifier(min_samples_leaf=1, random_state=15)
    DeepLearningTr = DeepLearningTr.fit(DLF, Y_train)
    
    # 決定木による識別境界可視化
    Xmin, Ymin, Xmax, Ymax = 0, 0, 1, 1 # 空間の最小最大値
    resolution = 0.01
    x_mesh, y_mesh = np.meshgrid(np.arange(Xmin, Xmax, resolution),
                             np.arange(Ymin, Ymax, resolution))
    MeshDat = np.array([x_mesh.ravel(), y_mesh.ravel()]).T
    z = DeepLearningTr.predict(MeshDat)
    z = np.reshape(z, (len(x_mesh), len(y_mesh)))
    im1 = plt.scatter(x_mesh[z=='N'], y_mesh[z=='N'], s=5, alpha=0.3, c='blue')
    im2 = plt.scatter(x_mesh[z=='A1'], y_mesh[z=='A1'], s=5, alpha=0.3, c='orange')
    im3 = plt.scatter(x_mesh[z=='A2'], y_mesh[z=='A2'], s=5, alpha=0.3, c='red')

    # 深層学習による特徴量空間プロットの可視化
    im4 = plt.scatter(DLF[0:30, 0], DLF[0:30, 1], s=30, c='blue', marker='s', linewidths=0.5, edgecolors='black')
    im5 = plt.scatter(DLF[30:60, 0], DLF[30:60, 1], s=30, c='orange', marker='s', linewidths=0.5, edgecolors='black')
    im6 = plt.scatter(DLF[60:90, 0], DLF[60:90, 1], s=30, c='red', marker='s', linewidths=0.5, edgecolors='black')
    
    plt.title("Deep Learning + Tree")
    plt.xlabel("Deep Feature 1")
    plt.ylabel("Deep Feature 2")
    plt.legend(["Class N", "Class A1", "Class A2"], loc="upper right")
    plt.grid(True)
    
    # フレームを積む
    DL_Frame.append([im1, im2, im3, im4, im5, im6]) 

plt.show()
#aniDLTr = animation.ArtistAnimation(fig, DL_Frame, interval=500)
#aniDLTr.save('DLTr_Movie.mp4', writer="ffmpeg")

4. おわりに

本日はお疲れ様でした。本日は、加速度・角速度という時系列信号から、異常状態を検出する様々な手法について説明しました。また、わかりやすさを優先にするため、単純な問題、単純なモデルを主体として解説しました。

  • 人工知能や深層学習は魔法のような方法だと思っていた方:
    • そんなにすごいものでもないことがわかったでしょうか。
  • 人工知能をブラックボックスで使ってい方:
    • 多次元だと訳がわからないままモデルが出来上がっていたと思います。
    • 2次元で考えると、面に対し、ただ単に、線を引く方法の集まりでしかないので、結構単純です。
  • 深層学習とそれ以外の人工知能は対立する手法だと思っていた方：
    • 基本的にはお互い仲良しです。
  • 深層学習は古典的な手法よりも有効だと思っていた方：
    • 条件をかなり揃えないと、そうなりません。。。

簡単なモデルを作る分には、今回の内容で十分です。 ぜひ、いろいろな問題解決に応用していただければと思います。

Appendix. 環境構築など

  • 環境構築
    • 以下に、Linux系OSでの環境構築方法を公開しています（手軽にやりたい場合はMacbookなどでokです）。
    • 自由自在にAIを構築したい場合、Windows系OSだと結構苦しいのであまりオススメしません。
    • http://int-info.com/index.php/env/