Think more, try less

Owen Zhangの言葉を大切にして,日々Kaggleに取り組んでいます.

データサイエンティスト養成読本 登竜門編(11章 機械学習超入門)を執筆しました。

データサイエンティスト養成読本 登竜門編がついに3/25に発売されます(一部の大型書店にてすでに先行販売されています)。これまでのデータサイエンティスト養成読本シリーズでは内容が難しかったという声が多かったため、ビギナー向けにやさしい本を提供するために、この登竜門編が刊行されることとなりました。本書は共著という形で執筆しており、さまざまな分野の現役の若手のデータサイエンティストが執筆者として参画しています。

データサイエンティスト養成読本 登竜門編 (Software Design plus)

データサイエンティスト養成読本 登竜門編 (Software Design plus)

対象読者

本書の対象読者としては、データサイエンスのビギナー向けとなっています。これからデータサイエンスの世界に入っていこうとする学生さんや、これから実務でデータサイエンスを活用していこうと検討している社会人などを想定しています。また、データサイエンスを活用している企業の新人やあまりデータサイエンスに詳しくない担当者向けの教育に利用できる内容となっているため、データサイエンスのプロにも有効活用してもらえる内容になっています。

データサイエンティスト養成読本 登竜門編の読みどころ

これまでのデータサイエンティスト養成読本と比較すると、各章が、前後の繫がりを意識した構成になっており、非常に読みやすい内容となっています。また、各章の著者が得意な分野を執筆したという内容ではなく*1、全体のバランスを意識して章を構成したため、無駄がなくかつ中途半端な内容となっていないことも特徴です。きっと本書を通読すれば、データサイエンスに関する基礎が一通り身につくと思います。

また、個人的にお薦めしたい登竜門編の読みどころは、「10章 さまざまなデータの理解と表現(水上ひろき)」です。おそらく現存する一般化線形モデルの解説のなかで最もやさしく、すっきり書かれている内容だと思います。統計モデリングの初歩の初歩として本章を読んだ後、みどり本などにステップアップしていくと良いと思います。

11章 機械学習超入門の読みどころ

私は「11章 機械学習超入門」を担当しました。これだけ世の中に機械学習の入門書が溢れているなかで、入門的な内容を執筆することは苦労しました。しかし、さまざまな機械学習の解説書籍を参考にしましたが、どれも実務で機械学習を実践する上では、何か部分的に欠けている解説が多いなと感じる部分があったため、私が執筆するこの機械学習の入門記事は、とにかく必要な基礎基本はすべて網羅することを意識しました。本章で、意外にもこんなことが載ってて嬉しいかなと個人的に思っている所を、以下にあげます。

  • 汎用人工知能と特化型人工知能と違いが解説されている
  • 機械学習モデリングにおける基本用語が解説されている
  • ハイパーパラメータの探索方法として、ランダムサーチも載っている
  • 評価尺度の解説で、AUC*2対数損失についても触れられている
  • ジニ不純度の計算方法が解説されている
  • 過学習の直観的イメージを解説している
  • mlrパッケージについて解説している

特に最後の、mlrパッケージについては、これまでのRによる機械学習のプログラミングよりも非常に理解しやすいパッケージなので、ぜひ使用してもらいたいです。他の書籍では、mlrパッケージについてはあまり解説がないため、本章の良い点と思っています。まだ本書を手にしていない方は、mlrパッケージの使い方について、TokyoRの勉強会で発表した資料があるので、そちらを参照してイメージを掴んでみてください。

また、「11-3 Rで機械学習を試してみよう」のソースコードのリンクも掲載しておきます。

データサイエンティスト養成読本 登竜門編 「11-3 Rで機械学習を試してみよう」のソースコード · GitHub

ちなみに機械学習といえば、Pythonのscikit-learnが有名であり、本来的には、scikit-learnで解説するべきと思われるかもしれませんが、scikit-learnの使い方については、改訂2版 データサイエンティスト養成読本の「第3章 Pythonによる機械学習(早川敦士(@gepuro))」に詳しく基礎基本が解説されているため、こちらを参照される方が良いと思います。この改訂2版の機械学習の解説との重複を避けるためにも、Rによる機械学習のコーディング方法を解説した次第です。

改訂2版 データサイエンティスト養成読本 [プロになるためのデータ分析力が身につく! ] (Software Design plus)

改訂2版 データサイエンティスト養成読本 [プロになるためのデータ分析力が身につく! ] (Software Design plus)

11章 機械学習超入門では載せられなかった内容

完全に網羅することを意識しましたが、それでも紙面の都合上、載せられなかった内容がいくつかあります。

