Pythonの制御構文周りの備忘録です。
随時追加・修正を行っていきます。

比較演算子

演算子	説明
A == B	AとBが等しい
A != B	AとBが等しくない
A < B	AがBよりも小さい
A > B	AがBよりも大きい
A <= B	AがB以下
A >= B	AがB以上
A in [LIST]	[LIST]の中にAがある
A not in [LIST]	[LIST]の中にAがない
A is None	AがNoneである
A is not None	AがNoneでない
条件A and 条件B	条件Aと条件Bのどちらも満たす
条件A or 条件B	条件Aまたは条件Bのどちらかを満たす

if

num = 5
if num == 0:
    print('数値は0')
elif num < 0:
    print('数値は0より小さい')
else:
    print('数値は0より大きい')

while

limit = input('Enter:') # 入力を受け付ける
count = 0
while True:
    if count >= limit:
        break # countが10以上ならループを抜ける

    if count == 5:
        count += 1
        continue # countが5なら次のループへ

    print(count)
    count += 1
else: # breakされずにループが終了したら実行する
    print('Done')

for

for i in [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]:
    if i == 5: # iが5なら次のループへ
        continue:

    if i == 8: # iが8ならループを抜ける
        break:

    print(i)
# breakされずにループが終了したら実行する
else:
    print('Done')

# 指定回数処理したいが値を取り出す必要が無い時に_を利用する
for _ in range(10):
    print('hello')

# 2から3つ飛ばしで10を超えるまで処理を行う
for i in range(2, 10, 3):
    print('hello')

# インデックスも取得したい場合
for i, animal in enumerate(['dog', 'cat', 'bird']):
    print(i, animal)

# 複数のリストを同時に展開して取得したい場合
animals = ['dog', 'cat', 'bird']
foods = ['meat', 'fish', 'bean']
for animal, food in zip(animals, foods):
    print(animal, food)

# 辞書のループ処理
data = {'x': 10, 'y': 20}
for k, v in d.items():
    print(k, ':', v)

Pythonの制御構文周りの備忘録です。
随時追加・修正を行っていきます。

演算子

演算子	説明
A == B	AとBが等しい
A != B	AとBが等しくない
A < B	AがBよりも小さい
A > B	AがBよりも大きい
A <= B	AがB以下
A >= B	AがB以上
A in [LIST]	[LIST]の中にAがある
A not in [LIST]	[LIST]の中にAがない
A is None	AがNoneである
A is not None	AがNoneである
条件A and 条件B	条件Aと条件Bのどちらも満たす
条件A or 条件B	条件Aまたは条件Bのどちらかを満たす

if

num = 5
if num == 0:
    print('数値は0')
elif num < 0:
    print('数値は0より小さい')
else:
    print('数値は0より大きい')

while

limit = input('Enter:') # 入力を受け付ける
count = 0
while True:
    if count >= limit:
        break # countが10以上ならループを抜ける

    if count == 5:
        count += 1
        continue # countが5なら次のループへ

    print(count)
    count += 1
else: # breakされずにループが終了したら実行する
    print('Done')

for

for i in [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]:
    if i == 5: # iが5なら次のループへ
        continue:

    if i == 8: # iが8ならループを抜ける
        break:

    print(i)
# breakされずにループが終了したら実行する
else:
    print('Done')

# 指定回数処理したいが値を取り出す必要が無い時に_を利用する
for _ in range(10):
    print('hello')

# 2から3つ飛ばしで10を超えるまで処理を行う
for i in range(2, 10, 3):
    print('hello')

# インデックスも取得したい場合
for i, animal in enumerate(['dog', 'cat', 'bird']):
    print(i, animal)

# 複数のリストを同時に展開して取得したい場合
animals = ['dog', 'cat', 'bird']
foods = ['meat', 'fish', 'bean']
for animal, food in zip(animals, foods):
    print(animal, food)

# 辞書のループ処理
data = {'x': 10, 'y': 20}
for k, v in d.items():
    print(k, ':', v)

概要

Python2のサポート終了からそろそろ1年経とうとしていますが、Python3に切り替えたことによってdict.valuesの挙動が変わりましたので、改めて紹介したいと思います。
というのも、書き方によってはPython2時代と同じ挙動するので、そもそも以下で説明する問題に気付いていない可能性もあります。

ざっくり公式ページの内容を要約しますと、
Python2で辞書型の変数をvalues()、items()、keys()メソッドで取得すると結果はリスト型で取得できていたが、Python3では、辞書ビューオブジェクトが返されるようになりました。
なので、listフィルタを追加することでPython2と同じ結果が得られるようになります。

さっそくコード

Python3ではブロックスタイルと通常のスタイルで結果が異なるので、大丈夫だと思っていても、書き方(通常のスタイル)を変えただけで意図しない結果を招く可能性があります。
ちなみにansible_python_interpreterで挙動が変わるわけではなく、ansibleそのものがPython2.xかPython3.xによって挙動が変わります。

Python3実行環境

# 設定なし
> ansible-config view
ERROR! Invalid or no config file was supplied

# Python3版 ansible
> ansible --version

ansible 2.9.14
  config file = None
  configured module search path = ['/root/.ansible/plugins/modules', '/usr/share/ansible/plugins/modules']
  ansible python module location = /usr/local/lib/python3.9/site-packages/ansible
  executable location = /usr/local/bin/ansible
  python version = 3.9.0 (default, Nov 18 2020, 13:28:38) [GCC 8.3.0]

Python2実行環境

# 設定なし
> ansible-config view
ERROR! Invalid or no config file was supplied

# Python2版 ansible
> ansible --version

ansible 2.9.14
  config file = None
  configured module search path = [u'/root/.ansible/plugins/modules', u'/usr/share/ansible/plugins/modules']
  ansible python module location = /usr/local/lib/python2.7/site-packages/ansible
  executable location = /usr/local/bin/ansible
  python version = 2.7.18 (default, Apr 20 2020, 19:27:10) [GCC 8.3.0]

dict_values.yml

---
- name: Dict values 
  hosts: localhost
  connection: local
  gather_facts: False
  vars:
    mydict: {"a": {"aa": 1 }, "b": {"bb": 2 }, "c": {"cc": 3 }}
  tasks:
    - name: block style debug
      debug:
        msg: >-
          {%- for d in mydict.values() -%}
            {{ d }}
          {%- endfor -%}

    - name: normal style debug
      debug:
        msg: "{{ item }}"
      with_items: "{{ mydict.values() }}"

ブロックスタイルと通常の結果が全く同じ

実行結果(python2)

TASK [debug] *************************************************************************
ok: [localhost] => {
    "msg": "{u'aa': 1}{u'cc': 3}{u'bb': 2}"
}

TASK [debug] *************************************************************************
ok: [localhost] => (item={u'aa': 1}) => {
    "msg": {
        "aa": 1
    }
}
ok: [localhost] => (item={u'cc': 3}) => {
    "msg": {
        "cc": 3
    }
}
ok: [localhost] => (item={u'bb': 2}) => {
    "msg": {
        "bb": 2
    }
}

通常の結果がdict_valuesで出力されています。
なので、この結果をループしたりすると動作が変わるので、要注意です。

実行結果(python3)

TASK [debug] *************************************************************************
ok: [localhost] => {
    "msg": "{'aa': 1}{'bb': 2}{'cc': 3}"
}

TASK [debug] *************************************************************************
ok: [localhost] => (item=dict_values([{'aa': 1}, {'bb': 2}, {'cc': 3}])) => {
    "msg": "dict_values([{'aa': 1}, {'bb': 2}, {'cc': 3}])"
}

対策

listフィルタを追加することで解決します。
全体のコードの整合性を保つためにブロックスタイルにもlistフィルタを追加するのも良いかと思います。（出力結果は変わりません）

dict_values.yml

---
- name: Dict values 
  hosts: localhost
  connection: local
  gather_facts: False
  vars:
    mydict: {"a": {"aa": 1 }, "b": {"bb": 2 }, "c": {"cc": 3 }}
  tasks:
    - name: block style debug
      debug:
        msg: >-
          {%- for d in mydict.values() | list -%}
            {{ d }}
          {%- endfor -%}

    - name: normal style debug
      debug:
        msg: "{{ item }}"
      with_items: "{{ mydict.values() | list }}"

print文

print文

%s # 文字列
%r・・・
%d # 整数
%f # 固定小数点表記
%1.5f # 固定小数点表記(,5のところは小数以下桁)
%e # 指数表記

print("数値1=%f, 数値2=%.3f" % (1/3, 1/3))
数値1=0.333333, 数値2=0.333

引数"end"に値を入れることで、末尾に何をつけるかを設定

print("a", end=",")
print("b", end=",")
print("c")

対話モード時に最後に表示した式を格納している変数

>>> tax = 12.5 / 100
>>> price= 100.5
>>> price * tax
12.5625
>>> price + _
113.0625
>>> round(_,2)
113.06

モジュールの検索パス

1.ビルトインモジュールの中
2.sys.path変数で得られるディレクトリのリスト
3.入力スクリプトのあるディレクトリ
4.PYTHONPATH
5.インストールごとのデフォルト

コーディングスタイル

Pythonのコーディングスタイルを PEP8 という

コーディングスタイル

・インデントはスペース４つ
・タブ禁止
・４スペースは狭いインデントと広いインデントのちょうどよい妥協点。
・コメント行の独立。
・docstring
・エンコーディングは、UTF-8(default)・ASCII
・演算子の周囲やカンマの後ろにスペースを入れる
・カッコのすぐ内側にはスペースを入れない

シーケンスオブジェクトの比較

・どちらかのシーケンスがなくなるまで比較
・同じ要素は比較しない
・違う要素の比較結果

クラス変数とインスタンス変数

・クラス変数とは・・・すべてのインスタンスが共有するメモリ
・インスタンス変数とは・・・各インスタンスに固有のメモリ 

以下は間違いの例・・・クラス変数にmutable(変更可)を使ってはいけない

クラス変数とインスタンス変数

class Sample:
  c_list = []・・・クラス変数の使い方を間違えた例
  def add_c_list(self,data):
    self.c_list.append(data)

print("出力結果:", end=" ")
sample1 = Sample()
sample1.add_c_list("データ1")

sample2 = Sample()
sample2.add_c_list("データ2")

for item_data in sample1.c_list:
  print(item_data, end=" ")
=============================
出力結果: データ1 データ2

エスケープシーケンスの無効化

len関数

エスケープシーケンスは１文字でカウント
改行の\nは１文字としてカウント

range関数

range関数

>>>print(range(5))
range(0,5)

・range関数は反復可能体（イテラブル）
・range関数はオブジェクト

