アプリとサービスのすすめ

アプリやIT系のサービスを中心に書いていきます。たまに副業やビジネス関係の情報なども気ままにつづります

クラスタリング(k-means)で画像から色の検出(機械学習、opencv)

今回はクラスタリング手法で、画像の重要な色を検出するタスクをやった。ニューラルネットワークならより正確な検出が可能だけど、データセット作成もろもろコストがでかい。なので、昔からあるクラスタリング手法で手軽に、かつ精度よく色を検出してみた。そのロジックと工程を書いてく。

目次
・rgb値から色を特定する方法
・用意したもの&ライブラリ
クラスタリングで色検出
・色同士のノルム計算で色の特定
・実験


全体コードはこちら

rgb値から色を特定する方法



色はrgb値の3つの値で表される。赤なら255:0:0だし、紫なら128:0:128。各3つの値のとる範囲は0〜255。rgb値はよく数学であるxyz空間と同じで、rgb空間上の点で表せる。

rgb空間は図で表すと下のような、6角形の形。


数学では三平方の定理とかでよく出てくるけど、「ベクトル間の距離(ノルム)が近いもの同士は似てる」という性質がある。


色でも同じ性質が使えて、rgb値のノルムが近いほど、互いの色が似ていることを表す。
この性質を使って、

・rgb値のノルムが近いもの同士=似た色

として色を検出する。



ただし、rgb空間は六角形なので、xyz空間のように3次元に変換してやる必要がある(座標変換)。rgb空間の3次元空間はlab空間というらしい。

lab値の各値のとりうる範囲は

L: [0, 100]
a: [-127, 127]
b: [-127, 127]

これでノルムを計算する。




色を特定する手順はこんな感じ

①色のlab値のリストの作成

②rgbで画像の読み込み

③rgb値のままクラスタリング

④skimage.color.rgb2lab()でrgbをlabに変換

⑤リストのlab値と、読み込んだ画像のlabのノルムを計算。リスト内の色と類似した色のtop4を特定



lab値の変換について
ライブラリのopencvでrgbをlabに変換すると、値は0〜255の範囲のまま。それだと困るので、ここではskimageというライブラリを使う。skimageならlabの値の範囲内で正確に変換してくれる。

参考サイト





用意したもの&ライブラリ


ライブラリ
opencv=> rgb値での画像読み込み用

・skimage => rgb値をlab値に変換用

・k-means => クラスタリング



用意したもの
・色のlab値のリスト

=> lab値の取得は下記のサイトを使った
https://syncer.jp/color-converter

・試す画像複数枚 =>上手く検出できるか試す画像






クラスタリングで色検出



まず、下の画像をopencvでrgbで読み込む。

def load_image(image_file):
    # cv2 load images as BGR
    image_bgr=[cv2.imread(image_file+'/'+i) for i in os.listdir(image_file)]
    image_rgb = [cv2.cvtColor(i, cv2.COLOR_BGR2RGB) for i in image_bgr]
    image_rgb = [cv2.resize(i, (150, 150)) for i in image_rgb]
    return image_rgb


img = load_image("/Users/d/clust_img")
img=np.reshape(img, (len(img),150,150,3))
hstack=np.hstack(img)
plt.imshow(hstack)
plt.show()

そのままこの画像をクラスタリング、して色の抽出。最後は、5色のrgb値を算出してる。
メインの関数はコード参照


img = load_image("/Users/tatsuyahagiwara/d/clust_img")
img=np.reshape(img, (len(img),150,150,3))
b,w,h,c=img.shape
N=7
sample_img=[i[int(w/10):int(N*w/10), int(h/10):int(N*h/10)] for i in img]
sample_img=[i.reshape(-1,3) for i in sample_img]
color=[]
for i in sample_img:
    clt = KMeans(n_clusters = 5)
    clt.fit(i)
    hist = centroid_histogram(clt)
    bar = plot_colors(hist, clt.cluster_centers_)
    plt.figure()
    plt.axis("off")
    plt.imshow(bar)
    plt.show()

k-meansはお手頃だし、割と正確に検出できたので使った。クラスタリングした結果の一部はこんな感じ。








色同士のノルム計算で色の特定



算出したrgb値をskimage.color.rgb2lab()でlab値に変換。
あとは関数でリストの色と読み込んだ色のlab値のノルムを計算して、一番近い(似ている)色を出した。

for i in sample_img:
    clt = KMeans(n_clusters = 5)
    clt.fit(i)
    hist = centroid_histogram(clt)
    bar = plot_colors(hist, clt.cluster_centers_)
    plt.figure()
    plt.axis("off")
    plt.imshow(bar)
    plt.show()
    bar_lab=rgb2lab(bar) # LAB値に変換
    c1, c2, c3, c4, c5=get_lab_from_list(bar_lab[0].tolist())
    dic1, dic2, dic3, dic4, dic5=lab_distance_dic(color_name, lab_list, c1, c2, c3, c4, c5)
    t1,t2,t3,t4,t5=get_topN_color(dic1,dic2,dic3,dic4,dic5)
    color.append([t1,t2,t3,t4,t5])

main_color=[]
for i in color:
    for n in i:
        main_color.append(n[0][0])
        
N=sorted(list(collections.Counter(main_color).values()),reverse=True)
for key, value in collections.Counter(main_color).items():
    if value in N[:4]:
        print('top_color:{}'.format(key))

色のtop4算出結果はこんな感じ。割と画像のメインの色を特定できてる感じなので、十分使えると思う。

ブラック系
ピンク系
ブルー系
グレー系



精度を上げるためにしたこと


メインのカラーが全部で13 。なので、色ごとに一つのlab値(青なら青一色のlab値)しか用意してない場合、単一のlab値でしか色を判断できない。
しかし、色一つとっても、その色に近い色はたくさんある。(例えば青でも、青っぽい色はたくさんある)
http://irononamae.web.fc2.com/colorlist/wa.html




なので、精度を上げる方法としては、

①増やしたい代表色(列の一番目)を決める

②代表色に似ている色とそのlab値を、列の2番目以降に入れる。




こうすれば、代表色をいろんなlab値で表せるので、精度は上がる。今回は青、イエロー、赤を増やしてみた結果、精度が上がった。

(ex:青っぽい色を全部「青」と区分すれば、「青っぽい色=全て青」となるので、いろんなlab値でも青と特定できる)。





実験



一応実験として、下の画像群でも色の特定をしてみた

やることは上の作業を別の画像でやるだけ


クラスタリング結果と検出した色名は以下の通り


ブラウン系
ピンク系
ブラック系
グレー系
イエロー系

赤とブラウンの区別が難しいので、主観でlab値を選んだので、検出結果はよくわからない。
総括では、割といい感じにできて、実用できるレベルだと思う




ガチのところは、ニューラルネットワークを使ってかなり精度高く色検出やってる(メルカリとか)。ただ、管理、データセット作成、計算もろもろコストがでかすぎるのが難点。その分、k-meansみたいな昔からある手法ても、工程さえしっかりやれば、簡単かつ高い精度でできるケースのが多い。アルゴリズムは手段にすぎないので使えさえすれば、k-meansでも十分精度の良いものができた。

AWSのGPU環境下でkerasを使った百万単位(ビックデータ)の画像分類の訓練、テスト、予測までの過程まとめ

今回はkerasを使って、AWSGPU環境下で5百万枚の画像を訓練してみた。ラベル数は200ラベル。おそらくビックデータと呼ばれる規模だと思う。エラーとか、障壁が多々あったので、備忘録もかねて工程を一通りまとめてく

目次
・EC2にGPU適用&jupyter環境構築
・tfrecordから画像の読み込み(input image)
・訓練(training)
・学習率
GPU環境下での訓練
・モデルの保存
・テスト・予測
・最後に

コードはこちら





