機械学習(pcaとkmeans)による画像のグルーピング

はじめに

自分で集めた画像でCNNするために、TwitterAPIの検索機能を使って画像を集めています。
集めている画像は特定の作品のキャラクターだったりするわけですが、CNNで分類モデルを作る上で画像にラベルを付けなくてはなりません。
このラベル付けは言うまでもなく面倒で、時間がかかります。そこである程度自動化できないかと思ったわけです。
したがって、今回はタイトルにある通り、pca(主成分分析)とkmeans(クラスタリング)を用いて画像のグルーピングを行いたいと思います。

データ

キャラクターの画像は顔画像だけにあらかじめトリミングしてあるものを用います。
その手法を知りたい方は以下を参照ください。
古い記事ですが、lbpcascade_animeface.xmlGoogleで検索した限りかなり広く使われています。
検索すればいくらでも出てくるかと思いますので、オリジナルのものを今回は添付させていただきました。

  • ブログ記事

ultraist.hatenablog.com

github.com

問題

では、具体的に画像のラベル付けで何が面倒かを以下に挙げます。

  1. キャラクターの顔認識精度に限界があり、誤検知が発生するため、顔画像以外の画像がデータ内に混在する

  2. 特定のキャラクター1人のみラベルをつける場合、1枚1枚人が判定するのは時間がかかる

  3. Twitterから集めているデータのため、ほぼ同じ画像がいくつか集まってしまう

1番の問題

lbpcascade_animeface.xmlは非常に良くできていますが、さすがに完璧な検知精度はありません。
opencv側のパラメータをいじることで、ある程度は見落とさないよう検知したり過検知することは防げますが、それにも限界はあります。
そのため、顔画像ではない画像が一つのクラスタとしてグルーピングされれば非常に手間が省けます。
結論を先に述べると、残念ながら本記事ではこの問題は解決できておりませんのでご注意ください

2番の問題

今回の場合、2021年2月現在、話題のゲーム「原神」のキャラクターの画像をTwitterから集めているわけですが、この中で筆者が欲しい画像はパイモンだけです。
なぜパイモンの画像ばかり集めているかは別の記事で書くとして、パイモン以外の画像はすべて不要ということになります。
クラスタリングをすることで、あるクラスタに一人もパイモンがいなければ、そのクラスタディレクトリ名で保存されている画像は不必要な画像として扱えます。
特定のキャラクターの画像を集める場合、こうすることで非常に効率が良くなることがわかるかと思います。

3番の問題

TwitterAPIの検索機能を使って、例えば「パイモン 原神」のように画像を集めるわけですが、どうしても同じ画像ばかり収集してしまうことがあります。
これを問題視するかどうかはCNNの学習に関わるかと思いますが、削除したいモチベーションがあるなら、これもクラスタリングで解決できます。 あるクラスタに同じ画像が集まれば、その画像の内、一枚を除いて削除すれば画像の重複がなくなるからです。

ps

これらの問題はクラスタリングすることでラベル付けの手間を省くことに注力しているため、最終的には人の判断が介入します。

プログラム

全体

import glob as gb
import shutil
import cv2
import os
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import seaborn as sns
from mpl_toolkits.mplot3d import Axes3D


# 画像の読み込み
def load_image(path):
    image_list = []
    npy_image_list = []
    image_path_list = gb.glob(path)

    # 画像データを一枚ずつ読み込む
    for i, image in enumerate(image_path_list):
        image_path = image.replace(chr(92), '/') # \を/に置換(windows特有)->macはchr(165)
        if i % 100 == 0: # 雑に進行状況出力
            print(i)
        
        img_npy = cv2.imread(image_path, cv2.IMREAD_COLOR) # デフォルトカラー読み込み
        img_npy = cv2.cvtColor(img_npy, cv2.COLOR_BGR2RGB) # RGB変換
        img_npy = cv2.resize(img_npy, (64, 64)) # リサイズ64x64

        # plt.imshow(img_npy)
        # plt.show()
        # break
        img_npy = img_npy.flatten() # 一次元化
        npy_image_list.append(img_npy/255) # 0~1に正規化

    return npy_image_list


# kmeansのモデル構築
def build_kmeans(df, cluster_num):
    kmeans = KMeans(n_clusters=cluster_num, random_state=2021)
    kmeans.fit(df)

    return kmeans


# 主成分分析のモデル構築
def build_pca(df):
    pca = PCA()
    pca.fit(df)

    return pca
    

# 主成分分析の累積寄与率を可視化(この結果をもとに特徴ベクトルを決める)
def plot_contribution_rate(pca):
    fig = plt.figure()
    plt.gca().get_xaxis().set_major_locator(ticker.MaxNLocator(integer=True))
    plt.plot([0] + list(np.cumsum(pca.explained_variance_ratio_)), "-o")
    plt.xlabel("Number of principal components")
    plt.ylabel("Cumulative contribution rate")
    plt.grid()
    plt.show()
    # plt.savefig('../figure/pca_contribution_rate.png') 


# 主成分分析の第一主成分と第二主成分で散布図による可視化
def plot_scatter2d(df):
    fig = plt.figure()
    sns.scatterplot(data=df, x='PC1', y='PC2', hue='label', palette='bright', legend='full')
    plt.show()
    # plt.savefig('../figure/pca_scatter2d.png')


# 主成分分析の第一主成分と第二主成分と第三主成分で散布図による可視化
def plot_scatter3d(df):
    # https://qiita.com/maskot1977/items/082557fcda78c4cdb41f
    fig = plt.figure()
    ax = Axes3D(fig)

    #軸にラベルを付けたいときは書く
    ax.set_xlabel("PC1")
    ax.set_ylabel("PC2")
    ax.set_zlabel("PC3")

    #.plotで描画
    for label in df['label'].values:
        ax.plot(df.loc[df['label']==label, 'PC1'], 
                df.loc[df['label']==label, 'PC2'], 
                df.loc[df['label']==label, 'PC3'], 
                alpha=0.8, marker=".", linestyle='None')
    plt.show()
    # plt.savefig('../figure/pca_scatter3d.png')


# 結果をクラスタごとにディレクトリに保存
def make_cluster_dir(load_path, save_path, kmeans):
    # 保存先のディレクトリを空にして作成
    shutil.rmtree(save_path)
    os.mkdir(save_path)

    # クラスタごとのディレクトリ作成
    for i in range(kmeans.n_clusters):
        cluster_dir = save_path + "cluster{}".format(i)
        if os.path.exists(cluster_dir):
            shutil.rmtree(cluster_dir)
        os.makedirs(cluster_dir)

    # 各ディレクトリにコピー保存
    image_path_list = gb.glob(load_path)
    for label, path in zip(kmeans.labels_, image_path_list):
        shutil.copyfile(path, save_path + 'cluster{}/{}'.format(label, os.path.basename(path)))

    print('クラスタごとにファイル作成完了')


def main():
    LOAD_PATH = 'D:/Illust/Paimon/interim/face_only/*'            # 画像データの読込先
    SAVE_PATH = 'D:/Illust/Paimon/interim/face_only_clustering/'  # 画像データをクラスタリングした結果の保存先
    CSV_PATH = 'D:/Illust/Paimon/interim/face_only_pca.csv'       # 画像データを主成分分析した結果の保存先
    
    try:
        # すでに画像データを主成分分析した結果のCSVファイルがあれば読み込む、なければexceptへ
        pca_df = pd.read_csv(CSV_PATH)
    except FileNotFoundError:
        # 画像読み込み
        npy_image_list = load_image(LOAD_PATH)
        df = pd.DataFrame(npy_image_list)
        print(df.shape)
        
        # 主成分分析の実行
        pca = build_pca(df)
        pca_df = pd.DataFrame(pca.transform(df), columns=["PC{}".format(x + 1) for x in range(len(df))])
        plot_contribution_rate(pca)          # 累積寄与率可視化
        pca_df.to_csv(CSV_PATH, index=False) # 保存
    
    # kmeansによるクラスタリング
    train_df = pca_df.iloc[:, :1200]               # 学習データ
    cluster_num = int(input('cluster_num >'))      # クラスタ数を入力
    kmeans = build_kmeans(train_df, cluster_num)   # kmeansモデル構築
    make_cluster_dir(LOAD_PATH, SAVE_PATH, kmeans) # クラスタリング結果からディレクトリ作成

    # 可視化
    pca_df['label'] = kmeans.labels_ 
    plot_scatter2d(pca_df)              # 二次元散布図
    plot_scatter3d(pca_df)              # 三次元散布図

    
if __name__ == "__main__":
    main()

補足説明

全体の流れ

  1. 画像読み込み(一次元化済み)

  2. 一次元化した画像を主成分分析にかけて、kmeans用の特徴ベクトルに

  3. kmeansを実行してクラスタリング

  4. クラスタリングした結果に基づいて、ディレクトリを作成し、各クラスタディレクトリ名で画像をコピー保存

  5. 最後に クラスタリング結果を可視化

mainの例外処理

main関数の例外処理はCSV_PATHにあらかじめファイルを用意しない限り、原則except側に入ります。
画像を1枚1枚読み込んで主成分分析するのには、そこそこ時間がかかります。
何度もこの処理をするのは効率が悪いため、主成分分析した結果を一度実行したら保存しているわけです。二回目以降に実行時には保存した結果を読みに行きます。
主成分分析した結果がクラスタリングするときの特徴ベクトルになります。

