Keras for R
RstudioがR上でKerasによるディープラーニングのモデルを構築するためのライブラリ{keras}を公開した。
以前から{tensorflow}を使えばtensorflow::import(module = "keras")でKerasを導入することができたようだが、{keras}を先にインストールすることでpythonさえ入っていればtensorflowのインストールもできるらしい。随分と楽にtensorflowを動かせることができそうなので早速試してみた。なお類似のプロジェクトがRstudio以外でも行われているようで、以下でもKerasをRから動かすためのライブラリを提供している。
実行環境はAWS上のWindowsマシン。Rのバージョンは3.3.2だけれど、Microsoft R Clientがインストールされている。
実行したスクリプトは上で紹介したページの内容なのだけれど、サンプルデータであるx_trainやx_testなどは自前で用意する必要があったのでirisを使って適当に作成した。なおこちらのページはKerasの公式ドキュメントをR用に書き換えただけの様子。
{Keras}のインストール
早速{keras}のインストールとtensorflowのインストールを実行してみる。なお私の環境では{keras}のインストールの際に{reticulate}が一緒に入らず、{keras}のインストールに失敗したため{reticulate}だけ先にインストールしている。
## library install & load devtools::install_github("rstudio/reticulate") devtools::install_github("rstudio/keras") library(keras) library(dplyr) # for data handling ## install tensorflow install_tensorflow() Error: Installing TensorFlow requires a 64-bit version of Python 3.5 Please install 64-bit Python 3.5 to continue, supported versions include: - Anaconda Python (Recommended): https://www.continuum.io/downloads#windows - Python Software Foundation : https://www.python.org/downloads/release/python-353/ Note that if you install from Python Software Foundation you must install exactly Python 3.5 (as opposed to 3.6 or higher).
{keras}のインストールがうまくいけばinstall_tensorflow()という関数一つでtensorflowをインストールすることができる。とはいえpythonのインストールまではやってくれないようで私の環境ではエラーが出てしまった。指示に従い、Anaconda経由でpythonをインストールする。なお現時点でのpythonのバージョンは3.6だった。
Anacondaのインストールが済んだので再度tensorflowのインストールを試してみる。
## install tensorflow install_tensorflow()
今度はうまくいった。
サンプルデータの作成
サンプルスクリプトでは分析用のデータが提供されていなかったため、irisを使ってサンプルデータを作成する。特別なことはしておらず、sampleによって乱数を振った上で学習用とテスト用にデータを分割した。なおデータはmatrixにする必要がある。
set.seed(123) learn_id <- sample(1:nrow(iris), 100, replace = FALSE) train_data <- iris %>% filter(row.names(.) %in% learn_id) test_data <- iris %>% filter(!row.names(.) %in% learn_id) ## training data x_train <- train_data %>% select(-Species) %>% as.matrix() y_train <- as.matrix(model.matrix(~ train_data$Species - 1)) ## test data x_test <- test_data %>% select(-Species) %>% as.matrix() y_test <- as.matrix(model.matrix(~ test_data$Species - 1))
Kerasによるmodelオブジェクトの生成
サンプルデータが出来上がったので、{keras}によるフィッティングを行う。まずはmodelオブジェクトの生成から。
## generate model object model <- keras_model_sequential() > Model Model _______________________________________________________________________________ Layer (type) Output Shape Param # =============================================================================== Total params: 0 Trainable params: 0 Non-trainable params: 0 _______________________________________________________________________________
入力層 → 隠れ層 → 出力層という簡単なニューラルネットワークを構築する。layer_denseで簡単にレイヤーを積み上げることができるようで、入力層のみinput_shapeを指定する必要がある。これはpython上でKerasを実行する際のmodel.addに相当するものと思われる。
model %>% layer_dense(units = 10, input_shape = 4) %>% # Output & Input dimension layer_activation(activation = 'relu') %>% layer_dense(units = 3) %>% layer_activation(activation = 'softmax') > model Model _______________________________________________________________________________ Layer (type) Output Shape Param # =============================================================================== dense_3 (Dense) (None, 10) 50 _______________________________________________________________________________ activation_3 (Activation) (None, 10) 0 _______________________________________________________________________________ dense_4 (Dense) (None, 3) 33 _______________________________________________________________________________ activation_4 (Activation) (None, 3) 0 =============================================================================== Total params: 83.0 Trainable params: 83.0 Non-trainable params: 0.0 _______________________________________________________________________________
パラメータの数は層ごとに(入力データの次元数+1) * 出力データの次元で決まる。+1は切片。
ここで出来上がったmodelオブジェクトの属性や構造などを見てみようとstrを実行してみたのだが、print(model)と結果は変わらなかった。やはりR的なオブジェクトではない様子。
ニューラルネットワークの構造が決まれば、最適化の際の条件を指定する。
## model parameter model %>% compile( loss = 'categorical_crossentropy', optimizer = optimizer_sgd(lr = 0.02), metrics = c('accuracy'))
フィッティング
以上で準備が整ったため、フィッティングを行う。
## fitting model %>% fit(x_train, y_train, epochs = 1000, batch_size = 32) ## 省略 # Epoch 990/1000 # 100/100 [==============================] - 0s - loss: 0.0432 - acc: 0.9800 # Epoch 991/1000 # 100/100 [==============================] - 0s - loss: 0.0442 - acc: 0.9800 # Epoch 992/1000 # 100/100 [==============================] - 0s - loss: 0.0455 - acc: 0.9800 # Epoch 993/1000 # 100/100 [==============================] - 0s - loss: 0.0508 - acc: 0.9800 # Epoch 994/1000 # 100/100 [==============================] - 0s - loss: 0.0602 - acc: 0.9800 # Epoch 995/1000 # 100/100 [==============================] - 0s - loss: 0.0437 - acc: 0.9800 # Epoch 996/1000 # 100/100 [==============================] - 0s - loss: 0.0433 - acc: 0.9800 # Epoch 997/1000 # 100/100 [==============================] - 0s - loss: 0.0432 - acc: 0.9800 # Epoch 998/1000 # 100/100 [==============================] - 0s - loss: 0.1129 - acc: 0.9600 # Epoch 999/1000 # 100/100 [==============================] - 0s - loss: 0.0529 - acc: 0.9700 # Epoch 1000/1000 # 100/100 [==============================] - 0s - loss: 0.0647 - acc: 0.9600
サンプルデータがものの100行程度しかないため一瞬で終わる。Accuracyが0.96とまずまず。続いてこちらのモデルを用いて、テストデータに対する精度評価を実行する。
## test loss_and_metrics <- model %>% evaluate(x_test, y_test, batch_size = 128) > loss_and_metrics [[1]] [1] 0.1853937 [[2]] [1] 0.96
テストデータに対しても精度が0.96と高く、過学習を起こしている様子はない。
テストデータに対するラベルの予測は以下のように実行する。
## predict classes <- model %>% predict(x_test, batch_size = 128) > classes [,1] [,2] [,3] [1,] 9.983064e-01 1.693639e-03 6.461022e-20 [2,] 9.971897e-01 2.810343e-03 7.951551e-19 [3,] 9.983315e-01 1.668578e-03 3.816257e-20 [4,] 9.984826e-01 1.517428e-03 4.521292e-20 [5,] 9.994913e-01 5.087039e-04 3.988174e-23 ## 省略
サンプルスクリプトは以上である。
注意点として、パイプ演算子(%>%)が使えるからといってdplyrが読み込まれている訳ではないので、データ加工を行うのであれば別途読み込む必要がある。
Rでディープラーニングを簡単に実行できる環境が整ってきているので、これを機会にRユーザーが増えると良いな。
帰無仮説は採択できない
統計的仮説検定は通常、以下のような手順に従って行われる：
- 帰無仮説を設定する
- 対立仮説を設定する
- 検定統計量を設定する
- 有意水準を設定する
- 実験やデータ解析によって検定統計量を求める
- 帰無仮説を真と仮定した時の検定統計量の得られる確率を求める
- 設定した有意水準に従い、判断を下す
このうち７番で行われる判断において、帰無仮説が棄却された場合は「対立仮説を採択する」という表現がなされる。では対立仮説が棄却された場合は「帰無仮説を採択する」と言うのかといえばそうではない。帰無仮説は採択できないからだ。
そもそも統計的仮説検定は「帰無仮説が真であるとき、この統計量が偶然に得られる確率はどの程度であろうか？」を定量的に評価するための手続きである。大事なポイントとしては「帰無仮説が真である」という点で、この前提の上でどの程度の検定統計量が得られれば帰無仮説を棄却できるか、ということを見ている。帰無仮説と対立仮説は排反であるため、帰無仮説が棄却された場合は自動的に対立仮説が採択される。
では帰無仮説が棄却されなかった場合はと言うと、これがややこしいところで、「帰無仮説を採択する」ことはできないのである。「棄却できない」と「採択する」は異なるからだ。
先にも述べたように、統計的仮説検定はあくまで「帰無仮説を棄却できるか」を検証するための手続きである。この手続きから得られる結果は「帰無仮説を棄却できる」「できない」の２つだけであり、採択に関する考慮は一切ない。そのため、例えば回帰分析を行った際にβ = 1.5, P < 0.01という回帰係数が得られたとしても、ここで言えることは「回帰係数が0であること（=帰無仮説）が棄却されたので、回帰係数は0ではない」ということだけであり、1.5という数値そのものに対しては何も言えないのである。
ちなみに検定統計量（t値やχ2値など）の大小をもって「極めて有意」「わずかに有意」などといった定量的は判断が下されることがあるが、厳密にはこれはよろしいものではない。伝統的な手続きに則った場合、棄却できるか否か（1か0か）の判断だけが可能なのである。
（念のため付け加えると、もちろん「t値が大きく、めったに起こらないことが発生した」といった表現は問題なく、「棄却できるか否か」について定量的な表現がよくない、ということである。最初に設定した有意水準を超えるか否かだけが焦点であり、5%と設定したのであれば、P値が0.04であろうと0.000001であろうと等しく「棄却できる」のである。）
Ad-Stock効果を推定しつつ回帰を回したい
最近ずっとMarketing Mix Modeling（MMM）をやっている。その中で広告効果（いわゆるROI）を推定しているのだけれど、広告の効果というものは出稿した時点だけでなく将来に渡って影響を及ぼすため、過去の広告の累積による影響（いわゆる残存効果・Ad-Stock効果）を取り上げる必要がある。ところが、このAd-Stock効果をモデルにより推定しようとすると少し難しい。
Ad-Stock効果として以下のような分布ラグ型：
を仮定するのが最も単純なのだけれど、考慮するラグの次数mによってはモデルに投入する変数の数が多くなり、また多重共線性も強くなりがち。時系列データを用いているためデータポイントが少なく、このような方法は避けたい。
そこで「m期前の広告効果はλのm乗に従って減衰する」という仮定をおき（いわゆるKoyck型）、推定すべきパラメータの数を減らしている。この場合変数の数は減らせるが、非線形な効果であるλをどのように推定するかが鍵となる。
考えられる方法として、
といったものが思いついた。しかし1は季節性や他の広告の効果などを同時に考慮できないため✖️、2は使えそうな数値がなかった。そこで今回は回帰の中でλを同時に推定できないか、以下のように調査してみた。
set.seed(123) # λに従って広告の効果が減衰するデータを生成する関数 simulate_y_01 <- function(pars) { X1 <- rnorm(100, 100, 2) * rep(c(0, rep(1, 7), 0, 1), 10) X2 <- filter(X1, pars[1], "recursive") Z1 <- rnorm(100, 5, 2) * rep(c(rep(0, 5), rep(1, 5)), 10) dat <- data.frame( "Y" = 50 + pars[2] * X2 + pars[3] * Z1 + rnorm(100, 0, pars[4]), "X1" = X1, "Z1" = Z1) return(dat) } # λを与えるとAd-Stock変数を作成して回帰を実施し、AICを返す関数 return_AIC_01 <- function(a) { dat$Tmp <- filter(dat$X1, a, "recursive") AIC(lm(Y ~ Tmp + Z1, dat)) }
Koyck型の変数の作成にはstats::filterが便利だが、使う際には{dplyr}を読み込んでいないか注意しよう。これでreturn_AIC_01を最小化した時に得られるλが、simulate_y_01で用いた数値を一致するか確認する。実際にやってみると以下の結果が得られる。
## λを小さな値とした場合 pars <- c(0.1, 0.5, 4, 5.0) n <- 500 res_01_01 <- matrix(0, n, 4) for (i in 1:n) { dat <- simulate_y_01(pars) a <- optimise(return_AIC_01, c(0, 1))$minimum X2 <- filter(dat$X1, a, "recursive") res_01_01[i, ] <- c(a, coef(lm(Y ~ X2 + Z1, data = dat))) } ## λを大きな値とした場合 pars <- c(0.9, 0.5, 4, 5.0) n <- 500 res_01_02 <- matrix(0, n, 4) for (i in 1:n) { dat <- simulate_y_01(pars) a <- optimise(return_AIC_01, c(0, 1))$minimum X2 <- filter(dat$X1, a, "recursive") res_01_02[i, ] <- c(a, coef(lm(Y ~ X2 + Z1, data = dat))) }
> summary(res_01_01) V1 V2 V3 V4 Min. :0.01545 Min. :43.81 Min. :0.4633 Min. :3.486 1st Qu.:0.08780 1st Qu.:48.73 1st Qu.:0.4918 1st Qu.:3.877 Median :0.10450 Median :49.90 Median :0.4998 Median :3.999 Mean :0.10324 Mean :49.88 Mean :0.4996 Mean :3.997 3rd Qu.:0.11992 3rd Qu.:51.02 3rd Qu.:0.5081 3rd Qu.:4.109 Max. :0.17782 Max. :55.08 Max. :0.5408 Max. :4.474 > summary(res_01_02) V1 V2 V3 V4 Min. :0.8972 Min. :40.15 Min. :0.4768 Min. :3.414 1st Qu.:0.8992 1st Qu.:47.93 1st Qu.:0.4946 1st Qu.:3.869 Median :0.9000 Median :50.04 Median :0.4995 Median :3.994 Mean :0.9000 Mean :49.94 Mean :0.5002 Mean :3.999 3rd Qu.:0.9007 3rd Qu.:52.10 3rd Qu.:0.5060 3rd Qu.:4.133 Max. :0.9031 Max. :58.68 Max. :0.5236 Max. :4.526
λが小さい時にはある程度バラついているものの、この誤差なら十分に実用的に思える。λが大きい時は全く問題がなさそうだ。
少し問題を難しくするため、2つの変数を共にAd-Stockに変更すると：
simulate_y_02 <- function(pars) { X1 <- rnorm(100, 100, 2) * rep(c(0, rep(1, 7), 0, 1), 10) X2 <- filter(X1, pars[1], "recursive") Z1 <- rnorm(100, 5, 2) * rep(c(rep(0, 5), rep(1, 5)), 10) Z2 <- filter(Z1, pars[2], "recursive") dat <- data.frame( "Y" = 50 + pars[3] * X2 + pars[4] * Z2 + rnorm(100, 0, pars[5]), "X1" = X1, "Z1" = Z1) return(dat) } return_AIC_02 <- function(param) { a <- param[1] b <- param[2] dat$TmpX <- filter(dat$X1, a, "recursive") dat$TmpZ <- filter(dat$Z1, b, "recursive") AIC(lm(Y ~ TmpX + TmpZ, dat)) } ## λを小さな値とした場合 pars <- c(0.1, 0.1, 0.5, 4, 5.0) n <- 100 res_02_01 <- matrix(0, n, 5) for (i in 1:n) { dat <- simulate_y(pars) res <- optim(par = optim(par = init_par, fn = return_AIC)$par, fn = return_AIC)$par X2 <- filter(dat$X1, res[1], "recursive") Z2 <- filter(dat$Z1, res[2], "recursive") res_02_01[i, ] <- c(res, coef(lm(Y ~ X2 + Z2, data = dat))) } ## λを大きな値とした場合 pars <- c(0.8, 0.7, 0.5, 4, 5.0) n <- 100 res_02_02 <- matrix(0, n, 5) for (i in 1:n) { dat <- simulate_y(pars) res <- optim(par = optim(par = init_par, fn = return_AIC)$par, fn = return_AIC)$par X2 <- filter(dat$X1, res[1], "recursive") Z2 <- filter(dat$Z1, res[2], "recursive") res_02_02[i, ] <- c(res, coef(lm(Y ~ X2 + Z2, data = dat))) }
計算に少し時間がかかったので、反復は100回とした。結果は以下の通り。
> summary(res_02_01) V1 V2 V3 V4 V5 Min. :0.03261 Min. :-0.02245 Min. :44.54 Min. :0.4685 Min. :3.597 1st Qu.:0.08483 1st Qu.: 0.06107 1st Qu.:48.62 1st Qu.:0.4899 1st Qu.:3.902 Median :0.09598 Median : 0.09694 Median :50.23 Median :0.4997 Median :3.997 Mean :0.09810 Mean : 0.09459 Mean :50.06 Mean :0.4997 Mean :4.017 3rd Qu.:0.11280 3rd Qu.: 0.11921 3rd Qu.:51.30 3rd Qu.:0.5081 3rd Qu.:4.122 Max. :0.14426 Max. : 0.24930 Max. :55.33 Max. :0.5419 Max. :4.606 > summary(res_02_02) V1 V2 V3 V4 V5 Min. :0.7894 Min. :0.6499 Min. :41.16 Min. :0.4630 Min. :3.701 1st Qu.:0.7961 1st Qu.:0.6833 1st Qu.:48.19 1st Qu.:0.4905 1st Qu.:3.906 Median :0.7990 Median :0.7003 Median :50.51 Median :0.5016 Median :3.984 Mean :0.7997 Mean :0.7000 Mean :49.99 Mean :0.5006 Mean :3.994 3rd Qu.:0.8030 3rd Qu.:0.7140 3rd Qu.:52.09 3rd Qu.:0.5093 3rd Qu.:4.065 Max. :0.8137 Max. :0.7496 Max. :57.14 Max. :0.5257 Max. :4.287
先ほどと変わらず、λが小さいと推定の精度が落ちるようだ。また一部のデータセットではλが負となってしまった。
ついでなので、{rBayesianOptimization}パッケージによる最適化を通常のoptimと比較してみた。
library(rBayesianOptimization) # BayesianOptimization関数では最大化を行うので、返り値となるAICに負の符号を加えてある return_AIC_BO <- function(a, b) { dat$TmpX <- filter(dat$X1, a, "recursive") dat$TmpZ <- filter(dat$Z1, b, "recursive") out <- list(Score = -AIC(lm(Y ~ TmpX + TmpZ, dat)), Pred = -AIC(lm(Y ~ TmpX + TmpZ, dat))) return(out) } pars <- c(0.1, 0.1, 0.5, 4, 5.0) dat <- simulate_y(pars) res_BO <- BayesianOptimization(return_AIC_BO, bounds = list(a = c(0, 1), b = c(0, 1)), init_points = 20, n_iter = 1, acq = "ucb", kappa = 2.576, eps = 0, verbose = TRUE) init_par <- array(c(0.5, 0.5)) res_Opt <- optim(par = optim(par = init_par, fn = return_AIC)$par, fn = return_AIC)$par
> res_BO$Best_Par a b 0.8274328 0.4757574 > res_Opt [1] 0.7959288 0.7930622
BayesianOptimizationの結果はあまり良くなく、通常のoptimの方が精度よく推定できた。以前の記事（optimってあんまり信用できないなぁ、って話 - 統計コンサルの議事メモ）ではoptimによるパラメータの推定はうまくいかなかったのだけど、lmと組み合わせているのが良いのかもしれない。
この方法でならAd-Stock効果の推定できそうだ。optimでは"L-BFGS-B"を指定すれば上限・下限を指定できたと思うので、もう少し検証してみよう。
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"と呼ばれるアルゴリズムを実装した関数はかなり強力な感触を受けた。こちらの検証結果も後日書いてみようと思う。