ビルトイン関数

ビルトイン関数

>>>dir(モジュール名) #dirがビルトイン関数。モジュールが定義している名前をすべて表示。

1章 食欲をそそってみようか

・データ型の一般性が高く、問題領域は、Awk・Perlより広く、他言語と同等以上

2章　インタープリタの使い方

インタープリタの起動

>>> python -c コマンド
>>> python -m モジュール名
>>> 

3章　気楽な入門編

演算子

・べき乗演算子は、他の演算子より優先順位が高いため例外的に右から左へと評価される
・演算対象の型が混同していた場合(int,float)、整数は浮動小数点に変換される

文字列の特徴

>>> word[10000] # 大きすぎるインデックスを指定するとエラー
IndexError  Traceback (most recent call last)
<ipython-input-4-47f442646512> in <module>
----> 1 Zen[50]

IndexError: string index out of range

>>> word[10000:20000] # 範囲外のスライシングはうまく処理
''

複数行にまたがる改行文字を取り消す

print("""\
Usage:thingy[options]
    -h                Display this usage message
    -H hostname       Hostname to connect to
""")

対話モードの電卓。アンダースコア

最後に表示した式を変数「_」（アンダースコア）に代入してある。
>>> tax = 12.5/100
>>> price = 100.50
>>> price * tax
12.5625
>>> price + _
113.0625
>>> round(_,2)
113.006

対話モード

列挙された文字リテラルは自動的に連結される

>>> 'Py' 'thon'
'Python'

4章　制御構造ツール

リストのアンパック

>>> list(range(3,6)) # 一般的な個別の引数を使ったコール
[3,4,5]

>>> args = [3,6] # ここからは特殊な方法
>>> list(range(*arg))　# *arg を使うとアンパックされる。ただの　3,6　になる。
[3,4,5] # range(3,6) と同じ

>>> 

ドキュメンテーション(docstring)

・１行目：簡潔な要約。大文字始まる、ピリオド終わり。
・２行目：空白
・３行目：

関数アノテーション

def 関数名(arg1: 'arg1の説明', arg2: 'arg2の説明', , ,)->'戻り値の説明':
    処理

docstringと関数アノテーションの例

def my_func(n: 'この値から足し始める', m:'この値まで足す') -> 'nからmの合計値':
    """ nからmまでの合計を返す関数 """
    ret = 0
    for i in range(n, m+1):
        ret += i 
    return ret

5章　データ構造

リストをキューとして使う

>>> from collections import deque
>>> queue = deque(["A","B","C"]) # キューを作った
>>> queue.append("D") # Dを追加
>>> queue.popleft() # 最初を取出
>>> queue.pop() # 最後を取出
>>> queue.pop(idx) # idx番目を取出

ディクショナリ

・キー：変更不能型
・値：変更可能型
・キーkeyの存在を確認、取得（検索）: in演算子
・値valueの存在を確認、取得（検索）: in演算子, values()
・キーkeyと値valueの組み合わせの存在を確認: in演算子, items()

多次元リストのソート

演算子

比較演算子 < <= == != is  is not  in  not in

条件についての補足

・比較演算子 in および not in ・・・シーケンスの値の有無
・演算子 is および is not ・・・オブジェクトの比較
・ブール演算子 and および or・・・短絡演算子
・比較の組み合わせ（複数条件） if x<y and x>z

演算子の優先順位

数値演算子　＞　比較演算子

シーケンスの比較、その他の型の比較

・2つが基本的に同じシーケンスで、片方の長さが短い時は、この短い方が小となる。
・文字列の辞書的順序には、個々の文字のUnicodeコードポイント番号で比較

6章　モジュールとパッケージ

dir関数

>>> import sys,fibo
>>> dir(fibo)

・モジュールが定義している名前を確認するのに使う。

モジュールの検索は、○○モジュールの場合、○○.pyを以下の順番で検索する

モジュールの検索パス

１．ビルトインモジュール内で探す
２．sys.path変数で得られるディレクトリのリストを使ってspam.pyを検索
　　２－１．入力スクリプトのあるディレクトリ
　　２－２．PYTHONPATH
　　２－３．インストールごとのデフォルト

module

・モジュールは、ファイルである
・モジュールのファイルは、「.py」

package

・パッケージはフォルダーである

>>> from パッケージ名 import モジュール名 #これを使うとモジュールを参照する際に短くてフルネームにする必要がない
>>> import パッケージ名.モジュール名 #これを使うとモジュールを参照する際に長くてフルネームにする必要がある。[ドット区切りモジュール名]と呼ぶ。

package実装例

>>> import sound.effects.echo
>>> sound.effects.echo.echofilter()#サブモジュールのロード。参照はフルネーム。長くてダメ！！

>>> from sound.effects import echo
>>> echo.echofilter()#参照が短くて良い！！

コンパイル済みPythonファイル

・Pythonは .py 以外に、Pythonコードをコンパイルした .pyc というファイルも実行できる
・インタプリタなので１行ごとにバイナリファイルに変換
・コンパイラはまとめてバイナリファイルに変換

8章 エラーと例外

概要

・”エラー”は、「構文エラー」と「例外」に大別される
・"構文エラー"は「パース上のエラー」「構文解釈エラー」と呼ばれる
・”例外”は「文や式が正しくても、実行すると起こるエラー」である

例外

ZeroDivisionError
NameError
TypeError
KeyboardInterrupt # キーボード割込例外 [Ctrl]+[c]

9章　クラス

10章　標準ライブラリめぐり（モジュール）

モジュール名を参照短縮

>>> from パッケージ名.サブモジュール名 import モジュール名 # これを使うとモジュールを参照する際に短くてフルネームにする必要がない
>>> import パッケージ名.モジュール名 #これを使うとモジュールを参照する際に長くてフルネームにする必要がある。[ドット区切りモジュール名]と呼ぶ。

module

import os # オペレーティングシステムとやり取りする関数
import glob # ファイルをワイルドカード検索
import sys # コマンドライン引数を処理する
import re # 正規表現
import math # 浮動小数点数学
from struct import * # バイナリー
import random # ランダム
import collections # リスト
import logging # ログ取り

コマンドライン引数

>>> import sys
>>> print(sys.argv)

randomモジュール・・・無作為抽出のツール

>>> import random
>>> random.choice(['apple','banana','lemon']) # choice はリスト内から選択
'apple'
>>> random.sample(range(100),10) # sample は、第１引数から第２引数個を重複なく抽出
>>> random.random() # ランダムな浮動小数点

>>> radom.randrange(6) # range(6)からランダムに選んだ整数

11章　標準ライブラリめぐり

ログ取り

・ログ出力なのか、ユーザ(プログラム実行者)に伝えたい情報としての出力なのかを切り分けられる
・ログの種類をErrorやDebugのように、レベル分けが可能
・フォーマットを指定すれば、簡単に統一化された出力が可能

ログの優先順位（左から優先順位が低い）

低い＜－　　　　　　　　　　　　－＞高い
DEBUG、INFO、WARNING、ERROR、CRITICAL

12章　仮想環境とパッケージ

pipによるパッケージ管理

>>> pip install パッケージ名 # 最新バージョンのパッケージをインストール 
>>> pip install パッケージ名 ==2.6.0 # 特定バージョンのパッケージをインストール
>>> pip install --upgrade パッケージ名 # 最新バージョンにアップグレード
>>> pip uninstall パッケージ名 # パッケージをアンインストール
>>> pip list # インストール済み確認
>>> pip freeze # インストール済み確認(出力形式がpip install)
>>> pip show パッケージ名 # パッケージの詳細表示　version,author,summar,説明hp

仮想環境

>>> deactivate # 仮想環境終了

14章　対話環境での入力行編集とヒストリ置換

ソースコード・エンコーディング

・Pythonソースコードのエンコーディング：UTF-8
・敢えてエンコードを変える時は以下になる
#-*- coding:エンコーディング名 -*-

インタープリタの終了

ctrl+d
>>>exit()
>>>quit()

タブ補完とヒストリ編集

対話型インタープリタ

bpython
IPython

アクティベート状態から抜ける

対話モードのプライマリプロンプト／セカンダリプロンプト

>>> 
...

1. メソッド（関数）を宣言する

Python3

def sample_func_1(text : str) -> str:
    add_text = "追加した文字列"
    return (text + add_text)


def sample_func_2(text : str) -> str:
    add_text2 = "追加した文字列"
    return (text + add_text2)

2. 辞書型（dict）オブジェクトのValueに、関数を括弧を付けずに記述する

Python3

func_dict = {
    'label_1' : sample_func_1,
    'label_2' : sample_func_2
    }

3. 辞書型オブジェクトKeyラベル名で、関数を呼び出して実行する

辞書（dict）のvalue値として、格納されているのは、関数（メソッド）オブジェクトそのものです。
そのため、関数（メソッド）を実行するには、辞書（func_dict）のvalue値に、引数を受け取る丸括弧（"( )"）を付ける必要があります。

丸括弧（"( )"）を付けないで、辞書（dict）のvalue値として格納されている関数（メソッド）「それ自体」を呼び出した場合の結果は、以下になります。

Python3

func_dict['label_1']
<function sample_func_1 at 0x100bafb80>

上記のコードで辞書から呼出した関数（メソッド）を実行するには、引数を受け取る丸括弧（"( )"）をつける必要があります。

Python3

func_dict['label_1']()

なお、funct_dictが引数を1つ以上、受け取るメソッドである場合は、引数を与えないで関数呼出しをすると、エラーが出ます。

定義したメソッド（関数）は、str型のオブジェクトを引数を受け取る関数でした。
引数に適当な文字列を渡して、func_dict内に格納されているメソッド（関数）を実行してみます。

Python3

func_dict['label_1']("vxw")
# 実行結果
'vxw追加した文字列'

（ 参考 ）

Python3

func_dict
# 実行結果
{'label_1': <function sample_func_1 at 0x100bafaf0>, 'label_2': <function get_text_from_ppt at 0x100bafa60>}

func_dict['label_1']
<function sample_func_1 at 0x100bafb80>

type(func_dict['label_1'])
# 実行結果
<class 'function'>

type(func_dict['label_1']())
# 実行結果
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: sample_func_1() missing 1 required positional argument: 'text'

type(func_dict['label_1']("giho"))
# 実行結果
<class 'str'>

func_dict['label_1']("giho")
# 実行結果
'giho追加した文字列'

returned_value = func_dict['label_1']("giho")
returned_value
# 実行結果
'giho追加した文字列'

（ 参考にしたWebページ ）

1. @satoichiさんのQiita記事 「【Python】辞書を利用して関数をまとめて実行する」
2. 「[Python] **を用いたキーワード引数の辞書化」