try:
        # すでに画像データを主成分分析した結果のCSVファイルがあれば読み込む、なければexceptへ
        pca_df = pd.read_csv(CSV_PATH)
    except FileNotFoundError:
        # 画像読み込み
        npy_image_list = load_image(LOAD_PATH)
        df = pd.DataFrame(npy_image_list)
        print(df.shape)

        # 主成分分析の実行
        pca = build_pca(df)
        pca_df = pd.DataFrame(pca.transform(df), columns=["PC{}".format(x + 1) for x in range(len(df))])
        plot_contribution_rate(pca)          # 累積寄与率可視化
        pca_df.to_csv(CSV_PATH, index=False) # 保存

train_dfの1200

今回の場合、主成分分析した結果は約5000x5000の行列です。行が各主成分(つまり第5000主成分まである)、列が画像の枚数です。 主成分分析は次元圧縮のアルゴリズムで、5000ある特徴ベクトルを減らして(必要な情報だけに圧縮して)クラスタリング精度向上のために用います。 つまり、下記のコードで意味する1200は5000を1200まで減らしたということになります。

train_df = pca_df.iloc[:, :1200]               # 学習データ

1200にした理由は、主成分分析の累積寄与率というものを用います。plot_contribution_rate(pca)関数はそれを可視化したものになります。
縦軸が1に近づけば近づくほど、その主成分だけでデータを説明できていることになります。
1200で限りなく1に近いと判断し、今回は1200にしました。600でも、1800でも試す価値はあると思います(今回はやっていません)。
なお、あらかじめ指定した累積寄与率を閾値として設定しておき、その値を超えたときの数を与えてやることもできます。

f:id:Noleff:20210221230731p:plain

kmeansのクラスタ

クラスタの数は標準入力で指定します。これをいくつにするかは正直適当です。
今回は50にしました。

cluster_num = int(input('cluster_num >'))      # クラスタ数を入力
kmeans = build_kmeans(train_df, cluster_num)   # kmeansモデル構築

参考までにクラスタリングした可視化結果を載せます。クラスタ数が50だと見えにくいので、クラスタ数を10で可視化しています。

  • 二次元散布図

f:id:Noleff:20210221231657p:plain

  • 三次元散布図

f:id:Noleff:20210221231706p:plain

結果

ディレクトリ全体

f:id:Noleff:20210221235048p:plain

ディレクト

うまくいった例

2番の問題を解決

このクラスタ数のディレクトリは特に触る必要がなくなります。

f:id:Noleff:20210221233359p:plain
パイモンしかいない

左上4番目にだけ凝光様がぽつんといます。白髪のため紛れ込んでしまったのでしょう。

f:id:Noleff:20210221233749p:plain
凝光様がいる

今度は金髪のキャラクターが集まっているクラスタがありました。主人公がいっぱいいますね。一部クレーや凝光様もいるのが見てとれると思います(原神わからない方ごめんなさい)。
なお、今回は髪色で同じクラスタになっているケースは他にも多くありました。

f:id:Noleff:20210221234019p:plain
金髪キャラ

3番の問題を解決

ほぼ同じ画像のパイモンばかりしかいません。お好みで画像を削除できます。

f:id:Noleff:20210221233529p:plain
ほぼ同じパイモン1
f:id:Noleff:20210221235426p:plain
ほぼ同じパイモン2

うまくいっていない例

しいて共通点を挙げるなら全体的に暗い画像となっている点です。しかし、明るい画像も一部含まれてもいます。
このような場合、髪色のクラスタにグルーピングされず、結果として共通点のないクラスタになっていましまいた。
今回うまくいっていないクラスタの大半がこのような暗めの画像でした。RGBだけでなくHSVの特徴ベクトルも必要かもしれませんね。
パイモンがいなければディレクトリごと除去できますが、実は真ん中あたりにしれっといます。

f:id:Noleff:20210221234423p:plain
暗い画像??

さて、これは本当に共通点がないクラスタです。カラー画像や白黒画像が紛れ込んでます。
同じ画像は同じクラスタには入っているようですが、それくらいしか共通点が見つかりませんでした。

f:id:Noleff:20210221234948p:plain
共通点不明

まとめ

今回はpcaとkmeansを用いて画像のグルーピングをしました。結果はそこそこ実用的ではありますが、まだまだ精度向上はできそうな気がします。
今後の課題としては

  • 1番の問題を解決する
  • 画像の明るさを踏まえた特徴ベクトルを作成する

などがありますね。
あとは別の次元圧縮アルゴリズムクラスタリングアルゴリズムを使ってみると良いかもしれません。正直pcaとkmeansは王道ワンパターンなので(ただ王道を馬鹿にできないのもまた事実)。
なお、最適なクラスタ数については考えません。kmeansでは正直不毛です。
クラスタリングアルゴリズムではxmeans使うか階層型クラスタリングを使うか。
次元圧縮アルゴリズムではt-SNEを使うか。
といったところでしょうか。

参考文献

【機械学習】ブログのサムネ画像をクラスタリングしてみる! - 株式会社ライトコード

主成分分析を Python で理解する - Qiita

Python: seaborn を使った可視化を試してみる - CUBE SUGAR CONTAINER

matplotlib 3次元散布図 | Python学習講座

【Matplotlib】3D散布図の作成 | 西住工房

[Python入門]shutilモジュールによる高水準ファイル操作:Python入門(1/3 ページ) - @IT

Pythonでファイル・ディレクトリを削除するos.remove, shutil.rmtreeなど | note.nkmk.me

カテゴリ変数(質的データ)の前処理の違いまとめ

はじめに

一般的に機械学習においてカテゴリ変数は、前処理として数値化する必要があります。
本記事ではその前処理の方法と違いについてまとめです。

データの種類と意味

下図のように変数は4つの尺度に分けられます。
今回説明するのは名義尺度と順序尺度に関する前処理の方法です。

f:id:Noleff:20210102193614p:plain

引用:人工知能プログラミングのための数学がわかる本

データ

まず、簡単なデータを以下のコードで準備します。あるユーザがある飲み物を買ったときの値段、サイズ、飲み物の種類が記載されています。
このデータでは、sizeデータはL>M>Sの関係にあるため順序尺度、drinkデータは大小関係がないため名義尺度になります。

import numpy as np
import pandas as pd

df = pd.DataFrame([
        ['Aさん', 100, 'S', 'cola'],
        ['Bさん', 150, 'M', 'tea'],
        ['Cさん', 200, 'L', 'coffee'],
        ['Dさん', 100, 'S', 'tea'],
        ['Eさん', 200, 'L', 'coffee'],
        ['Fさん', 200, 'L', 'tea'],
        ['Gさん', 150, 'M', 'tea'],
        ['Hさん', 200, 'L', 'coffee'],
        ['Iさん', 100, 'S', 'cola'],
        ['Jさん', 200, 'L', 'tea']],
        columns=['user', 'price', 'size', 'drink'])
df

f:id:Noleff:20210102203642p:plain

前処理の種類と違い

一般に2つのアプローチがあります。
1つ目は順番にラベリングしていく方法です。colaを0、teaを1、coffeeを2のようにラベルをつけます。
2つ目はダミー変数化(OneHotエンコーディング)する方法です。cola、tea、coffeeのカラムを新しく追加し、任意のカラムを1、それ以外のカラムを0にします。

1つ目の方法は決定木ベースのモデルには効果はありますが、線形モデルやNNに使う場合は注意が必要です。「teaはcoffeeと2倍の関係にある」といった解釈をされかねません。
2つ目の方法は1つ目より使われる気がします。ただし、カテゴリ変数が1000個あるといったふうに、量が多ければ多いほどカラムも増えます。膨大なカテゴリ変数は適宜集約するといった処理が必要になるかと思います。

名義尺度

ラベリング

factorize

pandasのfactorizeを使えば、簡単にカテゴリ変数をラベリングしてくれます。

ft_array, ft_index = pd.factorize(df['drink']) # tupple型で返却される
df_ft = pd.DataFrame(ft_array, columns=['ft_drink'])
df_factrize = pd.concat([df, df_ft], axis=1) # 元のデータフレームと連結
df_factrize

f:id:Noleff:20210102205057p:plain

LabelEncoder

sklearnのLabelEncoderを使えば、簡単にカテゴリ変数をラベリングしてくれます。

from sklearn.preprocessing import LabelEncoder

lenc = LabelEncoder()
lenc.fit(df['drink'])
lenc_vec = lenc.transform(df['drink'])
df_le = pd.DataFrame(lenc_vec, columns=['le_drink'])
df_lenc = pd.concat([df, df_le], axis=1) # 元のデータフレームと連結
df_lenc

f:id:Noleff:20210102205138p:plain

factorizeとLabelEncoderの違い

以下のように学習データとテストデータに分けられたデータがあるとします。
Nさんが学習データにはないサイズと飲み物の種類(LLとcider)が記載されていることに注意してください。

train = df.copy()
test = pd.DataFrame([
        ['Kさん', 200, 'L', 'cola'],
        ['Lさん', 100, 'S', 'tea'],
        ['Mさん', 150, 'M', 'coffee'],
        ['Nさん', 250, 'LL', 'cider']],
        columns=['user', 'price', 'size', 'drink'])

