トポロジカルデータアナリシス：TDAパッケージを使ってみる

導入

とあることがきっかけで、とっても久しぶりにRでTDA（トポロジカルデータアナリシス）をしてみました。だいぶいろんなことを忘れていたので、単純な例を使ったメモを残しておきます。

トポロジカルデータアナリシスとは

とてもざっくりいうと、位相幾何学という数学の知見をつかって、データから「形」の情報を抽出するような手法になっています。

導入は、過去の記事にも載っていますので、こちらも参照していただけるとうれしいです。
tekenuko.hatenablog.com

「形」の情報の抽出のしかたですが、データ点のまわりにある半径の球をとり、その球たちの重なりによって点の間に線を引く、といった方法でデータ点から図形を見立てます。ただし、データに付随するノイズなどにロバストな形で図形的な情報を抜き出すために、半径を連続的に変化させた場合の構造の変化を見るといったテクニックを用います。

これらの簡単な解説は、一応過去の記事でも紹介しています（わかりづらかったらすいません）。
tekenuko.hatenablog.com

半径を連続的に変化させた場合、図形が生まれて、亡くなる…という動きが起こります。その一連を動きをプロットしたものが、パーシステント図と呼ばれます。今回、具体例を用いてパーシステント図がどうなるかを説明したいと思います。

使用パッケージ

ここでは、Rでトポロジカルデータアナリシスを行うためのパッケージである{TDA}の関数を使います。以下のようにしてインストールします。Rcppなど、いくつか依存パッケージがあり、それらのバージョンなどでエラーが出ることがあります。その場合は、依存パッケージを入れ直すなどするなどの対処をする必要があります（うろ覚え）。

# TDAパッケージ
if(!require(TDA)){
  install.packages("TDA", quiet = TRUE)
}
require(TDA)

使い方（マニュアル）はこちら：
https://cran.r-project.org/web/packages/TDA/vignettes/article.pdf

日本語の解説だと、ALBERTさんの記事が参考になります。
blog.albert2005.co.jp

サンプルデータ

大きな円と小さな円が繋がっている＋微小なノイズがのっているような人工データを生成します。ここでは、{TDA}パッケージにある、半径1の円上に沿ってランダムサンプリングを行う関数circleUnif()をうまく組み合わせます。

X1 <- circleUnif(100) * 2 + rnorm(100, 0, 0.03)
X2 <- circleUnif(50) * 0.3 + c(1.65,1.65) + + rnorm(50, 0, 0.02)
X <- rbind(X1, X2)
X <- data.frame(x = X[, 1], y = X[, 2])

可視化すると以下のようになります。

ggplot(X) + 
  geom_point(aes(x = x, y = y))

パーシステント図

パーシステント図を描きます。図形の構築法は、ヴィートリス・リップス複体を利用します。ヴィートリス・リップス複体を含めた構築法の概要は、過去記事を参照いただけると幸いです。
tekenuko.hatenablog.com

{TDA}パッケージでは、以下のような関数を使って実現します。

# パーシステント図の大元（パーシステントホモロジー）を計算
Diag_Cir <- ripsDiag(X = X, maxdimension = 1, maxscale = 5)

Xはデータ点、maxdimensionは何次元の三角形の移り変わりを見るかを指定するオプションです。今回は1としており、穴が空いているかの移り変わりを見ています。maxscaleは図を描くときの軸の最大値の値を指定するオプションです。パーシステント図を作るために動かす半径のパラメータは、本来無限大まで取れますが、変化は実質どこかで止まるはずで、その止まるパラメータ点をmaxscaleとしてスケーリングします。

可視化すると以下のようになります。{TDA}に内蔵されているplot関数を使うのが基本ですが、ここでは{ggplot2}で可視化を試みました。

# ggplot2で描くためにripsDiagの結果を整形
dimension_Cir <- Diag_Cir$diagram[,1]
birth_Cir <- Diag_Cir$diagram[,2]
death_Cir <- Diag_Cir$diagram[,3]
Diag_Cir_DF <- data.frame(cbind(Dimension = dimension_Cir, 
                              Birth = birth_Cir,
                              Death = death_Cir))