EC2にGPU適用&jupyter環境構築



EC2インスタンスGPUの適用

Deep learning ubuntu(linux)で検索、インスタンス作成、作成時に「EPS-optimize( EPS最適化)」を選択

②あとはチュートリアル通りに最後までやれば、GPU環境を構築できる。
(https://docs.aws.amazon.com/ja_jp/AWSEC2/latest/UserGuide/install-nvidia-driver.html )





**補足**

チュートリアルの下のコマンドはlinuxパッケージのみ必要らしい。

sudo yum erase nvidia cuda


GPUの型番を確認コマンド ([ ]の中がGPUの型番に該当)

sudo update-pciids 
lspci | grep -i nvidia 




jupyter notebook 設定について

① 以下のサイトの$ipython以降を実行し、パスワードとか設定https://qiita.com/Salinger/items/c7b87d7000e48be6ebe2

② EC2インスタンスでjupyter用のIPアドレスの設定

EC2のインスタンス画面→対象のインスタンスのセキュリティグループを選択→インバウンドから編集→カスタムTCPルールのポート範囲/送信元を設定(8888/カスタム)

③必要なパッケージpipでインストール

④次のコマンドでjupyterの起動

$jupyter notebook & 

⑤jupyter上でGPU適用の確認 (「nvidia-smi」コマンドはGPU動作中のみ有効)

from tensorflow.python.client import device_lib
device_lib.list_local_devices()

上のコマンド打って、「name: "/device:GPU:0"〜」みたいなメッセージが出ればOK







tfrecordから画像の読み込み(input image)



だいたい、cifar-10のチュートリアル通り。けど、バッチを作るメソッドの

tf.train.shuffle_batch()

は画像枚数が大きすぎると、メモリを食い過ぎてショートする。tfrecordを複数にしても、ラベルが上手くシャッフルできないせいか、損失関数が上手く減少しない。ここら辺のロジックは、なかなかセオリー通りにいかないので、以下のコードで読み込ませた。

def distorted_input(data_files, batch_size, train=True, num_readers = 60):
    num_class=200
    if train:
        filename_queue = tf.train.string_input_producer(data_files, shuffle=True, capacity=16)
    else:
        filename_queue = tf.train.string_input_producer(data_files, shuffle=False, capacity=1)
    num_preprocess_threads = 60
    examples_per_shard = 1024
    min_queue_examples = examples_per_shard * 16
    if train:
        examples_queue = tf.RandomShuffleQueue(capacity=min_queue_examples + 3 * batch_size,
                min_after_dequeue=min_queue_examples, dtypes=[tf.string])
    else:
        examples_queue = tf.FIFOQueue(capacity=examples_per_shard + 3 * batch_size, 
                                      dtypes=[tf.string])

    # Create queue
    if num_readers > 1:
        enqueue_ops = []
        for _ in range(num_readers):
            reader = tf.TFRecordReader()
            _, value = reader.read(filename_queue)
            enqueue_ops.append(examples_queue.enqueue([value]))
        tf.train.queue_runner.add_queue_runner(
            tf.train.queue_runner.QueueRunner(examples_queue, enqueue_ops))
        example_serialized = examples_queue.dequeue()
    else:
        reader = tf.TFRecordReader()
        _, example_serialized = reader.read(filename_queue)
        
    images_and_labels = []
    for thread_id in range(num_preprocess_threads):
        image, label_index = parse_example_proto(example_serialized)
        images_and_labels.append([image, label_index])
    
    images, label_index_batch = tf.train.batch_join(images_and_labels,
             batch_size=batch_size, capacity=2 * num_preprocess_threads * batch_size)
    
    height = 150
    width = 150
    images = tf.reshape(images, shape=[batch_size, height, width, 3])
    
    return tf.subtract(tf.div(images,127.5), 1.0), tf.one_hot(tf.reshape(label_index_batch, [batch_size]), num_class)


def parse_example_proto(serialized_example):
    height = 150
    width = 150
    features = tf.parse_single_example(serialized_example,
                        features={"label": tf.FixedLenFeature([], tf.int64),
                                  "image": tf.FixedLenFeature([], tf.string)})
    label = tf.cast(features["label"], tf.int32)
    imgin = tf.reshape(tf.decode_raw(features["image"], tf.uint8), tf.stack([height, width, 3]))
    image = tf.cast(imgin, tf.float32)
    distorted_image = tf.image.random_flip_left_right(image)
    distorted_image = tf.image.random_brightness(distorted_image, max_delta=63)
    distorted_image = tf.image.random_contrast(distorted_image, lower=0.2, upper=1.8)
    distorted_image.set_shape([height, width, 3])
    return distorted_image, label

num_readers と num_preprocess_threadsの数を増やしてやれば、大容量でも上手く読み込めた。枚数が億単位になるとできるか不明。

batch数は32×GPU数(今回は32×8=256)が一番安定して、損失関数が減少した。






tfrecordファイルの作成について

tfrecordはprotocol buffer 形式で画像とラベルがペアで入ってる。けど、tfrecordにこだわらず、pickleファイルとか他のファイルでも問題ない。
tfrecordはたまたま使いやすいから、個人的に使っているだけ。枚数が増えてメルカリみたいに転移学習を用いる環境下ではむしろtfrecordは向いてないので、適材適所でファイルは使うべき。

ちなみに1ラベルごとに1つのtfrecordを作成したので、合計200ファイル。ラベルごとの作成過程は、
urlを呼び出してcsvに保存→s3からディレクトリに保存し→tfrecord作成

まで一連の過程を1ラベルごとにした。


label_0.csv → label_0(ディレクトリ)→train_label_0.tfrecords

こんな感じでやると以下の利点があって、一番効率がよく作成できた。



・問題が起こってもすぐ作り直せる

・for文で連鎖的に呼び出せる

path='/home/ubuntu/train_tf'
filenames = [os.path.join(path, 'record_%d.tfrecords' % i) for i in range(0, 200)]


・ミスなく作れる


ちなみに全tfrecordの中身の画像枚数は下のコードでカウント

# pathの読み込み
path=[os.path.join('/home/ubuntu/train_tf/record_%d.tfrecords' % i) for i in range(0, 200)]

# データ数カウント
cnts=0
for i, p in enumerate(path):
    cnt = len(list(tf.python_io.tf_record_iterator(p)))
    print("NO:{}, データ件数:{}".format(i, cnt))
    cnts+=cnt
print("合計データ件数:{}".format(cnts))





訓練(training)



モデルはinception_ resnet_v2を使用。




訓練のメイン部分

sess = tf.Session(config=tf.ConfigProto(allow_soft_placement=True, log_device_placement=False))
K.set_session(sess)

# lr = tf.train.exponential_decay(0.1, global_step, decay_steps, learning_rate_decay_factor, staircase=True)

def step_decay(epoch):
    initial_lrate = 0.01
    decay_rate = 0.5
    decay_steps = 8.0
    lrate = initial_lrate * math.pow(decay_rate,  
           math.floor((1+epoch)/decay_steps))
    return lrate


NUMN_GPUS = 8
batch_size=32*NUMN_GPUS

train_image, train_labels=distorted_input(filenames, batch_size, train=True)

with tf.device('/cpu:0'):
    train_model= InceptionResNetV2(train_image)
pmodel= multi_gpu_model(train_model, gpus=NUMN_GPUS)
pmodel.compile(optimizer=SGD(lr=0.01, momentum=0.9, decay=0.001, nesterov=True),
                    loss='categorical_crossentropy', 
                    metrics=['accuracy'], target_tensors=[train_labels])

# callback
history = History()
callback=[]
callback.append(history)
callback.append(LearningRateScheduler(step_decay))
  

tf.train.start_queue_runners(sess=sess)
coord = tf.train.Coordinator()
threads = tf.train.start_queue_runners(sess, coord)
try:
    pmodel.fit(steps_per_epoch=int(5101031/batch_size), epochs=3, callbacks=callback)
finally:
    train_model.save('/home/ubuntu/check_dir/inception_model_ep8.h5')
coord.request_stop()
coord.join(threads)

K.clear_session()


損失関数はこんな感じで順調に減少。3エポックくらいから減少しなくなった。





学習率
cifar-10とかで使ってる学習率の調整用メソッド
tf.train.exponential_decay()はkerasの場合、次で代用できる。

(tf.train.exponential_decayの引数「global_step」は kerasではepoch)

def step_decay(epoch):
     learning_rate = 0.01
    decay_rate = 0.5
    decay_steps = 8.0
    lrate = learning_rate  * math.pow(decay_rate,
           math.floor((1+epoch)/decay_steps))
    return lrate

これをcallback.append(LearningRateScheduler(step_decay))
でcallbackに入れれば、epochごとに最適な学習率に調節してくれる。パラメータは自分で調整した。

resnet系のoptimizerのベストプラクティスに関しての詳細はここら辺のサイトに詳しくのってる。→参考サイト

SGD+Momentumとnesterov=Trueがベストプラクティスらしい。パラメータはチュートリアルとか、論文読みあって総括して決めた






GPU環境下での訓練
kerasの場合、複数のGPUを使用するなら、modelをmulti_gpu_model()の引数に代入する必要がある。インスタンスタイプのp2_8xlargeはGPUを8こ搭載してるので、gpus=8にしてる。

またモデルの読み込みはCPU上でやらないと、OOMエラーが出る

with tf.device('/cpu:0'):
    train_model= InceptionResNetV2(train_image)

またOOMエラー対策かわからないけど、GPU使うときはtf.Session内で、次のおまじないが必要。

sess = tf.Session(config=tf.ConfigProto(allow_soft_placement=True,
log_device_placement=False))




モデルの保存
訓練はmulti_gpu_model()メソッドで作成したモデルでやってるけど、
保存はmulti_gpu_model()メソッドに入れたモデル(ここではtrain_model)を保存しないとダメ。このときモデルの保存にcallbackの
ModelCheckpointはmulti_gpu_model()メソッドのモデルを保存してしまうので使えない。

またモデルの保存は訓練終了時に指定。

coord = tf.train.Coordinator()
threads = tf.train.start_queue_runners(sess, coord)
try:
    pmodel.fit(steps_per_epoch=int(5101031/batch_size), epochs=3, callbacks=callback)
finally:
    train_model.save('/home/ubuntu/check_dir/inception_model_ep8.h5')
coord.request_stop()
coord.join(threads)








テスト・予測



テスト、予測のメイン部分のコード

sess = tf.Session(config=tf.ConfigProto(allow_soft_placement=True, log_device_placement=False))
K.set_session(sess)

NUM_GPUS=8
batch_size=125*NUM_GPUS

test_image, test_labels= distorted_input(filenames, batch_size)
with tf.device('/cpu:0'):
    test_model = InceptionResNetV2(test_image)
test_model.load_weights('/home/ubuntu/check_dir/inception_model_ep8.h5')
test_model= multi_gpu_model(test_model, gpus=NUM_GPUS)
test_model.compile(optimizer=SGD(lr=0.01, momentum=0.9, decay=0.001, nesterov=True),
                    loss='categorical_crossentropy', 
                    metrics=['accuracy'], target_tensors=[test_labels])


coord = tf.train.Coordinator()
threads = tf.train.start_queue_runners(sess, coord)

_, acc = test_model.evaluate(x=None, verbose=1, steps=100)
print('\nTest accuracy: {0}'.format(acc))

y_pred = test_model.predict(x=None,verbose=1,steps=100)
LABEL=[]
for i in range(100):
    LABEL.append(sess.run(test_labels))
top_k=sess.run(tf.nn.top_k(y_pred,k=3,sorted=False))
coord.request_stop()
coord.join(threads)

test_model.evaluate(x=None, verbose=1, steps=100)のstepsとバッチサイズは、読み込むテスト画像とぴったり合うようする必要がある。
今回はテスト画像枚数が10万枚なので、steps=100 、batch_size=125*8
にして、1000バッチを100回繰り返して、10万枚の画像を読み込んだ。

また評価用にラベルもtfrecordから取り出した。

LABEL=[]
for i in range(100):
    LABEL.append(sess.run(test_labels))

あとは混合行列と、f値で評価

# Accuracy , F-score.etc
print('acuracy:{}'.format(accuracy_score(label_batch,f_pred)))
label_string = ['{}'.format(i) for i in range(200)]
print(classification_report(label_batch, f_pred,target_names=label_string))





最後に


今回は、

・環境が整ってない
・5百万単位の大規模データでの訓練

の2つが要因でエラーが多々発生した。
反省点としては


・停滞したら、立ち止まって調査、考える時間にあてる
・思考時間とPC作業の割合は7:3くらいが自分にはベストプラクティス

これ以上枚数が増えると、一部の界隈(メルカリ)みたいに

転移学習+次元削減(PCAとか)+faiss

を使うような、効率性を追求した環境構築の必要あると思った。




EC2、S3関係のコマンドメモ



・キーの作成
# ec2インスタンスへのログインキーの作成・保存

cp path/to/フィル名.txt ~/.ssh/任意の名前

chmod 400  ~/.ssh/名前


sshでログイン

$ ssh -i ~/.ssh/任意の名前 root@ 

・scpコマンドで、フォルダorファイルコピー
# EC2からローカルにファイルのコピー

scp -i ~/.ssh/aws_develop_sec /Users/Downloads/test.tfrecords ubuntu@13.231.234.143:/home/ubuntu/test.tfrecords 

# ローカルからEC2へファルダのコピー

scp -r -i ~/.ssh/aws_develop_sec /Users/Downloads/prediction_img centos@13.231.108.156:/home/centos/prediction_img

# EC2からローカルへフォルダのコピー

scp -r -i ~/.ssh/aws_develop_sec centos@13.231.234.17:/home/centos/test /Users/Downloads/tfrecords/test 


# EC2からs3へフォルダのアップロード

aws s3 cp a s3://cnn-image --recursive --acl public-read --cache-control “max-age=604800"

# S3にディレクトリごとアップロード

$ aws s3 cp <ディレクトリ名> s3://<バケット名>/<バケット名>/ --recursive --acl public-read --cache-control "max-age=604800"


# sshのログインがキレてもpythonを動かし続けられるコマンド(バックグラウンド→ https://qiita.com/alancodvo/items/15dc36d243e842448d33

$ nohup python Inception_keras_gpu_train.py & 

xgboostの回帰モデルで精度検証から重要な特徴量選択までやってみた全行程まとめ(機械学習)

今回はkaggleでよくある特徴量エンジニアリングのテクを使って、精度向上から重要な特徴選択までをやった。普通は精度高ければ終わり的な感じだけど、今回は精度検証からさらに掘り下げて、特徴量の選択までやったので、その過程を書いてく。

目次
・プロジェクト紹介
・デーセット
・特徴量作成
・モデル構築
・正解率
・特徴量の探索(分析)
・最後に(反省点)

コードはこちら




プロジェクト紹介


一回しか購入しないユーザー(u1)と2回以上購入しているユーザー(u2)を分ける重要な特徴量を見つけるプロジェクト。目的は一回しか買わないユーザーに有効な施策をうつこと。

よくある特徴量エンジニアリングで、回帰モデルを使って精度を競うkaggleとかのコンペと似ている。けど、今回は精度検証後に重要な特徴量選択までしたので、特徴量選択の分析過程も含まれてる。全体の工程としては、特徴量作成が7、8割を占めた感じ。

特徴量エンジニアリングが、精度向上に重要な過程だと再認識。





データセット


データセットはデータベースに保存してあるテーブルのレコード(特徴量)から選択して作成。kaggleみたいにあらかじめ綺麗なのが用意されてないし、データベースとbigqueryに別々に保存してある。そんな中、無数にある特徴量から適切なのをSQLで取り出して、ベストな特徴量を作成。

最終的な特徴量は、2回目の検証結果から決めた。アルゴリズムでの検証は全体で3回。

まず1回目は商品に結びつきそうな特徴量を選んで検証。1回目はtouroku_dateという特徴量が一番重要らしく、次の特徴量を選ぶ指針になった。

2回目は施策から「u1が1回で購入をやめてしまう」ことに関する仮説を立てて、この仮説に直結する特徴を選んで使うことにした。
仮説の設定はAmazonとか、リピート率が鬼のように高いサイトレビューとかを参考にした。


一回で購入をやめてしまう or 継続して購入する理由の仮説

=> 評価の悪いユーザーから購入してる
=> レコメンド関係の施策で欲しい商品を購入してる
=> 欲しい商品が見つからない(時間をかけて探していない)



3回目は特徴量を2個追加して、若干の修正後に再検証。






***仮説を立てることで作業を効率化
無数にある特徴量の中から、適切なものを選択する過程で、施策に応じて仮説を立てる作業はかなりの効率化につながった。はじめに仮説を立てるのは、データサイエンティスト的な手法らしい。

さらに、選んだ特徴量から、u1とu2間の違いをあらかじめ調べる作業をする事で、さらに効率化。違いを見るのには「標準偏差・平均・中央値」とかの分析手法を使った。

特徴量は何を選べばいいかわからない状況下で、

・仮説の設定

・違いを見る作業

この二つをすることで、かなり効率化できた。

特徴量に必要な特徴量を選択するまでの、全行程はこんな感じ

1、まず分類に使えそうな特徴量を作成・検証

2、試作にあった仮説を立てる(重要)

3、仮説に直結する特徴量を選んで、標準偏差とかで違いをみる(重要)

4、特徴量を決めて、作成

2、3の作業が一番精度向上につながった。

最終的なデータセットは、
特徴量10こ
u1は250853件、u2は1169898件のデータを使った。




特徴量作成


特徴量の作成過程は上で書いた通りなので、ここではデータベースのレコードをSQLで取り出し、それを使いやすい形に変えて特徴量を作成する特徴量エンジニアリングのテクニックについて触れたい。

特徴量エンジニアリングでは主に質的変数(赤、青、白とかみたいな数字じゃない形)をダミー変数に、量的変数(体重とか身長みたいな数値)を使いやすく変形させるテクがある。

特にダミー変数は離散化とか、いろいろあるので、ここら辺の資料を参考にした。

機械学習レシピ#3

データの前処理だけで競馬は強くなる

特徴量の形を使いやすい形に変えるのも、めっちゃ大事。
今回は可視化した結果の分布が、「0が多く、0以上が少ない」場合、SQLでこんな感じに取り出すのが有効だった。

SUM(case when レコード= 値 then 1 else 0 end) 


逆に分布がバラバラな特徴量は、AVG() で平均をとるか、そのまま使うのが有効だった。

特徴量は検証前に、標準化しないと分類精度がすごい下がるので、データ作成後は必ず標準化した。






モデル構築


モデルは、kaggleで有名なランダムファレストの改良版のxgboostを使った。
特にGBDTの改良手法Dropouts meet Multiple Addtive Regression Trees(DART)とパラメーターチューニングのhyperoptの組み合わせが、xgboostとドンピシャで、かなりいい分析精度になった。

xgboostのパラメータは山ほどあるので、チューニングはベストプラクティスを載せてるサイトがを参考にした。

Parameter_Tuning_XGBoost_with_Example

パラメータチューニングはhyperoptとgridseachを使い分けて、手早くチューニング。

精度検証で、交差検証はもはや当たり前で、さらに複数のアルゴリズムでの検証や、複数の精度指標での試行錯誤も王道らしい。

def objective(params):

    skf = cross_validation.StratifiedKFold(
        train_y, # Samples to split in K folds
        n_folds=5, # Number of folds. Must be at least 2.
        shuffle=True, # Whether to shuffle each stratification of the data before splitting into batches.
        random_state=30 # pseudo-random number generator state used for shuffling
    )

    boost_rounds = []
    score = []

    for train, test in skf:
        _train_x, _test_x, _train_y, _test_y = \
            train_x[train], train_x[test], train_y[train], train_y[test]

        train_xd = xgb.DMatrix(_train_x, label=_train_y)
        test_xd = xgb.DMatrix(_test_x, label=_test_y)
        watchlist = [(train_xd, 'train'),(test_xd, 'eval')]

        model = xgb.train(
            params,
            train_xd,
            num_boost_round=100,
            evals=watchlist,
            early_stopping_rounds=30
        )

        boost_rounds.append(model.best_iteration)
        score.append(model.best_score)

    print('average of best iteration:', np.average(boost_rounds))
    return {'loss': np.average(score), 'status': STATUS_OK}

def optimize(trials):
    space = {'booster':'dart',
         'learning_rate':0.1,
         'n_estimators':1000,
         'sample_type':'uniform',
         'normalize_type': 'tree',
         'objective':'binary:logistic',
         'min_child_weight':1,
         'max_depth':9,
         'gamma':0.0,
         'subsample':0.6,
         'colsample_bytree':0.9,
         'reg_alpha':1e-05,
         'nthread':4,
         'scale_pos_weight':1,
         'seed':27,}
    
    # minimize the objective over the space
    best_params = fmin(
        fn=objective,
        space=space,
        algo=tpe.suggest,
        trials=trials,
        max_evals=10
    )

    return best_params

# parameter tuning
trials = Trials()
best_params = optimize(trials)
print(best_params)

# 損失関数(loss)の計算
print(objective(best_params))






正解率


正解率は


AUCで89.33%



MCCで0.548(-1~1の範囲をとり、1で100%の精度)



特徴量の重要度


評価指標はAUCがよく使われるけど、不均衡データとかより正確な精度判定にはMCCの方がいいっぽい。

from sklearn.metrics import matthews_corrcoef
thresholds = np.linspace(0.01, 0.99, 50)
mcc = np.array([matthews_corrcoef(test_y, pred_y>thr) for thr in thresholds])
plt.plot(thresholds, mcc)
best_threshold = thresholds[mcc.argmax()]
print(mcc.max())

この他にも混合行列とかF値を使ってもいいけど、今回はMCCがF値と同じくらい精密なので、MCCで代用。




特徴量の探索(分析)



特徴量の重要度、そのツリー、xgboostの結果から、u1とu2に有効な特徴量選択。検証からさらに特徴量を探すまで深掘りするケースはサイトでは見つからなかったので、このフェーズが一番大変だった。
とりあえず、一般的な統計手法とか、マーケット分析事例で使われてるっぽい分析手法を漁った。

あとは知識、ひらめき、そしてひたすら根気の作業。
結果、ヒストグラムクラスタリング、次元削減、重回帰分析とかが使えるっぽくて、一番合理的な「重回帰分析」を使った。

重回帰分析は、ある一つの特徴量と、他の複数の特徴量の相関関係を計算してくれる

結果的に、

u1=>海外商品系の特徴量

u2=>満足度が高い特徴量

との相関が高く、これらが重要な特徴量と判断。

結論とその理由とかのロジック的な部分は、他の特徴量との相関関係の数値、仮説の消去法で導いた感じ。ちゃんと根拠つきで




最後に(反省点)



kaggleとかの特徴量エンジニアリングを使った分析は、あらかじめ綺麗な特徴量が用意されていて、分析精度が良ければ終わり、みたいなケースが多い。
けど今回はさらに掘り下げて、重要な特徴量の選択までやった。前例がないとほぼ手探りなので、かなり大変だった。

とりあえず、反省点として

・試行錯誤より頭を使って、作業量を減らす(think more, do less)

・行き詰まったり、停滞したら、一回止まって問題の分解とか本質の問題探しをする。絶対にただタスクをこなすだけのループにはまらないこと
(今回はマッキンゼーの問題解決法の書籍が本質の問題探しにかなり使えた
イシューからはじめよ―知的生産の「シンプルな本質」
世界一やさしい問題解決の授業―自分で考え、行動する力が身につく


・チームで動く以上、進歩確認は随時やる


・データ量が多ければ、少ないデータ量(1万件くらい)で試す


・やることの意思疎通にすれ違いがないようにして、やり直しは極力やらない


・段取り8割


・わからないことは、詳しい人に積極的に、遠慮なく聞くこと

あえてxgboostのDARTみたいなハイリスク・ハイリターンのスキルを使うことで、かなりのスキルが身についたなと実感した。




kaggleのテクニック系の参考資料
top2%の私が教えるKaggleの極意

Kaggleで使われた特徴量エンジニアリングとアルゴリズムまとめ

Djangoでアップロードした画像をCNNで予測し、結果を返すアプリを作ってみた(画像認識、機械学習)

kerasで作った画像分類器に画像を読み込ませ、予測したラベルのidを返すアプリ作った。以前、rubyで作ったことがあるけど、今回はpython専用のフレームワークDjangoを使って作成。

画像分類器にはCNNを使ったので、GPUとか学習のところは割愛して、アプリ作成過程についてだけ書いてみる。

目次
1.アプリについて
2.メインのファイル



1.アプリについて


一通りの動作
デフォルトのホーム画面はこんな感じでシンプル

写真をアップロードして、クリックすると予測した確率が高いトップ3のラベルのidが表示される仕組み


アプリ名と構成ファイル
project名はimage_pred

アプリ名はmyapp

メインの構成ファイルは
・setting.py
・forms.py
・views.py
・index.html
・main.py
・urls.py

いろんなファイルを保存するディレクトリは、必要に応じて作成した。




MVCフレームワークについて

DjangoでのMVCフレームワーク

・Mはmodels.pyでデータベースの操作

・Vはviews.pyで画面の表示を操作

・Cはurls.pyでアクセス関係の操作

今回はデータベースにデータを保存しなくていいので、models.pyは使わない。使うときはマイグレーションが必要で、使えるデータベースは任意に指定可能

便利なパッケージでdjango-cleanupがあった。データベースを使うならpipでインストールして、setting.pyのINSTALLED_APPSに入れておきたい





2.メインのファイル


まずパッケージは
・Pillow
・keras
・tensorflow
opencv

とかをpipでインストール





setting.py
myapp作成後に、INSTALLED_APPSにmyappを追加。

画像をアップロードするので以下のディレクトリを設定
・MEDIA_ROOT
=>サーバから見たメディアルートの絶対パス

・MEDIA_URL
=>メディアファイル公開時のURLのプレフィクス。 メディアファイルのURLは「http://アプリのドメイン+MEDIA_URL+メディアファイル名」

# 一部掲載
import os

# Build paths inside the project like this: os.path.join(BASE_DIR, ...)
BASE_DIR = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))