factorizeは値の出現順にラベルが振られるため、データフレームが別々にある場合、別のラベルが振らてしまう可能性があります。

train_ft, idx = pd.factorize(train['drink'])
test_ft, idx = pd.factorize(test['drink'])

train_df = pd.concat([train, pd.DataFrame(train_ft, columns=['ft_drink'])], axis=1)
test_df = pd.concat([test, pd.DataFrame(test_ft, columns=['ft_drink'])], axis=1)

display(train_df)
display(test_df)

f:id:Noleff:20210102213152p:plain

一度trainデータとtestデータを連結させ、1つのデータフレームとしてラベルを振れば回避することができます。

all_df = pd.concat([train, test], axis=0).reset_index(drop=True)
all_ft, idx = pd.factorize(all_df['drink'])
train_test_df = pd.concat([all_df, pd.DataFrame(all_ft, columns=['ft_drink'])], axis=1)

しかし、kaggleのようにあらかじめ学習データとテストデータがわけられている場合、わざわざ連結するのは面倒です。LabelEncoderを使えば連結せずに済みます。
なお、ラベリングの順番はfactorizeは値の出現順ですが、LabelEncoderはアルファベット順です。名義尺度のデータでは、ここはあまり気にする必要はないかと思います。

lenc = LabelEncoder()
train_df = lenc.fit(['cola', 'tea', 'coffee', 'cider']) # ここでカテゴリ変数の種類を指定
train_lenc = lenc.transform(train[['drink']])
test_lenc = lenc.transform(test[['drink']])

train_df = pd.concat([train, pd.DataFrame(train_lenc, columns=['le_drink'])], axis=1)
test_df = pd.concat([test, pd.DataFrame(test_lenc, columns=['le_drink'])], axis=1)

display(train_df)
display(test_df)

f:id:Noleff:20210102214232p:plain

ダミー変数化

get_dummies

pandasのget_dummiesを使えば、簡単にカテゴリ変数をダミー変数化してくれます。

df_gd = pd.get_dummies(df['drink'])
df_dummy = pd.concat([df, df_gd], axis=1) # 元のデータフレームと連結
df_dummy

f:id:Noleff:20210102214359p:plain

OneHotEncoder

sklearnのOneHotEncoderを使えば、同様にダミー変数化してくれます。

from sklearn.preprocessing import OneHotEncoder

oenc =  OneHotEncoder(sparse=False, dtype=int)
oenc.fit(df[['drink']]) # pandas.core.frame.DataFrame型もしくは二次元のnumpy.ndarray型が引数でないとエラー
oenc_vec = oenc.transform(df[['drink']]) # numpy.ndarray型で返却される

df_oenc = pd.DataFrame(oenc_vec, columns=['coffee', 'cola', 'tea'])
df_oht =  pd.concat([df, df_oenc], axis=1) # 元のデータフレームと連結
df_oht

f:id:Noleff:20210102214533p:plain

LabelBinarizer

sklearnのLabelBinarizerを使えば、同様にダミー変数化してくれます。

from sklearn.preprocessing import LabelBinarizer

lbnr = LabelBinarizer()
lbnr.fit(df[['drink']])
df_lbnr =  pd.concat([df, pd.DataFrame(lbnr.transform(df[['drink']]), columns=['coffee', 'cola', 'tea'])], axis=1) # OneHotEncoderとほぼ同じコードなためワンライナーで記述
df_lbnr

f:id:Noleff:20210102214533p:plain

get_dummiesとOneHotEncoder、LabelBinarizerの違い

factorizeとLabelEncoderの違い同様、 学習データとテストデータにデータがわけられているとします。

get_dummisを用いた場合、学習データとテストデータで作成されるカラムの数が異なってしまいます。factorizeとLabelEncoderの違いと同じように学習データとテストデータを連結させる方法もありますが、OneHotEncoderとLabelBinarizerを使えば回避できます。
以下のコードではLabelBinarizerの例です。

lbnr = LabelBinarizer()
lbnr.fit(train['drink'])
display(pd.DataFrame(lbnr.transform(train[['drink']]), columns=['coffee', 'cola', 'tea']))
display(pd.DataFrame(lbnr.transform(test[['drink']]), columns=['coffee', 'cola', 'tea']))

f:id:Noleff:20210102215708p:plain

ただし、データの中にnanもしくはinfが含まれている場合、get_dummisはエラーが出ませんが、OneHotEncoderとLabelBinarizerはエラーが出ます。

oenc = OneHotEncoder(sparse=False)
df['nan_and_inf'] = ['A', 'A', 'A', 'A', 'A', 'A', np.nan, np.inf, 'B', 'B']
df

f:id:Noleff:20210102225157p:plain

# エラーが出る
oenc.fit_transform(df[['nan_and_inf']]) 
# エラーが出ない
pd.get_dummies(df['nan_and_inf'])

f:id:Noleff:20210102225227p:plain

nanはスルーされますが、infはされないことに注意です。

OneHotEncoderとLabelBinarizerの違い

OneHotEncoderとLabelBinarizerの違いは複数のカラムをまとめてダミー変数化できるかどうかです。順序尺度データではありますが、sizeデータもまとめてダミー変数化してみます。

# エラーが出る
pd.DataFrame(lbnr.fit_transform(train[['drink', 'size']]), columns=['coffee', 'cola', 'tea', 'L', 'M', 'S'])
# エラーが出ない
pd.DataFrame(oenc.fit_transform(train[['drink', 'size']]), columns=['coffee', 'cola', 'tea', 'L', 'M', 'S'])

f:id:Noleff:20210102220004p:plain

順序尺度

名義尺度とは異なりラベリングが一般的かと思います。しかし、factorizeやLabelEncoderでは任意の順番にラベリングできません。また、色々調べてみましたが、pandas、sckit-learnともに任意の順番にラベリングできるメソッドはなさそうです。

そのため自分で関数なり作る必要があります。以下に例を示します。 

apply

lambda式を使えば一行で書けます。

df['ordinal_size'] = df['size'].apply(lambda x: ['S', 'M', 'L', 'LL'].index(x))

map

mapでもlambda式を使えば一行で書けます。ただmapの場合は辞書型もいける口です。ここは個人の好みでしょう。

# df['ordinal_size'] = df['size'].map(lambda x: ['S', 'M', 'L', 'LL'].index(x)) # lambdaもいける
df['ordinal_size'] = df['size'].map({'S': 0, 'M': 1, 'L': 2, 'LL':3})

また、自分でカテゴリ変数のリストを作るのが面倒な場合(もしくは多すぎる)は、学習データとテストデータ含めたすべてのデータからユニークな値を取るしかない思います。

unique_size_list = list(all_df['size'].unique())

train['ordinal_size'] = train['size'].apply(lambda x: unique_size_list.index(x))
test['ordinal_size'] = test['size'].apply(lambda x: unique_size_list.index(x))

display(train)
display(test)

f:id:Noleff:20210103012832p:plain

factorize

一応factorizeとLabelEncoderの実行結果も示します。

sizeの出現順がS、M、Lだったためたまたまうまくいっています。順番があらかじめラベリングしたい順序にソートされていればfactorizeは使えなくもないです。

ft_array, ft_index = pd.factorize(df['size']) # tupple型で返却される
df_ft = pd.DataFrame(ft_array, columns=['ft_size'])
df_factrize = pd.concat([df, df_ft], axis=1) # 元のデータフレームと連結
df_factrize

f:id:Noleff:20210102203834p:plain

LabelEncoder

アルファベット順=ラベリングしたい順序であれば使えます。sizeの例では無理です。

from sklearn.preprocessing import LabelEncoder

lenc = LabelEncoder()
lenc.fit(df['size'])
lenc_vec = lenc.transform(df['size'])
df_le = pd.DataFrame(lenc_vec, columns=['le_size'])
df_lenc = pd.concat([df, df_le], axis=1) # 元のデータフレームと連結
df_lenc

f:id:Noleff:20210102204314p:plain

まとめ

まとめると以下の表になります。

  • ○:この手法がベター、もしくはこの手法しかない
  • △:この手法よりも良い手法がある
  • ✕:この手法は一般的ではない
名義尺度 順序尺度 感想
get_dummis 学習データとテストデータのカテゴリ変数に注意
OneHotEncoder nanとinfに注意
LabelBinarizer nanとinfに注意
factorize 楽だが中途半端
LabelEncoder 楽だが中途半端
オリジナル関数 オレオレ関数は最強

表から分かる通り、名義尺度はOneHotEncoderもしくはLabelBinarizerを使う、順序尺度はオリジナル関数を定義する、というのが結論です。

参考文献

人工知能プログラミングのための数学がわかる本 | 石川 聡彦 |本 | 通販 | Amazon

Pythonでのカテゴリ変数(名義尺度・順序尺度)のエンコード(数値化)方法 ~順序のマッピング、LabelEncoderとOne Hot Encoder~ - Qiita

One-HotエンコーディングならPandasのget_dummiesを使おう | Shikoan's ML Blog

【python】機械学習でpandas.get_dummiesを使ってはいけない - 静かなる名辞

【python】sklearnでのカテゴリ変数の取り扱いまとめ LabelEncoder, OneHotEncoderなど - 静かなる名辞

pandasでカテゴリー変数を数値に変換する | 分析ノート

モデリングのための特徴量の前処理について整理した - For Your ISHIO Blog

python - LabelEncoderパンダdfの適合順序