# plot
ggplot(Diag_Cir_DF) + 
  geom_point(aes(x = Birth, y = Death, colour = as.factor(Dimension)), size = 4) + 
  xlim(0,5) + ylim(0,5) + 
  geom_abline(slope = 1) + 
  labs(color = "Dimension") + 
  scale_color_manual(name="Dimension",
                     values=c("0"="black","1"="#F8766D"))

赤い点が穴に対応しており、対角線から距離が離れているほど半径の変化で安定的な穴であることを表しています。今回、赤い点は相対的に安定的なものが2点、不安定なものが1点出現しています。前者はデータ点を大きな円（対角線から大きく離れている点）と小さな円（対角線からちょっと離れている点）に対応しており、後者はノイズ起源でたまたま穴の空いた図形が生まれてすぐ消えていったものだと考えられます。このように、パーシステント図を利用すると、データ点の形に付随した量を抽出することができます。

Next Step

今回は、非常に簡単な例ではありますが、トポロジカルデータアナリシスのRでの実行例を紹介しました。このパッケージには他にも色々な機能が搭載されていますので、気になる方はマニュアルを見つつ色々お試しください。本ブログでも利用例を定期的に紹介していきたいと思っています（たぶん）。

Pythonでデータ分析：Prophetを使ってビットコインの予測（笑）をやってみる

導入

直近、これといって緊急の業務がなく、「自分の時間だ何勉強しようかなー」とPyStanとかをいじっていた矢先、「暇なら技術調査やってよ、Deep Learning的な何かとか」というお達しがきました。あいにく私は天邪鬼なので、２つ返事をして気になっていたけど触っていなかったProphetを調べることにしたのでした。

注：仕事はちゃんとしました（Seq2Seqの論文や書籍見て簡単な実装をしました）。

Prophet

Facebookが出した時系列予測のツールです。
facebook.github.io
すでに様々な方が紹介をしたり、Contributeしていたりするので、釈迦に説法感がありますが、このツールの良い点は、簡単に（分析の専門知識がなくても）ある程度それらしい予測値を出してくれるところです。ビジネス側でデータを活用したい場合や、分析者でもいったん簡単にデータから言えることを見てみる、といった場合に便利そうです。

というわけで、分析者としても抑えておきたいと前々から思っていたので、いい機会だと思って少し動かしてみました。一通りサンプルを見たあと、題材として何かあるかなと思って探していたところ、ビットコインのデータをPythonから取得できるライブラリを見つけたので、ビットコインを題材とすることにしました。

参考

• Giihub（メソッドの引数とか見るのに使った）

github.com

• ホクソエムさんの資料（最初に超見た）

d.hatena.ne.jp

• Pythonでビットコイン（US）のデータを取得するために参考

darden.hatenablog.com

データ取得

以下のコードを使います。
github.com
インストール方法はサイトに書かれています。pipが使えるなら

$ pip install https://github.com/s4w3d0ff/python-poloniex/archive/v0.4.6.zip

とすれば使用できる状態になるはずです。

データですが、以下のコードを実行すると、ビットコイン（US）のデータ（日次）を過去500日分取得してくれます。ざくっと紹介すると、polo.returnChartDataでデータを辞書型で取得しています。periodでどれくらいの間隔のデータかを指定（5分足なども指定できます）しています。ただ、時間がUNIX timeからの経過時間（秒）になっているので、後でよしなに変換しています。

# numpyやpandas
import numpy as np
import pandas as pd
from pandas import DataFrame, Series

# APIに関係するライブラリ
import poloniex
import time
import datetime

def getDataPoloniex():
    polo = poloniex.Poloniex()
    polo.timeout = 2
    chartUSDT_BTC = polo.returnChartData('USDT_BTC', period=polo.DAY, start=time.time() - polo.DAY * 500, end=time.time())
    tmpDate = [chartUSDT_BTC[i]['date'] for i in range(len(chartUSDT_BTC))]
    date = [datetime.datetime.fromtimestamp(tmpDate[i]).date() for i in range(len(tmpDate))]
    data = [float(chartUSDT_BTC[i]['open']) for i in range(len(chartUSDT_BTC))]
    return [date, data]