（ 使いみち ）

今回、紹介した方法を用いることで、辞書型オブジェクトに定義されたKeyの値に応じて、そのKeyに対応する関数（メソッド）を呼び出すことができます。

データを処理する流れのなかで、条件に応じて、呼び出すべき適切な関数（メソッド）が切り替わるロジックが必要な場面で、役に立ちます。

基本的な例として、読み込んだファイルの拡張子の内容に応じて、そのファイルに所要の処理を行う関数（メソッド）を動的に選び出して呼出して実行する場面、などが考えられます。

PDFファイルやWordファイル、Excelファイルから、テキスト文字列のみを切り抜いてデータ取得するコードは、以下のWebページで紹介されています。

1. deecode blog 「[Python] Word/Excel/PowerPoint/PDFからテキスト抽出するライブラリ・サンプルコード」
2. @butadaさんのQiita記事 「pdf/docxファイルからのテキストマイニング」

Python3

def get_text_from_txt(filename : str) -> str:
    filenameとして受け取ったパスに格納されているテキストファイルをreadして、テキストを*str*オブジェクトとして取得して、返り値として出力する処理

    return text

def get_text_from_pdf(filename : str) -> str:
    filenameとして受け取ったパスに格納されているPDFファイルをreadして、テキスト文字列を抜き出して、ひとつの*str*オブジェクトに連結して、返り値として出力する処理
    return text

def get_text_from_word(filename : str) -> str:
    filenameとして受け取ったパスに格納されているWordファイルをreadして、テキスト文字列を抜き出して、ひとつの*str*オブジェクトに連結して、返り値として出力する処理
    return text

def get_text_from_excel(filename : str) -> str:
    filenameとして受け取ったパスに格納されているExcelファイルをreadして、テキスト文字列を抜き出して、ひとつの*str*オブジェクトに連結して、返り値として出力する処理
    return text

func_dict = {
    'txt' : get_text_from_txt,
    'pdf' : get_text_from_pdf,
    'word' : get_text_from_word,
    'excel' : get_text_from_excel
         }

（ 備考 ）

辞書のValueに関数（メソッド）を格納するときに、関数に括弧を付けると、以下のエラーが発生する。
関数の実行を行おうとしていると、解釈されてしまうため。

PZython3

func_dict = {
         'label_1' : sample_func_1(),
         'label_2' : sample_func_2()
         }
# 実行結果
Traceback (most recent call last):
  File "<stdin>", line 2, in <module>
TypeError: sample_func_1() missing 1 required positional argument: 'text'

辞書を定義する上のコードのなかで、sample_func_1()の括弧のなかに特定の引数を与えると、辞書は定義できる（エラーは起きない）。

しかし、辞書のValueに格納された関数は、辞書を定義する際に受け取った引数しか、受け取ることのできない関数になってしまう。

（ 辞書の定義 ）

Python3

func_dict = {
    'label_1' : sample_func_1("abc"),
    'label_2' : sample_func_2("def")
}

（ 辞書内の関数の呼出し実行 ）

func_dict辞書内で、Key"label_1"に対応するValueとして定義されている、関数sample_func_1("abc")が、引数"abc"を与えられた状態で、実行された。

Python3

func_dict['label_1']
# 実行結果
'abc追加した文字列'

func_dict['label_1']（を評価すると）は、すでに返り値として、str型のインスタンス'abc追加した文字列'を返している。そのため、func_dict['label_1']にさらに括弧を付けると、TypeError: 'str' object is not callableエラーが起きる。

Python3

func_dict['label_1']()
# 実行結果
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'str' object is not callable

Python3

func_dict['label_1']("abc")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'str' object is not callable

Python3

unc_dict['label_1']('abc')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'str' object is not callable

---

（ 備考２ ）

辞書のValにテキスト文字列を格納して、evalとexecを用いて、関数をよびだせる。

しかし、キーボードのストローク数が増えてしまい、エレガントではない。

Python3

func_dict = {
    'label_1' : "sample_func_1({arg})",
    'label_2' : "sample_func_2({arg})"
    }

Python3

func_dict
# 実行結果
{'label_1': 'sample_func_1({arg})', 'label_2': 'sample_func_2({arg})'}

func_dict['label_1']
# 実行結果
'sample_func_1({arg})'

expression = func_dict['label_1'].format(arg="'123'")
expression
# 実行結果
"sample_func_1('123')"


eval(expression)
# 実行結果
'123追加した文字列'

result = eval(expression)
result
# 実行結果
'123追加した文字列'

add_text = "result2="
text = add_text + expression
text
# 実行結果
"result2=sample_func_1('123')"

exec(text)
result2
# 実行結果
'123追加した文字列'

Pythonで利用可能な四則演算一覧

数値型で使用可能な四則演算の一覧まとめです。数値型と言って急に出てきました。もう急に難しい感じです。要は「各データは型に分類されるよ、データ取り扱いの際は注意してね」ってことで納得しました。詳しくは下記を参照してください。最初の数行で挫折しそうですが目を通しました。脳みそをすべっていく感じがします。
データ型とは

ためしに

足し算のことです。これを記載する前にtest.pyにひたすら足し算の式を記入しターミナルで実行してもなにも表示されなくて途方にくれましたが、きちんとprintでくくらないとターミナルに表示されないんですね。もしくはターミナルに「python」て入力してインタラクティブモードで行う必要があります。

一覧

出力しないのでprintは記載しません。ターミナルに出力した場合はprintで囲ってください。普通の四則演算はできます。

test.py

x + y #xとyの和
x - y #xとyの差
x * y #xとyの積
x / y #xとyの商
x % y #xをyで割った余り
x ** y #xのy乗

型の異なる演算

これらは何も数値型だけのものではなく、文字列型でも適用可能です。ためしに"apple"と"banana"を足してみるとターミナルにはapplebananaと結果が表示されます。
しかし数値型と文字列型はそのままでは演算できません。
試しに文字列型appleと数値型2をそのまま足し算するとエラーがでました。
数値型の2を文字列型の"2"とするとapple2とエラーなく計算されました。pythonのような動的型付け言語の場合はそこまで意識することはありませんが、意識して考えていくとその他の言語を勉強する際にも役に立つかもしれません。
Javaの本とかを読んでいると「そんなのみればわかるだろうよ！！」という部分にも型を指定しなければならないので大変そうでした。

はじめに

ターミナルでpythonを実行するとモジュールのインポートでエラーが出てしまいました。今回はその備忘録です。

環境

mac：Big Sur
エディタ：VScode
Python：3.8

トラブル内容

VScodeで書いたpythonファイルをターミナルで実行するとそんなモジュールねぇよって怒られる。
でもインポートしてるのnumpyだし、ないわけないじゃん！！てことでひたすらググったら
どうやらターミナルから使用するpythonとVScodeで設定したpythonが異なるとインポートエラーが起こる模様

解決方法

ターミナルで

pip show numpy

と入力してLocationのところを見る。自分の場合は

anaconda3/lib/python3.8/site-packages

となっていてanaconda3の中のpythonを使用していた。なのにVScodeで/usr/bin/python3を指定
していたのでエラーが出てしまった。わかれば簡単なことだったなぁ。

後日談

このトラブルを解決した数日後にまたしても同じエラーが起きた。今回は使用するpythonも一致しているのにどしてと思ってまたひたすらググる。結果としては、モジュールのバージョンが古かったらしい。
なので

pip install -U numpy

とすることで解決。numpyのところをエラーの出ているモジュール名に変えれば汎用できます。
ちなみに

pip list -o

とするとバージョンの古いやつが一覧で表示されます。

ポケモンの環境とは

ポケモン対戦で使われる環境という言葉は，ポケットモンスター ソード・シールドのランクバトルにおいてどんなポケモン・戦略・構築が使われているかを指す言葉です．環境を知ることで，それに対するメタを考察することで現環境で勝ちやすい構築を作っていくのがポケモン対戦の醍醐味と言えます．

環境の中心

ガチ勢向けポケモン対戦情報共有サイトのポケモンソルジャーさんが運営するキャラランクメーカーではS ~ Eランクまでの評価を個人的に設定することができ，有名な実況者や解説動画投稿者が主観的に評価を下したキャラランクが参考にされています．

ポケモンソルジャーさん自身も定期的にキャラランクを公開しています．
https://youtu.be/TMXQ497foRg

このキャラランクのSランクについて，ポケモンソルジャーさんでは「環境の中心」と動画内で定義しています．
一方でその環境の中心の評価については主観的であることが問題であると考え，定量的な評価手法が必要であると考えました．

目標

目標はずばり，「環境の中心を定量的に評価する」ことです．この記事では，可能な限り環境の中心をデータに基づいて考察したいと思います．

ポケモンホーム

ポケモンホームはスマホとSwitchで利用可能なポケモンの管理アプリです．ソフト内に保管しておけるポケモンの数に限りがあるため，こちらに不要なポケモンを一時的に保管するなどして使うことができます．

このポケモンホームにはそれ以外には交換機能やランクバトルでのポケモンの環境情報・プレイヤーのランキング機能なども備わっていて，ガチ勢の間ではこまめに環境調査に利用されています．

API

コチラの記事で紹介されている通り，ポケモンホームはUser-Agentを偽装することでアプリ内で利用されているAPIを直接叩いてjsonデータを取得することができます．今回使用するのは以下のAPIです．詳細についてはリンク先からご確認ください．

# シーズン・タームの情報の取得
https://api.battle.pokemon-home.com/cbd/competition/rankmatch/list
# 該当シーズン・タームの対戦データの取得
https://resource.pokemon-home.com/battledata/ranking/$id/$rst/$ts2/pokemon
# 該当シーズン・タームにおけるポケモンのデータ
https://resource.pokemon-home.com/battledata/ranking/$id/$rst/$ts2/pdetail-$j

評価手法

環境を視覚的にもわかりやすく，かつ分析がしやすい形に加工するために，今回はグラフを用いた評価を検討します．
グラフについては過去に書いた記事があるので，コチラを参照してください．

この記事では，使用率ランキング上位のポケモンとの同時採用率の高さを評価することにしました．
環境の中心は，単体の性能ではなく他のポケモンに与える影響やパーティへの組み込みやすさを重視すべきだと考えたためです．
ポケモンホームの使用率ランキング上位30匹までの同時採用率が高いポケモンTOP10までを考慮し，その関係性を無向グラフで可視化することにします．

ノードとエッジ

無向グラフではノードとエッジを定義する必要があります．

今回はノードをポケモン,エッジを同時採用率上位の有無としてグラフを作成します．