機械学習のおおまかな流れを理解する

この記事は SLP KBIT Advent Calendar 2020 12日目の記事です。(じゃあなぜ 13日に公開しているのか) adventar.org

はじめに

機械学習を勉強する上で、重要なことはいくつもあると思います。
その中でも本記事は全体像の部分、機械学習のおおまかな流れに焦点を当てます。

機械学習を勉強し続けてはや2年。機械学習の分野は学ぶことが多く、勉強したてのときは気にも止めていなかったおおまかな流れですが、これを理解しているのとしていないとのでは今後の理解度に大きく響くと感じています。

ですので、本記事では〇〇モデルの具体的なアルゴリズムや数学的内容(理論)だったり、精度向上の前処理だったりは省いてますのでご了承を。

機械学習勉強したての後輩の誰かの目に入れば幸いです。

機械学習のおおまかな流れ

下図が機械学習のおおまかな流れです。大きく二つにわかれます。

一つ目は学習です。すでに手元にあるデータを学習データとして入力し、モデルを出力します。
モデルとはデータからある目的を解くために機械が分かる形に数値化(数式化)したものです。書いてて思いますが、非常にわかりにくい。
ある目的には大きく分類と回帰があります。分類は例えば犬と猫の画像を分類する。回帰は株価を予測するなどがあります。

二つ目は推論です。モデルと新しく入手したデータを現在データとしてマッチングさせることで、結果を出力します。
推論と書いてますが、これは分類や回帰など目的によって変わるため、一般的に推論と呼ばれています。そのため出力結果は二値分類なら0or1、回帰なら予測値になります。

長くなりましたが、この学習推論がごっちゃにならないようにしようぜ、というのが言いたいことです。

f:id:Noleff:20201213014755p:plain

線形回帰

機械学習を学ぶとき線形回帰から触れることが多いんじゃないでしょうか。数学的な理論を理解するには大学レベルの数学はいりますが、直感的なイメージは中学生でも理解できるからです。
Excelの散布図で近似曲線を引いたことがある人はすぐに理解できると思います。 そうです、y=ax+bです。

例えば人の身長を予測したいとします。身長が高ければ、体重も重くなるでしょう。これは正の相関があると言えます。正の相関があるということは体重からある程度身長がいくつか予測できることはなんとなくわかるかと思います。

ここで、身長(予測対象の変数)を一般的に目的変数と呼びます。対して体重(予測するために使う変数)を説明変数と言います。前者は一つ、後者は一つ以上と考えてください。
説明変数は今は体重だけです。人には痩せ型の人もいれば、太っている人もいるでしょう。そのため体重からだけで身長を正確に予測することは不可能に近いです。

では別に年齢、性別、手足の大きさ、腕の長さ、足の長さなどなど、人に関する情報がより詳しくデータとしてあるとしましょう。体重だけよりも正確に予測できることは想像できると思います。

線形回帰の話はこれくらいにしますが、より詳しく知りたい人は参考文献を参照してみてください。

データ

ボストン市郊外の地域別住宅価格のデータを使わさせていただきます。線形回帰を使う機械学習のデータセットとして一番有名だと思います。

プログラム

Python機械学習関連のライブラリが豊富ですが、ライブラリがあまりにできすぎていて、ほとんど数行足らずで書き終わります。そのため実装は簡単でも本質的に理解しにくく、非常に混乱してしまいます。

その結果、世の中にある機械学習のコードはほとんど、学習と推論を一括で書かれているケースが多いです。数行で書けるゆえの弊害ですかね。

とは言ったものの、世の中の実際のデータは機械学習のデータセットのようにきれい整形されていないため、地獄の前処理に苦しむわけですが、ここでは割愛。

学習と推論を強調するためにモデルをファイルとして一度保存してます。
ライブラリを使い倒してますが、とりあえずそこは一旦横に置いておいて。
図と照らし合わせながら読むと良いかと思います。おおまかな流れだけを意識してみてください。

  • X_train: 学習データ(説明変数)
  • y_train: 学習データ(目的変数)
  • X_test: 現在データ(説明変数)
  • y_test: 現在データ(目的変数)

gistb11b42b7ac823d54430f5cacf418e4d2

まとめ

機械学習のおおまかな流れについてまとめました。
繰り返しになりますが学習と推論です。可視化などしてますが、説明変数を決めるためにしているだけです。
線形回帰の精度評価は本記事では行いません。当たり前ですが、何の前処理もしていないので、精度は低いことは明らかです。
ボストン市郊外の地域別住宅価格を扱った記事は多くあると思いますので、そちらをご参照を。
キーワードとしては標準化、重回帰あたりを調べると精度が向上するかと思います。

参考文献

http://archive.ics.uci.edu/ml/index.php https://qiita.com/0NE_shoT_/items/08376b08783cd554b02e https://momonoki2017.blogspot.com/2018/01/scikit-learn_28.html https://localab.jp/blog/save-and-load-machine-learning-models-in-python-with-scikit-learn/

画像が見えにくいところもあるので、GitHubのリンクも載せておきます。

github.com

オートエンコーダによる時系列データの教師なし異常検知

はじめに

深層学習を用いた異常検知手法では有名なオートエンコーダを使ってプログラミングしたことをまとめます。オートエンコーダによる再構成誤差とLSTMによる予測誤差などとも比較予定です。

対象データ

今回もこちらの心電図のデータを使わせていただきました。
心電図のデータですので比較的周期性のあるデータですね。では実装していきます。

なお、データに関しては以前異常検知に関するブログでまとめてあるのでご参照ください。 noleff.hatenablog.com

プログラム

環境

  • OS: macOS Catalina
  • numpy: 1.18.5
  • pandas: 1.0.5
  • matplotlib: 3.3.0
  • scikit-learn: 0.23.1
  • tensorflow: 2.3.0
  • Keras: 2.4.3

関数

関数定義します。意味はコメントアウトを見てください。

# センサデータ取得
def read_df(path):
    df = pd.read_csv(path)
    return df

# Timeから時間の差を追加
def differ_time(df):
    # Timeから時間の差を割り当てる
    df['dif_sec'] = df['time'].diff().fillna(0)
    df['cum_sec'] = df['dif_sec'].cumsum()
    return df

# 正規化
def act_minxmax_scaler(df):
    scaler = MinMaxScaler()
    scaler.fit(df)
    # 正規化したデータを新規のデータフレームに
    # df_mc = pd.DataFrame(scaler.transform(df), columns=df.columns) 
    # 正規化したデータをnumpyリストに
    mc_list = scaler.fit_transform(df)
    return mc_list

# 部分時系列の作成
def split_part_recurrent_data(data_list, window):
    data_vec = []
    for i in range(len(data_list)-window+1):
        data_vec.append(data_list[i:i+window])
    
    return data_vec

# オートエンコーダの層を作成
def create_auto_encorder(in_out_shape):
    model = Sequential()
    # エンコード
    model.add(Dense(units=200, activation='relu', input_shape=(in_out_shape,)))
    model.add(Dense(units=100, activation='relu'))
    model.add(Dense(units=50, activation='relu'))
    # デコード
    model.add(Dense(units=100, activation='relu'))
    model.add(Dense(units=200, activation='relu'))
    # 出力層
    model.add(Dense(in_out_shape, activation='sigmoid'))

    # 作成したネットワークの確認
    model.summary()

    return model

# オートエンコーダ実行
def act_auto_encoder(model, train_vec, batch_size, epochs):
    # 学習条件の設定 誤差関数=平均二乗誤差、最適化手法=Adam法
    model.compile(loss='mse', optimizer='adam', metrics=['acc'])

    # 学習の実行
    hist = model.fit(x=train_vec, y=train_vec, batch_size=batch_size, verbose=1, epochs=epochs, validation_split=0.2)

    return hist

# 評価系のグラフをプロット
def plot_evaluation(eval_dict, col1, col2, xlabel, ylabel):
    plt.plot(eval_dict[col1], label=col1)
    plt.plot(eval_dict[col2], label=col2)
    plt.ylabel(ylabel)
    plt.xlabel(xlabel)
    plt.legend()
    plt.show()

以下ベタ書き

前処理

まずデータを読み込んで、学習データとテストデータにわけます。

# データ読み込み
df = read_df('./ecg.csv')
# 使うデータだけ抽出
test_df = df.iloc[:5000, :] 
train_df = df.iloc[5000:, :]

グラフにするとこんな感じ。

plt.figure(figsize=(16,4))
plt.plot(train_df['signal2'], label='Train')
plt.plot(test_df['signal2'], label='Test')
plt.legend()
plt.show()

f:id:Noleff:20201019182827p:plain

秒数を計算して、時系列情報を秒として持ちます。 差分をとったdif_secを加算して、cum_secカラムを新規追加します。

# 秒数を計算
test_df = differ_time(test_df)  
train_df = differ_time(train_df) 

display(test_df)
display(train_df)

f:id:Noleff:20201019182949p:plain

教師なしで異常検知する上で最終的に評価したいので正解のラベルを振ります。
今回は完全にグラフを目で見て、勝手に決めました。
そのため、これが実際に異常かどうかはわかりません。注意してください。
テストデータのグラフを15秒から20秒の間でクローズアップしました。
赤色に示すとおり、17秒から17.5秒の間を異常としました。