# データ取得
dat = getDataPoloniex()
# DataFrameへ格納
df = DataFrame(dat).T
df.columns = ['ds', 'y']

ちなみに、プロットすると以下のようになっています。

import matplotlib.pyplot as plt
%matplotlib inline
df.plot(x = df.columns[0], figsize = (15, 8))

途中で下落したりすれど、急成長しているような系列になっています。

Prophetによるモデル作成

ProphetはPythonの場合はpipで導入することができます。

pip install fbprophet

PyStanとか色々依存ライブラリがあるようなんですが、そういったものたちは先に入れとくといいのかもしれません。自分はすべて先に入っていたので比較ができませんが。

モデル構築はsklearnライクにできます。

mod = Prophet()
mod.fit(df)

非線形トレンドを入れたり、変化点を考慮したりといろいろオプションはあるのですが、ここではよしなにできる感を演出するのに何も指定しないでおきます。線形トレンド以上の成長をしていそうな系列なので、ホントはちゃんと処理したいですけどね。

予測に関しては、make_future_dataframeとpredictで新たに生成した系列（100日分）に予測値を格納します。

future_df = mod.make_future_dataframe(periods = 100)
forecast_df = mod.predict(future_df)

簡単なプロットはProphetじたいにplotというメソッドがあり、それでよしなな図が見れます。

mod.plot(forecast_df)
plt.show()

まあ、線形トレンド的な予測がされています。もっと長期まで予測しようと思ったら大まかな傾向（指数的）とは異なってくるとは思いますが、おおかなに上昇していきそうというのが捉えられている、としましょう。

もう少し細かい成分でプロットする機能もあります。

mod.plot_components(forecast_df)

週の上下している感を見ると

• 火曜日は上がり調子
• 水曜、木曜にかけて下がる

というホントかな、、という傾向が見られます。ちょっとこれは期間を変えて検証を繰り返す、とかしないと正確なことは何も言えません。ちなみに、他の記事では学習期間を変えると週の振る舞いは変わると言っているものもありますので、やはりこれから何かを主張するのは難しいのかもしれません。

qiita.com

感想

データがあるときに、チャチャッと分析をするのに便利です。これからは割りと使用するかもしれない、そう思われるツールでした。ただ、やはり適用限界はありそうなので、100%頼るというよりはできそうなところまで使ってみる、という付き合い方になるかなと思います。

Pythonでデータ分析：主成分分析（PCA）による異常検知

導入

データ分析の種類の一つとして、教師なし学習による異常検知というものがあります。ほとんどが正常なデータでまれに異常なデータが混じっている、その異常発生のパターンや異常と他の要因との紐付きがいまいちつかみきれていないというような場合、教師あり学習による2値分類がうまくワークしない、といった状況がありえます。そういった場合には、正常パターンを教師なし学習で学び、その正常パターンから外れているものを異常とする、という方法が有効です。

異常検知の方法ですが、正常パターンを正規分布でフィットするようなものから状態空間モデルに即したものまで多くのバリエーションがあります。その中の一つに、主成分分析（PCA）をベースにした異常検知の方法があります。正常データがある程度相関を持っているような多変数正規分布に従っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています（正規分布から外れている場合は、Box.Cox変換などで無理やり正規分布に近づける、他の手法を使うなどを検討します）。この方法ですが、Rでは実行例は井出さんの異常検知の本
www.coronasha.co.jp
にありますが、Pythonでの実行例というのが書籍やネットを探してもいまいち見つかりません。そのため、メモがてらこちらに簡単に残しておこうと思います。

多変数正規分布による人工データ

異常検知のデモ用に、人工データを用意します。PCAによる異常検知の場合、データは正規分布に従い、かついくつかの変数で相関が高いような場合にうまくはまります。というわけで、多変数正規分布で変数間にいくらか相関のあるようなデータを生成します。

# 必要なライブラリ
import numpy as np
import pandas as pd
from pandas import DataFrame, Series

# 多変数正規分布に従う乱数を生成
mean = (4, 2, 3)
sigma = [[1, 0.7, 0.75], [0.7, 1, 0.3], [0.75, 0.3, 1]]
df = np.random.multivariate_normal(mean, sigma, 500)