例えば，エースバーンと同時採用率が高い10匹を考えるとすると，エースバーンのノードから10本のエッジをそれぞれのポケモンのノードに伸ばすことになります．

これを上位30匹のポケモンに対して繰り返し，次数が多いノードに君臨するポケモンが環境上位のポケモンと相互に関係が強く影響を与えるポケモンであると仮定しました．

データ収集

必要なデータは以下の通りです．

• ポケモン図鑑（図鑑番号と名前の対応づけ）
• シーズン・タームのID情報（rst・ts1・ts2）
• 当該シーズン・タームの対戦環境情報（使用率上位のポケモンリスト）
• 当該シーズン・タームのポケモンに関する情報（同時採用率上位のポケモンリスト）

ポケモン図鑑

8世代に対応したポケモン図鑑を探したところ，jsonファイルで配布してくださっている方がいたのでコチラをお借りしました．まだ鎧の孤島以降に追加されたウーラオス・ブリザポス等のポケモンの情報が更新されていないため，最新の環境で試すと図鑑にないポケモンのせいで少し結果にノイズが入ります．

> git clone https://github.com/dayu282/pokemon-data.json

import json

pokedex = {}
for i in range(8):
  # すべての世代のポケモンを追加
  j = i + 1
  with open("pokemon-data.json/en/gen" + str(j) + "-en.json", "r", encoding = "utf-8") as f:
    pokemons = json.loads(f.read())
    for pokemon in pokemons:
      pokedex[pokemon["no"]] = pokemon["name"]

シーズン・タームのID情報

# シーズンの情報を丸ごと全取得してjson形式で保存
> curl https://api.battle.pokemon-home.com/cbd/competition/rankmatch/list \
  -H 'accept: application/json, text/javascript, */*; q=0.01' \
  -H 'countrycode: 304' \
  -H 'authorization: Bearer' \
  -H 'langcode: 1' \
  -H 'user-agent: Mozilla/5.0 (Linux; Android 8.0; Pixel 2 Build/OPD3.170816.012) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/83.0.4103.61 Mobile Safari/537.36' \
  -H 'content-type: application/json' \
  -d '{"soft":"Sw"}' \
  -o season.json

with open("season.json", "r", encoding = "utf-8") as f:
  data = json.load(f)["list"]

# タームごとに必要な情報だけを配列にまとめる
terms = []
for season in data:
  for id in data[season]:
    terms.append({"id" : id, "season" : data[season][id]["season"], "rule" : data[season][id]["rule"], "rst" : data[season][id]["rst"], "ts1" : data[season][id]["ts1"], "ts2" : data[season][id]["ts2"]})

当該シーズン・タームの対戦環境情報・ポケモンに関する情報

Jupyterを使ってコマンドに変数を埋め込むことでループを使って取得していく．ポケモン情報については5分割されているので，pdetail-{1から5}すべてを取得して連結する

for term in terms:
  if term["rule"] == 0:
    # 対戦環境の取得
    id = term["id"]
    rst = term["rst"]
    ts1 = term["ts1"]
    ts2 = term["ts2"]
    pokemons_file = str(id) + "-pokemons.json"
    !curl -XGET https://resource.pokemon-home.com/battledata/ranking/$id/$rst/$ts2/pokemon -H 'user-agent: Mozilla/5.0 (Linux; Android 8.0; Pixel 2 Build/OPD3.170816.012) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/83.0.4103.61 Mobile Safari/537.36'  -H 'accept: application/json' -o $pokemons_file
    for i in range(5):
      # 当該タームのポケモン情報を取得
      j = i + 1
      pokeinfo_file = str(id) + "-pokeinfo-" + str(j) + ".json"
      !curl -XGET https://resource.pokemon-home.com/battledata/ranking/$id/$rst/$ts2/pdetail-$j -H 'user-agent: Mozilla/5.0 (Linux; Android 8.0; Pixel 2 Build/OPD3.170816.012) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/83.0.4103.61 Mobile Safari/537.36'  -H 'accept: application/json' -o $pokeinfo_file

id = 10011 # シーズン1のシングルバトル

# 上位のポケモンリスト
with open(str(id) + "-pokemons.json", "r", encoding = "utf-8") as f: 
  top_pokemons = json.load(f)[:30]

# ポケモンの情報
pokemons_info = {}
for i in range(5):
  j = i + 1
  with open(str(id) + "-pokeinfo-" + str(j) + ".json") as f:
    data = json.load(f)
    for index in data.keys():
      pokemons_info[index] = data[index]

グラフの描画と次数ランキングの作成

必要なデータが集まったので，グラフを描画していきます．

nodes = [] # ノードリスト
edges = [] # エッジリスト

# トップ30匹を考慮
for pokemon in top_pokemons:
  # 当該ポケモンをノードに追加
  nodes.append(pokedex[pokemon["id"]])
  # 一緒に採用されるポケモンを追加
  with_poke_list = pokemons_info[str(pokemon["id"])][str(pokemon["form"])]["temoti"]["pokemon"]
  for with_poke in with_poke_list:
    nodes.append(pokedex[with_poke["id"]])
    if (pokedex[with_poke["id"]], pokedex[pokemon["id"]]) not in edges: # 順番による重複をなくす
      edges.append((pokedex[pokemon["id"]], pokedex[with_poke["id"]]))

nodes = list(set(nodes)) # 一度集合に変換して重複なしにする

import networkx as nx
import matplotlib.pyplot as plt
from matplotlib import font_manager

# 無効グラフ作成
G = nx.Graph()
G.add_nodes_from(nodes)
G.add_edges_from(edges)

# ネットワーク全体の次数の平均値を計算
average_deg = sum(d for n, d in G.degree()) / G.number_of_nodes()
# ノードの次数に比例するようにサイズを設定
sizes = [2500*deg/average_deg for node, deg in G.degree()]

plt.figure(figsize=(15, 15))
nx.draw(G, with_labels=True, node_color = "lightblue", node_size = sizes, edge_color = "black")
plt.axis('off')
plt.savefig(str(id) + '-graph.png')
plt.show()

# 次数の多い順番でnodeを並び替える
nodes_sorted_by_degree = sorted(G.degree(), key=lambda x: x[1], reverse=True)
print(nodes_sorted_by_degree)

# 次数ランキングの出力
x, height = list(zip(*nodes_sorted_by_degree))
plt.figure(figsize=(15, 15))
plt.xlabel("Number of degrees")
plt.ylabel("Name of node")
plt.barh(x, height)
plt.savefig(str(id) + '-degs.png')
plt.show()

結果

シーズン1から順番に示していきます．鎧の孤島解禁以降はノイズが乗っています．

シーズン1

剣盾環境開幕．次数的にはミミッキュ・ロトム・ドラパルト・ドリュウズなどが圧倒的に環境の中心に．グラフで見ると中心に9体のポケモンが固まっていてその周りを囲むようにほかのポケモンが描画されているので，その9体が環境の中心でSランクと考えてもよいと思います．

シーズン2

環境の変化が始まる．端から迫ってくる新しい並びが面白い．ヌルアント・ブリムオンドサイドン・ヌオードヒドなどが環境入りする．

シーズン3

キョダイマックスカビゴンとキョダイマックスラプラスが解禁．すぐさま環境に入ってくる．この段階だとラプラスは採用率は高いものの，環境上位のポケモンとの同時採用率はそこまで高くないため，実はこのグラフで考えるとSランクではなくAかBくらいの評価になりますね．

シーズン4

新たにキョダイマックスポケモンが何匹か解禁されるもさほど変化なし．環境にコータス・リザードンが入ってくる．

シーズン5

特に変化なし．

シーズン6

若干ナットレイが流行り始める．そのほかに変化特になし．

シーズン7

鎧の孤島解禁．エースバーン・ゴリランダーの夢特性の解禁とキョダイマックスの解禁により環境の中心に一気に食い込む．カビゴン・アーマーガアが環境の中心から外れ始める．ホルードが少し流行る．

シーズン8

ポケモンホームが解禁．エアームド・ハピナス・ナットレイ・ドヒドイデなどの受けループ気質のポケモンが環境入り．カビゴンは受けきられてしまうため，環境から外れる．ポリゴン2が猛威を振るい環境の中心に．起点作成のためのカバルドンの価値が上がり，カバルドンが環境の中心に入る．ちなみにCharizard-1は図鑑の問題でばぐっているウーラオスです．

シーズン9

環境の中心の境界がくっきりとし始める．ギャラドスとドリュウズが少し下火になる．受けループも環境の外側に出始める．

シーズン10

上位ポケモン禁止ルールになり，環境が激変．ウーラオス・パッチラゴン・リザードン・ナットレイ・ロトム・アシレーヌが環境に中心に．カビゴンとアイアントが環境に戻ってくる．モロバレル・ピクシーが環境に顔を出すようになる．

シーズン11

ゲンガーが環境の中心入り．その他特に変化なし．

シーズン12

冠の雪原解禁．ランドロス・テッカグヤ・カプ・レヒレ・サンダー・ウツロイドなどの準伝説が環境のほとんどを占める．
シーズン12現在のデータからは，現在の環境の中心はミミッキュ・エースバーン・ポリゴン2・ゴリランダー・テッカグヤ・カプレヒレ・ランドロス・サンダー・ウツロイドあたりだといえそう．

最後に

シーズン12現在のキャラランクを次数をもとに以下のように設定してみました．是非参考にしてみてください．理論に突っ込みどころ満載かもしれませんが，一つの指標として皆さんの役にたつとうれしいです．それではよきポケモンライフを！

追記: 似たようなことやってる方がいました
https://qiita.com/b_aka/items/9020e3237ff1a3e676e4

始め

彡(ﾟ)(ﾟ)「学外からアクセスする以外の方法で、Gaussianを利用したいなぁ」

彡(ﾟ)(ﾟ)「せや、計算機の所定のフォルダーにファイルを放り込んだら自動でかけさせたらええんや!」

彡(ﾟ)(ﾟ)「それならわいでもできるやろ()」

利用した物

CentOSにあるcron

python3およびそのライブラリー

実物

auto.py

import os
import subprocess

path = '******'
A=os.listdir(path)
C=[i for i in A if '.com' in i] #計算したいファイルのリストを作成
L=[i.replace('.log', '.com') for i in A if '.log' in i] #計算の終了したファイルのリストの作成
G=[i for i in C if i not in L] #計算していないファイルを抽出
command1 = 'ps ax |grep g09 | grep -v grep | wc -l'
command2 = 'cd ****** | g09 '
res = subprocess.check_output(command1, shell=True)
try:
    if res==b'0\n':
        print("計算可能です")
        subprocess.Popen(command2 + G[0] + '&', shell=True)
        print("計算を実行しました")

    else:
        print("計算中です")