anomaly_df = test_df[(test_df['cum_sec']>15)&(test_df['cum_sec']<20)]
plt.figure(figsize=(16,4))
plt.plot(anomaly_df['cum_sec'], anomaly_df['signal2'], color='tab:orange', label='Test')

anomaly_df = test_df[(test_df['cum_sec']>17)&(test_df['cum_sec']<17.5)] 
plt.plot(anomaly_df['cum_sec'], anomaly_df['signal2'], color='tab:red', label='Test')
plt.xlabel('sec')

plt.show()

f:id:Noleff:20201019183505p:plain

赤色の部分にラベル1、それ以外をラベル0にして、新規カラムlabelを定義します。

test_df.loc[(test_df['cum_sec']>17)&(test_df['cum_sec']<17.5), 'label'] = 1
test_df['label'] = test_df['label'].fillna(0)

次に部分時系列を作ってきます。 異常検知には「外れ値検知」「変化点検知」「異常部位検知」の大きく3種類ががあります。 3種類の異常に関してはこちらを参照ください。

www.albert2005.co.jp

今回の心電図データは「異常部位検知」に該当します。
異常部位検知において、言葉の通り、ある時系列の部位を異常として検知したいです。そこで時系列を部分ごとに分割した部分時系列を作ります。
生データそのままだとうまく学習できないため、正規化してデータを0~1の範囲に値を収めます。

window = 250 # 部分窓

# 正規化
train_vec = act_minxmax_scaler(train_df[['signal2']])
test_vec = act_minxmax_scaler(test_df[['signal2']])

# 一次元リスト
train_vec = train_vec.flatten()
test_vec = test_vec.flatten()

# 部分時系列作成
train_vec = split_part_recurrent_data(train_vec, window)
test_vec = split_part_recurrent_data(test_vec, window)

学習

続いてメインのオートエンコーダです。学習するための層を構築していきます。
in_out_shapeはオートエンコーダの入力と出力のデータ数です。今回の場合はwindowと同じです。
オートエンコーダの中間層は適当に決めています。ここをチューニングすることこそが、ニューラルネットワークの醍醐味といえるでしょう。

in_out_shape = window
# オートエンコーダ作成
model = create_auto_encorder(in_out_shape)
Model: "sequential"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
dense (Dense)                (None, 200)               50200     
_________________________________________________________________
dense_1 (Dense)              (None, 100)               20100     
_________________________________________________________________
dense_2 (Dense)              (None, 50)                5050      
_________________________________________________________________
dense_3 (Dense)              (None, 100)               5100      
_________________________________________________________________
dense_4 (Dense)              (None, 200)               20200     
_________________________________________________________________
dense_5 (Dense)              (None, 250)               50250     
=================================================================
Total params: 150,900
Trainable params: 150,900
Non-trainable params: 0
_________________________________________________________________

学習を実行させていきます。オートエンコーダは学習データとして入力したデータを再現するように学習するニューラルネットワークです。
そのため、説明変数と目的変数は同じtrain_vecが入ります(act_auto_encoder関数の話)。

batch_size = 100
epochs = 20

train_vec = np.array(train_vec)
test_vec = np.array(test_vec)
hist = act_auto_encoder(model, train_vec, batch_size, epochs)
Epoch 1/20
318/318 [==============================] - 1s 3ms/step - loss: 0.0090 - acc: 0.0584 - val_loss: 0.0014 - val_acc: 0.1011
Epoch 2/20
318/318 [==============================] - 1s 3ms/step - loss: 0.0011 - acc: 0.1096 - val_loss: 7.4762e-04 - val_acc: 0.0962
Epoch 3/20
318/318 [==============================] - 1s 3ms/step - loss: 6.4285e-04 - acc: 0.1120 - val_loss: 5.4386e-04 - val_acc: 0.0996
Epoch 4/20
318/318 [==============================] - 1s 3ms/step - loss: 5.2786e-04 - acc: 0.1153 - val_loss: 4.9082e-04 - val_acc: 0.1028
Epoch 5/20
318/318 [==============================] - 1s 3ms/step - loss: 4.7851e-04 - acc: 0.1195 - val_loss: 4.2493e-04 - val_acc: 0.1082
Epoch 6/20
318/318 [==============================] - 1s 3ms/step - loss: 4.3090e-04 - acc: 0.1236 - val_loss: 3.7944e-04 - val_acc: 0.1099
Epoch 7/20
318/318 [==============================] - 1s 3ms/step - loss: 4.0277e-04 - acc: 0.1275 - val_loss: 3.7550e-04 - val_acc: 0.1146
Epoch 8/20
318/318 [==============================] - 1s 3ms/step - loss: 3.8248e-04 - acc: 0.1302 - val_loss: 3.4081e-04 - val_acc: 0.1204
Epoch 9/20
318/318 [==============================] - 1s 3ms/step - loss: 3.7038e-04 - acc: 0.1336 - val_loss: 3.3807e-04 - val_acc: 0.1229
Epoch 10/20
318/318 [==============================] - 1s 3ms/step - loss: 3.5362e-04 - acc: 0.1362 - val_loss: 3.4684e-04 - val_acc: 0.1269
Epoch 11/20
318/318 [==============================] - 1s 3ms/step - loss: 3.4186e-04 - acc: 0.1391 - val_loss: 3.1406e-04 - val_acc: 0.1294
Epoch 12/20
318/318 [==============================] - 1s 3ms/step - loss: 3.2798e-04 - acc: 0.1446 - val_loss: 3.1423e-04 - val_acc: 0.1332
Epoch 13/20
318/318 [==============================] - 1s 3ms/step - loss: 3.2826e-04 - acc: 0.1473 - val_loss: 2.9715e-04 - val_acc: 0.1377
Epoch 14/20
318/318 [==============================] - 1s 3ms/step - loss: 3.1464e-04 - acc: 0.1475 - val_loss: 3.0181e-04 - val_acc: 0.1439
Epoch 15/20
318/318 [==============================] - 1s 3ms/step - loss: 3.0473e-04 - acc: 0.1517 - val_loss: 2.9136e-04 - val_acc: 0.1446
Epoch 16/20
318/318 [==============================] - 1s 4ms/step - loss: 2.9575e-04 - acc: 0.1564 - val_loss: 2.9008e-04 - val_acc: 0.1465
Epoch 17/20
318/318 [==============================] - 1s 3ms/step - loss: 2.9276e-04 - acc: 0.1569 - val_loss: 2.8658e-04 - val_acc: 0.1494
Epoch 18/20
318/318 [==============================] - 1s 3ms/step - loss: 2.8440e-04 - acc: 0.1574 - val_loss: 2.7905e-04 - val_acc: 0.1521
Epoch 19/20
318/318 [==============================] - 1s 3ms/step - loss: 2.7920e-04 - acc: 0.1611 - val_loss: 2.6296e-04 - val_acc: 0.1561
Epoch 20/20
318/318 [==============================] - 1s 4ms/step - loss: 2.7432e-04 - acc: 0.1648 - val_loss: 2.5073e-04 - val_acc: 0.1508

戻り値histから、誤差値と評価値の収束具合を描画します。
val_lossよりもlossが大きければ未学習の可能性が高いです。

# 誤差の収束具合を描画
plot_evaluation(hist.history, 'loss', 'val_loss', 'epoch', 'loss')

f:id:Noleff:20201019184351p:plain

val_accよりもaccが大きければ過学習の可能性が高いです。

# 評価の収束具合を描画
plot_evaluation(hist.history, 'acc', 'val_acc', 'epoch', 'acc')

f:id:Noleff:20201019184404p:plain

異常検知

作成したモデルから予測してみます。
異常部位の予測が変になってますね。

# モデルをテストデータに適用
pred_vec = model.predict(test_vec)

# テストデータと出力データの可視化
plt.plot(test_vec[:,0], label='test')
plt.plot(pred_vec[:,0], alpha=0.5, label='pred')
plt.legend()
plt.show()

f:id:Noleff:20201019184456p:plain

テストデータとテストデータを再現したデータの再構成誤差を計算します。シンプルな距離の計算です。
グラフをみると、まあまあうまくできているように見えます。

# 再構成誤差(異常スコア)の計算
dist = test_vec[:,0] - pred_vec[:,0]
dist = pow(dist, 2)
dist = dist / np.max(dist)

# 異常スコアの可視化
plt.plot(dist)
plt.show()

f:id:Noleff:20201019184709p:plain

もとのデータと比べてみるとこんな感じです。

plt.plot(test_vec[3800:,0], label='test')
plt.plot(pred_vec[3800:,0], label='pred')
plt.plot(dist[3800:], label='dist')
plt.legend()
plt.ylabel('anomaly score')
plt.show()

f:id:Noleff:20201019190850p:plain

評価

実際に評価してみます。 まず、異常度からラベルが0か1かを決めます。異常度は0から1に収まりますので、試しに0.8以上を異常としてみます。

pred_label_list = [1 if d >= 0.8 else 0 for d in dist]

pred_df = pd.DataFrame(pred_label_list, columns=['pred_label'])

new_test_df = pd.concat([test_df, pred_df], axis=1)
new_test_df = new_test_df.fillna(0)

0と1の離散型のデータをグラフにしているため、カクカクしてます。
青色が初めに定義した異常のラベル。オレンジ色が異常度が0.8以上を異常としたラベルです。