主成分分析（PCA）

※以下のコードは、無駄にDataFrameが使われていますが、もっとすっきり書ける可能性があります。

主成分分析は、次元圧縮の方法の一つです。データセットに含まれる変数が非常に多いがいくらか相関が強いものがあるとき、適切に基底の回転・取り替えを行うことによってより少ない軸でデータを近似することが期待できます。主成分分析（PCA）は今のデータセットから新しい基底（主成分）を求め、その基底のうちデータの情報を失わない（軸に射影しても分散が小さくならない）ものでいくらか代表して近似してしまおう、といった手法になっています。

PCAは、sklearnのメソッドで簡単に行うことができます。

# 異常度算出の都合上、正規化を行う
std = StandardScaler()
X = DataFrame(std.fit_transform(df))

from sklearn.decomposition import PCA
# インスタンス化
pca = PCA()
# 最初の1/5を学習に用いる
pca.fit(X[0:round(X.shape[0] / 5)])

pca.fitでレコード数×特徴量数の行列に関して特異値分解（もしくはこの行列と転置をかけた共分散行列というものの固有値分解）を行い、新しい基底を求めています。今回、変換前の軸は3つあったのですが、変換後は2つの軸でデータを近似することにします（こちらに関しては、本来は軸の寄与率などでどれくらい軸を採用するかを決めるなどします）。新しい基底はcomponents_というメソッドで引き出すことができ、その行方向を制限します。

Um = DataFrame(pca.components_[0:2, :]).T

ここで、後の都合上転置をしています。また、共分散行列というもの（）を以下で定義しておきます。

S2 = np.identity(X.shape[1]) - Um.dot(Um.T)

異常度算出

PCAで次元圧縮した際に用いたデータは、正常なデータのみ、もしくはほとんど異常なデータは存在しないとします。すると、正常データは近似した軸で貼られる部分空間上でうまく近似でき、逆に異常なデータはその部分空間から大きく離れるのではないかと期待されます。この部分空間は正常部分空間と呼ばれます。異常度を評価する一つの方法が、データと正常部分空間からの距離を用いる方法です。今、正規化（中心化でもよい）されたデータを（次元は特徴量の数）とすると、異常度は次のように書かれます。
$$a(\vec{x}) = \vec{x}^T S^2 \vec{x}^T$$
導出は井出さんの異常検知の本などを参照してください。

Pythonでは、おそらくライブラリで異常度を出力してくれるメソッドはないので、手で書く必要があります。ここでは、1レコードに対して異常度を計算する関数を作っておき、それを検証用データに対してループさせて使う、としています。

1レコードに対して異常度算出
def CalcAnomaly(x, S2):
    return x.T.dot(S2).dot(x)

# 検証データ（全データ）に対して異常度算出
a = [CalcAnomaly(X.iloc[i, :], S2) for i in range(X.shape[0])]   

検証・可視化

各レコードに対し、異常度を可視化してみます。

import matplotlib.pyplot as plt
%matplotlib inline

plt.figure(figsize = (15, 8))
plt.plot(a)
plt.xlabel('n_step')
plt.ylabel('anomaly score')

概ね1より小さい値になっています。まれに1程度の相対的に異常度が高い値が出力されていますが、こういった値は乱数の関係上発生しうるものだと思われます。本当に異常なデータがある場合は、今のスケール感でいうと数十といった値が出てくることもあるので、今回はほぼ正常データしかないだろうと思われるデータになっています。

また、もとの空間での散布図と異常度の関係も可視化しておきましょう。3軸なので、散布図も3次元でプロットしています。

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize = (12, 10))
ax = Axes3D(fig)
p = ax.scatter(df[:, 0], df[:, 1], df[:, 2], marker = 'o', c = a, cmap = 'jet', s = 50)
fig.colorbar(p)

新しい軸はデータの相関のある方向（図で言うと斜め方向）などになっているのですが、その軸たちから離れているデータほど異常度が大きくなっている、といった様子が図から読み取れます。このように、PCAによる異常検知は正常データがある程度相関を持っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています。

Next Step