確率的勾配降下法(Stochastic Gradient Descent:SGD

これは、水上さんとも話して、10章、11章で載せられなかったことを悔やみました。確率的勾配降下法とは、多くの機械学習アルゴリズムにおける誤差関数の最適化に利用される方法で、10章だとロジスティック回帰において最尤推定するときに用いられます。11章では、これは登竜門編としては、難しい内容となると判断して、載せないことにしました。確率的勾配降下法については、下記の本が参考になると思います。

ITエンジニアのための機械学習理論入門

ITエンジニアのための機械学習理論入門

深層学習(Deep Learning)

現在、これだけ注目されている深層学習(Deep Learning)ですが、こちらもデータサイエンスのビギナーには必要ないと判断して、載せないこととしました。深層学習が必要になるのは、現状では、これを活用してビジネスを検討している研究機関や一部の先進的なWeb系企業の方が多いかなと思います。初心者でもすぐにできるような手書き文字認識のようなタスクを理解したところで何も活用できません。それでも理解したいということであれば、下記の本が参考になると思います。

深層学習 (機械学習プロフェッショナルシリーズ)

深層学習 (機械学習プロフェッショナルシリーズ)

サポートページ

一部、最後の入稿間際でバタバタしていたこともあり、初版については、第11章のP.222の右下段(図8 ROC曲線の例)が誤っています。正しくは、サポートページの図をご参照ください。大変申し訳ありませんでした。

サポートページ:データサイエンティスト養成読本 登竜門編:|技術評論社

おわりに

2016年4月にキックオフで企画がスタートして、1年近くずっと執筆してきましたが、技術評論社のデータサイエンティスト養成読本 登竜門編のご担当の高屋様、また共著で執筆されていた皆様には大変ご迷惑をおかけしました。執筆中は、内容構成などでかなり迷走して、本当に入稿間際までどうなることやらという感じでしたが、なんとか形になってほっとしています。辛抱強くご支援頂きましてありがとうございました。そしてこの本書を執筆するにあたり、R勉強会やKaggle Tokyo Meetupに参加されている多くの方々から日頃学んでいることが、このような書籍の執筆に大きく影響していると感じています。日頃お世話になっている皆様にも心より御礼申し上げます。最後に、本書が、これからデータサイエンスを始めようとする多くの入門者に届き、最初の登竜門として意義のある書籍になることを願います。

*1:そのため、執筆者側がかなり苦労しました。

*2:初版では、図の差し替えが間に合わず、P.222 第11章 右下段(図8 ROC曲線の例)に間違った図が挿入されています。大変申し訳ありません。正しい図については、サポートページをご参照ください。

Python機械学習プログラミング/福島真太朗(監訳)のご紹介

これは技術書献本大感謝 Advent Calendar 2016の6日目の記事です。

qiita.com

福島さんよりPython機械学習プログラミング/福島真太朗(監訳)をご恵贈頂きました。2016/06/30に発売され、すでに多くの方に知れ渡っているかと思いますが、この場を借りて書評という形で宣伝させていただきます。

対象読者

これまでに機械学習について学んで来なかったエンジニアで、一通り機械学習の基礎を身につけたい方、もしくは、基本的な機械学習は習得しているが、Pythonによる機械学習の実践方法を学びたい方に非常に有益な書籍となっています。

本書の特徴と良い点

よく他の書籍でありがちな機械学習の各種手法を試しましたというようなチュートリアルで終わるような書籍ではなく、基本的な機械学習の理論とその実装まで詳しく書かれているため、しっかり読むことでロジックがよくわからないまま機械学習を使うようなことも避けることができます。

また、本書はscikit-learnによる実装が詳しく解説されているため、scikit-learnを一通り学ぶことできます。本家のscikit-learnのチュートリアルはやや難しいところや説明不足のところもあり初学者には中々難しいため、scikit-learnの入門としては本書から始められると良いと思います。

機械学習やscikit-learnにすでに熟練している読者の場合、本書のiPython Notebookがレファレンスとしても役立つと思います。

github.com

本書で触れられていない内容

本書は本当によくこの厚さでここまで盛り込んだなと思わせるほど密度の濃い内容となっていて全く問題ない内容ですが、わたしが思いつく範囲で触れられていない内容を確認すると次のような内容があります。

例えば、Gradient Boosting Decision Tree(GBDT)のライブラリのXGBoostや、Deep LearningのライブラリのChainerといった有名なライブラリの解説はありません。そのため最近の理解しておくべきライブラリについては別途、勉強する必要があります。

また機械学習の内容の中でも発展的な内容として、トピックモデルに関する機能を備えたgensimや、またFactorization Machinesなどをもちいたレコメンデーションの解説もありません。またCTR(Click Through Rate)の問題を解くためのオンライン学習のアルゴリズムであるFTRL-Proximalなども載っていません。そのため、興味のある個々の内容については各自で勉強する必要があります。あくまで本書は基本的な機械学習の内容をカバーしているものになります。

特におすすめの章

個人的にここは読むべきと思われる章を紹介します。記事の文面上おすすめの内容を2つの章に絞っておすすめします。

第4章 データ前処理 -- よりよいトレーニングセットの構築

この章は、特に機械学習初学者にとって必読の内容であると思います。機械学習によるモデル構築方法などは色々な解説がありますが、前処理について丁寧に解説している書籍は中々ありません。しかもそれをscikit-learnpandasを使った処理で解説しているため、非常に実用的です。特に初学者が陥りがちな誤りとして、カテゴリカルデータの取扱方法がありますが、しっかりと4.2.2. クラスラベルのエンコーディング4.2.3 名義特徴量でのone-hotエンコーディングの節で解説されています。

from sklearn.preprocessing import LabelEncoder

class_le = LabelEncoder()
y = class_le.fit_transform(df['classlabel'].values)
y
from sklearn.preprocessing import OneHotEncoder

ohe = OneHotEncoder(categorical_features=[0])
ohe.fit_transform(X).toarray()

このデータ前処理に関しては、この内容だけで、1冊の書籍が書けるほど奥深い分野です。より詳しい解説としては、本書の翻訳者である福島さんが執筆されたデータ分析プロセスを読まれることをおすすめします。

データ分析プロセス (シリーズ Useful R 2)

データ分析プロセス (シリーズ Useful R 2)

7章 アンサンブル学習 —異なるモデルの組み合わせ

この章はアンサンブル学習について知っていると思われる読者にも一読されたい内容でした。たいていのアンサンブル学習の解説は、最初は、決定木をもちいたアンサンブル手法であるRandom Forestの解説から始まり、Random Forestの実装をRandomForestClassifierをもちいた実装の解説で終わる場合が多いと思います。

しかし、本書では丁寧にBaggingの解説が載っています。DecisionTreeClassifierのEstimatorをBaggingClassifierに通す実装が紹介されています。

from sklearn.ensemble import BaggingClassifier
from sklearn.tree import DecisionTreeClassifier

tree = DecisionTreeClassifier(criterion='entropy', 
                              max_depth=None,
                              random_state=1)

bag = BaggingClassifier(base_estimator=tree,
                        n_estimators=500, 
                        max_samples=1.0, 
                        max_features=1.0, 
                        bootstrap=True, 
                        bootstrap_features=False, 
                        n_jobs=1, 
                        random_state=1)

これに対応する形として、Boostingを解説するために、DecisionTreeClassifierのEstimatorをAdaBoostClassifierに通す実装も紹介されているため、単なるアンサンブル学習の手法を学ぶための実装の紹介がなされているわけではなく、アンサンブル学習の本質を学ぶことができる構成になっています。

from sklearn.ensemble import AdaBoostClassifier

tree = DecisionTreeClassifier(criterion='entropy', 
                              max_depth=1,
                              random_state=0)

ada = AdaBoostClassifier(base_estimator=tree,
                         n_estimators=500, 
                         learning_rate=0.1,
                         random_state=0)

この章の最後に、Netflixコンペティションのエピソードでアンサンブル手法の注意点にも触れられています。このあたりの補足についても非常に読み応えがあります。

おわりに

以上、簡単ですが、Python機械学習プログラミングのご紹介をしました。この本書ほど機械学習の理論的な基礎を触れつつ、Pythonによる実装を解説した書籍は他にないと思います。ぜひ本書を手にとって頂き、実務で機械学習をもちいる際に役立ててもらえたら幸いです。

追記(2016/12/09)

よくscikit-learnを使ったソースコードを見ていると、X_train、x_train、train_x、Y_train、y_train、train_y、target、Yなどさまざまな書き方が散見されますが、本書では、P.10で解説されているように、特徴量のベクトルを小文字のx、行列を大文字のX、ラベルの列ベクトルを小文字のyと表して、X_train、y_train、X_test、y_testのように統一して記述されています。このコーディング規則に従うことを推奨します。

Deep Learningライブラリ「MXNet」のR版をKaggle Otto Challengeで実践してみた

最近、Kagglerの間で話題となっていたDeep Learningのライブラリ「MXNet」のRパッケージ{mxnet}についてTJOさんよりご紹介がありました。 tjo.hatenablog.com

多くの方がこのブログで初めて{mxnet}を知った、もしくはあまり触れてこなかったが久々に使ってみようという段階かもしれませんが、私をはじめ多くのKagglerはすでに{mxnet}を多少なり使用していることもあり、下記のように解説記事の書くようにと無茶ぶりを頂きましたので、私からもご紹介させて頂きたいと思います。

TJOさんの解説記事では、Deep Learningの代表的な事例でもある「Convolutional Neural Network」について解説されていますが、多くのデータサイエンスに関わる方、皆がコンピュータビジョン(画像認識等)を業務としているわけではないと思います。逆にコンピュータビジョンを業務とされている方々は、最近ではPythonから利用できるChainer、KerasといったDeep Learningライブラリをすでに実践しているというのが現状です。そこでより多くの方が気軽にRの{mxnet}を利用できるようにKaggleでの事例をご紹介したいと思います。

Kaggle「Otto Group Product Classification Challenge」

今回は、Otto Challengeで実践することにしました。 www.kaggle.com

このコンペティションは、Otto Group(20カ国以上で子会社を持つ世界でも最大級のECカンパニーのひとつ)から提供されたデータセットを題材としており、20万以上の商品を特徴付けた93変数(これら変数名は匿名化された数値変数)とその商品についての9種類の商品カテゴリ(この目的変数も匿名化されており9種類のカテゴリカル変数)のラベルを持ったデータセットとなっています。本コンペ参加者は、未知のラベルのデータセットに対してこの9種類の商品カテゴリを予測し、その精度を競いました。

このOtto Challengeは過去のKaggleの歴史上、最も多くの参加者を集めたコンペティションであり、Kaggleの中でも代表的なコンペのひとつです。マーケテイングにおける商品カテゴリ分類はセグメンテーションの下地にもなる重要なマーケテイング課題のひとつであると思います。その問題設定の取っ付きやすさもあり多くの参加者を集めた理由のひとつと考えられます。

また、最近では、Kaggle等のコンペティションにてGBDTのライブラリであるXGBoostが猛威をふるっているとたびたび話題にあがりますが、Kagglerから正確に現状を伝えると、XGBoostとDeep Learningのアンサンブルが猛威をふるっているというのが今の状況です。そしてこのOtto Challengeは、XGBoostだけでなく、Deep Learningの威力をまざまざと見せつけられた代表的な事例であると言えるでしょう。このコンペティションの上位者は絶対と言っていいほど、XGBoostとDeep Learningのアンサンブルを行っています。XGBoostでほとんどの参加者が良い精度を出すのは当たり前となっているため、ゆえにDeep Learningがいかにできるかが勝敗を分ける鍵であったと思います。

今回は、その勝負の鍵であったDeep LearningをRの{mxnet}をもちいて実践してみることとしました。

Multi Layer Perceptron

今回の問題設定は、9種類の商品カテゴリ分類となるため、下記の図のようなニューラルネットワークを組むことが基本となります。入力層はデータセット93変数を入力とし、隠れ層ではユニットを512、活性化関数としてRectified Linear Unit (ReLU)、さらにBatch NormalizationとDropoutを設定しています。この隠れ層をさらに同じアーキテクチャーの層を追加し、最終的に出力層としてSoftmax関数で9つの出力をしています。

このアーキテクチャーは、Kerasの開発者で現在Googleにお勤めのfcholletさんから提供されたサンプルコードを参考にしています。詳しくは下記リンクを参照してください。

Achieve 0.48 in 5 min with a deep net (feat. BatchNorm, PReLU) - Otto Group Product Classification Challenge | Kaggle

またMulti Layer PerceptronやReLU、Batch Normalizationの教科書的な解説は書籍に任せたいと思います。

f:id:Keiku:20160330184900p:plain

Customized mlogloss metricとEarly Stopping Ruleの設定

Otto Challengeで競う評価指標(Metric)は、multi-class logarithmic lossという指標であり、下記のように表されます。

{ \displaystyle
log loss = -\frac{1}{N}\sum_{i=1}^N\sum_{j=1}^My_{ij}\log(p_{ij})
}

Rの{mxnet}では評価指標としてAccuracyやRMSEは標準で実装されていますが、このmulti-class loglossという指標は、標準でサポートされていなかったため、mx.metric.customを用いて自作する必要があります。以下のようになります。

# Customized mlogloss eval_metric
mx.metric.mlogloss <- mx.metric.custom("mlogloss", function(label, pred){
  label_mat <- data.frame(i = 1:length(label),
                          j = label + 1,
                          x = rep(1, length(label)))
  label_mat <- sparseMatrix(i = label_mat$i,
                            j = label_mat$j,
                            x = label_mat$x)
  label_mat <- as.matrix(label_mat)
  return(MultiLogLoss(label_mat, t(pred)))
})

また、XGBoostには標準でサポートされているEarly Stopping Ruleの設定は{mxnet}ではmx.callback functionを使って自分で実装する必要があります。以下のようになります。

# Customized mx.callback function for early stopping
mx.callback.early.stop <- function (period, nrounds, logger = NULL, maximize=TRUE, verbose=FALSE)
{
  function(iteration, nbatch, env, verbose) {
    if (nbatch%%period == 0 && !is.null(env$metric)) {
      result <- env$metric$get(env$train.metric)
      if (nbatch != 0) 
        cat(paste0("Batch [", nbatch, "] Train-", result$name, 
                   "=", result$value, "\n"))
      if (!is.null(logger)) {
        if (class(logger) != "mx.metric.logger") {
          stop("Invalid mx.metric.logger.")
        }
        logger$train <- c(logger$train, result$value)
        if (!is.null(env$eval.metric)) {
          result <- env$metric$get(env$eval.metric)
          if (nbatch != 0) 
            cat(paste0("Batch [", nbatch, "] Validation-", 
                       result$name, "=", result$value, "\n"))
          logger$eval <- c(logger$eval, result$value)
          if (iteration > nrounds) {
            if(maximize) {
              if(sum(result$value > tail(logger$eval,nrounds)) == 0 ) { 
                bst.round <- which.max(logger$eval)
                if (verbose) {
                  print(paste0("Best round: ", bst.round, "; Best metric: ", max(logger$eval)))
                  print(paste0("This round: ", result$value))
                }
                return(FALSE) 
              }  
            } else if(sum(result$value < tail(logger$eval,nrounds)) == 0 ) { 
              bst.round <- which.min(logger$eval)
              if (verbose) {
                print(paste0("Best round: ", bst.round, "; Best metric: ", min(logger$eval)))
                print(paste0("This round: ", result$value))
              }
              return(FALSE)                           
            } 
          }
        }
      }
      return(TRUE)
    }
  }
}

実際にKaggleにSubmitしてみた

以下のコードを実行することで、submissionファイルが生成されます。次の手順でモデリングを行っています。

  1. データを読み込み、読み込んだデータをscaleを使って標準化する。
  2. caret::createFoldsを使ってValidation datasetを用意する。
  3. Validation datasetに対して予測を行う。
  4. Validation datasetの結果、最適なnum.roundを見つけて、Test datasetに対して予測を行う。

結果は、Batch Normalizationありのモデルで0.50705、Batch Normalizationなしのモデルで0.47950となりました。

options(scipen = 100)
options(dplyr.width = Inf)
options(dplyr.print_max = Inf)

require(mxnet)
require(caret)
require(dplyr)
require(data.table)
require(MLmetrics) # to evaluate multi-class logarithmic loss
require(stringr)
require(Matrix)

setwd("kaggle/otto-group-product-classification-challenge")

train <- fread("data/train.csv", data.table = F)
test <- fread("data/test.csv", data.table = F)
sampleSubmission <- fread("data/sampleSubmission.csv", data.table = F)

y <- as.numeric(as.factor(train[["target"]])) - 1
# drop id and target
X <- train[, -c(1, ncol(train))]
dim(X)
# [1] 61878    93
X_test <- test[, -c(1)]
dim(X_test)
# [1] 144368     93

All <- bind_rows(X, X_test)
# Scaling
All <- All %>%
  dplyr::mutate_each(funs(scale))
X <- data.matrix(All[1:nrow(X),])
X_test <- data.matrix(All[(nrow(X)+1):nrow(All),])

# -----------------------------------------
# Evaluate validation dataset
# -----------------------------------------

# Extract validation dataset
set.seed(123)
cv_list <- createFolds(y, k = 5, list = TRUE, returnTrain = FALSE)
X_train <- X[-cv_list[[1]], ]
dim(X_train)
# [1] 49501    93
X_valid <- X[cv_list[[1]], ]
dim(X_valid)
# [1] 12377    93
y_train <- y[-cv_list[[1]]]
y_valid <- y[cv_list[[1]]]

# Customized mlogloss eval_metric
mx.metric.mlogloss <- mx.metric.custom("mlogloss", function(label, pred){
  label_mat <- data.frame(i = 1:length(label),
                          j = label + 1,
                          x = rep(1, length(label)))
  label_mat <- sparseMatrix(i = label_mat$i,
                            j = label_mat$j,
                            x = label_mat$x)
  label_mat <- as.matrix(label_mat)
  return(MultiLogLoss(label_mat, t(pred)))
})

# Customized mx.callback function for early stopping
mx.callback.early.stop <- function (period, nrounds, logger = NULL, maximize=TRUE, verbose=FALSE)
{
  function(iteration, nbatch, env, verbose) {
    if (nbatch%%period == 0 && !is.null(env$metric)) {
      result <- env$metric$get(env$train.metric)
      if (nbatch != 0) 
        cat(paste0("Batch [", nbatch, "] Train-", result$name, 
                   "=", result$value, "\n"))
      if (!is.null(logger)) {
        if (class(logger) != "mx.metric.logger") {
          stop("Invalid mx.metric.logger.")
        }
        logger$train <- c(logger$train, result$value)
        if (!is.null(env$eval.metric)) {
          result <- env$metric$get(env$eval.metric)
          if (nbatch != 0) 
            cat(paste0("Batch [", nbatch, "] Validation-", 
                       result$name, "=", result$value, "\n"))
          logger$eval <- c(logger$eval, result$value)
          if (iteration > nrounds) {
            if(maximize) {
              if(sum(result$value > tail(logger$eval,nrounds)) == 0 ) { 
                bst.round <- which.max(logger$eval)
                if (verbose) {
                  print(paste0("Best round: ", bst.round, "; Best metric: ", max(logger$eval)))
                  print(paste0("This round: ", result$value))
                }
                return(FALSE) 
              }  
            } else if(sum(result$value < tail(logger$eval,nrounds)) == 0 ) { 
              bst.round <- which.min(logger$eval)
              if (verbose) {
                print(paste0("Best round: ", bst.round, "; Best metric: ", min(logger$eval)))
                print(paste0("This round: ", result$value))
              }
              return(FALSE)                           
            } 
          }
        }
      }
      return(TRUE)
    }
  }
}

# 2 layers MLP (Multi layer perceptron)
data <- mx.symbol.Variable("data")
fc1 <- mx.symbol.FullyConnected(data, name="fc1", num_hidden=512)
act1 <- mx.symbol.Activation(fc1, name="relu1", act_type="relu")
act1 <- mx.symbol.BatchNorm(act1, eps=1e-06, momentum=0.9, fix.gamma=FALSE, name="bn1")
dp1 <- mx.symbol.Dropout(act1, p = 0.5, name="dp1")

fc2 <- mx.symbol.FullyConnected(dp1, name="fc2", num_hidden=512)
act2 <- mx.symbol.Activation(fc2, name="relu2", act_type="relu")
act2 <- mx.symbol.BatchNorm(act2, eps=1e-06, momentum=0.9, fix.gamma=FALSE, name="bn2")
dp2 <- mx.symbol.Dropout(act2, p = 0.5, name="dp2")

fc3 <- mx.symbol.FullyConnected(dp2, name="fc3", num_hidden=9)
softmax <- mx.symbol.SoftmaxOutput(fc3, name="sm")

devices <- mx.cpu()
mx.set.seed(123)

logger <- mx.metric.logger$new()

model <- mx.model.FeedForward.create(
  softmax,
  X=X_train, y=y_train, eval.data=list(data=X_valid, label=y_valid),
  ctx=mx.cpu(),
  num.round=100, array.batch.size=128,
  optimizer = "sgd", learning.rate=0.01, momentum=0.9,
  initializer=mx.init.uniform(0.01),
  eval.metric=mx.metric.mlogloss,
  array.layout="rowmajor",
  epoch.end.callback = mx.callback.early.stop(1, 5, logger=logger, maximize=FALSE, verbose = TRUE)
  )
# With Batch Normalization
# [1] "Best round: 19; Best metric: 0.499036037381019"
# [1] "This round: 0.505431101369955"
# Without Batch Normalization
# [1] "Best round: 51; Best metric: 0.480038855711653"
# [1] "This round: 0.481640159668552"

# Check network architecture
graph.viz(model$symbol$as.json())

# Predict validation dataset
preds = predict(model, X_valid)
preds = t(preds)

y_valid_mat <- data.frame(i = 1:length(y_valid),
                          j = y_valid + 1,
                          x = rep(1, length(y_valid)))
y_valid_mat <- sparseMatrix(i = y_valid_mat$i,
                            j = y_valid_mat$j,
                            x = y_valid_mat$x)
y_valid_mat <- as.matrix(y_valid_mat)

# Evaluate multi-class logarithmic loss
print(MultiLogLoss(y_true=y_valid_mat, y_pred=preds))
# With Batch Normalization
# [1] 0.5047899
# Without Batch Normalization
# [1] 0.4805768

# -----------------------------------------
# Predict test dataset
# -----------------------------------------

# 2 layers MLP (Multi layer perceptron)
data <- mx.symbol.Variable("data")
fc1 <- mx.symbol.FullyConnected(data, name="fc1", num_hidden=512)
act1 <- mx.symbol.Activation(fc1, name="relu1", act_type="relu")
act1 <- mx.symbol.BatchNorm(act1, eps=1e-06, momentum=0.9, fix.gamma=FALSE, name="bn1")
dp1 <- mx.symbol.Dropout(act1, p = 0.5, name="dp1")

fc2 <- mx.symbol.FullyConnected(dp1, name="fc2", num_hidden=512)
act2 <- mx.symbol.Activation(fc2, name="relu2", act_type="relu")
act2 <- mx.symbol.BatchNorm(act2, eps=1e-06, momentum=0.9, fix.gamma=FALSE, name="bn2")
dp2 <- mx.symbol.Dropout(act2, p = 0.5, name="dp2")

fc3 <- mx.symbol.FullyConnected(dp2, name="fc3", num_hidden=9)
softmax <- mx.symbol.SoftmaxOutput(fc3, name="sm")

devices <- mx.cpu()
mx.set.seed(123)

model <- mx.model.FeedForward.create(
  softmax,
  X=X, y=y,
  ctx=mx.cpu(),
  num.round=19, array.batch.size=128,
  optimizer = "sgd", learning.rate=0.01, momentum=0.9,
  initializer=mx.init.uniform(0.01),
  eval.metric=mx.metric.mlogloss
)
# Check network architecture
graph.viz(model$symbol$as.json())

# Predict test dataset
preds = predict(model, X_test)
preds = t(preds)
preds = data.frame(preds)
names(preds) <- str_c("Class_", 1:9)

# Make a submission
submission <- bind_cols(data_frame(id = sampleSubmission$id),
                        preds)
write.csv(submission,
          "submit/submission.csv",
          quote = F,
          row.names = F)
# Private Leaderboard Score
# With Batch Normalization
# 0.50705
# Without Batch Normalization
# 0.47950

モデリングに関するメモ

Batch Normalizationの有無での比較

上記のコードで、Batch Normalizationありなしで比較も行ってみましたが、Batch Normalizationを入れた場合の方が精度を落としてしまいました。元の変数が頻度を表すような変数であり、バッチごとに正規化されたことが悪さをしたのか?実のところよくわかっていません。

Test dataset評価時におけるnum.roundの設定

上記のコードでは、Early Stopping RuleによりValidation datasetにおいて最も精度の高いnum.roundを見つけてそれをTest datasetでの評価時に設定していますが、Validation datasetは80%のデータで作成したモデルであるため、Test datasetを予測するためのTraining dataset(100%)とは少し異なります。num.roundは少し多めに設定しても問題ないと思いますし、最終的には適当にsubmitしてLeaderboard Scoreが最も高くなるnum.roundにすべきとしか言いようがないですが、このあたりどのように設定することがいいのか常に悩むところになります。

今後どのDeep Learningフレームワークを使うべきか?

Rの{mxnet}についてほんの少しご紹介しましたが、実際これだけDeep Learningフレームワークがあるなかでどれを選べばいいのか困ると思います。私のTwitterのタイムラインでは多くのKagglerが「Keras」をあげていました。私も同じくKerasを使っていることがほとんどですが、整理をすると以下のように考えています。

  1. Webサービスに組み込む、研究開発などの用途、またKaggleコンペティションにおいて魔改造が必要となる場合:Theano、Tensorflow、Caffe、Torch7、Chainer、MXNet(C++)など

  2. 通常のKaggleコンペティションやアドホックな分析用途:Keras、Lasagne、MXNet(Python/R)、Neonなど

この2番目の用途において、Rの{mxnet}によって分析可能となったことは多くのデータ分析者にとって喜ばしいことです。また、私としてはまだDeep Learningのモデリングは非常に難しく、さまざまなフレームワークでのDemoやExampleの実装を見て勉強する必要があると考えているため、Kerasに絞らずできるだけ多くのフレームワークに慣れ親しんでおこうと思っています。

追記 (2016/04/01)

TOP Kagglerのid:ultraistさんから貴重なコメント頂きました。

「R言語徹底解説」をご恵贈いただきました

Hadley Wickhamの著書「Advanced R」の邦訳「R言語徹底解説」を訳者のみなさまよりご恵贈いただきました,ありがとうございました.

R言語徹底解説

R言語徹底解説

読んでみて気になった所

R言語徹底解説」の原著の「Advanced R」についてですが,これは無料で公開されており,2013年11月頃にはほぼ内容も充実したものでした.当時から多くのRユーザーにとって新しい上級な教科書としての位置づけのドキュメントとして参照されているものになります.

現在,2016年の今でも「Advanced R」は色褪せない内容ではありますが,2年前ということもあり一部,新しい機能の実装などについては反映されていません.当然その翻訳本である「R言語徹底解説」にも同様に載っていない内容があります.「Advanced R」から親しんできたユーザーであるような人なら最新のR事情はキャッチアップされていると思うので問題ないと思いますが,「R言語徹底解説」からはじめる新規のRユーザーである場合,現在の新しい機能の実装との差分がどこにあるのかわからないかもしれません.そこで,すべてを取り上げることはできませんが,私が気になった部分について内容を補足をしたいと思います.

13章「非標準評価」を一通り読まれた方に紹介したい記事として,

Lazyeval: a new approach to NSE (2015-01-02)

をご紹介したいと思います.

Lazyeval: a new approach to NSE (2015-01-02)

この記事では,{lazyeval}パッケージによる新しい非標準評価(Non-Standard Evaluation)に対するアプローチの解説をしています.要点は3つあります.

  • substitute()を使用する代わりに,表現式と環境の両方を捕捉することができるlazyeval::lazy()を使用する.(また...でプロミスを捕捉するlazyeval::lazy_dots(...)を使用する.)

  • 非標準評価(NSE)を使用するすべての関数は,実際の計算をする標準評価(SE)のエスケープハッチを持っているべきである.標準評価(SE)関数の名前の末尾は_とする.

  • 標準評価(SE)関数は,ユーザーがプログラミングしやすいように柔軟な入力仕様を持っている.

とあります.以下,どういうことなのか例を見て確認してみたいと思います.

lazy()lazy_eval()

lazy()は,関数の引数に関連した表現式と環境の両方を捕捉するsubstitute()と等価のものとあります.Vignetteにsubstitute()の振る舞いが載っていないので比較してみます.substitute()の方は,environmentの情報が得られませんが、一見、同様の振る舞いに見えます.

g <- function(x = a - b) {
  substitute(x)
}
g()
# a - b
g(a + b)
# a + b
library(lazyeval)

f <- function(x = a - b) {
  lazy(x)
}
f()
#> <lazy>
#>   expr: a - b
#>   env:  <environment: 0x7fd2f3a406b0>
f(a + b)
#> <lazy>
#>   expr: a + b
#>   env:  <environment: R_GlobalEnv>

eval()を補完するものとして,{lazyeval}パッケージは,遅延オブジェクト(lazy object)に関連した環境を使用するlazy_eval()を備えているとあります.eval()と比較してみると以下のようになります.

a <- 10
b <- 1
eval(g())
# [1] 9
eval(g(a + b))
# [1] 11
a <- 10
b <- 1
lazy_eval(f())
#> [1] 9
lazy_eval(f(a + b))
#> [1] 11

標準評価(Standard evaluation)

私たちは非標準評価をするある関数を必要とするときはいつでも,常に最初に標準評価版を書こうとあります.あとで非標準評価版を実装するために,標準評価版をlazy()lazy_eval()を用いて実装してみます.

subset2_ <- function(df, condition) {
  r <- lazy_eval(condition, df)
  r <- r & !is.na(r)
  df[r, , drop = FALSE]
} 
subset2_(mtcars, lazy(mpg > 31))
#>     mpg cyl disp hp drat   wt qsec vs am gear carb
#> 18 32.4   4 78.7 66 4.08 2.20 19.5  1  1    4    1
#> 20 33.9   4 71.1 65 4.22 1.83 19.9  1  1    4    1

以下でも同様に振る舞います.

subset2_(mtcars, ~mpg > 31)
#>     mpg cyl disp hp drat   wt qsec vs am gear carb
#> 18 32.4   4 78.7 66 4.08 2.20 19.5  1  1    4    1
#> 20 33.9   4 71.1 65 4.22 1.83 19.9  1  1    4    1
subset2_(mtcars, quote(mpg > 31))
#>     mpg cyl disp hp drat   wt qsec vs am gear carb
#> 18 32.4   4 78.7 66 4.08 2.20 19.5  1  1    4    1
#> 20 33.9   4 71.1 65 4.22 1.83 19.9  1  1    4    1
subset2_(mtcars, "mpg > 31")
#>     mpg cyl disp hp drat   wt qsec vs am gear carb
#> 18 32.4   4 78.7 66 4.08 2.20 19.5  1  1    4    1
#> 20 33.9   4 71.1 65 4.22 1.83 19.9  1  1    4    1

非標準評価(Non-standard evaluation)

手元にある標準評価版で,非標準評価版を書くことは簡単である.私たちは環境に関連する未評価の表現式を捕捉するlazy()を使用するだけとあります.ここで13章「非標準評価」のsubset2()と比較してみます.:

subset2 <- function(df, condition) {
  subset2_(df, lazy(condition))
}
subset2(mtcars, mpg > 31)
#>     mpg cyl disp hp drat   wt qsec vs am gear carb
#> 18 32.4   4 78.7 66 4.08 2.20 19.5  1  1    4    1
#> 20 33.9   4 71.1 65 4.22 1.83 19.9  1  1    4    1

別の書き方をすると,

subset2 <- function(df, condition) {
  condition_call <- lazy(condition)
  r <- lazy_eval(condition_call, df)
  r <- r & !is.na(r)
  df[r, , drop = FALSE]
}
subset2(mtcars, mpg > 31)
# mpg cyl disp hp drat    wt  qsec vs am gear carb
# Fiat 128       32.4   4 78.7 66 4.08 2.200 19.47  1  1    4    1
# Toyota Corolla 33.9   4 71.1 65 4.22 1.835 19.90  1  1    4    1

これは,13章「非標準評価」のsubset2()(P288.)と対応しています.

subset2 <- function(df, condition) {
  condition_call <- substitute(condition)
  r <- eval(condition_call, df)
  r <- r & !is.na(r)
  df[r, , drop = FALSE]
}
subset2(mtcars, mpg > 31)
# mpg cyl disp hp drat    wt  qsec vs am gear carb
# Fiat 128       32.4   4 78.7 66 4.08 2.200 19.47  1  1    4    1
# Toyota Corolla 33.9   4 71.1 65 4.22 1.835 19.90  1  1    4    1

このように,13章「非標準評価」ではいきなり非標準評価の実装を行っていますが,手続きとしては標準評価版の実装を行い,その後lazy()を用いて非標準評価版の実装を行うことが推奨されています.これから私たちが実装する際に常にこうするべきとは思いませんが、少なくともHadley Wickhamはそのように実装しているんだということを知っているだけで彼のパッケージをさらに使いやすくなると思います.また、非標準評価(NSE)については批判もありますが,標準評価版の実装を常に用意してくれているという配慮を考えるとこのような批判はなくなるのかなと思います.

interp()

{lazyeval}パッケージには,interp()という便利な関数が実装されています.Interpolate values into an expression.(表現式に値を挿入する)と説明されていますが,下記のように振る舞います.

var_ <- "mpg"
threshold <- 31
interp(~var > x, var = as.name(var_), x = threshold)
# ~mpg > 31
subset2_(mtcars, interp(~var > x, var = as.name(var_), x = threshold))
#                 mpg cyl disp hp drat    wt  qsec vs am gear carb
# Fiat 128       32.4   4 78.7 66 4.08 2.200 19.47  1  1    4    1
# Toyota Corolla 33.9   4 71.1 65 4.22 1.835 19.90  1  1    4    1

このinterp()を用いれば,標準評価版の関数は,私たちがプログラミングしやすいように柔軟な入力仕様を持つことができます.

var_ <- "mpg"
threshold <- 31
above_threshold_ <- function(df, var_, threshold) {
  subset2_(mtcars, interp(~var > x, var = as.name(var_), x = threshold))
}
above_threshold_(mtcars, var_, 31)
#                 mpg cyl disp hp drat    wt  qsec vs am gear carb
# Fiat 128       32.4   4 78.7 66 4.08 2.200 19.47  1  1    4    1
# Toyota Corolla 33.9   4 71.1 65 4.22 1.835 19.90  1  1    4    1

above_threshold <- function(df, var_, threshold) {
  cond <- interp(~ var > x, var = lazy(var_), x = threshold)
  subset2_(df, cond)
}
above_threshold(mtcars, mpg, 31)
#                 mpg cyl disp hp drat    wt  qsec vs am gear carb
# Fiat 128       32.4   4 78.7 66 4.08 2.200 19.47  1  1    4    1
# Toyota Corolla 33.9   4 71.1 65 4.22 1.835 19.90  1  1    4    1

スコーピング(Scoping)

またlazy()は関数の引数に関連した環境を捕捉するので,以下のようなバグを回避することができます.

x <- 31
f1 <- function(...) {
  x <- 30
  subset(mtcars, ...)
}
# Uses 30 instead of 31
f1(mpg > x)
#>     mpg cyl disp  hp drat   wt qsec vs am gear carb
#> 18 32.4   4 78.7  66 4.08 2.20 19.5  1  1    4    1
#> 19 30.4   4 75.7  52 4.93 1.61 18.5  1  1    4    2
#> 20 33.9   4 71.1  65 4.22 1.83 19.9  1  1    4    1
#> 28 30.4   4 95.1 113 3.77 1.51 16.9  1  1    5    2

f2 <- function(...) {
  x <- 30
  subset2(mtcars, ...)
}
# Correctly uses 31
f2(mpg > x)
#>     mpg cyl disp hp drat   wt qsec vs am gear carb
#> 18 32.4   4 78.7 66 4.08 2.20 19.5  1  1    4    1
#> 20 33.9   4 71.1 65 4.22 1.83 19.9  1  1    4    1

lazy()substitute()よりもう一つ利点を持っており,非標準評価(NSE)をよりカジュアルに使用することが可能になります.

x <- 31
g1 <- function(comp) {
  x <- 30
  subset(mtcars, comp)
}
g1(mpg > x)
#> Error: object 'mpg' not found
g2 <- function(comp) {
  x <- 30
  subset2(mtcars, comp)
}
g2(mpg > x)
#>     mpg cyl disp hp drat   wt qsec vs am gear carb
#> 18 32.4   4 78.7 66 4.08 2.20 19.5  1  1    4    1
#> 20 33.9   4 71.1 65 4.22 1.83 19.9  1  1    4    1

以上,簡単ではありますが、{lazyeval}パッケージのご紹介でした*1.この内容は,13章「非標準評価」にも記載されていない内容となるのでぜひキャッチアップしておきたいところです.(私もまだまだ勉強中です.間違い等あればご教授頂ければ幸いです.)

おわりに

上記内容は「Advanced R」の補足ですが,「R言語徹底解説」の方について触れますと,私自身で{lazyeval}パッケージの簡単なご紹介を行うことだけでも,英語の解釈で詰まったり,適切な訳語がわからず,説明に悩みました.こんな少しのVignetteで非常に時間のかかったことを考えると,「Advanced R」という書籍の分量を翻訳されたこと,そしてその翻訳が極めて精緻なレベルであることに敬服するばかりです.このような書籍をRユーザーに届けてくれたことに感謝しかありません.またこの素晴らしい書籍を多くの方に読んでもらいたいと思いました.ぜひ手にとって頂ければ幸いです.

*1:GPL-3の問題でややグレーな内容かもしれません...