except IndexError:
    print("計算するファイルがありません")

*/30 * * * * /usr/bin/tcsh ./autopy.sh >>&  /*******/analog-err.log

色々とあって直でpythonのプログラムを動かせなかったので、迂遠な方法を利用しています

詰まった所

gaussian等々のPathがcronでは通ってなかった事に気づくのに結構な時間を使った件について

一個一個Pathを設定するのが面倒だったので、

cronにprintenvの結果を張り付けるとかいう力技を行っています

だから、実際のcronの中身はもっと長いです

あと、pythonのps ax |grep g09 | grep -v grep | wc -lは自力で思いついたわけじゃないです

ネットでジョブの実行状況を知るにはどうしたらいいか、

的なので検索して引っかかった物を自分なりに弄っただけです

このプログラムにおける一番の難題はこの部分でしたね…

終わりに

私の環境ではこれで元気に動いています

ここをこうすればもっとスマートにいくのでは?みたいなのがあったら教えてください(乞食)

ただし、python使う意味ある?全部コマンドで良くね?みたいなことに気付いた勘の良いガキは嫌いだよ

どうして私の研究室にはジョブマネージャーが無いんですか(血涙)

次の目標として、ファイルの送信・転送について自動化できればいいなぁと思っています

転送は出来てるんですけど、送信がね……

最終的にはGUIを作成したいなぁと思っていたり…

それこそGo言語で何とかできないかなぁ

おしまい

異常検知について勉強したのでまとめておきます。

参考文献

下記文献を大いに参考にさせていただきました：
[1] Lukas Ruff, et al.：A Unifying Review of Deep and Shallow Anomaly Detection
[2] 井手,杉山：入門 機械学習による異常検知―Rによる実践ガイド
[3] 井手,杉山：異常検知と変化検知 (機械学習プロフェッショナルシリーズ)

§１ 異常検知の概要

異常検知の適応例

• 医療×異常検知
  医療用画像からの疾患部位の特定

出典：Thomas, et al.Unsupervised Anomaly Detection with Generative Adversarial Networks to Guide Marker Discovery

• 金融×異常検知
  マネーロンダリング検知

出典：Mark, et al.Anti-Money Laundering in Bitcoin: Experimenting with Graph Convolutional Networks for Financial Forensics

異常は多種多様にわたる

異常と一口に言っても多種多様な異常があり、大きく静的な異常と動的な異常に分けられます（私が勝手に作った切り口なのでおかしな点があればご指摘ください）。

静的な異常とは、各データ間に順序性等の関係がなく、データの値が異常の尺度になるものをいいます。

• 例えば，医療用画像において該当部位の濃淡が健康な場合のそれと大きく異なる場合に異常（疾患の存在）と判断されます。
• このような静的な異常は特に外れ値（Outlier）と呼んだりします。

動的な異常とは、時系列データのようにデータ間に順序性等の関係があり、データの振る舞いが異常の尺度になるものをいいます。

• 振る舞いが異常であるとは，同一区間での観測値の増減が激しかったり逆に変化が著しく小さくなったりと個々の問題に依存します。
• 静的な異常との大きな違いは，動的な異常は値そのものが異常であるとは限らない点です。
• 例えば，下グラフの中央部は異常のように見えますが、値自体は取り得る値です。このケースでは、「周波数の変化」という振る舞いに異常性を見出すことになります。

import numpy as np
import matplotlib.pyplot as plt

x1 = np.concatenate([np.linspace(-10, -3, 50), np.zeros(50), np.linspace(3, 9, 50)])
x2 = np.concatenate([np.zeros(50), np.linspace(-3, 3, 50), np.zeros(50)])
y = np.sin(x1) + np.sin(5*x2)
plt.plot(np.linspace(-10, 9, 150), y)

数式で表現すると、系列データ$\{x_1, ..., x_t\}$において各点$x_j$が$p(x)=\int p(x,t){\rm d}t$の意味で正常であっても、時間の条件付き確率$p(x_j \mid t) = \int \cdots \int p(x_1,...x_t \mid t){\rm d}t_1 \cdots {\rm d}t_{j-1}{\rm d}t_{j+1}\cdots{\rm d}t_t$においては異常となり得えます。

ただしこのケースにおいても、中央部が異常ではない場合もあります。 例えば、上グラフが工場の機械動作シグナルであり中央部だけオペレーションモードが切り替わっていたため異常ではないということも考えられます。

以上のように異常には色々な側面がありますが、文献[1]ではこのように異常の種類をまとめています：

出典:[1]

左上図のPoint Anomalyは観測エラー等による単発ノイズデータ（外れ値）に見えますが、Group Anomalyは外れ値というよりも別の分布から発生しているデータ群に見えます。両者では異なった異常検知アプローチが必要と想定されます。右下のSemantic Anomalyまで来ると古典的な異常検知手法では対応が難しく、ディープラーニングによるアプローチが必要になってきます。

このように「異常」に統一解はなく、ケースバイケースで「異常」を定義しなければいけません。

異常検知アルゴリズムは基本的に教師なし学習である

異常検知アルゴリズムは教師あり学習ではなく基本的に教師なし学習を用います。

理由は2つで：

1. 大概、異常データは少ないから
2. 異常データは多種多様であり網羅的にモデリングすることが難しいから

です。特に異常がクリティカルな現場（医療や工場）では、異常データを数百~数千も集められないはずです。また、理由２に関しては冒頭で説明したとおりです。「正常・異常」を当てるよりも、「正常・正常じゃない」を当てる方が簡単で筋が良い問題設定になります。後者は異常データに関して言及しておらず、正常データの構造理解のみで達成可能です。教師なし学習とは、サンプル$x_i(i=1,...,n)$から$p(x)$を推定することでデータの構造を理解する枠組みですから、異常検知は基本的に教師なし学習を用いることになります。ただし、教師なし学習を行う際のサンプルには異常データが含まれていない、もしくは圧倒的少数であると仮定できる必要があります。

余談ですが、文献[2]では以下のようなことを述べておりお気に入りです：

トルストイの『アンナ・カレーニナ』という小説に「幸せな家族は皆どこか似ているが，不幸せな家族はそれぞれ違う（Happy families are all alike; every unhappy family is unhappy in its own way）」という一節があり示唆的です.

異常検知アルゴリズムは4種類に大別される

それでは具体的に異常検知アルゴリズムの紹介をしていきます。

アルゴリズムは大きく4つのアプローチにより分類できます。それぞれ、Classification Approach, Probabilistic Approach, Reconstruction Approach, Distance Approachです。この4つの分類に加えて、ディープラーニングを用いるか否かで分類をした表がこちらになります：

出典：[1]

PCAやOC-SVM等のワードひとつひとつが過去提案されてきた異常検知アルゴリズムの名前です。各アルゴリズムの詳細は後続で紹介するので、ここでは4つのアプローチそれぞれの概要の紹介に留めます。（私自身、下記テーブルへの理解が十分ではないためご意見をいただけますと勉強になります。）

アプローチ	概要	Pros	Cons
Classification	データの正常・異常を分ける境界を直接推定するアプローチ	分布の形状を仮定する必要がない。また、SVMを用いた手法においては比較的少数データで学習可能。	正常データが複雑な構造をしていると上手く動作しない場合がある
Probabilistic	正常データの確率分布を推定し、生起確率が低い事象を異常とみなすアプローチ	データが仮定したモデルに従っていれば非常に高精度を達成	分布モデル等多くの仮定が必要。また高次元への対応も困難。
Reconstruction	正常データを再構成するようwell-trainされたロジックが上手く再構成できないデータは異常である、というアプローチ	近年盛んに研究されているAEやGAN等の強力な手法が適応できる。	あくまでも次元圧縮や再構成精度を高めるために目的関数がデザインされているため、必ずしも異常検知に特化した手法ではないこと。
Distance	データ間の距離に基づき異常検知を行うアプローチ。異常データは正常データから距離的に離れているという発想。	アイディアがシンプルで、出力への説明可能性も高い。	計算量が大きい（$O(n^2)$など）

上記テーブルは以下資料を参考に作成しました。
1. 比戸：Jubatus Casual Talks #2 異常検知入門
2. PANG, et al. Deep Learning for Anomaly Detection: A Review

§2 アルゴリズムの紹介と実験

紹介する手法まとめ

• Classification Approach
  • One-Class SVM
• Probability Approach
  • ホテリングのT2法
  • GMM
  • KDE
• Reconstruction Approach
  • AE
• Distance Approach
  • k-NN
  • LOF

Classification Approach

One-Class SVM

データ$D=\{x_1,...,x_N\}$が与えられているとします。ただし$D$には異常データが含まれていないか、含まれていたとしても圧倒的少数であるとします。

前章でClassification Approachは正常・異常を分ける境界を推定すると紹介しました。One-Class SVMは、正常データをすっぽり囲むような球を計算し境界とするアルゴリズムです。（球といっても三次元とは限りません。）

半径をとても大きくとれば$D$をすっぽり内包する球は作れますが、大きすぎる球ではFalse Negativeが1になってしまい嬉しくありません。そこで、

• 球の体積はできるだけ小さくする
• $D$のなかで、球からはみ出るデータがあることを許容する

の2点を許すことで汎用的な球を計算します。具体的には以下の最適化問題を解くことで最適な球の半径$R$と中心$b$が求まります：

\min_{R,b,u} \left\{R^2 + C\sum_{n=1}^{N}u_n \right\} \quad {\rm s.t.} \quad \|x_n - b\|^2 \le R^2 + u_n

ただし、$u$は球からはみ出した分のペナルティ項、$C$はそのペナルティをどれくらい重要視するかのパラメータです。SVMではカーネル関数によるデータの非線形変換が可能ですが、もちろんOne-Class SVMでも利用可能です。カーネル関数を適切に利用することで、適切な球が求まります。

数値実験をしてみます。データは以下のように準備しました。

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns

sns.set()
np.random.seed(123)

#平均と分散
mean1 = np.array([2, 2])
mean2 = np.array([-2,-2])
cov = np.array([[1, 0], [0, 1]])

# 正常データの生成(2つの正規分布から生成)
norm1 = np.random.multivariate_normal(mean1, cov, size=50)
norm2 = np.random.multivariate_normal(mean2, cov, size=50)
norm = np.vstack([norm1, norm2])

# 異常データの生成(一様分布から生成)
lower, upper = -10, 10
anom = (upper - lower)*np.random.rand(20,20) + lower

s=8
plt.scatter(norm[:,0], norm[:,1], s=s, color='green', label='normal')
plt.scatter(anom[:,0], anom[:,1], s=s, color='red', label='anomaly')
plt.legend()