今回は、主成分分析（PCA）をベースにした異常検知の方法があります。PCAによる異常検知は正常データがある程度相関を持っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています。異常検知の方法は他にも色々あるので、機会があればまた別の方法も紹介したいと思います。

Memo：Gluonの解説やコード紹介（海外）

導入

2017年10月12日（現地時間）に、MicrosoftとAWSがGluonというDeep Learningのライブラリを公開しました。

www.itmedia.co.jp

日本語の解説記事があまり見当たらなかったので、簡単なところは自分で試してみるなどし、いくつか記事にもしました。色々調べていたところ、Documentや国際会議などではすでにGluonじたいやコードの紹介があることがわかってきました。こういった情報を整理しておくのは後々の自分にとって有用だろうと思ったので、メモがてら残しておきます。

2017年の国際会議でのTalk

CVPR2017

mli.github.io

ICML2017, KDD2017

github.com

GithubのLecture Note

ICML2017, KDD2017での発表者のGithubにMXNetも含めたLecture Noteがあります。形式はJupyter Notebookです。MXNetのメソッドの使い方から、深層強化学習まで、広い範囲を扱っています。これを最初に見つけてたらどんなに楽だったか。
github.com

また、上記のLectureに対応するWeb Pageもあります。

• Deep Learning - The Straight Dope — The Straight Dope 0.1 documentation

所感

日本では、ついこの間登場した感が出てますが、海外の情報や国際会議のTalkではPreliminalなものも含め前々から紹介されていたりする、そんな当たり前のことを再確認しました。日本人だと日本語の記事を情報源にしたくなりますが、真面目に情報収集したいならソース元を探しにいかないと。論文のサーベイと同じですね。

とはいえGluonは海外も含めてもまだまだ情報が少ない状況なので、情報収集は定期的に続けていき、それを自分が紹介することで日本での紹介や実例がより活発になっていけばよいなと思います。

GluonでDeep Learning：CNNを組んでみる

導入

前回、MicrosoftとAWSが公開したライブラリであるGluonの紹介をしました。
tekenuko.hatenablog.com

前回紹介したのは、Tutorialの多層パーセプトロン（MLP）でしたが、Gluonは他のネットワークもサポートしています。今回は、畳み込みニューラルネットワーク（Convolutional Neural Network, CNN）の例を紹介しようと思います。

参考

参考記事などはまだ見当たらなかったので、四苦八苦してたのですが、MXNetのDocumentを探してみると、使用例がありました。
Convolutional Neural Networks in gluon — The Straight Dope 0.1 documentation
概ねこの使用例を踏襲して記述していけばよさそうです。

データセット：MNIST

前回と同じく、MNISTを使います。いったん、バッチサイズや出力層の次元数を変数と読み込んだデータ（MXNetのNDArray形式のデータ）を変換する関数を定義しておきます。nd.transpose()でうまく転置するのがポイントで、自分はこの処理をやってなくてエラーが出て対処に困りました（汗）。

import mxnet as mx
from mxnet import nd, autograd
from mxnet import gluon
import numpy as np

batch_size = 64
num_outputs = 10

def transform(data, label):
    return nd.transpose(data.astype(np.float32), (2,0,1))/255, label.astype(np.float32)

MNISTデータはMXNetのサイトからロードします。

train_data = mx.gluon.data.DataLoader(mx.gluon.data.vision.MNIST(train=True, transform=transform),
                                      batch_size, shuffle=True)
test_data = mx.gluon.data.DataLoader(mx.gluon.data.vision.MNIST(train=False, transform=transform),
                                     batch_size, shuffle=False)

ネットワーク構築

Gluonでは、Sequential()というメソッドを定義し、そこに積み木のように層を足していく感覚でネットワークのモデルを構築できます。今回は畳み込み層とプーリング層を2回挟み、その後全結合層を入れて最後に出力層、といった構成にしました。畳み込み層→プーリング層の組み合わせは、ある程度微小変換に対してロバストな特徴をデータから抽出するのによく用いられます。その後、全結合層を入れることで、ざっくりそういった抽出した特徴をどういった重みで組み合わせるのがよいかを表現しています。

# 全結合層の次元
num_fc = 512

