Keras で少ない画像データから学習したい

深層学習で画像分類をしてみたいときとかあると思います。でも画像を用意するのが面倒です。すると Keras ブログに「少ないデータから強力な画像分類」とあります。GitHubスクリプトもあります。これを使ってみます。
Building powerful image classification models using very little data
fchollet/classifier_from_little_data_script_1.py - GitHub

データの準備

Google 検索で出てきた順にカレーとラーメンの画像を保存しました。でも時間がなかったので訓練用とテスト用に各クラス10枚ずつしか保存しませんでした。Keras ブログに書いてあるように、「少ないデータ」とは少ないといえども "just a few hundred or thousand pictures" ということなので、真面目にやる場合はきちんと用意してください。これらの画像は Keras ブログの指示通り、data/train/curry/ のようにフォルダを切った下に置きます。

訓練

このカレーとラーメンの画像分類を、畳み込み層とプーリング層を積み上げたニューラルネットワークで訓練するわけですが、少ないデータを最大限に活用して学ぶために、画像を常にランダムに少しの回転や平行移動、シア変形(平行四辺形のように歪める変形)、拡大縮小、左右反転などをさせながらモデルに入力します(どれだけの変形まで許容するかは自分で指定します)。こうすることによって過学習が抑制されモデルの汎化性能が上がるそうです。

結果

訓練データが20枚というやる気のなさでもテストデータの9割を正解してくれました(下表)。何回か訓練を試行するとテスト正解率が100%になることも多かったです。

f:id:cookie-box:20170810000634p:plain:w480
なお「カレー」と誤判断されたラーメンは以下のサイトの「天日地鶏」というお店のもので、ニューラルネットはこのチャーシューがカレールーに、麺とネギがライスに見えたのかもしれません(?)。
キングオブ静岡ラーメン~人気ランキングTOP10|静岡新聞SBS-アットエス
「ラーメン」と誤判断されたカレーは以下でした。何がラーメン要素だったのかよくわかりません。
​新潟発。衝撃200円カレーが東京に初進出 - エキサイトニュース(1/2)
逆に以下のあまりカレーに見えない MOKUBAZA のチーズキーマカレーはカレーに判定されました。ラーメンにも見えませんが。
【保存版】2015年のカレーまとめ / もっとも印象に残ったお店20選 | ロケットニュース24

スクリプト

ほぼ GitHubスクリプトの通りですが、サンプル数を合わせてあるのと、最後に個別データの確率を csv 出力する部分を追加してあります。

from keras.preprocessing.image import ImageDataGenerator
from keras.models import Sequential
from keras.layers import Conv2D, MaxPooling2D
from keras.layers import Activation, Dropout, Flatten, Dense
from keras import backend as K
import pandas