from sklearn import svm
xx, yy = np.meshgrid(np.linspace(-10, 10, 500), np.linspace(-10, 10, 500))

# fit the model
clf = svm.OneClassSVM(nu=0.1, kernel="rbf", gamma=0.1)
clf.fit(norm)

def draw_boundary(norm, anom, clf):
    # plot the line, the points, and the nearest vectors to the plane
    Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)

    plt.title("One-Class SVM")
    plt.contourf(xx, yy, Z, levels=np.linspace(Z.min(), Z.max(), 7), cmap=plt.cm.PuBu)
    bd = plt.contour(xx, yy, Z, levels=[0], linewidths=2, colors='darkred')
    plt.contourf(xx, yy, Z, levels=[0, Z.max()], colors='palevioletred')

    nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)
    am = plt.scatter(anom[:, 0], anom[:, 1], c='red', s=s)

    plt.axis('tight')
    plt.xlim((-10, 10))
    plt.ylim((-10, 10))
    plt.legend([bd.collections[0], nm, am],
               ["learned frontier", "normal", "anomaly"],
               loc="upper left",
               prop=matplotlib.font_manager.FontProperties(size=11))
    plt.show()

draw_boundary(norm, anom, clf)

ちなみに、パラメータ$C$を小さくする、すなわち球からはみ出すペナルティを軽くみるようにするとこのように境界が変わります。（sklearnではパラメータ$\nu$で表していますが、どうやら分母にいるらしいので$\nu$を小さくする$=C$を大きくするのようです）

# fit the model
clf = svm.OneClassSVM(nu=0.5, kernel="rbf", gamma=0.1)
clf.fit(norm)
draw_boundary(norm, anom, clf)

他にもOne-Class NN等が提案されています。One-Class SVMではカーネル関数を調整して適切な非線形変換を探索していましたが、そこをNNで代用＆自動化しようというアイディアです。実装はめんどくさそうなのでやってません。

Probability Approach

ホテリングのT2法

ホテリングのT2法は古典的な手法であり、理論的な解析がなされています。 結論、正規分布を仮定することで異常度（Anomaly Score）が$\chi^2$分布に従う性質を利用したアルゴリズムになります。

ここでは簡単のため1次元データで説明します。 例によってデータ$D$が与えられているとします。ここで各データ$x$が正規分布

p(x\mid \mu, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{1}{2\sigma^2}(x-\mu)^2\right\}

に従うものと仮定します。この正規分布を最尤推定法によりデータへフィッティングすることで

\begin{align*}
    \hat{\mu} &= \frac{1}{N}\sum_{n=1}^{N}x_n \\
    \hat{\sigma}^2 &= \frac{1}{N}\sum_{n=1}^{N}(x_n - \hat{\mu})^2
\end{align*}

と各パラメータが推定されます。つまり、データ$D$は正規分布$p(x\mid \hat{\mu}, \hat{\sigma}^2)$に従うと推定したわけです。

さて、ここで異常度（Anomaly Score）を導入します。異常度とは、文字通り異常の度合いを示す統計量のことで、よく負の対数尤度が用いられます。

prob = np.linspace(0.01, 1, 100)
neg_log_likh = [-np.log(p) for p in prob]

plt.plot(prob, neg_log_likh)
plt.title('Negative Log Likelihood')
plt.ylabel('Anomaly Score')
plt.xlabel('Prob.')

グラフのように、負の対数尤度は確率が小さいところで大きな値をとる形状をしており、異常度として機能しそうです。今推定した正規分布の異常度を計算すると

-\log p(x\mid \hat{\mu}, \hat{\sigma}^2) = \frac{1}{2\hat{\sigma}^2}(x-\hat{\mu})^2 + \frac{1}{2}\log(2\pi\hat{\sigma}^2)

となります。異常度はデータ$x$に対して定義される量ですから、$x$に関係ない項は除いたうえで異常度$a(x)$は

a(x) = \left( \frac{x - \hat{\mu}}{\hat{\sigma}} \right)^2

と計算されます。異常度の分子は、「データが平均から離れているほど異常度が高くなる」ということを表しています。分母は、「どれくらいデータが散らばっているか？」を表しており、散らばっているデータでは平均から遠く離れていようと大目に見る、逆に密なデータでは少しのズレも許さないというお気持ちがあります。

実はこの異常度$a(x)$はサンプルサイズ$N$が十分大きいときに自由度$1$のカイ二乗分布$\chi^2_1$に従うことがわかっています。

from scipy.stats import chi2
x = np.linspace(0, 8, 1000)
for df in [1,2,3,4,5]:
    plt.plot(x, chi2.pdf(x ,df = df), label='deg. of freedom:'+str(df))
plt.legend()
plt.ylim(0, 0.5)
plt.title('Chi-Square Dist.')

あとは有意水準$5$%で閾値$a_{th}$を決めるなどしてあげれば異常検知が可能となります。まとめると、

1. データが正規分布に従うと仮定する。
2. 正規分布のパラメータ$\mu$, $\sigma$を最尤推定法により推定する
3. カイ二乗分布から、正常・異常を分ける点$a_{th}$を決める
4. 新たなデータ$x_{new}$が、$a(x_{new})>a_{th}$だったら異常とみなす

がホテリングのT2法です。今回は1次元で説明しましたが、多次元でも同様に異常度がカイ二乗分布に従います。ただし、自由度=次元数です。

では実験してみましょう。データは以下のように作成しました。

norm = np.random.normal(0, 1, 20) # 正常データ
anom = np.array([-3, 3]) # 異常データ
plt.scatter(norm, [0]*len(norm), color='green', label='normal')
plt.scatter(anom, [0]*len(anom), color='red', label='anomaly')

# 最尤推定
N = len(norm)
mu_hat = sum(norm)/N
s_hat = sum([(x-mu_hat)**2 for x in norm])/N

# カイ二乗分布の分位点を閾値とする
a_th_5 = chi2.ppf(q=0.95, df=1)
a_th_30 = chi2.ppf(q=0.7, df=1)

# 各データの異常度を計算
data = np.hstack([norm, anom])
anom_score = []
for d in data:
    score = ((d - mu_hat)/s_hat)**2
    anom_score.append(score)

#描画
norm_pred_5 = [d for d,a in zip(data, anom_score) if a<=a_th_5]
anom_pred_5 = [d for d,a in zip(data, anom_score) if a>a_th_5]
norm_pred_30 = [d for d,a in zip(data, anom_score) if a<=a_th_30]
anom_pred_30 = [d for d,a in zip(data, anom_score) if a>a_th_30]

fig, axes = plt.subplots(1,3, figsize=(15,5))

axes[0].scatter(norm, [0]*len(norm), color='green', label='normal')
axes[0].scatter(anom, [0]*len(anom), color='red', label='anomaly')
axes[0].set_title('Truth')
axes[1].scatter(norm_pred_5, [0]*len(norm_pred_5), color='green', label='normal pred')
axes[1].scatter(anom_pred_5, [0]*len(anom_pred_5), color='red', label='anomaly pred')
axes[1].set_title('Threshold:5%')
axes[2].scatter(norm_pred_30, [0]*len(norm_pred_30), color='green', label='normal pred')
axes[2].scatter(anom_pred_30, [0]*len(anom_pred_30), color='red', label='anomaly pred')
axes[2].set_title('Threshold:30%')
plt.legend()

異常度の閾値を小さくすれば、異常と見なされる区域が広がっていくことが見て取れます。

閾値は問題に応じて決める必要があります。例えば、医療の現場など異常の見逃しがクリティカルな場合は閾値は小さく設定すべきですし、データクレンジング等の場合はもう少し大きい閾値にしてもいいかもしれません。

GMM(Gaussian Mixture Model)

ホテリングのT2法ではデータに正規分布を仮定していましたが、これが不適切なケースもあるでしょう。例えば以下のようにデータが複数個所に集中している（単峰ではなく多峰）場合に正規分布でフィッティングすると

lower, upper = -3, -1
data1 = (upper - lower)*np.random.rand(20) + lower
lower, upper = 1, 3
data2 = (upper - lower)*np.random.rand(30) + lower
data = np.hstack([data1, data2])

# 最尤推定
N = len(data)
mu_hat = sum(data)/N
s_hat = sum([(x-mu_hat)**2 for x in data])/N

# プロット
from scipy.stats import norm
X = np.linspace(-7, 7, 1000)
Y = norm.pdf(X, mu_hat, s_hat)

plt.scatter(data, [0]*len(data))
plt.plot(X, Y)

のように、全くデータが観測されていない箇所に正規分布のピークが来るように推定されてしまいます。こういった場合は正規分布ではなく、もう少し複雑な分布を設定するのが良いでしょう。今回はその一例としてGMM(Gaussian Mixture Model, 混合正規分布)による異常検知を紹介します。

GMMとは正規分布をいくつか足し合わせることで多峰型を表現する分布です。足し合わせる個数はよく$K$で表現され、コンポーネント数と呼ばれます。GMMの数式とビジュアルは以下の通りです(1次元にしています)。

p(x\mid \pi, \mu, \sigma) = \sum_{k=1}^{K}\pi_k N(x\mid \mu_k, \sigma_k)

ただし、

0 \le\pi_k\le 1,\quad \sum_{k=1}^{K}\pi_k=1,\quad N(x\mid\mu,\sigma)={\text 正規分布}

です。

# プロット
from scipy.stats import norm
X = np.linspace(-7, 7, 1000)

Y1 = norm.pdf(X, 1, 1)
Y2 = norm.pdf(X, -3, 2)
Y = 0.5*Y1 + 0.5*Y2 

plt.plot(X, Y1, label='Component1')
plt.plot(X, Y2, label='Component2')
plt.plot(X, Y,  label='GMM')
plt.legend()

各正規分布の平均($\mu$)・分散($\sigma$)と、混合率$\pi$はデータから学習されるパラメータです。しかしコンポーネント数$K$は決め打ちする必要があります。このように分析者側で決め打ちにするパラメータをハイパーパラメータと呼びます。

ではGMMによる異常検知を実験してみます。コンポーネント数$K$によってどのように推定分布の形状が変わるのかも観察しましょう。データはOne-Class SVMと同様以下を使用します。