SECRET_KEY = '8bdf)xh2b4h!8#g$i0j_9pymjld*ov19!rpfgj2qsi6j)-$t9d'

DEBUG = True

ALLOWED_HOSTS = []


# Application definition

INSTALLED_APPS = [
    'django.contrib.admin',
    'django.contrib.auth',
    'django.contrib.contenttypes',
    'django.contrib.sessions',
    'django.contrib.messages',
    'django.contrib.staticfiles',
    'myapp' # 追加
]

STATIC_URL = '/static/'
MEDIA_ROOT = os.path.join(BASE_DIR, 'media')
MEDIA_URL = '/media/'






forms.py
ここは画像をアップロードする専用のフォームを作成するファイル
画像専用フォーム「forms.ImageField()」を使用

from django import forms

class PhotoForm(forms.Form):
    
    image = forms.ImageField()






index.html

ここはメインのホーム画面を表示するhtml。主にしてる処理は
・main.pyで予測したラベルのidをhtmlの変数に埋め込み、レンダリング
・フォームの設定
・静的ファイル(CSS)の読み込み
など。htmlにpythonの変数を直接、埋め込めるのがDjangoの強み

{% load static %}
<!doctype html>
<html lang="ja">
<head>
 <meta charset="utf-8">
 <title>img_prediction</title>
 <link rel="stylesheet" type="text/css"
      href="{% static 'myapp/css/style.css' %}" />