if __name__ == '__main__':
  img_width, img_height = 150, 150 # 訓練時の画像サイズ
  
  train_data_dir = 'data/train'
  validation_data_dir = 'data/validation'
  nb_train_samples = 20
  nb_validation_samples = 20
  epochs = 50
  batch_size = 5
  
  if K.image_data_format() == 'channels_first':
    input_shape = (3, img_width, img_height)
  else:
    input_shape = (img_width, img_height, 3)
  
  # モデル構築
  model = Sequential()
  model.add(Conv2D(32, (3, 3), activation='relu', input_shape=input_shape))
  model.add(MaxPooling2D(pool_size=(2, 2)))
  model.add(Conv2D(32, (3, 3), activation='relu'))
  model.add(MaxPooling2D(pool_size=(2, 2)))
  model.add(Conv2D(64, (3, 3), activation='relu'))
  model.add(MaxPooling2D(pool_size=(2, 2)))
  model.add(Flatten())
  model.add(Dense(64, activation='relu'))
  model.add(Dropout(0.5))
  model.add(Dense(1, activation='sigmoid'))
  
  model.compile(loss='binary_crossentropy', optimizer='rmsprop', metrics=['accuracy'])

  # 訓練データのプレ処理の設定: RGB値のスケーリングに加え、シア変形や拡大縮小によって訓練増強
  train_datagen = ImageDataGenerator(rescale=1. / 255, shear_range=0.2, zoom_range=0.2, horizontal_flip=True)
  # テストデータのプレ処理の設定: こちらはRGB値のスケーリングのみ
  test_datagen = ImageDataGenerator(rescale=1. / 255)
  
  # 訓練データをディレクトリ内のデータから随時作成するジェネレータ
  train_generator = train_datagen.flow_from_directory(train_data_dir,
    target_size=(img_width, img_height), batch_size=batch_size, class_mode='binary')
  # テストデータ作成をディレクトリ内のデータから随時作成するジェネレータ
  validation_generator = test_datagen.flow_from_directory(validation_data_dir,
    target_size=(img_width, img_height), batch_size=batch_size, class_mode='binary')
  
  # 訓練
  model.fit_generator(train_generator, steps_per_epoch=nb_train_samples//batch_size,
    epochs=epochs, validation_data=validation_generator, validation_steps=nb_validation_samples//batch_size)
  # 訓練済重みの保存
  model.save_weights('weight.h5')
  
  # 具体的に個々のデータに対してどのような確率分布になったかの出力
  train_generator = test_datagen.flow_from_directory(train_data_dir,
    target_size=(img_width, img_height), batch_size=1, class_mode='binary', shuffle=False)
  validation_generator = test_datagen.flow_from_directory(validation_data_dir,
    target_size=(img_width, img_height), batch_size=1, class_mode='binary', shuffle=False)
  
  predict = model.predict_generator(train_generator, steps=nb_train_samples)
  prob = pandas.DataFrame(predict, columns = ['curry'])
  prob['ramen'] = 1.0 - prob['curry']
  prob['predict'] = prob.idxmax(axis=1)
  prob.to_csv('train.csv', index=False, encoding='utf-8')
  
  predict = model.predict_generator(validation_generator, steps=nb_validation_samples)
  prob = pandas.DataFrame(predict, columns = ['curry'])
  prob['ramen'] = 1.0 - prob['curry']
  prob['predict'] = prob.idxmax(axis=1)
  prob.to_csv('test.csv', index=False, encoding='utf-8')

雑記: 最近のふぁぼより

後で読もうとふぁぼって結局読まないこととかあると思います。なので最近どんなことをふぁぼったか復習してみたんですが結局まだ読んでいません。

「A Distributional Perspective on Reinforcement Learning

 [1707.06887] A Distributional Perspective on Reinforcement Learning

強化学習において Bellman 最適方程式による行動価値の更新のとき、期待値に更新するのではなく分布のまま扱おうという話のようです。期待値で更新する方法でもリスク回避的な学習にはできるけどって言っていると思います。行動価値を分布として保持するのはイメージが湧くような気がしないでもないですがそれでどう方策を決めるのかはまだ全然読んでいないのでわかりません。確率的に行動するんだろうと思います。

「Reducing Dimensionality from Dimensionality Reduction Techniques」

 Reducing Dimensionality from Dimensionality Reduction Techniques

次元削減手法(PCA、t-SNE、自己符号化器)に関する解説のようです。

「深層学習とベイズ統計」

 深層学習とベイズ統計

@yutakashino 様のスライドです。決定論的モデルであるニューラルネットを確率論的モデルに modify したいため、ニューラルネットの各重みパラメータ  \theta を「ある分布からの実現値」として取り扱おうという話だと思います。
_人人人人人人人人人人人人人人人人_
> 学習の収束のふるまいがキモい < (スライド13頁)
  ̄Y^Y^Y^Y^Y^Y^Y^Y^Y^Y^Y^Y^Y^Y ̄

「Deep Reinforcement Learning, Decision Making, and Control」

 https://sites.google.com/view/icml17deeprl

深層強化学習のチュートリアルのようです。ICML2017 で講演があったのでしょうか。

「物理のいらない量子アニーリング入門」

 物理のいらない量子アニーリング入門 - Platinum Data Blog by BrainPad

量子アニーリング入門ということです。github でコードが公開されているということです。
GitHub - ohtaman/anneal

「Fundamentals of Nonparametric Bayesian Inference (Cambridge Series in Statistical and Probabilistic Mathematics)」

Fundamentals of Nonparametric Bayesian Inference (Cambridge Series in Statistical and Probabilistic Mathematics)Fundamentals of Nonparametric Bayesian Inference (Cambridge Series in Statistical and Probabilistic Mathematics)
Subhashis Ghosal Aad van der Vaart

Cambridge University Press 2017-06-26
売り上げランキング : 2101

Amazonで詳しく見る
by G-Tools

比較的最近出た本らしいです。670ページらしいです。購入したいですが置く場所がありません。

「On the State of the Art of Evaluation in Neural Language Models」

 [1707.05589] On the State of the Art of Evaluation in Neural Language Models

自然言語処理には結局 LSTM がよいという話のようです。なんかもう state-of-the-art がタイトルになってしまいました。

「~簡単!炎上のさせ方講座~」

@jagarikin 様のツイートより。炎上のさせ方が簡単にわかります。
ただ臨場感あふれる炎上にするにはコツが要るらしく、適当な実装では上手くいきませんでした(下図)。なんか地面が全体的にのっぺり燃えているようになってしまい、点々と炎を吹いているようにならなかったです。炎上のさせ方は奥が深いようです。また頑張ります。

f:id:cookie-box:20170808113251g:plain

あと PyTorch とか最近話題ですね。

岩波データサイエンス Vol.6: ノート1

6月下旬に以下の本が出ましたので、この本を読んで時系列解析の復習ポイントだと思ったことと、本の記述に対して自分で考えたことを書きます。誤っている可能性があります。
岩波データサイエンス Vol.6岩波データサイエンス Vol.6
岩波データサイエンス刊行委員会

岩波書店 2017-06-23
売り上げランキング : 944

Amazonで詳しく見る
by G-Tools

  • そもそもなぜ時系列解析には状態空間モデルなのかというと、時系列を予測したい場合は、本質的に、「観測できない真のダイナミクス」と「ノイズ」を分離する必要があるからだよという話。ノイズまで学習してしまったら予測にとっては害だから(「しばしば悲惨な結果をもたらす」5ページ)。
  • 直接観測できない真の姿を状態という。状態を観測するときにノイズが加わるし(観測ノイズ)、その状態が時間発展するときにもノイズが加わる(システムノイズ)(7ページ)。
    • 線形・ガウス型状態空間モデルとは発生する観測ノイズ  v_n とシステムノイズ w_n が平均ゼロのガウシアン(分散は時変でよい)、一期前の状態を現在の状態にうつす写像 F_n が線形(時変でよい)、システムノイズが状態に加わるときの変換  G_n も線形(時変でよい)、観測も線形 H_n(時変でよい;観測ノイズは平均ゼロのガウシアンのまま観測に加わるとする)の状態空間モデルを特にそういう(7~8ページ)。
       x_n = F_n(x_{n-1})+G_n(v_n)
      y_n = H_n(x_n) + w_n
      この式をみるとなんでシステムモデルの方のノイズだけ  G_n が付いてるんだと思うけど、平均ゼロのガウシアンを G_n で線形変換するというのは平均ゼロとは限らないガウシアンにするということで、例えば海の上で漂流していて陸地までの距離を知りたいとき(下図;2015-12-29の日記からリサイクル)沖の方に流される風が吹いていたら「大気によるノイズ」はゼロでないドリフト成分をもつだろう。「それって状態の時間発展 F_n の方に入らないの?」って思うかもしれないけど、ケースバイケースだと思う。状態空間モデルを描くとき現実のシステムを思い浮かべて描くわけで、「うーんこの要素はシステムの範疇」だったら  F_n に含めるし、システム外と考えるのだったら  G_n に含めると思う。他方観測ノイズはドリフトしないのは「観測ノイズがドリフトすんな」という要請だと思う。仮に何らかのドリフト要素があるなら  H_n に含めましょうということなのだと思います。
      f:id:cookie-box:20151229173717p:plain:w400
  • 「パラメータ推定に使われる尤度も、予測値だけでは計算できない(11ページ)」
    • 観測データ  y_0 を説明するモデルを  p(y | \theta) と考えているとする。パラメータ  \theta の尤度関数は  f(\theta) = p(y_0 | \theta) となる。もし  p(x | \theta)デルタ関数だったら(分布を考えず点のみ予測するのだったら)  p(y_0 | \hat{\theta})=1 を与える  \hat{\theta} を除いて尤度がゼロになってしまいパラメータが尤もらしいかどうかの比較ができない。尤度関数が微分できないので最尤法もしづらい。たぶん。なお、時系列の場合の尤度関数は  f(\theta) = p(y_{0:n} | \theta) = \prod_{i=0}^{n} p(y_{i} | \theta)
  • 季節調整法で、トレンド成分、季節成分に加えて定常AR成分を導入すると予測がよくなることがある(20ページ)。
    • 定常過程とは: 弱定常過程 - クッキーの日記
      • 時系列から任意の10時点の値を取ってきたときにそれらがしたがう同時確率分布と、さらにそれぞれの30ステップ後の10時点の値を取ってきたときにそれらがしたがう同時確率分布が同じ。時系列をプロットしていったとき、プロットはずっと水平な帯の中におさまる。
    • 「定常モデルと非定常モデルでは長期予測においてまったく働きが異なる(20ページ)」: ここが何を言いたいのかと思ったんだけど、次の文章と合わせると、定常モデルは「繰り返されるパターンが何か」を表現しさえすればよいけど、非定常モデルは「繰り返されるパターン」「あと時系列全体の成長トレンド」を同時に表現することになるので働きが異なるということなんだろうと思います。
    • ここでいう「トレンドが波打っている」というのは、季節調整で排除できない、「商品が売れた月の翌月は一時的に売れ行きが鈍る」のような短期変動なのだろうと思います。
  • 拡張カルマンフィルタ(言葉と参考文献のみ;20ページ):
  • 混合ガウス分布近似で「正規分布の項数が時間の進行とともに爆発的に増大する(21ページ)」:
    • 何が爆発するのかわからないんですが、毎ステップか数ステップ毎に何か情報量基準を最小化するように混合正規分布をその個数も含めてフィッティングし直すということなんでしょうか。
  • R の dlm はつかったことがあります。

雑記

Anaconda 環境を利用した Windows への Keras (TensorFlow backend) 導入

以下の記事の手順でできました。
qiita.com今回は Keras を利用するので Keras も導入しました。Anaconda インストール後、コマンドプロンプトを立ち上げて以下のように「仮想環境作成 → 仮想環境立ち上げ → TensorFlow導入 → Keras導入」しました。仮想環境を作成するとき Python のバージョンを指定できますが、3.5 にしないといけないようです。

C:\work> conda create -n tensorflow python=3.5 anaconda
C:\work> activate tensorflow
(tensorflow) C:\work> pip install --ignore-installed --upgrade https://storage.googleapis.com/tensorflow/windows/cpu/tensorflow-1.0.1-cp35-cp35m-win_amd64.whl
(tensorflow) C:\work> pip install keras

これで Keras (TensorFlow backend) がつかえるようになったことを確認するために以下のスクリプトを実行してみます。中間層が2層の全結合ニューラルネットワークで iris を分類します。

# -*- coding: utf-8 -*-
import numpy
import random
from keras.models import Sequential
from keras.layers.core import Dense, Activation
from keras.utils import np_utils
from sklearn import datasets

if __name__ == '__main__':
  # irisデータのロード
  iris = datasets.load_iris()
  x = iris.data
  y = np_utils.to_categorical(iris.target, num_classes=3)

  # シャッフルして訓練データとテストデータに分割
  n = iris.target.shape[0]
  indices = random.sample(range(n), n)
  x = x[indices]
  y = y[indices]
  n_train = int(n / 2)
  x_train = x[range(n_train)]
  y_train = y[range(n_train)]
  x_test = x[range(n_train, n)]
  y_test = y[range(n_train, n)]

  # モデルの構築と学習
  model = Sequential()
  model.add(Dense(10, input_shape=(x.shape[1],), activation='relu'))
  model.add(Dense(10, activation='relu'))
  model.add(Dense(3, activation='softmax'))
  model.compile(optimizer='adam', loss='categorical_crossentropy', metrics=['accuracy'])
  model.fit(x_train, y_train, epochs=100, batch_size=10, shuffle=True)

  # テスト結果表示
  scores = model.evaluate(x_train, y_train, verbose=0)
  print('Accuracy (train): %.2f%%' % (scores[1]*100))
  scores = model.evaluate(x_test,  y_test,  verbose=0)
  print('Accuracy (test) : %.2f%%' % (scores[1]*100))

以下のように学習できました。

(tensorflow) C:\work> python test.py
...
Epoch 99/100
75/75 [==============================] - 0s - loss: 0.3137 - acc: 0.9733
Epoch 100/100
75/75 [==============================] - 0s - loss: 0.3072 - acc: 0.9733
Accuracy (train): 97.33%
Accuracy (test) : 96.00%
さらに Anaconda 環境への OpenCV の導入

さらに画像処理ライブラリ OpenCV も導入します。
Anaconda で対象の仮想環境を立ち上げた状態で、以下のコマンドで導入できます。

(tensorflow) C:\work> conda install -c menpo opencv3

導入が上手くいったかどうか、以下のスクリプトで確認します。

# -*- coding: utf-8 -*-
import cv2
import numpy

if __name__ == '__main__':
  image = numpy.zeros((256, 256, 3), numpy.uint8)
  for i in range(256):
    for j in range(256):
      image[i, j, 0] = i
      image[i, j, 1] = j
  cv2.imwrite('hoge.png', image)

以下のような画像ファイル hoge.png が生成されます。

f:id:cookie-box:20170806220129p:plain

OpenCV で動画作成

OpenCV で動画を作成することができます。以下のスクリプトで movie.avi が生成されます。以下では各フレームを画素の行列から生成していますが、画像ファイルを読み込んでつなげて動画にすることもできます。

# -*- coding: utf-8 -*-
import cv2
import numpy

if __name__ == '__main__':
  vvw = cv2.VideoWriter('movie.avi', cv2.VideoWriter_fourcc('X','V','I','D'), 10, (200, 200))
  for i in range(10):
    for j in range(10):
      j_ = (9 - j) if (i % 2 == 1) else j
      image = numpy.zeros((10, 10, 3), numpy.uint8)
      image[i, j_, 0] = i  * 25
      image[i, j_, 1] = j_ * 25
      image = cv2.resize(image, (200, 200), interpolation=cv2.INTER_AREA)
      vvw.write(image)

以下のような動画ができます。.avi ははてなブログに貼り付けられないようなので、以下のページで .gif に変換しました。
AVI GIF 変換。オンライン フリー — Convertio

f:id:cookie-box:20170806222204g:plain

なお、cv2.VideoWriter の第2引数の4文字はコーデック指定です( Video Codecs by FOURCC - fourcc.org )。最適かどうかはわかりませんがとりあえず ('X','V','I','D') で生成できました。第3引数は1秒あたりのフレーム数、第4引数は画像サイズです。

雑記

データを分析するとき、取り込んだままのデータの列をそのままつかうのではなくて、元々あった列をもとに新しい列をつくることとかあると思います。でも Python とかよくわからないので Google 検索してどうやって加工するか調べると思います。

例えば、元々あった列のうち a, b, c 列の中央値をとり、新しい列 d に追加したいと思います。
以下の (1) (2) はどちらも望みを満たすようであることがわかります。でも (1) は遅いのでやってはいけないとわかります。
というか (2) のように一気にできるとわかれば (1) のように行ごとの処理にはしないと思いますが調べ方が難しいと思います。

# -*- coding: utf-8 -*-
import pandas
import numpy

if __name__ == "__main__":
  data = pandas.DataFrame({'name': numpy.repeat('Alice', 100000),
                           'a': range(0, 100000, 1),
                           'b': range(0, 200000, 2),
                           'c': range(0, 400000, 4)})

  # (1) 44.674 sec
  data['d'] = data.apply(lambda x: numpy.median(x[['a', 'b', 'c']]), axis=1) 
  # (2) 0.010 sec
  data['d'] = numpy.median(data[['a', 'b', 'c']].values, axis=1)

ともあれ、numpy.mean や numpy.median などで各行への処理を一気にできるのだったらそうすればいいことはわかります。
ただ、もっと自分で定義した処理をしたいときはベクトル演算が定義されているとは限らないので困ると思います。
そのような場合は numpy.vectorize で関数をベクトル化することにより処理速度を改善できると以下の記事にあります。
sinhrks.hatenablog.com

例えば、データに元々ある time 列をもとに、米国サマータイム期間かどうかを is_dst 列として付加したいとします。
apply をつかうと (1) のようにかけます。でも上の記事のアドバイス通りこの処理をベクトル化した (2) の方が速いです。

# -*- coding: utf-8 -*-
import pandas
import numpy

def is_dst(dt):
  return(pandas.Timestamp(dt).tz_localize('Asia/Tokyo').tz_convert('US/Eastern').dst().seconds > 0)

if __name__ == "__main__":
  data = pandas.DataFrame({
    'name': numpy.repeat('Alice', 100000),
    'time': numpy.repeat(['2017/03/12 15:59', '2017/03/12 16:00'], [50000, 50000])})
  data.time = pandas.to_datetime(data.time)

  # (1) 28.834 sec
  data['is_dst'] = data.apply(lambda x:
    x.time.tz_localize('Asia/Tokyo').tz_convert('US/Eastern').dst().seconds > 0, axis=1)
  # (2) 9.768 sec
  data['is_dst'] = numpy.vectorize(is_dst)(data.time)

ところで、上の実装で自作の is_dst 関数は受け取った日付時刻を pandas.Timestamp 型に変換していますが、time 列は元々 pandas.Timestamp 型です。なのになぜこの変換をしているのかというと、ベクトル化した関数を time 列に適用したときに time 列が勝手に numpy.datetime64 型になってしまったからです。なんでそうなるのかは Python とかよくわからないのでわかりません。