plt.plot(new_test_df['label'], label='true anomaly')
plt.plot(new_test_df['pred_label'], label='pred anomaly')
plt.ylabel('anomary label')
plt.legend()

f:id:Noleff:20201019193921p:plain

分類タスクで評価に使われる混同行列から、精度、適合率、再現率を計算します。

print(confusion_matrix(new_test_df['label'], new_test_df['pred_label']))
print('accuracy_score: ', accuracy_score(new_test_df['label'], new_test_df['pred_label']))
print('precision_score: ', precision_score(new_test_df['label'],new_test_df['pred_label']))
print('recall_score: ', recall_score(new_test_df['label'],new_test_df['pred_label']))

果たしてこれはうまくいっていると言えるでしょうか。
結果から分かる通り、異常の誤検知は一回もしていません。それは適合率から明らかです。
ですが、再現率が非常に低いことから、多くの異常を見逃していることがわかります。

[[4876    0]
 [ 110   14]]
accuracy_score:  0.978
precision_score:  1.0
recall_score:  0.11290322580645161

ここで異常度からラベルを振る閾値をまじめに考えてみます。
まずは異常度をヒストグラムで描画してみます。
ほとんど0近辺の値しかないことがわかります。
つまり、閾値0.8は本当にかけ離れた異常度のみを異常としてラベルをふっていることになります。

f:id:Noleff:20201019190928p:plain

閾値を0.1 以上にしてみましょう。
閾値を下げたのでもちろん誤検知の数は増えます。しかし、同時に異常の見逃しの数を減らすことができました。
誤検知と見逃しはトレードオフの関係なので、どういう異常検知をしたいかで閾値を調整する必要があります。

pred_label_list = [1 if d >= 0.1 else 0 for d in dist]

pred_df = pd.DataFrame(pred_label_list, columns=['pred_label'])

new_test_df = pd.concat([test_df, pred_df], axis=1)
new_test_df = new_test_df.fillna(0)

plt.plot(new_test_df['label'], label='true anomaly')
plt.plot(new_test_df['pred_label'], label='pred anomaly')
plt.ylabel('anomary label')
plt.legend()
plt.savefig('figure/11_0.1以上の異常ラベル.png')

print(confusion_matrix(new_test_df['label'], new_test_df['pred_label']))
print('accuracy_score: ', accuracy_score(new_test_df['label'], new_test_df['pred_label']))
print('precision_score: ', precision_score(new_test_df['label'],new_test_df['pred_label']))
print('recall_score: ', recall_score(new_test_df['label'],new_test_df['pred_label']))
[[4852   24]
 [  46   78]]
accuracy_score:  0.986
precision_score:  0.7647058823529411
recall_score:  0.6290322580645161

f:id:Noleff:20201019191710p:plain

まとめ

オートエンコーダによる時系列データの教師なし異常検知をプログラミングしました。
時系列データは有名な心電図のデータです。
簡単にプログラムしたため、異常検知の精度はいまいちです。心電図の生データを正規化しただけの前処理しかしていないため、もう少し工夫をする必要があるかもしれません。
また、オートエンコーダの層の次元数に関しても色々試してみる価値はあると思います。

参考文献

https://www.matlabexpo.com/content/dam/mathworks/mathworks-dot-com/images/events/matlabexpo/jp/2017/b5-neural-net-sensor-data.pdf

https://www.renom.jp/ja/notebooks/tutorial/time_series/lstm-ae-ad/notebook.html

https://note.nkmk.me/python-tensorflow-keras-basics

https://engineer-fumi.hatenablog.com/entry/2018/09/18/161914

https://intellectual-curiosity.tokyo/2019/06/29/keras%E3%81%AB%E3%81%8A%E3%81%91%E3%82%8Btrain%E3%80%81validation%E3%80%81test%E3%81%AB%E3%81%A4%E3%81%84%E3%81%A6/

https://www.albert2005.co.jp/knowledge/machine_learning/anomaly_detection_basics/anomaly_detection_time

https://noleff.hatenablog.com/entry/2020/07/07/153450

LightGBMを使って気温予測してみた

はじめに

今回の記事では、Kagglerランカー達がこぞって使ってるという「LightGBM」なるものを使ってプログラミングしたことをまとめていきます。

LightGBMとは

LightGBMとは決定木アルゴリズムを応用したアルゴリズムです。よくXGBoostと比較されるのを目にします。

詳しくは以下のサイトをご参考ください。決定木、アンサンブル学習、勾配ブースティング(LightGBMやXGBoost)の順に、非常にわかりやすくまとめられています。

www.codexa.net

プログラミング

データ

実際に手を動かしていきます。まず、肝心のデータですが、今回は高松市の気象データを気象庁の方からダウンロードさせていただきました。
気象データの期間は2019-01-01から2020-08-29までです。

目標

今回の目標は「一時間単位に高松市の気温を予測する」こととします。そのためダウンロードした気象データは一日単位ではなく、一時間単位の気象データを収集しました。

プログラム

先に関数を定義しておきます。

# データ読み込み
def read_df(path):
    df = pd.read_csv(path, sep=',')
    return df

# データ抽出
def extract_train_test_data(df):
    train_df = df[df['year'] == 2019]
    test_df = df[(df['year'] == 2020)]

    return train_df, test_df

# 学習、テストデータ作成
def make_train_test_data(df, col_list):
    train_df, test_df = extract_train_test_data(df)
    
    # 学習用、テスト用にデータを分割
    # X_train = train_df[col_list]
    X_train = train_df.drop(['datetime', 'temperature'], axis=1)
    y_train = train_df['temperature']
    
    # X_test = test_df[col_list]
    X_test = test_df.drop(['datetime', 'temperature'], axis=1)
    y_test = test_df['temperature']

    return X_train, y_train, X_test, y_test

まず、データを読み込んで表示します。read_df関数の引数のファイル名は任意です。

# データ読み込み
df = read_df('20190101-20200829-takamatsu.csv')
df

気象庁のデータをそのまま使うと、pandasを使う上で不都合が多いので、以下のように整形しています。約14500レコードになりました。

f:id:Noleff:20200901185420p:plain

カラムそれぞれの意味は以下の表にまとめます。

カラム名 意味
datetime 日付
temperature 気温
pressure 気圧
relative_humidity 相対湿度
wind_speed 風速
rainfall 降水量
sea_lavel_pressure 海面気圧
day_length 日照時間
year
month
day

次に、気象データを学習用とテスト用に分割します。今回は2019の気象データを学習データ、2020の気象データをテストデータとしました。
col_listは目的変数である気温を予測するために使用する説明変数のカラムです。今回は使用していません。 特徴量選択や特徴量エンジニアリングについて触れず、気温と日付以外のすべてのカラムを説明変数(特徴量)とします。

# 学習、検証データ作成
col_list = ['relative_humidity', 'pressure', 'day_length', 'month', 'day', 'hour']
X_train, y_train, X_test, y_test = make_train_test_data(df, col_list)

LightGBM用にデータをセットし、学習データで回帰モデルを作成、テストデータで予測をします。
実際の2020-01-01から2020-08-29の気温データと、予測した結果を連結し、predicted_dfにまとめます。
評価基準はrmse(Root Mean Square Error)としています。

# LightGBM用のデータセット
lgb_train = lgb.Dataset(X_train, y_train)
lgb.test = lgb.Dataset(X_test, y_test)

# 評価基準 
params = {'metric' : 'rmse'}

# 回帰モデル作成
gbm = lgb.train(params, lgb_train)

# 予測
test_predicted = gbm.predict(X_test)
predicted_df = pd.concat([y_test.reset_index(drop=True), pd.Series(test_predicted)], axis = 1)
predicted_df.columns = ['true', 'pred']
predicted_df

f:id:Noleff:20200901185451p:plain

予測値と実測値の散布図を作成してみます。散布図と一緒に赤色で y = x の一次関数も描画しました。赤色にデータが近いほど良い結果が出ていると言えます。

# 予測値と実測値の散布図
RMSE = np.sqrt(mean_squared_error(predicted_df['true'], predicted_df['pred']))

plt.figure(figsize = (4,3))
plt.scatter(predicted_df['true'], predicted_df['pred'])
plt.plot(np.linspace(0, 40), np.linspace(0, 40), 'r-')

plt.xlabel('True Temperature')
plt.ylabel('Predicted Temperature')
plt.text(0.1, 30, 'RMSE = {}'.format(str(round(RMSE,3))), fontsize=10)

plt.show()

f:id:Noleff:20200901190031p:plain

今度は予測値と実測値をそれぞれ折れ線図として作成してみます。横軸は2020-01-01 01:00:00 からの経過時間(一時間単位)です。

# 予測値と実測値の折れ線図
plt.figure(figsize=(16,4))
plt.plot(predicted_df['true'], label='true')
plt.plot(predicted_df['pred'], label='pred')

plt.xlabel('time')
plt.ylabel('temperature')
# plt.xlim(0,1000)
plt.legend()

plt.show()

f:id:Noleff:20200901190341p:plain

大体うまくいけているのではないでしょうか。
もう少しクローズアップしてデータを見てみます。上記のプログラムにplt.xlimメソッドを使って横軸に制限をかけたときの折れ線図を示します。
0-1000の図で、横軸が500から600にかけてのところで予測値と実測値がずれていることがわかります。

  • 0-100

f:id:Noleff:20200901190643p:plain

  • 0-1000