def draw_GMM(norm, anom, clf, K):
    # plot the line, the points, and the nearest vectors to the plane
    xx, yy = np.meshgrid(np.linspace(-10, 10, 1000), np.linspace(-10, 10, 1000))
    Z = np.exp(clf.score_samples(np.c_[xx.ravel(), yy.ravel()]))
    Z = Z.reshape(xx.shape)

    plt.title("GMM:K="+str(K))
    plt.contourf(xx, yy, Z, levels=np.linspace(Z.min(), Z.max(), 7), cmap=plt.cm.PuBu)

    nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)
    am = plt.scatter(anom[:, 0], anom[:, 1], c='red', s=s)

    plt.axis('tight')
    plt.xlim((-10, 10))
    plt.ylim((-10, 10))
    plt.legend([nm, am],
               ["normal", "anomaly"],
               loc="upper left",
               prop=matplotlib.font_manager.FontProperties(size=11))
    plt.show()

# GMMによるフィッティング
from sklearn.mixture import GaussianMixture

for K in [1,2,3,4]:
    clf = GaussianMixture(n_components=K, covariance_type='full')
    clf.fit(norm)
    draw_GMM(norm, anom, clf, K)

このようにコンポーネント数を多くすると、データの細部の傾向まで捉えていることがわかります。しかし、細部まで捉えている＝汎化性能が高いわけではないので注意が必要です。コンポーネント数はWAIC等の情報量規準を用いて決定することができます（とりあえずWAICを使ってみました：GitHub）

KDE(Kernel Density Estimation)

今までは正規分布や混合正規分布を仮定し、データへのフィッティングを通じて内部パラメータの推定を行ってきました。このようなアプローチをパラメトリックなアプローチといいます。しかし、パラメトリックな分布では上手く表現できないような場合にはどうしたら良いでしょうか？

ノンパラメトリックなアプローチはそのような場合に有効であり、KDEは代表的なノンパラメトリック推定手法です。KDEは、各データ点にカーネル（通常はGaussian）を重ねることで確率分布を推定する手法です。可視化すると以下のようになります。

from scipy.stats import norm
np.random.seed(123)
lower, upper = -1, 5
data = (upper - lower)*np.random.rand(3) + lower
for i, d in enumerate(data):
    X = np.linspace(lower,upper,500)
    Y = norm.pdf(X, loc=d, scale=1)
    knl = plt.plot(X,Y,color='green', label='kernel'+str(i+1))

# KDE
n = 3
h = 0.24
Y_pdf = []
for x in X:
    Y = 0
    for d in data:
        Y += norm.pdf((x-d)/h, loc=0, scale=1)
    Y /= n*h
    Y_pdf.append(Y)
kde = plt.plot(X,Y_pdf, label='KDE')
plt.scatter(data, [0]*len(data), marker='x')
plt.ylim([-0.05, 1])
plt.legend()

このようにカーネルの重ね合わせで確率分布を推定しています。カーネルに正規分布を用いる場合は分散（バンド幅と呼びます）を決め打ちする必要があります。後続の実験では、このバンド幅を変えた時の推定分布の変化も観察しましょう。

KDEの理論詳細および長所・短所は鈴木.データ解析第十回「ノンパラメトリック密度推定法」がわかりやすかったです：

推定分布

p(x) = \frac{1}{Nh}\sum_{n=1}^{N}K\left(\frac{x-X_n}{h}\right)

ただし、$h>0$をバンド幅、$K()$をカーネル関数、$X_n$は各データ点とする。
長所・短所


出典：鈴木.データ解析第十回「ノンパラメトリック密度推定法」

ではKDEによる異常検知を実験してみましょう。データは例によって

を使用します。バンド幅を変えた時の推定分布は以下のようになりました。

def draw_KDE(norm, anom, clf, K):
    # plot the line, the points, and the nearest vectors to the plane
    xx, yy = np.meshgrid(np.linspace(-10, 10, 1000), np.linspace(-10, 10, 1000))
    Z = np.exp(clf.score_samples(np.c_[xx.ravel(), yy.ravel()]))
    Z = Z.reshape(xx.shape)

    plt.title("KDE:BandWidth="+str(bw))
    plt.contourf(xx, yy, Z, levels=np.linspace(Z.min(), Z.max(), 7), cmap=plt.cm.PuBu)

    nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)
    am = plt.scatter(anom[:, 0], anom[:, 1], c='red', s=s)

    plt.axis('tight')
    plt.xlim((-10, 10))
    plt.ylim((-10, 10))
    plt.legend([nm, am],
               ["normal", "anomaly"],
               loc="upper left",
               prop=matplotlib.font_manager.FontProperties(size=11))
    plt.show()

from sklearn.neighbors import KernelDensity
for bw in [0.2,0.5,1.0,2.0]:
    clf = KernelDensity(bandwidth=bw, kernel='gaussian')
    clf.fit(norm)
    draw_KDE(norm, anom, clf, bw)

このように、バンド幅が大きくなると推定分布が滑らかになることがわかります。また、カーネル関数は正規分布に限らず他にもあるので、「バンド幅・カーネル関数」の両者を適切に設定してあげる必要があります。

Reconstruction Approach

Reconstruction Approachは画像データに用いられることが多いので、ここでも画像データを使用します。正常データをMNIST、異常データをFashion-MNISTとします。

# MNISTのロード
import torch
from torch.utils.data import DataLoader
from torchvision import datasets, transforms

# 学習データの読み込み
train_mnist = datasets.MNIST(
    root="./data", 
    train=True, 
    transform=transforms.Compose([
               transforms.ToTensor(),
               transforms.Normalize((0.5,), (0.5,)) # 区間[-1,1]内で平均、標準偏差 = 0.5, 0.5に正規化
               ]), 
    download=True
)

train_mnist_loader = torch.utils.data.DataLoader(
    train_mnist, batch_size=64, shuffle=True
)

# テストデータの読み込み
test_mnist = datasets.MNIST(
    root="./data", 
    train=False, 
    transform=transforms.Compose([
               transforms.ToTensor(),
               transforms.Normalize((0.5,), (0.5,)) # 区間[-1,1]内で平均、標準偏差 = 0.5, 0.5に正規化
               ]), 
    download=True
)

test_mnist_loader = torch.utils.data.DataLoader(
    test_mnist, batch_size=64, shuffle=True
)

# Fashion-MNISTのロード
# 学習データの読み込み
train_fashion = datasets.FashionMNIST(
    root="./data", 
    train=True, 
    transform=transforms.Compose([
               transforms.ToTensor(),
               transforms.Normalize((0.5,), (0.5,)) # 区間[-1,1]内で平均、標準偏差 = 0.5, 0.5に正規化
               ]), 
    download=True
)

train_fashion_loader = torch.utils.data.DataLoader(
    train_fashion, batch_size=64, shuffle=True
)

# テストデータの読み込み
test_fashion = datasets.FashionMNIST(
    root="./data", 
    train=False, 
    transform=transforms.Compose([
               transforms.ToTensor(),
               transforms.Normalize((0.5,), (0.5,)) # 区間[-1,1]内で平均、標準偏差 = 0.5, 0.5に正規化
               ]), 
    download=True
)

test_fashion_loader = torch.utils.data.DataLoader(
    test_fashion, batch_size=64, shuffle=True
)

import matplotlib.pyplot as plt

# 画像を一気に並べて可視化する関数
def imshow(image_set, nrows=2, ncols=10):
    plot_num = nrows * ncols
    fig, ax = plt.subplots(ncols=ncols, nrows=nrows, figsize=(15, 15*nrows/ncols))
    ax = ax.ravel()
    for i in range(plot_num):
        ax[i].imshow(image_set[i].reshape(28, 28), cmap='gray')
        ax[i].set_xticks([])
        ax[i].set_yticks([])

# MNISTのサンプル可視化
imshow(train_mnist.data)

# Fashion-MNISTのサンプル可視化
imshow(train_fashion.data)

AE(Auto Encoder)

AE(Auto Encoder, オートエンコーダ)とは、NN（Neural Network, ニューラルネットワーク）の一種であり、入力と出力が等しくなるように学習するNNです。ただの恒等関数ではつまらないので、通常は入力次元よりも小さい次元の中間層を挟みます。

これは何をしているのかというと、中間層で入力データを「要約」しています。数理的には、入力データを再構成できるだけの情報を保持しつつ低次元空間へマッピングしている、ということになります。

AEを用いた異常検知は、「正常データを再構成するよう訓練されたAEは、異常データを上手く再構成できないのではないか？」という思想に基づきます。

では、AEによる異常検知を実験しましょう。AEをMNISTで訓練した後に、MNISTおよびFashion-MNISTの再構成精度を異常度として検知を試みます。

import torch.nn as nn
import torch.nn.functional as F

# AEの定義
class AE(nn.Module):
    def __init__(self):
        super(AE, self).__init__()
        self.fc1 = nn.Linear(28*28, 256)
        self.fc2 = nn.Linear(256, 28*28)

    def forward(self, x):
        x = torch.tanh(self.fc1(x))
        x = torch.tanh(self.fc2(x))
        return x
ae = AE()

# 損失関数と最適化手法の定義
import torch.optim as optim

criterion = nn.MSELoss()
optimizer = optim.Adam(ae.parameters(), lr=0.0001)

# 学習ループ
N_train = train_mnist.data.shape[0]
N_test = test_mnist.data.shape[0]

train_loss = []
test_loss = []
for epoch in range(10):  
    # Train step
    running_train_loss = 0.0
    for data in train_mnist_loader:
        inputs, _ = data
        inputs = inputs.view(-1, 28*28)

        # zero the parameter gradients
        optimizer.zero_grad()

        # forward + backward + optimize
        outputs = ae(inputs)
        loss = criterion(outputs, inputs)
        loss.backward()
        optimizer.step()

        running_train_loss += loss.item()
    running_train_loss /= N_train
    train_loss.append(running_train_loss)

    # Test step
    with torch.no_grad():
        running_test_loss = 0
        for data in test_mnist_loader:
            inputs, _ = data
            inputs = inputs.view(-1, 28*28)
            outputs = ae(inputs)
            running_test_loss += criterion(outputs, inputs).item()
        running_test_loss /= N_test
        test_loss.append(running_test_loss)
    if (epoch+1)%1==0:    # print every 1 epoch
        print('epoch:{}, train_loss:{}, test_loss:{}'.format(epoch + 1, running_train_loss, running_test_loss))

plt.plot(train_loss, label='Train loss')
plt.plot(test_loss, label='Test loss')
plt.ylabel('Loss')
plt.xlabel('epoch')
plt.legend()

# 再構成の様子を可視化してみる
with torch.no_grad():
    for data in test_mnist_loader:
        inputs, _ = data
        inputs = inputs.view(-1, 28*28)
        outputs = ae(inputs)