net = gluon.nn.Sequential()
with net.name_scope():
    net.add(gluon.nn.Conv2D(channels=20, kernel_size=5, activation = 'relu'))
    net.add(gluon.nn.MaxPool2D(pool_size=2, strides=2))

    net.add(gluon.nn.Conv2D(channels=50, kernel_size=5, activation = 'relu'))
    net.add(gluon.nn.MaxPool2D(pool_size=2, strides=2))

    net.add(gluon.nn.Dense(num_fc, activation = 'relu'))

    # 出力層の次元は10
    net.add(gluon.nn.Dense(num_outputs))

畳み込み層、プーリング層を今回は入れています。どんな引数が使えるかは、以下のページを見るとわかります。
github.com

• Conv2D
  • 2次元の畳み込みを行う層です。フィルタの数（channel）やサイズ（kernel）が最低限指定する引数になります。
  • 他、画像データに対するフィルタのスライドのさせ方（strides）や、入出力の画像のサイズを変えないための処方（padding）のやり方などを指定できます。
    • ここはDocumentを見て、デフォルト以外を試したかったら変更するとよいです。
• MaxPool2D
  • 最大プーリングを行う層です。特定のWindowの値を最大値で代用し、集約します。
  • 何も指定しないとデフォルト（2×2のWindowで最大プーリング、ストライドはNoneなど）
    • pool_sizeは数値 or list/tupleで指定できます。今回はあえて数値で指定しています。

学習の設定

前回とほぼ同じ設定です。まず、CPUを使用することを宣言します。

ctx = mx.cpu()

初期値はXavierの初期値というものを今回は用いています。目的関数や最適化手法は変更していません。

net.collect_params().initialize(mx.init.Xavier(magnitude=2.24), ctx=ctx)
softmax_cross_entropy = gluon.loss.SoftmaxCrossEntropyLoss()
trainer = gluon.Trainer(net.collect_params(), 'sgd', {'learning_rate': .1})

初期化はDocumentがやっていたことを変えていないだけです（笑）。Xavierの初期化は、乱数で重みに初期値を与える際に、層のノードの数で重み付けを行う方法になっています。重みのばらつきをノードの数を考慮して均一化し、学習をうまく進めるよう工夫しているようです。解説に関しては、以下のようなページがあります。
qiita.com

学習

大まかな設定は前回と同じです。今回は、学習データとテストデータのAccuracyを出力できるようにしました。データとモデルを指定するとAccuracyを出力する関数を予め作っておきます。

def evaluate_accuracy(data_iterator, net):
    acc = mx.metric.Accuracy()
    for i, (data, label) in enumerate(data_iterator):
        data = data.as_in_context(ctx)
        label = label.as_in_context(ctx)
        output = net(data)
        predictions = nd.argmax(output, axis=1)
        acc.update(preds=predictions, labels=label)
    return acc.get()[1]

epoch数は10で学習を進めます。epochごとにlossと学習データ、テストデータのAccuracyを出力します（lossは移動平均も考えていたりしてますが、Documentに書いてあることをそのまま書いてあるだけです（笑））。

epochs = 10
smoothing_constant = .01

for e in range(epochs):
    for i, (data, label) in enumerate(train_data):
        data = data.as_in_context(ctx)
        label = label.as_in_context(ctx)
        with autograd.record():
            output = net(data)
            loss = softmax_cross_entropy(output, label)
        loss.backward()
        trainer.step(data.shape[0])

        ##########################
        #  損失関数の移動平均
        ##########################
        curr_loss = nd.mean(loss).asscalar()
        moving_loss = (curr_loss if ((i == 0) and (e == 0))
                       else (1 - smoothing_constant) * moving_loss + (smoothing_constant) * curr_loss)

    test_accuracy = evaluate_accuracy(test_data, net)
    train_accuracy = evaluate_accuracy(train_data, net)
    print("Epoch %s. Loss: %s, Train_acc %s, Test_acc %s" % (e, moving_loss, train_accuracy, test_accuracy))