f:id:Noleff:20200901190715p:plain

最後に有効な特徴量が何だったのか見てみます。
dayやmonthなどが上位に来ています。

# 特徴量の重要度
lgb.plot_importance(gbm, height = 0.5, figsize = (4,8))
plt.show()

f:id:Noleff:20200901192312p:plain

所感

思ったよりもうまくできたので正直驚きました。Lightと名前につくだけあってとても高速です。
ハイパーパラメータをチューニングしたり、特徴量を選択してないため実装も簡単です。
検討すれば精度はさらに向上すると思います。

参考文献

気象庁|過去の気象データ検索

Kaggler がよく使う「LightGBM」とは?【機械学習】 – 株式会社ライトコード

LightGBM 徹底入門 – LightGBMの使い方や仕組み、XGBoostとの違いについて

【機械学習実践】LightGBMで回帰してみる【Python】

【データサイエンティストへの道⑴】lightGBMの軽ーい実装とまとめ|すりふと/医学生|note

Appendix

今回のデータのペアプロットと相関係数を以下に添付します。
このように可視化することで、目的変数との相関を見比べれば、特徴量選択を決める一つの指標になると思います。

f:id:Noleff:20200907155908p:plain

  • 相関

f:id:Noleff:20200907155945p:plain

Pythonのクラス内におけるメソッドについてまとめてみた

目的

備忘録用です。普段Pythonでプログラミングするとき、めんどくさいのであまりクラスを使ってプログラミングしないことが多いです(jupyter notebookのベタ書き脳死コーディングのせい)。
ですが、最近はちょっとこのままでは良くないなと、あまり使う必要ないと思う場面でもわざわざクラスにしてPythonのコードを書いてます。

今回の記事は二番煎じ感満載かもしれませんが、ご了承ください。

メソッド

クラス内のメソッドには以下の3種類があります。

インスタンスメソッド

普通のメソッドがこれです。簡単に例題プログラムを書きます。

class Student:
    
    def __init__(self, name, s_id):
        self.name = name
        self.s_id = s_id
        
    def reply_name(self):
        print("私は" + self.name + "と言います")
        
    def reply_id(self):
        print("私の学籍番号は" + str(self.s_id) + "です")

if __name__ == '__main__':
    st = Student('山田太郎', '20A999')
    st.reply_name()
    st.reply_id()
山田太郎と言います
学籍番号は20A999です

学生を表すStudentクラスを作りました。結果は上記のようになります。
クラス内のインスタンスメソッドであるreply_nameとreply_idはどちらもpublicなメソッドになります。
privateなメソッドも追加してみましょう。

class Student:

    def __init__(self, name, s_id):
        self.name = name
        self.s_id = s_id
        # 成績は順番に「秀、優、良、可、不可」
        self.school_credit = {}

    def reply_name(self):
        print(self.name + "と言います")
        
    def reply_id(self):
        print("学籍番号は" + str(self.s_id) + "です")

    def reply_gpa(self):
        print("GPAは" + str(self.__calc_GPA()) + "です")

    def input_credit(self, credit_dict):
        self.school_credit = credit_dict

    def __calc_GPA(self):
        gpa = 0
        gpa_list = ['不可', '可', '良', '優', '秀']

        # 取得した単位はすべて1とする
        for sc in self.school_credit.values():
            if sc in gpa_list:
                gpa += gpa_list.index(sc)
               
        gpa /= len(self.school_credit)
        
        return gpa

if __name__ == '__main__':
    st = Student('山田太郎', '20A999')
    st.reply_name()
    st.reply_id()
    st.input_credit({'英語':'可', '実験':'優', 'プログラミング':'秀', 'DB':'良'})
    st.reply_gpa()
山田太郎と言います
学籍番号は20A999です
GPAは2.5です

GPAの計算式は雑ですがこんな感じで。計算式間違ってたらごめんなさい。
__calc_GPAメソッドがprivateなので、次のようにするとエラーを吐きます。

if __name__ == '__main__':
    st = Student('山田太郎', '20A999')
    st.reply_name()
    st.reply_id()
    st.input_credit({'英語':'可', '実験':'優', 'プログラミング':'秀', 'DB':'良'})
    print(st.__calc_GPA())
山田太郎と言います
学籍番号は20A999です
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
 in 
      4     st.reply_id()
      5     st.input_credit({'英語':'可', '実験':'優', 'プログラミング':'秀', 'DB':'良'})
----> 6     print(st.__calc_GPA())

AttributeError: 'Student' object has no attribute '__calc_GPA'

privateなので外からアクセスできません。
しかし、実は「インスタンス._クラス名__変数名」でアクセスできます。推奨されないので覚えなくてもいいですが……。ちなみに変数も同様です。

if __name__ == '__main__':
    st = Student('山田太郎', '20A999')
    st.reply_name()
    st.reply_id()
    st.input_credit({'英語':'可', '実験':'優', 'プログラミング':'秀', 'DB':'良'})
    print(st._Student__calc_GPA())
山田太郎と言います
学籍番号は20A999です
2.5

クラスメソッド

インスタンスメソッドとの違い

ここから個人的には本題です。
インスタンスメソッドとの違いを挙げるとすれば、まずはクラスから直接呼び出せることだと思います。メソッドの頭に@classmethodとデコレートし、慣例的に引数をclsにします。selfはインスタンスなのに対し、こちらはクラスを意味します。

class Student:

    def __init__(self, name, s_id):
        self.name = name
        self.s_id = s_id

    def reply_name(self):
        print(self.name + "と言います")
        
    @classmethod
    def reply_greeting(cls):
        print('こんにちは')

if __name__ == '__main__':
    st = Student('山田太郎', '20A999')
    st.reply_greeting()
    Student.reply_greeting()
こんにちは
こんにちは

クラスメソッドの使い所

これだけの例だと、いまいち利点や使い所がわかりません。
例えば外部ファイル(JSONなど)から学生の情報を取得するとします(普通はDBでしょうが)。今までと同様に書けば以下のようになります。

import pandas as pd

class Student:

    def __init__(self, name, s_id):
        self.name = name
        self.s_id = s_id

    def reply_name(self):
        print(self.name + "と言います")

if __name__ == '__main__':
    info_df = pd.read_json('student.json', encoding='UTF-8')
    st = Student(info_df['name'].values[0], info_df['id'].values[0])
    st.reply_name()

関数にするなら以下です。

import pandas as pd

class Student:

    def __init__(self, name, s_id):
        self.name = name
        self.s_id = s_id

    def reply_name(self):
        print(self.name + "と言います")

def get_student():
    info_df = pd.read_json('student.json', encoding='UTF-8')
    return Student(info_df['name'].values[0], info_df['id'].values[0])

if __name__ == '__main__':
    st = get_student()
    st.reply_name()

では、続いてクラスメソッドで書きます。

import pandas as pd

class Student:

    def __init__(self, name, s_id):
        self.name = name
        self.s_id = s_id

    def reply_name(self):
        print(self.name + "と言います")

    @classmethod
    def get_student(cls):
        info_df = pd.read_json('student.json', encoding='UTF-8')
        return cls(info_df['name'].values[0], info_df['id'].values[0])

if __name__ == '__main__':
    st = Student.get_student()
    st.reply_name()

違いはなんとなくわかるかもしれませんが、クラスメソッドの利点がわかりにくいので解説します。
最大のうまみはクラスの中にインスタンスを作るメソッドを書けることです。

関数にした場合とクラスメソッドにした場合を比較するとわかりやすいですが、get_studentをクラスメソッドではStudentクラス内でまとめて管理できます。
MainクラスからStudentクラスをimportするときなど、まとめられている方が使い勝手が良いです。
クラスに依存したメソッドはクラスメソッドで定義するほうが好ましいと思います。

また、実装方法によって一概には言えませんが、クラスメソッドを使わない場合、毎回

info_df = pd.read_json('student.json', encoding='UTF-8')

を書かなければなりません。非常に面倒です。

スタティックメソッド

最後にスタティックメソッドとクラスメソッドの違いをまとめます。
スタティックメソッドもクラスメソッド同様クラスからもインスタンスからも呼び出せます。
しかし、いくつか違いあり、1つ目は@staticmethodでデコレートすること。2つ目は引数になにも受け取らない実装が可能なことです。
インスタンスメソッドではself、クラスメソッドではclsを暗黙的に引数として書く必要がありますが、スタティックメソッドは書く必要はありません。そのため、クラスに依存しないメソッドと明示的に記述することができます。

class Student:

    def __init__(self, name, s_id):
        self.name = name
        self.s_id = s_id

    def reply_name(self):
        print(self.name + "と言います")

    @staticmethod
    def reply_greeting():
        print('こんにちは')

if __name__ == '__main__':
    Student.reply_greeting()
    st = Student('山田太郎', '20E999')
    st.reply_name()
こんにちは
山田太郎と言います

reply_greetingはこんにちはと出力するだけなため、まったくクラスに依存しません。このようなメソッドはクラスメソッドではなく、スタティックメソッドにすることが好ましいです。

しかし、クラスメソッドがあれば、スタティックメソッドがなくても実装上は問題ありません。
最初にクラスメソッドのプログラム例で書いたように、スタティックメソッドを使って書けるコードはクラスメソッドでも書けるからです(reply_greetingの話)。

Pythonにスタティックメソッドが必要なのかという記事もあるため、スタティックメソッドを使うのはある意味自己満足の世界な気もします。