</head>
<body>
<p>{{pred}}</p>
<table>
<form action="{% url 'index' %}" method="POST" enctype="multipart/form-data">
   {% csrf_token %}
   {{ form }}
   <tr><td></td><td><input type="submit" value="click" /></td></tr>
</form>
</table>
</body>
</html>


静的ファイルは{% load フォルダ名 %}を記述してCSSを読み込ませた。 staticファルダ内にstyle.cssが保存してある。javascriptとかも同じように読み込める

変数は{{pred}}の部分。

フォームはformタグ内に記述。画像の場合はenctype="multipart/form-data"を指定しなきゃだめ







main.py
コードはkerasで作成した学習済みの画像分類器(CNN)と予測結果を返す処理。学習済みの重み、id用のcsvファイルは専用フォルダを別途作成して、保存してる。

# pred関数のみ
def pred(img_path):
    sess = tf.Session()
    K.set_session(sess)
    cnt=pd.read_csv('/Users/d/image_pred/myapp/weight_dir/cnt.csv')
    cnt=cnt.drop('Unnamed: 0', axis=1)
    cnt=cnt.drop('cnt', axis=1)
    
    
    dic = {}
    for i, v in cnt.iterrows():
        dic.setdefault(v['index'], []).append([v['brand_id'], v['model_id'], v['cate_id']])
    
    x=np.asarray(Image.open(img_path))
    x = cv2.resize(x, (100, 100))
    x = np.expand_dims(x, axis=0)
    image = preprocess_input(x)

    test_model = InceptionResNetV2(include_top=True)
    test_model.load_weights('/Users/d/image_pred/myapp/weight_dir/incep_model.h5')
    test_model.compile(optimizer=SGD(lr=0.01, momentum=0.9, decay=0.001, nesterov=True),
                    loss='categorical_crossentropy', 
                    metrics=['accuracy'])

    y_pred = test_model.predict(image)
    top_k=sess.run(tf.nn.top_k(y_pred,k=3,sorted=True))
    idxs=list(np.reshape(top_k[1],(3,)))
    id=[]
    for idx, v in dic.items():
        if idx in idxs:
            id.append(['brand_id:{} model_id:{}, cate_id:{}'.format(v[0][0], v[0][1], v[0][2])])
    return id





