dplyr::filterの注意点
最近すっかりHadley信者になってしまいデータ加工にもdplyrをよく使っているのだけれど、filterで少し躓いてしまったのでメモ。
まずは{dplyr}と{dtplyr}を読み込む：
library(dplyr) library(dtplyr)
やりたかった処理とは以下のようなもので、irisを例とするとSpeciesごとにfilter内での条件を変えたかった。
> iris %>% + filter(Species == "setosa", Sepal.Length >= 3 | + Species == "versicolor", Sepal.Length >= 1 | + Species == "virginica", Sepal.Length >= 2) %>% tbl_dt() Source: local data table [50 x 5] # tbl_dt [50 × 5] Sepal.Length Sepal.Width Petal.Length Petal.Width Species <dbl> <dbl> <dbl> <dbl> <fctr> 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa 7 4.6 3.4 1.4 0.3 setosa 8 5.0 3.4 1.5 0.2 setosa 9 4.4 2.9 1.4 0.2 setosa 10 4.9 3.1 1.5 0.1 setosa # ... with 40 more rows
しかしこのOutputにはsetosaしか残っておらず、どうやらOR条件("|")以下が評価されていないようなのである。。
{dplyr}のfilterでは"&"もAND条件として使用可能なので、そちらで試してみるとうまく行く。
> iris %>% + filter(Species == "setosa" & Sepal.Length >= 3 | + Species == "versicolor" & Sepal.Length >= 1 | + Species == "virginica" & Sepal.Length >= 2) %>% tbl_dt() Source: local data table [150 x 5] # tbl_dt [150 × 5] Sepal.Length Sepal.Width Petal.Length Petal.Width Species <dbl> <dbl> <dbl> <dbl> <fctr> 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa 7 4.6 3.4 1.4 0.3 setosa 8 5.0 3.4 1.5 0.2 setosa 9 4.4 2.9 1.4 0.2 setosa 10 4.9 3.1 1.5 0.1 setosa # ... with 140 more rows
数値などは一切変えていないが、150行抽出されている。
しかし以下の2つは同じ結果を返すので、","と"&"だけの問題ではなさそう。
# ","でつなぐ > iris %>% + filter(Species == "setosa" , Sepal.Length >= 5) %>% tbl_dt() Source: local data table [30 x 5] # tbl_dt [30 × 5] Sepal.Length Sepal.Width Petal.Length Petal.Width Species <dbl> <dbl> <dbl> <dbl> <fctr> 1 5.1 3.5 1.4 0.2 setosa 2 5.0 3.6 1.4 0.2 setosa 3 5.4 3.9 1.7 0.4 setosa 4 5.0 3.4 1.5 0.2 setosa 5 5.4 3.7 1.5 0.2 setosa 6 5.8 4.0 1.2 0.2 setosa 7 5.7 4.4 1.5 0.4 setosa 8 5.4 3.9 1.3 0.4 setosa 9 5.1 3.5 1.4 0.3 setosa 10 5.7 3.8 1.7 0.3 setosa # ... with 20 more rows
# "&"でつなぐ > iris %>% + filter(Species == "setosa" & Sepal.Length >= 5) %>% tbl_dt() Source: local data table [30 x 5] # tbl_dt [30 × 5] Sepal.Length Sepal.Width Petal.Length Petal.Width Species <dbl> <dbl> <dbl> <dbl> <fctr> 1 5.1 3.5 1.4 0.2 setosa 2 5.0 3.6 1.4 0.2 setosa 3 5.4 3.9 1.7 0.4 setosa 4 5.0 3.4 1.5 0.2 setosa 5 5.4 3.7 1.5 0.2 setosa 6 5.8 4.0 1.2 0.2 setosa 7 5.7 4.4 1.5 0.4 setosa 8 5.4 3.9 1.3 0.4 setosa 9 5.1 3.5 1.4 0.3 setosa 10 5.7 3.8 1.7 0.3 setosa # ... with 20 more rows
"|"も組み合わせて使うと異なる結果が返るのだろうか。
> iris %>% + filter(Species == "setosa" , Sepal.Length >= 5 | Species == "virginica") %>% + tbl_dt() Source: local data table [30 x 5] # tbl_dt [30 × 5] Sepal.Length Sepal.Width Petal.Length Petal.Width Species <dbl> <dbl> <dbl> <dbl> <fctr> 1 5.1 3.5 1.4 0.2 setosa 2 5.0 3.6 1.4 0.2 setosa 3 5.4 3.9 1.7 0.4 setosa 4 5.0 3.4 1.5 0.2 setosa 5 5.4 3.7 1.5 0.2 setosa 6 5.8 4.0 1.2 0.2 setosa 7 5.7 4.4 1.5 0.4 setosa 8 5.4 3.9 1.3 0.4 setosa 9 5.1 3.5 1.4 0.3 setosa 10 5.7 3.8 1.7 0.3 setosa # ... with 20 more rows
> iris %>% + filter(Species == "setosa" & Sepal.Length >= 5 | Species == "virginica") %>% + tbl_dt() Source: local data table [80 x 5] # tbl_dt [80 × 5] Sepal.Length Sepal.Width Petal.Length Petal.Width Species <dbl> <dbl> <dbl> <dbl> <fctr> 1 5.1 3.5 1.4 0.2 setosa 2 5.0 3.6 1.4 0.2 setosa 3 5.4 3.9 1.7 0.4 setosa 4 5.0 3.4 1.5 0.2 setosa 5 5.4 3.7 1.5 0.2 setosa 6 5.8 4.0 1.2 0.2 setosa 7 5.7 4.4 1.5 0.4 setosa 8 5.4 3.9 1.3 0.4 setosa 9 5.1 3.5 1.4 0.3 setosa 10 5.7 3.8 1.7 0.3 setosa # ... with 70 more rows
ちなみにfilter内ではAND/OR条件の優先順位を付けるために括弧を使うことができるが、以下のように","を括弧で括る書き方はエラーとなる：
# ","を括弧でくくるとエラーになる > iris %>% + filter((Species == "setosa" , Sepal.Length >= 5)) %>% Error: unexpected ',' in: "iris %>% filter((Species == "setosa" ,"
# "&"や"|"なら大丈夫 > iris %>% + filter((Species == "setosa" & Sepal.Length >= 5)) %>% + tbl_dt() Source: local data table [30 x 5] # tbl_dt [30 × 5] Sepal.Length Sepal.Width Petal.Length Petal.Width Species <dbl> <dbl> <dbl> <dbl> <fctr> 1 5.1 3.5 1.4 0.2 setosa 2 5.0 3.6 1.4 0.2 setosa 3 5.4 3.9 1.7 0.4 setosa 4 5.0 3.4 1.5 0.2 setosa 5 5.4 3.7 1.5 0.2 setosa 6 5.8 4.0 1.2 0.2 setosa 7 5.7 4.4 1.5 0.4 setosa 8 5.4 3.9 1.3 0.4 setosa 9 5.1 3.5 1.4 0.3 setosa 10 5.7 3.8 1.7 0.3 setosa # ... with 20 more rows > > iris %>% + filter((Species == "setosa" | Sepal.Length >= 5)) %>% + tbl_dt() Source: local data table [148 x 5] # tbl_dt [148 × 5] Sepal.Length Sepal.Width Petal.Length Petal.Width Species <dbl> <dbl> <dbl> <dbl> <fctr> 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa 7 4.6 3.4 1.4 0.3 setosa 8 5.0 3.4 1.5 0.2 setosa 9 4.4 2.9 1.4 0.2 setosa 10 4.9 3.1 1.5 0.1 setosa # ... with 138 more rows
原因が明確になったわけではないが、ひとまずANDの指定には"&"を使ったほうが良さそうだ。
summaryの罠
年月を６桁の数値（YYYYMM）で表すために以下のように書いて何気なくsummaryを実行したところ、思わぬ挙動となった。
> Year <- rep(2012:2016, each=12) > Month <- rep(1:12, 5) > YM <- Year * 100 + Month > summary(YM) Min. 1st Qu. Median Mean 3rd Qu. Max. 201200 201300 201400 201400 201500 201600
上記の通り、月は1〜12の範囲なのにいずれも下2桁が00となっている。
念のため他の関数を使ってみると、ちゃんとした値が返る。
> fivenum(YM) [1] 201201.0 201303.5 201406.5 201509.5 201612.0 > min(YM) [1] 201201 > max(YM) [1] 201612
summaryはgeneric functionなので、今回の例のようなnumeric型に対するsummaryの挙動を確認するためにmethodを調べてみると：
> methods("summary") [1] summary.aov summary.aovlist* summary.aspell* [4] summary.check_packages_in_dir* summary.connection summary.data.frame [7] summary.Date summary.default summary.ecdf* [10] summary.factor summary.glm summary.infl* [13] summary.lm summary.loess* summary.manova [16] summary.matrix summary.mlm* summary.nls* [19] summary.packageStatus* summary.PDF_Dictionary* summary.PDF_Stream* [22] summary.POSIXct summary.POSIXlt summary.ppr* [25] summary.prcomp* summary.princomp* summary.proc_time [28] summary.srcfile summary.srcref summary.stepfun [31] summary.stl* summary.table summary.tukeysmooth*
色々と出てきたが、どうやらsummary.defaultがそれっぽい。
summary.defaultの中身を調べてみる。
> summary.default … 中略 … }) else if (is.numeric(object)) { nas <- is.na(object) object <- object[!nas] qq <- stats::quantile(object) qq <- signif(c(qq[1L:3L], mean(object), qq[4L:5L]), digits) names(qq) <- c("Min.", "1st Qu.", "Median", "Mean", "3rd Qu.", "Max.") if (any(nas)) c(qq, `NA's` = sum(nas)) else qq } … 中略 …
どうやらsummary.defaultの中身は{stats}のquantileで、間にmeanを挟んでいるだけの様子。そこでquantileを直接実行してみると：
> quantile(YM) 0% 25% 50% 75% 100% 201201.0 201303.8 201406.5 201509.2 201612.0
それっぽい値が返ってきている。どうやら次の処理のsignifの引数のdigitsが怪しい。
そこで今度はsummaryのdigitsを指定して実行してみる。
> summary(YM, digits = 8) Min. 1st Qu. Median Mean 3rd Qu. Max. 201201.0 201303.8 201406.5 201406.5 201509.2 201612.0
狙った通りの挙動となった。どうやら6桁の整数では大きすぎて小数点以下が丸められてしまっていたようだ。
めでたしめでたし・・・・と思ったら、よく見るとsummary（つまりquantile）の返り値とfivenumの返り値が微妙に違っている。
> summary(YM, digits = 8) Min. 1st Qu. Median Mean 3rd Qu. Max. 201201.0 201303.8 201406.5 201406.5 201509.2 201612.0 > fivenum(YM) [1] 201201.0 201303.5 201406.5 201509.5 201612.0
25パーセンタイルと75パーセンタイルの値が違う。調べてみると、fivenumは四分位数ではなくヒンジを返すのだそう。
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/59.html
なのでデータが偶数個の場合に先ほどのような挙動となる様子。
こんな感じ。
> hoge <- 1:10 > summary(hoge) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.00 3.25 5.50 5.50 7.75 10.00 > quantile(hoge) 0% 25% 50% 75% 100% 1.00 3.25 5.50 7.75 10.00 > fivenum(hoge) [1] 1.0 3.0 5.5 8.0 10.0 # 合わない！
> foo <- c(0, hoge) > summary(foo) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0 2.5 5.0 5.0 7.5 10.0 > quantile(foo) 0% 25% 50% 75% 100% 0.0 2.5 5.0 7.5 10.0 > fivenum(foo) [1] 0.0 2.5 5.0 7.5 10.0 # 一致した！
SparklyrによるApache Sparkのインストールとロジスティック回帰の実行
{sparklyr}というパッケージを使うことでWindowsであってもApache Sparkのインストールが簡単にできる。また{sparklyr}にはSpark MLlibの機械学習用の関数がラップされており、それを使ってみた結果を記しておく。
基本的にはRstudioの{sparklyr}の紹介ページをなぞってみただけだが、せっかくなので前回の記事で紹介した{MicrosoftML}のrxGLMと比較してみた。
なお当然ながら環境はWIndows。
まずは{sparklyr}のインストール。普通にCRANからインストールできる。
install.packages("sparklyr") library(sparklyr)
その他、必要なライブラリも合わせて読み込んでおく。
library(dplyr) library(readr) library(MicrosoftML) library(RevoScaleR)
spark_installという関数でSparkがあっさりとインストールできる。なお私の環境（AWSに立ち上げたばかりのWindows Server）ではC++ SP1とJAVAが入っていなかったため以下のようなエラーがでた。
> spark_install(version = "1.6.2") Error: Running Spark on Windows requires the Microsoft Visual C++ 2010 SP1 Redistributable Package. Please download and install from: https://www.microsoft.com/download/en/details.aspx?id=13523 Restart your rsession after installation completes
もしくは
Error in shell_connection(master = master, spark_home = spark_home, app_name = app_name, : Java is required to connect to Spark. Please download and install Java from https://www.java.com/en/
エラーメッセージにしたがってC++ SP1とJAVAをそれぞれインストールする。インストール後はRのリスタートだけでなくWindowsの再起動をかけた方が良さそう。
なお{sparklyr}はまだまだ開発途上なので、なるべく最新版を用いた方が変なエラーに悩まされなくてすむようだ。以下のように{sparklyr}を最新版に更新する。なぜか{shiny}とRtoolsが一緒にインストールされた。
devtools::install_github("rstudio/sparklyr")
以上で準備が整ったので、Sparkに接続する。なお私はSparkに触ったことがないため、以降のコマンドは正直よくわからない。
sc <- spark_connect(master = "local")
copy_toによりRのオブジェクト（主にdata.frameが対象か？）をSparkに入れることができる。
install.packages(c("nycflights13", "Lahman")) iris_tbl <- copy_to(sc, iris) flights_tbl <- copy_to(sc, nycflights13::flights, "flights") batting_tbl <- copy_to(sc, Lahman::Batting, "batting")
Sparkに入れたテーブルはRのdata.frameと同様に、以下のように{dplyr}の関数が適用できる。最後のcollectは{dplyr}の関数でSpark上のテーブルをRに読み込むのだが、この辺りは{dplyr}のLazinessに関係しており、以下を参照すると良い。
delay <- flights_tbl %>% group_by(tailnum) %>% summarise(count = n(), dist = mean(distance), delay = mean(arr_delay)) %>% filter(count > 20, dist < 2000, !is.na(delay)) %>% collect
ローカルのCSVなどのファイルを読み込む時には、以下のようにspark_read_csvを使う。前回同様、SUSYをサンプルデータとして使用したが、ロードにかかった時間は4分ほどだった。
system.time( datTmp <- spark_read_csv(sc, "SUSY", path = "../MicrosoftML/Data/SUSY.csv", header = FALSE) )
Sparkに入れるのと、Rに直接読み込むのでは比較しても仕方ないかもしれないが、参考用に{readr}のread_csvで同様にSUSYを読み込んでみた。その結果こちらは10分程度かかった。
なお本来なら上記のようにcollectすれば良いのだが、RStudioがクラッシュしてしまったので仕方なくread_csvしている。もしかしたらas.data.frameの方が早いかもしれない。
system.time( datTmpDF <- read_csv("../MicrosoftML/Data/SUSY.csv", col_names = FALSE) ) ユーザ システム 経過 56.14 32.89 699.00
読み込んだデータを使い、ロジスティック回帰を実行する。本来ならGLMで比較するべきなのだが、ml_generalized_linear_regressionがSpark2.0以降しか使えないため、やむなくロジスティック回帰を用いている。
system.time( spark_logit <- datTmp %>% ml_logistic_regression(response = "V1", features = c("V2", "V3", "V4", "V5")) ) ユーザ システム 経過 0.07 0.01 270.53
結果は4分半で終了した。500万件の分析であることを考えると十分な速度と思われる。
system.time( res_rxGLM <- rxGlm(X1 ~ X2 + X3 + X4 + X5, family = binomial("logit"), data = datTmpDF) ) Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.271 seconds Starting values (iteration 1) time: 1.456 secs. Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.180 seconds Iteration 2 time: 1.229 secs. Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.556 seconds Iteration 3 time: 1.598 secs. Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.229 seconds Iteration 4 time: 1.272 secs. Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.291 seconds Iteration 5 time: 1.331 secs. Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.263 seconds Iteration 6 time: 1.304 secs. Elapsed computation time: 8.572 secs. ユーザ システム 経過 0.89 1.02 27.23
一方、rxGlmは27秒と、相変わらず凄まじい速度で分析が実行できる。
それぞれの結果を、回帰係数で比較してみるとほぼほぼ同じ結果となった。
> cbind(spark_logit$coefficients, res_rxGLM$coefficients) [,1] [,2] (Intercept) -1.7100304090 -1.7100301776 V2 2.6339536980 2.6339542674 V3 -0.0004056923 -0.0004060596 V4 -0.0009712714 -0.0009717918 V5 -0.9571219696 -0.9571225996
比較は以上である。前回の記事で書いた通りMRCを使えばかなりの規模のデータに対してRからモデリングを実行できるため、現時点でSparkあるいは{sparklyr}にはそれほど興味は引かれない。しかしRに読み込めないような規模のデータを、Sparkに格納しつつml_generalized_linear_regressionで分析が可能なのであれば、データの規模に応じて分析環境を変えることで分析の幅が広がると感じた。よって次の検証内容としては{sparklyr}を用いてのデータサイズの限界を探るようなものが必要となりそうだ。
なおSparkはH2Oに接続することでH2Oの機械学習用のライブラリが使え、これをH2O Sparkling Waterと言うが、Rでも{sparkling}というパッケージで対応しているのでそこも調査する必要がありそうだ。
最後に、{sparklyr}のチートシートを発見したのでついでにリンクを貼っておく。
Microsoft R Clientによる大規模データの分析
Revolution Analyticsを買収したMicrosoftが、Revolution R Openに代わりMicrosoft R OpenというRのラッパーのようなものを出している。それに更に大規模データ分析用の独自開発パッケージを追加したMicrosoft R Client（MRC）というツールがあり、MRCの独自パッケージ{MicrosoftML}に実装されているrxGLMを用いて簡単な動作テストを実行したのでその結果を記しておく。MRCのインストールについてはこちらを参照：
まずは大規模データを読み込むために{readr}パッケージをインストール。なお今回はAWSに新しく立てたWindows Serverを使うので、{tidyverse}を使って必要そうなパッケージをまとめてインストールすることにする。
install.packages("tidyverse")
続いてライブラリの読み込みと作業ディレクトリの変更を行う。
library(readr) library(MicrosoftML) library(RevoScaleR) library(lattice) wk_dir <-"C:/Users/Administrator/Desktop/MicrosoftML" setwd(wk_dir)
データを読み込む。{MicrosoftML}がハイパフォーマンスやスケーラブルを謳っているので、そこそこの大規模データ（500万行×18列、2GB程度）を使用してみた。サンプルデータはこちらから取得した。
UCI Machine Learning Repository
### 大規模データの読み込みにreadr::read_csvを使用 ### ヘッダーがないのでcol_namesはFALSEにする datTmp <- read_csv("Data/SUSY.csv", col_names = FALSE) ### X1が目的変数、X2〜X5を説明変数として使用した ### 全カラムを使用した時の結果は後日 ### まずは結果の比較用にglmを実行。まぁ、動かないだろうと予想 start_time <- Sys.time() res_GLM <- glm(X1 ~ X2 + X3 + X4 + X5, family = binomial("logit"), data = datTmp) end_time <- Sys.time() print(end_time - start_time) # 予想通りメモリの割り当てに失敗した旨のエラーメッセージを受け取る。。。 ### MicrosoftMLではどうか。rxGLMが{MicrosoftML}におけるGLMに相当する関数。 start_time <- Sys.time() res_rxGLM <- rxGlm(X1 ~ X2 + X3 + X4 + X5, family = binomial("logit"), data = datTmp) end_time <- Sys.time()
> print(end_time - start_time) Time difference of 9.806952 secs
なんと、500万行のデータのglmが10秒足らずで実行できた。なおWindowsインスタンスはAWSでのt2.medium（64bit、RAM4GB）を使用しており、特別スペックが良いという訳ではない。
このままではglmとの比較ができないので、もう少しデータを減らして再挑戦してみる。
### gcを使ってメモリを解放しておく gc() gc() ### 使用データを100万行にする datTmp_1MM <- datTmp[1:1000000, ] ### 再実行 start_time <- Sys.time() res_rxGLM_02 <- rxGlm(X1 ~ X2 + X3 + X4 + X5, family = binomial("logit"), data = datTmp_1MM) end_time <- Sys.time() print(end_time - start_time) # Time difference of 1.839989 secs gc() gc() start_time <- Sys.time() res_GLM_02 <- glm(X1 ~ X2 + X3 + X4 + X5, family = binomial("logit"), data = datTmp_1MM) end_time <- Sys.time() >print(end_time - start_time) # Time difference of 7.471947 secs
今度はglmでも10秒かからず実行できた。
そこで両者の結果が一致するかを確認する。
coef(res_GLM_02) (Intercept) X2 X3 X4 X5 -1.7128408655 2.6384939818 0.0002822391 -0.0012571389 -0.9573449310 coef(res_rxGLM_02) (Intercept) X2 X3 X4 X5 -1.7128408655 2.6384939818 0.0002822391 -0.0012571389 -0.9573449310
回帰係数は両者で一致しており、問題なさそうだ。
続いてsummaryを実行する。
### まずはrxGLM > summary(res_rxGLM_02) Call: rxGlm(formula = X1 ~ X2 + X3 + X4 + X5, data = datTmp_1MM, family = binomial("logit")) Generalized Linear Model Results for: X1 ~ X2 + X3 + X4 + X5 Data: datTmp_1MM Dependent variable(s): X1 Total independent variables: 5 Number of valid observations: 1e+06 Number of missing observations: 0 Family-link: binomial-logit Residual deviance: 1160233.9234 (on 999995 degrees of freedom) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.7128409 0.0052066 -328.976 2.22e-16 *** X2 2.6384940 0.0079644 331.285 2.22e-16 *** X3 0.0002822 0.0022039 0.128 0.898 X4 -0.0012571 0.0022356 -0.562 0.574 X5 -0.9573449 0.0062846 -152.332 2.22e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Condition number of final variance-covariance matrix: 7.0668 Number of iterations: 6 ### 続いてglm > summary(res_GLM_02) Call: glm(formula = X1 ~ X2 + X3 + X4 + X5, family = binomial("logit"), data = datTmp_1MM) Deviance Residuals: Min 1Q Median 3Q Max -8.4904 -0.8894 -0.7140 1.0654 3.0185 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.7128409 0.0052066 -328.977 <2e-16 *** X2 2.6384940 0.0079644 331.287 <2e-16 *** X3 0.0002822 0.0022039 0.128 0.898 X4 -0.0012571 0.0022356 -0.562 0.574 X5 -0.9573449 0.0062846 -152.333 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1379205 on 999999 degrees of freedom Residual deviance: 1160234 on 999995 degrees of freedom AIC: 1160244 Number of Fisher Scoring iterations: 5
summaryの結果も両者で大きく異なる様子はない。rxGLMの方でもう少し情報を追加したぐらいか。
一旦動作テストは以上である。{MicrosoftML}は少なくともオリジナルのRのglmと比較して大規模データに強いというのは間違いなさそうで、Rの弱点をカバーできる作りになっているようだ。従来MicrosoftはSQL ServerとRを連携させることで大規模データへの対応を行ってきたイメージがあったが、R単体（と言えないかもしれないが）でこれぐらいの規模のデータを処理できるのであれば、ウェブログを触るような状況とならない限りR一本で大抵の処理はできるように思う。RユーザーとしてはMRCを含めて近年のMicrosoftのRへの貢献はかなり大きいと感じている。
また今回の記事では割愛したが、{MicrosoftML}ではrxGLM以外にも様々な関数が用意されており、特に二項分類のための”Fast Tree"と呼ばれるアルゴリズムを実装した関数はかなり強力な感触を受けた。こちらの検証結果も後日書いてみようと思う。
人工知能と機械学習とDeep Learning
「人工知能（AI）による業務改革！」「Deep Learningを搭載した○○！」みたいな記事は巷に溢れかえっているが、そういった記事を見るにつけ思うのが「人工知能と機械学習、Deep Learningなどのアルゴリズムの区別ってどうなっているんだろう？」ということである。自分で機械学習のアルゴリズムを触ったことがあればわかるかと思うのだが、現在世間で言われているところの「人工知能」とは、「自律的な思考を行う」だとか「知性を獲得した」だとかそういったものでは全くなく、Deep Learningもそのための特殊なアルゴリズムというわけではない。
この辺りの関係性をもう少しスッキリと説明できないかなーと考えていたところ、以下の記事を発見した：
またDeep Learningの大家であるモントリオール大学のBengioも同様に整理している：
http://www.deeplearningbook.org/contents/intro.html
私の人工知能に対する認識も上記と同じで（権威に阿ることにする）、Deep Learningは機械学習の中の1アルゴリズム（いわゆる「弱いAI」）として理解している。そのためDeep Learningを指して「人工知能」と称することにはかなりの抵抗がある。もしそれが成り立つなら、同じく機械学習の1アルゴリズムである（例えば）Random Forestを指して「人工知能」と呼んだって良いわけで、それならば決定木だってアリということになるだろうし、さらには線形回帰だって機械学習の1アルゴリズムとして捉えれば「人工知能」となってしまうだろう（もしかしたら実際に線形回帰程度のものを人工知能と謳っているところがあるかもしれない）。
「人工知能による〜〜」という言葉に惑わされず、内容を吟味した上で評価することを心がけたい。
RでTensorFlow
知らない間にRでTensorFlowが使えるようになっていたので触ってみました。それにしてもRStudioは相変わらずイイ仕事をしますね。
まずは以下の通り、RStudioのGitHubからTensorFlowのライブラリをインストールします。なおこのライブラリはあくまでPCに事前にインストールされたTensorFlowをRから呼び出すためのものなので、TensorFlow自体は別途インストールする必要があります。
## TensorFlowのライブラリをインストールする devtools::install_github("rstudio/tensorflow") ## ライブラリの読み込み ## データ加工用にdplyrとcaret、分析結果の比較用にglmnetとrandom forestのライブラリも読み込む library(tensorflow) library(dplyr) library(caret) library(glmnet) library(ranger)
無事にインストールが完了したら、一旦TensorFlowのバージョンを確認しておきましょう。ちなみに私の環境は0.11です。
## TensorFlowのバージョン確認 tf$VERSION [1] "0.11.0rc0"
では早速、TensorFlowによる「Hello,World!」です。この辺りは普通にマニュアルにある通りです。
## hello, tensorflow sess <- tf$Session() hello <- tf$constant("Hello, TensorFlow!") sess$run(hello) [1] "Hello, TensorFlow!"
単回帰の実行
続いてTensorFlowによる単回帰を実行してみます。この書き方を理解できれば重回帰への拡張は容易でしょう。
## seedを固定しておく set.seed(123) ## データ作成。回帰のパラメータは切片が-12、傾きが1.5（適当です） x_data <- runif(100, min=0, max=1) y_data <- x_data * 1.5 - 12 ## TensorFlowによる変数の宣言 W <- tf$Variable(tf$random_uniform(shape(1L), -1.0, 1.0)) b <- tf$Variable(tf$zeros(shape(1L))) y <- W * x_data + b ## 損失関数の定義と最適化のパラメータ設定 loss <- tf$reduce_mean((y-y_data)^2) opt <- tf$train$GradientDescentOptimizer(0.2) train <- opt$minimize(loss)
上記で分析環境が整ったので以下の通り実行しますが、その際の注意点として、TensorFlowのバージョンによって初期化に用いる関数が異なることに気をつけましょう。私の環境（0.11）では初期化の関数は「initialize_all_variables()」なのですが、新しいバージョン（0.12）では「global_variables_initializer()」が用いられるようです。
sess <- tf$Session() ## 変数の初期化。TensorFlowのバージョンによって、初期化に用いる関数が異なる。 # sess$run(tf$global_variables_initializer()) sess$run(tf$initialize_all_variables()) ## 学習の実行 for (step in 1:1000){ sess$run(train) if (step %% 20 == 0) { cat(step, "-", sess$run(W), sess$run(b), "\n") } } 20 - -2.501043 -9.869881 40 - -0.8908119 -10.72716 60 - 0.0713779 -11.23942 80 - 0.646331 -11.54552 100 - 0.9898927 -11.72843 ： 920 - 1.499992 -12 940 - 1.499992 -12 960 - 1.499992 -12 980 - 1.499992 -12 1000 - 1.499992 -12
LearningRateは小さめに設定してありますが、ちゃんと収束しているようですね。
多項分類
続けてIrisのデータを用いてSpeciesの分類を行います。
## Irisデータの多項分類 W <- tf$Variable(tf$zeros(shape(4L, 3L))) x <- tf$placeholder(tf$float32, shape(NULL, 4L)) b <- tf$Variable(tf$zeros(3L)) y <- tf$nn$softmax(tf$matmul(x, W)+b) y_ <- tf$placeholder(tf$float32, shape(NULL, 3L)) ## 損失関数の定義 loss <- tf$reduce_mean(-tf$reduce_sum(y_*y, reduction_indices = 1L)) opt <- tf$train$GradientDescentOptimizer(0.9) train <- opt$minimize(loss) ## 分析用データ（Iris） trn_x <- iris %>% select(-Species) %>% as.matrix() trn_y <- iris %>% dummyVars(formula = ~ Species, sep = NULL) %>% predict(object = ., newdata = iris) %>% as.matrix() sess <- tf$Session() sess$run(tf$initialize_all_variables()) ## 学習の実行 N <- 20000 each <- 1000 resAll <- matrix(NA, N, 15) for (i in 1:N){ sess$run(train, dict(x=trn_x, y_=trn_y)) resAll[i, ] <- c(sess$run(W), sess$run(b)) if (i %% each == 0) { cat(i, "-", sess$run(W), sess$run(b), "\n") print(table( predict = sess$run(fetches = y, feed_dict = dict(x = trn_x, y_ = trn_y)) %>% apply(margen=1, FUN=which.max), true = trn_y %>% apply(margen=1, FUN=which.max))) } } 1000 - 1.820591 3.826579 -5.199587 -2.53192 1.656648 0.04738447 -0.9168726 -2.692287 -3.47724 -3.873963 6.116461 5.224204 0.8330443 1.588295 -2.421341 true predict 1 2 3 1 50 0 0 2 0 46 0 3 0 4 50 2000 - 2.170959 4.40071 -6.013396 -2.967016 1.705534 0.348204 -1.100592 -3.198786 -3.876493 -4.748915 7.113993 6.165799 0.9800377 2.52502 -3.505058 true predict 1 2 3 1 50 0 0 2 0 47 0 3 0 3 50 ： 19000 - 4.02061 6.796506 -9.696184 -4.990173 2.587524 0.9785848 -2.744421 -5.431025 -6.608002 -7.775012 12.4406 10.42116 1.755378 8.18106 -9.936396 true predict 1 2 3 1 50 0 0 2 0 47 0 3 0 3 50 20000 - 4.081785 6.860996 -9.807202 -5.050861 2.61721 0.9981934 -2.800498 -5.49942 -6.698854 -7.859092 12.60772 10.55024 1.781228 8.357506 -10.13868 true predict 1 2 3 1 50 0 0 2 0 47 0 3 0 3 50
20,000回ほど回してみたのですが、意外と分類が難しいのかパラメータは収束していません。また分類精度も最初の方で現在のレベルまで到達して以降、全く改善はありませんでした。この辺りはLearningRateや初期値の設定などに依存するのかもしれません。
最後に、以上の結果を他の分類器と比較してみます。比較対象にはglmnetとrandom forestを用いました。
## glmnetによる結果と比較してみる resGLM <- glmnet(trn_x, iris$Species, "multinomial") table(predict = predict(resGLM, newx = trn_x, type = "class", s = 0.01)[, 1], true = iris$Species) true predict setosa versicolor virginica setosa 50 0 0 versicolor 0 47 1 virginica 0 3 49 ## random forestによる結果と比較してみる resRF <- ranger(Species ~., data=iris, mtry=2, num.trees = 500, write.forest=TRUE) table(predict = predict(resRF, data=iris)$predictions, true=iris$Species) true predict setosa versicolor virginica setosa 50 0 0 versicolor 0 50 0 virginica 0 0 50
random forestは流石の精度ですね。Irisは分類が容易なデータなので、TensorFlow、glmnetもそれなりの精度で分類できています。
雑感
以上、RによるTensorFlowの流れを追ってみたのですが、正直に言えば「コレジャナイ」感が強いなーというのが私の感想です。一番引っ掛かりを感じる部分が「事前にTensorFlowのインストールが必要なこと」で、Rをメインに使用しているユーザーからはPythonの環境整備に対するハードルが高そうだなと思いました。他のライブラリのようにTensorFlow自体もinstall.packages()でインストールできると良いな、と。
もう一つ、これはPython用のTensorFlowをRから呼び出しているためなのでしょうが、書き方がRっぽくないなーと感じました。特に「tf$Session()」の部分などいかにもPythonっぽくてちょっと取っつきにくいなぁと思ってしまいます。
もちろんこれは私がそのように感じるというだけですので、Pythonにも馴染みのある方なら全く問題とはならないでしょう。ただ、私としてはmxnetの方が良さそうだなと思いましたので、もう少し検討を続けようと思います。
optimってあんまり信用できないなぁ、って話
タイトルの通りです。
仕事で使うことがあったのでRのoptimを使って回帰を解いてみたのですが、これが意外に安定しません。変数の数なのか、ダミー変数が含まれるとダメなのか、原因についてはよくわかりませんが想定以上に解がバラついてしまいました。
実行したコードは以下の通りです。まずは以下のようにデータとoptimに渡す目的関数を定義します：
## サンプルデータとしてirisを使用 dat <- iris datMdl <- as.data.frame(model.matrix(~.-1, data = dat)) ## Speciesを含まないデータ dat01 <- datMdl[, 1:4] ## Speceisを含むデータ（virginicaを除外） dat02 <- datMdl[, -7] ## 最小二乗法を以下のように定義 myLM <- function(b, dat) { Y <- as.matrix(dat[, 1]) X <- as.matrix(cbind(1, dat[, -1])) b <- as.matrix(b) prd <- X %*% b e <- Y - prd return(t(e) %*% e) }
以上の条件で、まずはSpeciesを含まないデータを用いてlmとoptimをそれぞれ実行します。
## Speciesを含まないデータでSepal.Lengthをyとして回帰をかける resLM_01 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data=dat01) ## optimによる回帰係数の取得 init_coef <- rep(1, 4) resOpt_01 <- optim(par=init_coef, fn=myLM, dat=dat01[, 1:4]) ## 係数の確認。だいたい同じ cbind( summary(resLM_01)$coef[, 1], resOpt_01$par) [,1] [,2] (Intercept) 1.8559975 1.8589403 Sepal.Width 0.6508372 0.6501880 Petal.Length 0.7091320 0.7087105 Petal.Width -0.5564827 -0.5559860
コメントでも書いてありますが、係数は概ね一致します。続けてSpeciesを含めたデータによりlmおよびoptimをそれぞれ再度実行します。
## 次にSpeciesを含むデータで同様に回帰を実行 resLM_02 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Speciessetosa + Speciesversicolor, data=dat02) ## 同じく、optimによる回帰係数の取得 init_coef <- rep(1, 6) resOpt_02 <- optim(par=init_coef, fn=myLM, dat=dat02) ## 係数の確認。結構違う cbind( summary(resLM_02)$coef[, 1], resOpt_02$par) [,1] [,2] (Intercept) 1.1477685 0.6123026 Sepal.Width 0.4958889 0.4621860 Petal.Length 0.8292439 0.8993939 Petal.Width -0.3151552 -0.1967130 Speciessetosa 1.0234978 1.5512063 Speciesversicolor 0.2999359 0.4695160
Rの最適化ソルバーはoptim以外にもありますし(optimizeとか)、optimを使うにしてもパラメータや用いるアルゴリズム(Nelder-Mead、BFGSとか)など設定できる箇所は色々とあります。なのでこの結果だけを元にoptimを信用できないとは言えないのですが、このくらいのデータ・モデルでこれほどlmとの差が出てしまうと複雑な問題への適用はちょっとためらってしまいますね。
なお当初の目的としてはMarketing Mix Modelによる広告効果の推定を行う際に、残存効果の減衰率をoptimで推定してみようというものでした。optimの結果があまり芳しくないためこういった検証を行っていたのですが、この辺はいずれ記事にまとめようと思います。それからoptimはオプションでHessian(ヘッセ行列)を出力することにより標準誤差を計算できますので、その辺もいずれ書いてみようと思います。