参考文献

【Python】インスタンスメソッド、staticmethod、classmethodの違いと使い方 - Djangoの学習ができるチュートリアルサイトDjangoBrothers

[Python] クラスメソッド・スタティックメソッドの違い - Qiita

Pythonのクラスメソッド(@classmethod)とは?使いどころとメソッドとの違いを解説 - Python学習チャンネル by PyQ

Python クラスについて - Qiita

機械学習を使わずに時系列データの異常検知

はじめに

研究で時系列データの異常検知に関する研究を行っています。そのため、機械学習による異常検知のアルゴリズムについて普段から文献調査などしているわけですが、機械学習を使わずとも異常検知できるという記事を見つけました。

大変興味深かったので、記事の内容に少しだけアレンジを加えて実装してみたいと思います。

対象データ

自分で異常検知用のデータを作成しても良いのですが、それだとどうしても作業的になってしまい、面白くないのでこちらの心電図のデータを使わせていただきました。
ダウンロードしたデータをそのままグラフのさせると以下のようになります。ただし、そのまま出力させるとデータが多過ぎて潰れたグラフになって見にくいので、データの頭から32秒付近までの間のデータを抽出しています。

f:id:Noleff:20200707152820p:plain

異常検知

見つけたい異常

上のグラフからわかるように、17秒付近で他とは振動をしていることがわかると思います。今回はこの異常を検知したいと思います。

特徴量

1. 単純移動平均

一つ目の特徴量は移動平均です。移動平均とはある区間を決めてその区間をずらしていき、各区間内で平均値を求めていく計算になります。
平均は一般的にあまり有効な特徴量となり得ませんが、今回は二つ目の特徴量を際立たせるためにあえて平均を採用しています。(といっても今回の場合この特徴量を使っていないとうまくいってないところもありますが)
以下に移動平均のグラフを載せます。移動平均なため、データの頭は0で埋めていることに注意です。

f:id:Noleff:20200707143036p:plain

2. ローパスフィルタとの差分

フーリエ変換

ローパスフィルタを説明するために、まずフーリエ変換について簡単に触れます。

フーリエ変換とは実変数の複素または実数値関数を別の同種の関数に写す変換である。変換後の関数はもとの関数に含まれる周波数を記述し、しばしばもとの関数の周波数領域表現 (frequency domain representation) と呼ばれる。実質的に、フーリエ変換は関数を振動関数に分解する。(Wikipediaより) f:id:Noleff:20200707011937g:plain https://ja.wikipedia.org/wiki/%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B

言葉にするとわかりにくいですが、同じくWikipediaの載ってあるGIF画像を見ると非常にわかりやすいです。
複雑な振動波形でも、複数の振動波形が合成されてできています。この一つ一つの振動波形を取り出し、その振動波形からどの周期にどのぐらいの振れ幅を持っているかわかります。この振動波形から周波数成分に分解する変換のことをフーリエ変換といいます。
詳しくはこちらのサイトが大変わかりやすかったです。

ローパスフィルタ

続いてローパスフィルタの説明をします。フーリエ変換をすることで、振動波形を周波数成分に変換しました。今回の心電図をデータをフーリエ変換すると以下のようになります。ただし、0から10Hzの範囲のみ出力させています。なお、縦軸は各周期数成分の大きさを表す振幅スペクトルになります。

f:id:Noleff:20200707140545p:plain

大きな振幅スペクトルが二つ、中くらいの振幅スペクトルが一つ、小さな振幅スペクトルが三つできているのがわかると思います。
ローパスフィルタはこれらの周期数の内、Lowな振幅スペクトルの周期数のみ抽出します。つまり、振幅スペクトルが小さい高周波を除去しているイメージです。 ローパスフィルタ後の心電図とそれをフーリエ変換したしたものを以下に載せます。

f:id:Noleff:20200707140640p:plain

f:id:Noleff:20200707140709p:plain

今回ローパスフィルタは2.5Hzよりも大きな値を除去したので、二つ目の振幅スペクトルがやや鋭角になってますが、今回はこれでいきたいと思います。

差分の計算

長くなりましたが、ここまでが前置きです。
2つ目の特徴量はローパスフィルタとの差分でした。つまり、ローパスフィルタと元のデータの差分をとり、それを異常度として算出します。上のグラフでいうと、橙がローパスフィルタ後の心電図、青が元のデータの心電図なため、橙から青を引いた絶対値を異常度とします。なお、以下のグラフは0から1の範囲に値が収まるように正規化してます。

f:id:Noleff:20200707135208p:plain

思ったよりきれいに異常度を算出できませんでした(汗)。
ですが、心電図に異常が起きている17秒付近で異常度が大きくなっていることは見てわかると思います。

プログラム

環境

  • OS:windows10
  • python:3.7.6
  • numpy:1.18.1
  • pandas:1.0.1
  • scikit-learn:0.22.1

コード

上の描画のコードなど一部省略します。ほとんど参考から流用したので、気になる方はそちらをご覧ください。

生データ

最初に読み込むecg.csvはこんな感じです。timeが200スタートだったり、signal1とsignal2の違いは詳しく調べてないのでわかりません。今回はsignal2の方を使用しました。データはtimeからわかるように0.004秒ごとに記録されています。

f:id:Noleff:20200707141231p:plain

関数

# データ取得
def read_df(path):
    df = pd.read_csv(path)
    return df

# Timeから時間の差を追加
def differ_time(df):
    # Timeから時間の差を割り当てる
    df['dif_sec'] = df['time'].diff().fillna(0)
    df['cum_sec'] = df['dif_sec'].cumsum()
    return df

# 正規化
def minxmax_scaler(df):
    scaler = MinMaxScaler()
    scaler.fit(df)
    # 正規化したデータを新規のデータフレームに
    df_mc = pd.DataFrame(scaler.transform(df), columns=df.columns) 
    # 正規化したデータをリストに
    # mc_list = scaler.fit_transform(df)
    return df_mc

以下べた書き

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler

# データ読み込み
df = read_df('./ecg.csv')
dif_df = differ_time(df) # 秒数を計算 
test_df = dif_df.iloc[0:8192, :] # 使うデータだけ抽出

N = len(test_df) # FFTのサンプル数 
dt = 0.004 # サンプリング周波数 
freq = np.fft.fftfreq(N, d=dt) # 周波数

# 高速フーリエ変換(FFT)
F = np.fft.fft(test_df['signal2'].values) 
F_abs = np.abs(F/(N/2))

# ローパスフィルタ
threshold_period = 0.4
threshold_freq = 1 / threshold_period
F_lowpass = np.where(abs(freq) > threshold_freq, 0, F)
lowpass = np.fft.ifftn(F_lowpass).real

# 特徴量
sma = test_df['signal2'].rolling(25).mean().fillna(0) # 0.1秒ごとの移動平均
diff = abs(lowpass - test_df['signal2']) # ローパスフィルタとの差分

# 結果をまとめるデータフレーム
result_df = pd.DataFrame(diff, columns=['sma', 'diff', 'anomaly', 'label'])
result_df['sma'] = sma
result_df['diff'] = diff

 # 異常度を正規化
result_df = minxmax_scaler(result_df)

# 移動平均とローパスフィルタとの差分を乗算
result_df['anomaly'] = result_df['sma'] * result_df['diff'] 

# ラベル振り
result_df.loc[result_df['sma'] == 0, 'label'] = 2
result_df.loc[result_df['anomaly'] > 0.6, 'label'] = 1
result_df['label'] = result_df['label'].fillna(0)

これで、結果をまとめたresult_dfができました。

結果

移動平均とローパスフィルタを乗算した結果を示します。

plt.figure(figsize=(16,4))
plt.plot(test_df['cum_sec'], minxmax_scaler(test_df)['signal2'], color='blue', label='row')
plt.plot(test_df['cum_sec'], result_df['anomaly'], color='red', label='anomaly')
plt.legend()
plt.xlabel('sec')
plt.ylabel('anomaly')
plt.show()

f:id:Noleff:20200707144458p:plain

異常度が0から1の範囲に収まるようにしているので、閾値を設定すれば異常検知できます。今回は異常度が0.6より大きい箇所を異常としています。
グラフを見ると、心電図の17秒付近のみが0.6を越えていることがわかります。

今度は散布図として出力させます。

c_list = ['blue', 'red', 'green']

for i, c in enumerate(c_list):
    plt.scatter(result_df.loc[result_df['label']==i, 'diff'], result_df.loc[result_df['label']==i, 'sma'], color=c, label=i)

plt.xlabel('diff')
plt.ylabel('sma')
plt.legend()
plt.show()

f:id:Noleff:20200707144934p:plain

青色が正常、赤色が異常、緑色が移動平均時に0を埋めた値です。 異常なデータのみだけをグラフの右上付近に寄せることができていると思います。

参考

【データ分析入門】機械学習未使用!Pythonでゼロから始める振動解析|はやぶさの技術ノート

FFT を使った時系列データ解析 - nykergoto’s blog

<NumPy> 高速フーリエ変換による周波数解析 - Helve’s Python memo

https://cpp-learning.com/hampel-filter/

https://qiita.com/hoto17296/items/d337fe0215907432d754

https://ja.wikipedia.org/wiki/%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B

http://www.cs.ucr.edu/~eamonn/discords/