views.py
ここは、index.htmlへのレンダリングなど、表示関連の処理を書くとこ。

主に、get時の処理とpost時の処理を関数に分けて書いてある。get時の処理はアクセス時のデフォルトの画面表示。

post時の処理はアップロードして画像を受け取って読み込み、予測したラベルのidを返す処理を書いてる。
フォームからアップロードされた画像を読み込み予測ラベルを返す処理はmain.pyのpred関数をimportで呼び出す。そして、変数を代入して、その結果をhtmlにレンダリング。画像じゃなきゃエラーを返す仕組み

from django.shortcuts import render, redirect
from django.http import HttpResponse
from django.views.generic import TemplateView
from .forms import PhotoForm
from myapp.main import pred

class MyappView(TemplateView):
   def __init__(self):
       self.params={'pred': 'idx',
                    'form': PhotoForm()}
   def get(self, req):
       return render(req, 'myapp/index.html', self.params)

   def post(self, req):
       form = PhotoForm(req.POST, req.FILES)
       if not form.is_valid():
           raise ValueError('invalid form')

       image = form.cleaned_data['image']
       self.params['pred'] = pred(image)
       return render(req, 'myapp/index.html', self.params)





urls.py

Djangoではプロジェクトの下にいくつもアプリを作れる。Djangoでは「1アプリ=1 url」が基本。つまりurls.pyはアプリを作るたびに、プロジェクトの階層以外にも。アプリ内に必ず一つurls.pyを配置する仕組み