# 出力
Epoch 0. Loss: 0.0838328978149, Train_acc 0.980033333333, Test_acc 0.9826
Epoch 1. Loss: 0.0540922844512, Train_acc 0.985, Test_acc 0.9856
Epoch 2. Loss: 0.0414136623197, Train_acc 0.991266666667, Test_acc 0.9897
Epoch 3. Loss: 0.0309304529114, Train_acc 0.989783333333, Test_acc 0.9868
Epoch 4. Loss: 0.0209868983092, Train_acc 0.995583333333, Test_acc 0.9917
Epoch 5. Loss: 0.0190144998394, Train_acc 0.988766666667, Test_acc 0.9848
Epoch 6. Loss: 0.0142381958556, Train_acc 0.997033333333, Test_acc 0.9915
Epoch 7. Loss: 0.0122593285998, Train_acc 0.997483333333, Test_acc 0.9912
Epoch 8. Loss: 0.0127596473961, Train_acc 0.997866666667, Test_acc 0.9915
Epoch 9. Loss: 0.00813980522647, Train_acc 0.9982, Test_acc 0.9913

ローカルPCだと計算に20分くらいかかります。ですが、学習、テストデータともにAccuracyが99%を超えています。感覚的には、中間層1層のMLPだとAccuracyが90%前半程度になので、CNNにしたことで大幅な精度向上になっていると期待できます。

別のネットワーク例

自作でCNNを作っていたとき、どうしてもエラーが解消できなくてネットサーフィンしてたときに見つけた例を紹介しておきます。

num_fc = 512
net = gluon.nn.Sequential()
with net.name_scope():
    net.add(gluon.nn.Conv2D(channels=20, kernel_size=5))
    net.add(gluon.nn.BatchNorm(axis=1, center=True, scale=True))
    net.add(gluon.nn.Activation(activation='relu'))
    net.add(gluon.nn.MaxPool2D(pool_size=2, strides=2))

    net.add(gluon.nn.Conv2D(channels=50, kernel_size=5))
    net.add(gluon.nn.BatchNorm(axis=1, center=True, scale=True))
    net.add(gluon.nn.Activation(activation='relu'))
    net.add(gluon.nn.MaxPool2D(pool_size=2, strides=2))

    net.add(gluon.nn.Flatten())

    net.add(gluon.nn.Dense(num_fc))
    net.add(gluon.nn.BatchNorm(axis=1, center=True, scale=True))
    net.add(gluon.nn.Activation(activation='relu'))

    net.add(gluon.nn.Dense(num_outputs))

途中でBatchNormalization（ざっくり学習を効率的にすすめるための工夫です）を施したりしています。上で自分が作ったネットワークは、スタンダードなCNNの要素に絞りたかったので、細かい工夫はカットしてネットワークを構築しました。

このネットワークを学習させると、Accuracyは以下のようになります。

# 学習と検証結果
Epoch 0. Loss: 0.0459737773102, Train_acc 0.992066666667, Test_acc 0.9892
Epoch 1. Loss: 0.0292876055092, Train_acc 0.992816666667, Test_acc 0.9882
Epoch 2. Loss: 0.0204365993603, Train_acc 0.995516666667, Test_acc 0.9908
Epoch 3. Loss: 0.015462121763, Train_acc 0.998183333333, Test_acc 0.9926
Epoch 4. Loss: 0.0106946259829, Train_acc 0.999033333333, Test_acc 0.9926
Epoch 5. Loss: 0.00973990939113, Train_acc 0.999466666667, Test_acc 0.9926
Epoch 6. Loss: 0.00662801339947, Train_acc 0.999466666667, Test_acc 0.993
Epoch 7. Loss: 0.00540244358498, Train_acc 0.9998, Test_acc 0.9927
Epoch 8. Loss: 0.00381258537248, Train_acc 0.999916666667, Test_acc 0.993
Epoch 9. Loss: 0.00272173117527, Train_acc 0.99995, Test_acc 0.9934

工夫しているだけあって、普通にCNN組むよりも性能が上がるようですね。

Next Step

今回は、Gluonで畳み込みニューラルネットワーク（Convolutional Neural Network, CNN）の例を紹介しました。ネットワークの構築は、CNNの場合も同様に非常に簡単に行うことができます。次回以降は、異なるデータを使う、CNN以外のネットワークを試す、などを検討しようと思います。