n = 5 # num to viz
fig, axes = plt.subplots(ncols=n, nrows=2, figsize=(15, 15*2/n))
for i in range(n):
    axes[0, i].imshow(inputs[i].reshape(28, 28), cmap='gray')
    axes[1, i].imshow(outputs[i].reshape(28, 28), cmap='gray')
    axes[0, i].set_title('Input')
    axes[1, i].set_title('Reconstruct')
    axes[0, i].set_xticks([])
    axes[0, i].set_yticks([])
    axes[1, i].set_xticks([])
    axes[1, i].set_yticks([])

きれいに再構成できていそうです。

それでは、MNISTとFashion-MNISTの再構成精度（今回は二乗誤差）の分布を描いてみます。AEによる異常検知が成功していれば、Fashion-MNISTの精度は低くなるはずです。

# 再構成精度の算出
mnist_loss = []
with torch.no_grad():
    for data in test_mnist_loader:
        inputs, _ = data
        for data_i in inputs:
            data_i = data_i.view(-1, 28*28)
            outputs = ae(data_i)
            loss = criterion(outputs, data_i)
            mnist_loss.append(loss.item())

fashion_loss = []
with torch.no_grad():
    for data in test_fashion_loader:
        inputs, _ = data
        for data_i in inputs:
            data_i = data_i.view(-1, 28*28)
            outputs = ae(data_i)
            loss = criterion(outputs, data_i)
            fashion_loss.append(loss.item())

plt.hist([mnist_loss, fashion_loss], bins=25, label=['MNIST', 'Fashion-MNIST'])
plt.xlabel('Loss')
plt.ylabel('data count')
plt.legend()
plt.show()

確かにFashion-MNISTの再構成は上手くいってなさそうなので、再構成精度による異常検知は成立しそうです。

最後に、Fashion-MNISTの中でもLossが低いもの（＝異常と判定し辛いもの）と高いもの（＝簡単に異常と判定できるもの）をピックアップして描画してみます。

lower_loss, upper_loss = 0.1, 0.8
data_low = [] # Lossが低いデータを格納
data_up = []  # Lossが高いデータを格納

with torch.no_grad():
    for data in test_fashion_loader:
        inputs, _ = data
        for data_i in inputs:
            data_i = data_i.view(-1, 28*28)
            outputs = ae(data_i)
            loss = criterion(outputs, data_i).item()
            if loss<lower_loss:
                data_low.append([data_i, outputs, loss])
            if loss>upper_loss:
                data_up.append([data_i, outputs, loss])

# ランダムに描画する用のインデックス
np.random.seed(0)

n = 3 #描画データ数
lower_indices = np.random.choice(np.arange(len(data_low)), n)
upper_indices = np.random.choice(np.arange(len(data_up)), n)
indices = np.hstack([lower_indices, upper_indices])

fig, axes = plt.subplots(ncols=n*2, nrows=2, figsize=(15, 15*2/(n*2)))

# Lossが小さいデータの描画
for i in range(n):
    inputs = data_low[indices[i]][0]
    outputs = data_low[indices[i]][1]
    loss = data_low[indices[i]][2]
    axes[0, i].imshow(inputs.reshape(28, 28), cmap='gray')
    axes[1, i].imshow(outputs.reshape(28, 28), cmap='gray')
    axes[0, i].set_title('Input')
    axes[1, i].set_title('Small Loss={:.2f}'.format(loss))
    axes[0, i].set_xticks([])
    axes[0, i].set_yticks([])
    axes[1, i].set_xticks([])
    axes[1, i].set_yticks([])

# Lossが大きいデータの描画
for i in range(n, 2*n):
    inputs = data_up[indices[i]][0]
    outputs = data_up[indices[i]][1]
    loss = data_up[indices[i]][2]
    axes[0, i].imshow(inputs.reshape(28, 28), cmap='gray')
    axes[1, i].imshow(outputs.reshape(28, 28), cmap='gray')
    axes[0, i].set_title('Input')
    axes[1, i].set_title('Large Loss={:.2f}'.format(loss))
    axes[0, i].set_xticks([])
    axes[0, i].set_yticks([])
    axes[1, i].set_xticks([])
    axes[1, i].set_yticks([])
plt.tight_layout()

Lossが低い画像は黒い面積が大きい傾向が見られる一方で、Lossが高い画像は白い面積が大きい傾向が見られます。MNISTは基本的に黒の面積が大きい画像ですので、白の面積が大きいとAEにとっては「異質」になり再構成が上手くいっていないと思われます。

Distance Approach

k-NN(k-Nearest Neighbor)

k-NNはデータ間の距離に基づき異常度を測る手法です。k-NNでは近傍のk個のデータ点を含むような円を考えます。下図は$k=5$の場合です。

図のように、「異常データが描く円は正常データよりも大きくなるはず」とう思想に基づきk-NNでは円の半径を異常度として定義します。ただし、近傍点の個数$k$と正常/異常を分ける半径の閾値は適切にチューニングする必要があります。

それでは数値実験を行います。閾値は、正常データの半径値分布における90%分位点とします。

# データの作成
np.random.seed(123)

#平均と分散
mean = np.array([2, 2])
cov = np.array([[3, 0], [0, 3]])

# 正常データの生成
norm = np.random.multivariate_normal(mean1, cov, size=50)

# 異常データの生成(一様分布から生成)
lower, upper = -10, 10
anom = (upper - lower)*np.random.rand(20,20) + lower

s=8
plt.scatter(norm[:,0], norm[:,1], s=s, color='green', label='normal')
plt.scatter(anom[:,0], anom[:,1], s=s, color='red', label='anomaly')
plt.legend()

k-NNによる異常検知を、どのようにライブラリで完結できるか分からなかったので自作します。。。

# 全データのk近傍円の半径を計算する関数
def kNN(data, k=5):
    """
    input: data:list, k:int
    output: 全データのk近傍円の半径:list

    距離はユークリッド距離を採用します。
    """
    output = []
    for d in data:
        distances = []
        for p in data:
            distances.append(np.linalg.norm(d - p))
        output.append(sorted(distances)[k])
    return output

# 閾値より大きければ1, 小さければ0を返す関数
def kNN_scores(point, data, thres, K):
    distances = []
    for p in data:
        distances.append(np.linalg.norm(point - p))
    dist = sorted(distances)[K]
    if dist > thres:
        return 1
    else:
        return 0

def draw_KNN(norm, anom, thres, K, only_norm=False):
    # plot the line, the points, and the nearest vectors to the plane
    xx, yy = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
    Z = [kNN_scores(point, norm, thres, K) for point in np.c_[xx.ravel(), yy.ravel()]]
    Z = np.array(Z).reshape(xx.shape)

    plt.title("KNN:K="+str(K))
    plt.contourf(xx, yy, Z)

    if not only_norm:
        nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)
        am = plt.scatter(anom[:, 0], anom[:, 1], c='red', s=s)

        plt.axis('tight')
        plt.xlim((-10, 10))
        plt.ylim((-10, 10))
        plt.legend([nm, am],
                   ["normal", "anomaly"],
                   loc="upper left",
                   prop=matplotlib.font_manager.FontProperties(size=11))

    if only_norm:
        nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)

        plt.axis('tight')
        plt.xlim((-10, 10))
        plt.ylim((-10, 10))
        plt.legend([nm],
                   ["normal"],
                   loc="upper left",
                   prop=matplotlib.font_manager.FontProperties(size=11))

    plt.show()

for k in [1,5,30]:
    anomaly_scores = kNN(norm, k=k)
    thres = pd.Series(anomaly_scores).quantile(0.90)
    draw_KNN(norm, anom, thres, k)

kが大きくなるにつれて、正常/異常の境界が滑らかになっていることがわかります。

さて、ここで少しいじわるな例で実験してみましょう。以下のように、正常データが密度の違う2つの正規分布から生成されているとします。

np.random.seed(123)

#平均と分散
mean1 = np.array([3, 3])
mean2 = np.array([-5,-5])
cov1 = np.array([[0.5, 0], [0, 0.5]])
cov2 = np.array([[3, 0], [0, 3]])

# 正常データの生成(2つの正規分布から生成)
norm1 = np.random.multivariate_normal(mean1, cov1, size=50)
norm2 = np.random.multivariate_normal(mean2, cov2, size=10)
norm = np.vstack([norm1, norm2])

s=8
plt.scatter(norm[:,0], norm[:,1], s=s, color='green', label='normal')
plt.legend()

for k in [1,5,30]:
    anomaly_scores = kNN(norm, k=k)
    thres = pd.Series(anomaly_scores).quantile(0.90)
    draw_KNN(norm, [], thres, k, only_norm=True)

特にk=1のケースに注目すると、正常データが密度の異なる別クラスターから発生している場合kNNが上手く動作していないことがわかります。「密なクラスターからのズレは重要視して、疎なクラスターからのズレは大目に見る」という処理を加えれば対処できそうだと想像できます。

実はそれが次節で紹介するLOFという手法になります。

LOF(Local Outlier Factor, 局所外れ値因子法)

前節で紹介したk-NNが上手く動作しないケースを図で表すと

となります。もちろん、左側の異常度を高くしたいのですが、LOFではデータの密度という概念を導入しそれを実現します。数理的詳細は割愛しますが、以下のようなことをしています。

自身のスカスカ度合いと近傍点のスカスカ度合いの比で異常度を定義しています。

それではLOFの実験を行いましょう。

def draw_LOF(norm, anom, K, only_norm=False):
    # plot the line, the points, and the nearest vectors to the plane
    xx, yy = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
    Z = - model.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = np.array(Z).reshape(xx.shape)

    plt.title("LOF:K="+str(K))
    plt.contourf(xx, yy, Z)

    if not only_norm:
        nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)
        am = plt.scatter(anom[:, 0], anom[:, 1], c='red', s=s)

        plt.axis('tight')
        plt.xlim((-10, 10))
        plt.ylim((-10, 10))
        plt.legend([nm, am],
                   ["normal", "anomaly"],
                   loc="upper left",
                   prop=matplotlib.font_manager.FontProperties(size=11))

    if only_norm:
        nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)

        plt.axis('tight')
        plt.xlim((-10, 10))
        plt.ylim((-10, 10))
        plt.legend([nm],
                   ["normal"],
                   loc="upper left",
                   prop=matplotlib.font_manager.FontProperties(size=11))

    plt.show()

from sklearn.neighbors import LocalOutlierFactor

for k in [1, 5, 30]:
    model = LocalOutlierFactor(n_neighbors=k, novelty=True)
    model.fit(norm)
    draw_LOF(norm, [], k, only_norm=True)

このように、クラスターの密度を考慮したうえで異常検知ができていそうです。

発展的な異常検知問題

今後異常検知の勉強をしていくにあたって、以下のようなトピックを勉強してきたいと思っています：

• 時系列データの異常検知
• グラフデータの異常検知
  • マネーロンダリング検知に使えそう