・プロジェクトの階層のurls.py
プロジェクトの階層では全てのアプリのurls.pyを管理するpathを指定する。「アプリのことは各アプリ内のurls.pyに聞け」って意味のinclude()を使う

urlpatterns = [
   path('admin/', admin.site.urls),
   path('myapp/', include('myapp.urls')),]


・各アプリ内のurls.py
各アプリ内ではurlspatternにurlを指定する。引数は3つあって、1つめはhttp://~/myapp/の後に続くアドレス、2つめはアクセスするファイル、3つめは指定したurlの名前


from django.conf.urls import url
from .views import MyappView
urlpatterns = [
              url(r'', MyappView.as_view(), name='index'),
              ]


画像をアップロードするたびにファルダにurlと画像が保存されるようにするには以下を追加。

if settings.DEBUG:
urlpatterns += static(settings.MEDIA_URL, document_root=settings.MEDIA_ROOT)




動作確認

画像をアップロードしてクリックを押す

すると、予測ラベルのidが表示される。

idを返す超単純なアプリを作成した。Djangoの勉強とアウトプットには、やっぱり基礎を抑えて自分でアプリとか作るのが一番っぽい


参考記事

Django - 画像ファイルのアップロード処理

Django-ファイルアップロード機能の使い方 (基本設定編)

Djangoで、けものフレンズキャラの顔を認識させる(Deep Learning)

LSTMを使って異常検知手法で急上昇ワードをやってみた(機械学習-変化点検知)

急上昇ワード(バースト検知)とえば、googleyahoo!でも似たようのがある。今回はあれほど高性能じゃないけど、急上昇ワードと同じ仕組みのものを異常検知手法でやってみた。

目次
・訓練
閾値の設定
・テスト
・実際に急上昇ワードをやってみる

今回はBigQuery(BQ)から抜き出した、一定期間内の検索ワードの急上昇を検知してみる。みんな大好きニューラルネットワークのLSTMで訓練し、バースト検知に利用した。訓練から、実際にBQから取り出した検索ワードの急上昇検知実行までの流れを一通り書いて行く。

実際の全体コードはgithubにあげてあります。




訓練


まず訓練データ、「mizugi」と「bts」というワードを2つ、8ヶ月分ほど抽出。btsは2桁(10台)の検索数が多く、mizugiは検索数が多く、3桁(100〜1000)のワードだ。
まずLSTMでbtsを8ヶ月分くらい、急上昇部分があっても構わず学習させた。


btsとmizugiの急上昇部分は特異変換スペクトルで検出した。結果がこれ

bts

mizugi

アルゴリズムはLSTMを使った。QRNNとかいう発展版も試したが、LSTMの方がkerasで簡単に使えるし、十分なので途中で断念。

データをLSTMで使えるようにする正規化とかの前処理と、モデル作成はこんな感じ

def get_data(data, time_steps):
    docX, docY = [], []
    for i in range(len(data)-time_steps):
        docX.append(data[i:i+time_steps])
        docY.append(data[i+time_steps])
    alsX = np.array(docX)
    alsY = np.array(docY)
    return alsX, alsY

def transform_data(data, inverse_option, scaler):
    data_shape = data.shape
    if inverse_option is True:
        data = scaler.inverse_transform(data)
    else:
        data = scaler.fit_transform(data)
    data = data.reshape(data_shape)
    return data, scaler

def prepare_data(original_data, time_steps):
    copy_data = original_data.copy()
    scaler = MinMaxScaler(feature_range=(0, 1), copy=False)
    data, scaler = transform_data(data=copy_data, 
                              inverse_option=False, scaler=scaler)
    
    x,y = get_data(data, time_steps=time_steps)
    return x, y, scaler

time_steps = 12

x, y, scaler = prepare_data(X_train, time_steps)

input_dim=x.shape[-1]
timesteps=x.shape[1]

model = Sequential()
model.add(LSTM(300, input_shape=(timesteps, input_dim),
         stateful=False,return_sequences=True))
model.add(Flatten())
model.add(Dense(1, kernel_initializer='lecun_uniform'))
model.add(Activation("linear"))
# model.load_weights('/Users/d/model.ep200.h5')
model.compile(loss="mean_squared_error", optimizer="adam",)

以下が訓練後、btsの学習期間を予測した図。赤がリアルの検索数値で青が予測値だ。btsは2桁の値なのか知らないけど、損失関数はあまり減少しなかった。学習結果としてはいまいちだったけど、それなりによく予測できてたかな。

btsの予測


閾値の設定


ここからbtsとmizugiから以下のように検証・テスト1・テスト2の3つにデータを分けた

・検証データ→btsの急上昇のない部分(異常なしデータ)

・テストデータ1→btsの急上昇のある部分(異常ありデータ)

・テストデータ2→mizugiの急上昇のある部分(異常ありデータ)


バースト検知は異常検知の中で「変化点検出」に当たる。

それには少し変わった閾値設定をする必要がある。以下の流れで閾値を設定した。

①検証データはオリジナルの閾値の設定に使用

②検知したいワードの平均値をオリジナルの閾値にかけて、その検知したいワードに最適な閾値を再設定する


変化点検知の閾値設定には、回帰(予測)モデルと、mseの組み合わせから求めるのが、結構メジャーらしい(参考記事)。




①検証データはオリジナルの閾値の設定に使用

今回の閾値の設定には、検証データを使って最初の(オリジナルの)閾値を求めた。コードは以下の通り。

def calculate_mse(value, predict_value, variance=1.0):
    value = value[:, 0]
    predict_value = predict_value[:, 0]
    mse_value = [(v - p_v)**2 for v, p_v in zip(value, predict_value)]
    return np.array(mse_value)


mse_value_valid = calculate_mse(x_scale_valid, predict_valid)
threshold = np.max(mse_value_valid)
# threshold:249.16269761799867

検証データ(異常なしデータ)の全期間から、MSEを求めて、その最大値を閾値とした(つまり、検証データの異常値の最大値)。




②検知したいワードの平均値をオリジナルの閾値にかけて、その検知したいワードに最適な閾値を再設定する


ここで問題なのは。btsは2桁の値なので閾値はかなり低い。これを水着のような3桁、4桁以上の検索数のワードに適用すると閾値が小さすぎて、全部の異常値を検出してしまう。


そこで検知するワードに対して最適な閾値を再設定してやる必要がある。以下のやり方で設定した。


1.検証データから求めたMSEの最大値をオリジナルの閾値とする

2.検知したいワードの全期間の検索数の平均値割る2の値を四捨五入して整数にする

3.それをオリジナルの閾値にかけて、検知したいワードの閾値(high_threshold)として再設定

わかりにくいので、コードで表すとこんな感じ

mse_value_normal = calculate_mse(x_scale_train, predict_train)
N=int(np.mean(x_scale_train)/2)
high_threshold = threshold*N 

なんで2で割ったのかは、いろいろ試して一番検出に適した値になるからなので、別に整数(int)じゃなくて、小数(float)のままでもよかった。





テスト


再設定した閾値で、テストデータ1、2のバースト(異常)を上手く検知できてるか試す。

上からテストデータ1(bts)、テストデータ2(mizugi)の時系列データと、再設定した閾値をプロットした

bts

mizugi

元の閾値よりも高い閾値で、程よく急上昇部分が検知されてる。一日、1週間、一月の期間でも試したが、上手く検出できてた。





実際に急上昇ワードをやってみる


いよいよ、バースト検知を試してみる。今回はjupyter上から直接、検索ワードをBigQueryから抽出して、検出してみた。


急上昇検知のメイン部分はこんな感じ

csvs=[]
original_threshold = 270.2200 
train_data_mean=4.964
num_words=len(words)
for i in range(num_words):
    cnt=pd.DataFrame(list(cnt_listed[i][0])).dropna()
    hour=pd.DataFrame(list(hour_listed[i][0])).dropna()
    hours=np.array(hour[window:])
    wo=cnt_listed[i][1] # each word
    
    x_predict, y_ ,scaler = prepare_data(cnt, time_steps)
    predicted, x_scale, mse_value = predict_model_show_graph(hour[window:], x_predict, y_, scaler, model, time_steps)
    mse_value = calculate_mse(x_scale, predicted)
    
    N=int(np.mean(x_scale)/2)
    high_threshold = original_threshold*N  # 閾値
    if high_threshold<original_threshold:
        #print('high_threshold: {}'.format(original_threshold))
        #show_graph_threshold(hour[window:], mse_value, original_threshold, 'Anomaly Score test', "r")
        # 急上昇ワードの時間を表示
        for mse, hour in zip(mse_value, hours):
            if mse >=original_threshold:
                csvs.append([hour, wo])
    else:
        #print('high_threshold: {}'.format(high_threshold))
        #show_graph_threshold(hour[window:], mse_value, high_threshold, 'Anomaly Score test', "r")
        # 急上昇ワードの時間を表示
        for mse, hour in zip(mse_value, hours):
            if mse >=high_threshold:
                csvs.append([hour, wo])

その前の前処理では、取り出したワードを時系列データに変換したり、LSTM用に正規化したりした。

今回は3月3日の分だけ、検索ワードを取り出して、検知。図をplotとするといい具合に検出できてるっぽい。

最終的には、ワードと時間を抽出する形にした。


閾値の設定は前に書いた動的時間伸縮法(DTW)でもできるけど、計算量がバカにならない。なので、MSEとかのローコストな計算式を臨機応変に使い分けるのが最適だと思う




今回はLSTMを用いた、変化点検出手法で急上昇ワードのバースト検知をやってみた。
昔はchange finderとかsmartfilterとか、ゴツい検出アルゴリズムがあったけど、ニューラルネットワークの登場で異常検知もかなりやりやすく、精度も高くなった。

異常行動検知もLSTMを使ってできる。機械学習シリーズにある密度比推定とか非構造学習みたなガチもんの専門的な内容には、機会があったら触れてみたい。

列名とかデータ加工に使うpandasの便利な機能まとめ (python・機械学習)

今回はデータ加工に使えるpandasの機能を紹介する。kaggleを含め、機械学習のデータ加工はpandasでの加工が多い。

理由は単純にpandasはデータ加工において、扱いやすいから。今回はxgboostの特徴量を加工する機会があった。そのときに使ったり、調査したkaggleで人気なpandasのメソッドの中で使えそうなやつをまとめた。


機能1〜5まではここではirisデータセットで試した。

iris = datasets.load_iris()
df_feature_iris = pd.DataFrame(iris['data'],columns=['sepal_length','sepal_width','petal_length','petal_width'])
df_species_iris = pd.DataFrame(iris['target'],columns=['species'])
df=df_feature_iris
df.head()

用語について

pandasで使う用語で
・行
・列

がある。pandasは行列のmatrix形式で扱ので、ググればすぐわかるが、
行=>縦
列=>横
と考えると分かりやすい

行が横方向にならんだものを示し、列が縦方向にならんだものを示す。表計算ソフトでも、横方向を行と言い、縦方向を列と言う(参考Wikipedia


pandas機能(1~19)


機能1
==行の複数条件指定==

これは列の内容で2つ以上の条件を指定する方法。

・純粋に&で指定するもの
SQLのクエリのように指定する方法
の2つがある。クエリライクのやり方はkaggleで人気

df[(df['sepal_length']>=7)&(df['sepal_width']>=3)]

#またはクエリのように指定する方法
df.query("sepal_length>=7 and sepal_width>=3")


機能2
==特定の型の列のみ取得==

列方向で特定のデータ型の値を取得する
pandasの場合、データ型は「整数('int')・小数('float')・真偽値('bool')」以外は全て「’object'」で指定

df.select_dtypes(['float']).head()


機能3
==列Aごとに列Bを集計する==

列Aのグループを作り、その値内で、列Bの値の数を集計する。

→irisデータで試す
'sepal_length'の列でグループを作り(4.3〜7.9)、'sepal_width'の値(2.0~4.4)を
値ごとに集計

df.groupby('sepal_length')['sepal_width'].value_counts().unstack().head()

‘sepal_length'が4.5で'sepal_width'が2.3の数=>1




機能4
==2つの列A、Bに着目し、3つ目の列Cの値を集計する==

これはわかりにくいので、irisデータで試した。


'sepal_length'でグループを作り(4.3〜7.9)、'sepal_width'の値(2.0~4.4)に
該当する'petal_length'の平均値を集計

df.pivot_table(index='sepal_length', columns='sepal_width',values='petal_length', aggfunc=np.mean).head()

‘sepal_length'が4.5で'sepal_width'が2.3の'petal_length'の平均値は1.3



機能5
==特定の型の削除==

「整数('int')・小数('float')・真偽値('bool')、それ以外(’object')」を指定して削除

df.select_dtypes(exclude=np.int)


機能6
==全列の欠損値カウント==

欠損値(Nanのデータ)を数える

print(df.isnull().sum())

機能7
==欠損値の値を置換==

欠損値を別の値に置換する。

例:欠損値を0に置換

df=df.fillna(0)

機能8
==標準化==

データを機械学習のモデルに読み込ませるとき、正規化(標準化)することが多いので、その手法の一つが標準化
numpyとかでもできるが、今回はScikit-learnを使った。

正規化=>データ等々を一定のルール(規則)に基づいて変形し、利用しやすくすること。

標準化=>元データを平均0、標準偏差が1のものに変換する正規化法の一つ

# 標準化
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
dff = sc.fit_transform(df)
dff=pd.DataFrame(dff)

# 列名を元に戻す
listed=list(df.columns)
dff.columns = listed

機能9
==行列内の文字列の末尾のN個削除 (’iraidate'は列名)==

列の文字列の末尾のN個の文字を削除する

b=[]
for i,v in f2.iterrows():
    b.append(str(v['iraidate'])[:3])

# bのDataFrame作成
f=pd.DataFrame(b, columns=['iraidate'])

# 前の同じ行削除
df=df.drop("iraidate", axis=1)

# 横に連結してもとに戻す
df=pd.concat([df, f], axis=1)


機能10
==データの散布図==

Dataframe(df)とplotしたい列名2つを指定

sns.jointplot('kaiin_id','safety_tesuryo', data=df)


pandasデータの散布図の可視化に便利





あとはそのまんまの機能(機能11〜19)

・全列の名前表示

df.columns

・列名の変更

df=df.rename(columns={'last_login_year_cnts1': 'last_login_year_cnt1'})

・列の削除

df=df.drop('列名', axis=1)


・列名の一括変更(df_colum_listに全列名を入れておく)

df.columns=df_colum_list

・地味に全列名の一括変更

df=df.ix[:,['A','B','C']]

・縦にDataFrame(df)を連結

df=pd.concat([df, df2])

・横にDataFrame(df)を連結

df=pd.concat([df, df2], axis=1)

・list(listed)から列名付きで、DataFarame(df)作成

df=pd.DataFrame(listed, columns=['列名'])


・DataFrame(df)をcsvに保存/読み込み

# 保存
df.to_csv('feature.csv')

# 読み込み
df=pd.read_csv('feature.csv')


備忘録もかねてpandasの機能をまとめてみた。使えるものも多いと思うので、気に入ったらストックしておきたい。

異常検知(変化検知)の詳細と特異変換スペクトルと動的時間伸縮法まとめ(機械学習)

異常検知とは、機械学習の一手法で、普通の値のデータの中から極端に大きかったり、小さい値の「異常」なデータを見つけ出すものだ。


異常検知の用途で有名どころは巷では以下のようなものがメインらしい  

 

マーケティング =>流行のブレイクの検出
コンプライアンス  =>情報漏洩/不正検知
・製造系   =>不具合発生検知
・機械系  =>故障検出
・システム運用系  =>障害検出
・ネットワーク運用系  =>トラフィック異常検出
・WEB系  =>トピック潮流変化検出
・セキュリティ系  =>攻撃・侵入検出(コンピュータウィルスやDOS攻撃
・サイバー犯罪に使うやつ=>情報潮流、なりすまし


いずれも普通じゃないデータとか、珍しいデータを検出するものだ。


異常検知の3種類
異常検知は主に3つのタイプに分かれる。機械学習シリーズの書籍とかに詳しく載ってるので、ここでは簡単に触れとく。


まず触れておきたいのは異常検知で扱うデータには「定常データ」と「非定常データ」がある。異常検知ではこの2つで使う手法が異なる。
定義は以下の通り


定常データ:同じ波形が続くような時系列データや、統計的にデータ分布図が正規分布(山が一つの単純な曲線)になる


非定常データ:同じ波形が続かない時系列データや、統計的なデータ分布図が正規分布にならない(山が2つ以上の曲線)


・外れ値検知
すでに持ってる(既存の)データから異常なデータや異常箇所を割り出す手法。
正常データと異常データの割合→99:1


くらいのデータ、つまり正常データがほとんどで異常データ少ししかないデータが対象になることが多い。


外れ値検出のほとんどのアルゴリズムは定常データを扱う。
その時は、マハラノビス距離で求められるアルゴリズムに限定させることが多い。代表的なものに、k-means、ランダムフォレスト、one class SVMなどがある

 


・変化検知
これは主に定常/非定常な時系列データなどに多く使うもので、データが大きく変化する点を検知する。googleである「急上昇ワード」が好例だ。扱うデータは非定常データであることが多い。


変化検知で、非定常データを扱うときのアルゴリズムには、いろいろある。メジャーなのは
・特異変換スペクトル
・動的伸縮法
・LSTMなどの時系列を予測するアルゴリズム


などを使うケースが多い。


・異常行動検知
これは多くのユーザーと違う行動をとるユーザーを検出するときに使う。
手法としては、同じ行動をとる多くのユーザーの行動を統計的に確率分布でパターン化し正常データとして扱う。逆に、その確率分布のパターンと異なる分布のデータになるユーザーを異常行動をしたユーザーとみなし、検出する。
ググるとLSTMとかがよく使われてる。

 


特異変換スペクトル(Singular Spectrum Transformation :SST)について

特異変換スペクトルは主に非定常データにおける異常部位の検出に使う(もちろん定常データにも使える)


特異変換スペクトルは、非定常データでもデータが大きく変化してる点(異常箇所)を検出できる。訓練データやテストデータはいらず、検出したいデータだけ用意すれば、numpyで簡単に実装できる
詳しい数式や論理とかは専門書や論文に任せてここではコードで実装みる。


まず検出したいデータは約半年分の検索ワードを2つ(AとB)を取り出して、特異変換スペクトルにかけてみた。

ワードA、Bの可視化

ワードA


ワードB


 
特異変換スペクトルのメイン実装コードは以下のようになる。ついでに結果を可視化した。

 

from sklearn.preprocessing import MinMaxScaler
 
mss=MinMaxScaler()
train1_frame=pd.DataFrame(mss.fit_transform(df_train1))
train2_frame=pd.DataFrame(mss.fit_transform(df_train2))
 
defembed(lst,dim):
emb=np.empty*1
emb=np.append(emb,tmp,axis=0)
returnemb
 
train1=np.array(train1_frame[0],dtype='float64')
 
w=168# 部分時系列の要素数
m=2# 類似度算出に使用するSVDの次元数
k=int(w/2)# SVD算出に用いる部分時系列数
L=int(k/2)# # 類似度を比較する2つの部分時系列群間のラグ
Tt=train1.size
anomaly_score=np.zeros(Tt)
 
# 異常値のスコアを算出するメソッド
score_list=[]
fortinrange(w+k,Tt-L+1+1):
tstart=int(t-w-k+1)
tend=t-1
# t以前の部分時系列群
X1=embed(train1[tstart:tend],w).T[::-1,:]
# 異常度算出対象の部分時系列群(test matrix)
X2=embed(train1[(tstart+L):(tend+L)],w).T[::-1,:]
 
# X1にSVDを行い上位m次元を抜き出す
U1,s1,V1=np.linalg.svd(X1,full_matrices=True)
U1=U1[:,0:m]
# X2にSVDを行い上位m次元を抜き出す
U2,s2,V2=np.linalg.svd(X2,full_matrices=True)
U2=U2[:,0:m]
 
# U1とU2の最大特異値を取得
U,s,V=np.linalg.svd(U1.T.dot(U2),full_matrices=True)
sig1=s[0]
 
# 最大特異値の2ノルムを1から引くことで異常値を得る
anomaly_score[t]=1-np.square(sig1)
score_list.append(anomaly_score[t])
 
 
score_=np.reshape(score_list,(len(score_list),))
 
# 変化度をmax1にするデータ整形
mx=np.max(score_)
score_=score_/mx
 
# trainデータの異常部位plot
train_for_plot=np.array(train1_frame[0],dtype='float64')
fig=plt.figure()
ax1=fig.add_subplot(111)
ax2=ax1.twinx()
 
p1,=ax1.plot(score_,'-b')
ax1.set_ylabel('degree of change')
ax1.set_ylim(0,1.2)
p2,=ax2.plot(train_for_plot,'-g')
ax2.set_ylabel('original')
ax2.set_ylim(0,12.0)
plt.title("Singular Spectrum Transformation")
ax1.legend([p1,p2],["degree of change","original"])
plt.show()

 

ワードA

ワードB


 
ご覧の通りAでは1カ所で急上昇し、Bは数カ所でデータが変化してる。いずれも非定常データだが上手く検出できてる。


数値も大きく変化してる点で

・Aが10→80

・Bが800→1800

まで上昇しているから、適切に異常個所を検出できてるっぽい。



特異変換スペクトルの活用方法は個人的に

・異常部位の検出
・訓練データやテストデータのを分ける

のに使えた。

 


閾値の設定について


反対に変化点検出のときに、異常値を判断する「閾値」の設定は特異変換スペクトルでは難しかった。
スペクトル変換しているので、リアルの値とかなり変わっている。そのため、閾値設定には工夫が必要になる。


非定常データの閾値設定には
・動的時間伸縮法(Dynamic time warping:DTW)

・LSTMのような回帰アルゴリズムとMSEの組み合わせ

とかの方が一般的だ。


動的時間伸縮法(DTW)は複数の時系列データ同士の距離(類似度)を計算し、最も高い類似度 or 最も低い類似度の時系列をk個選び、その平均値から閾値を求める。kが大きく設定すると閾値も低くなる。閾値の設定はケースバイケースで設定する必要がある。


DTWのコードはここを参照してほしい。


結論、特異変換スペクトルは定常/非定常に関わらず、どんなデータでも、訓練データとテストデータを分けたり、異常部位を見つけるのに役に立つ。かつ簡単に実装できるので、異常検知では色々な用途に使えるので、オススメの手法だ。