大規模データマイニングでのモデル探索手法：K-sample plot
巨大地震が日本を襲い、皆不安を感じながら生活していると思います。
そんな中せめて自分に出来ることをしようと思い、研究してきた内容をブログに記します。
サンプル数が大規模なデータでニューラルネットワークとかサポートベクターマシンとかをしたくても、
時間がかかってしょうがない！ってときに参考にしてみて下さい。
近年、特にweb関係の業界ではデータデータをいくらでも記録できるようになったため、データの規模が非常に大きくなっています。
大規模データに機械学習、マシーンラーニングを適用したいという要望は高まっています。
そういうときはデータからサンプリングして性能を確かめることが多いと思います。
ですがそんな時は、
- 「サンプル数はどれくらいがいいの？」
- 「一回やっただけじゃ真の性能は分からないよね？」
- 「しかもクロスバリデーションしなきゃいけないし。。」
などのような事を疑問に思うでしょう。
今回紹介するK-sample plot（以下、K's plot)のアイディアはとても単純で、「サンプル数を少しずつ増やしながら性能をプロットする」というものです。
これから事例を紹介します。
あるデータがあって、結果変数yと説明変数x1とx2が100万サンプル（人）であったとします。
yとx1は線形な関係、yとx2は2次的な関係があったとし、次のように乱数発生させます。
set.seed(1) x1 <- rnorm(1000000) set.seed(2) x2 <- rnorm(1000000) set.seed(3) y <- 2*x1 + x2**2 + rnorm(1000000) Data <- data.frame(x1 = x1, x2 = x2, y = y)
このデータに対して以下の3種類のモデルを仮定して適用してみます。
モデル1とモデル2は線形モデルなので、100万サンプル全体を使ってもそんなに時間はかかりません。
だいたい6〜7秒くらいでした。
AllModel1 <- lm(y~., data = Data) #6.24second summary(AllModel1)
しかし、モデル3はSVMなのでとても時間がかかります。
AllModel3 <- svm(y~., data = Data)
今30分くらい経ってますがまだ終わってませんorz
大規模データに対して不用意に機械学習を実行すると悲惨です。。。
そこで、K's plotを行ってみます。
天下り式ですが、結果からお見せします。
100万サンプルから10〜1,000サンプルを抜き出し、モデル1〜3を当てはめてR2乗を計算しています。
そしてR2乗とサンプル数の関係をプロットしています。
library(e1071) #SVM X1 <- data.frame(x1 = x1, x2 = x2) X2 <- data.frame(x1 = x1, x2 = x2, x3 = x2**2) y <- y Ksample1 <- c(seq(10,100,10), seq(100,1000,100)) KsResult1 <- KsamplePlot(X1, y, Ksample1) #2.39second KsResult2 <- KsamplePlot(X2, y, Ksample1) #2.39second KsResult3 <- KsamplePlot(X1, y, Ksample1, Method = "svm") #2.45second
3つのモデルのK's plotの結果です。
横軸がサンプル数で、縦軸が予測性能です（plateauが推定値）。
プロットの数だけ予測を行っていますので、SVMも何十回とやっているのです。
1. 線形モデル ：y=x1+x2
2. 2次線形モデル：y=x1+x2+x2^2
どのグラフも計算に2〜3秒しかかかってません！
サンプル数は10〜1000の間で変化させています。10〜100までは10ずつ、100〜1000までは100ずつ変化させています。
予測性能は1-R二乗なので、低いほど良いモデルになります。
このプロットは検出力とサンプル数の関係を擬似的に表したものであり、サンプルが増えるとどのように予測性能が上がるかを視覚的にチェックできます。
この3枚のプロットを見ると、モデル2が最も良いことが分かります。
モデル2はサンプルが200あたりで予測性能が一定になっているように思います。
しかしモデル3はまだ一定ではなく、まだ下がっている途中のようにも見えます。
そこでK-sampleの数を2000まで増やしてモデル3のK's plotを描いてみます。
Ksample2 <- c(seq(10,100,10), seq(100,1000,100), seq(1000,2000,100)) KsResult3 <- KsamplePlot(X1, y, Ksample2, Method = "svm") #7.54second
予想通り予測性能が少し上がりました。
しかし、モデル2には届かないようです。
このデータの場合はSVMより2次項を入れたモデルの方が良いということが示唆できます。
もう一つ分かる事は、このデータではSVMの方が予測性能の収束が遅いという事。
予測性能がいつ収束するかというポイントは、モデル探索する時に非常に興味がある点です。
場合によっては"elbow point"の推定にも使えるのではないかと思います。
さて、ここからK's plotのプログラムを紹介します。
KsamplePlot <- function(X, y, Ksample = c(seq(40,100,10), 150, seq(200,1000,100)),
Method = "lm", Caret = "No", size = 5, Type = "numeric"){
library(e1071)
library(caret)
library(nnet)
library(MASS)
library(caTools)
library(mda)
library(glmnet)
library(randomForest)
Rsq <- rep(NA, length(Ksample)*3)
for(j in 1:3){
for(i in 1:length(Ksample)){
ksample <- Ksample[i]
Sample <- sample(1:nrow(X), ksample)
KSamplex <- X[Sample, ]
KSampley <- y[Sample]
#---Half CV
Train <- sample(1:ksample, ksample / 2)
Trainx <- KSamplex[Train, ]
Trainy <- KSampley[Train]
TrainData <- data.frame(Trainx, y = Trainy)
Varidx <- KSamplex[-Train, ]
Varidy <- KSampley[-Train]
VaridData <- data.frame(Varidx, y = Varidy)
#------numeric
if(Type == "numeric"){
#---develop model
if(Caret == "No"){
if(Method == "lm") Model <- lm(y~., data = TrainData); Varidx <- data.frame(Varidx)
if(Method == "svm") Model <- svm(y~., data = TrainData)
if(Method == "nn") Model <- nnet(y~., data = TrainData, size = size, linout = T)
if(Method == "rf") Model <- randomForest(y~., data = TrainData)
if(Method == "mars") Model <- mars(Trainx, Trainy)
#if(Method == "bruto") Model <- bruto(Trainx, Trainy); Varidx <- do.call(cbind, Varidx)
if(Method == "cart") Model <- rpart(y~., data = TrainData)
if(Method == "lasso"){
cvob1 <- cv.glmnet(do.call(cbind, Trainx), Trainy)
CvLambda <- cvob1$lambda[cvob1$cvm == min(cvob1$cvm)]
Model <- glmnet(do.call(cbind, Trainx), Trainy, lambda = CvLambda)
Varidx <- do.call(cbind, Varidx)
}
}
if(Caret == "Yes"){
modeltxt <- parse(text=paste("train(y~., data=TrainData, method= Method)", sep=""))
Model <- eval(modeltxt)
}
Err <- as.numeric(Varidy - predict(Model, Varidx))
MSE <- Err%*%Err
Var <- (Varidy-mean(Varidy))%*%(Varidy-mean(Varidy))
Rsq[i + length(Ksample)*(j-1)] <- (Var - MSE)/Var
}
#------class
if(Type == "binary"){
#---develop model
if(Caret == "No"){
if(Method == "lm"){
Model <- lda(y~., data = TrainData)
Post <- predict(Model, VaridData)$posterior[, 2]
ROC <- caTools::colAUC(Post, Varidy)
Rsq[i + length(Ksample)*(j-1)] <- as.numeric(ROC)
}
if(Method == "svm"){
Model <- svm(y ~ ., data = TrainData, kernel="polynomial", degree=3, probability = T)
Post <- predict(Model, VaridData, probability = T)
ROC <- caTools::colAUC(Post, Varidy)
Rsq[i + length(Ksample)*(j-1)] <- as.numeric(ROC)
}
}
if(Caret == "Yes"){
}
}
print(paste("Ksample", j, " ", ksample, sep = ""))
}
}
Ks <- cbind(Ksample, 1 - Rsq)
Ks <- Ks[(0 <= Ks[, 2] & Ks[, 2] <= 1), ]
EstPar <- KsPlot(Ks, Method = Method, Type = Type)
list(OneExpPar = EstPar, Ksample = cbind(Ksample, 1 - Rsq))
}
KsPlot <- function(UseData, Type = "numeric", Method){
Par <- data.frame(k = rep(seq(0.01, 1, 0.01), 100), plateau = rep(seq(0.01, 1, 0.01), rep(100, 100)))
Object <- rep(NA, 10000)
for(i in 1:10000){
k <- Par[i, 1]
plateau <- Par[i, 2]
x <- UseData[, 1]
y <- UseData[, 2]
if(Type == "numeric") yhat <- (1 - plateau)*exp(-k*x) + plateau
if(Type == "binary") yhat <- (0.5 - plateau)*exp(-k*x) + plateau
Object[i] <- sum*1, lwd = 2, ylim=c(0, 1), xlab="K-Sample", ylab="1 - Explained Var", main = paste("Method = ", Method, sep = ""))
if(Type == "binary") plot(UseData, xlim = c(0, max(UseData[, 1])), lwd = 2, ylim=c(0, 1), xlab="K-Sample", ylab="1 - AUC", main = paste("Method = ", Method, sep = ""))
k <- as.numeric(EstPar[1])
plateau <- as.numeric(EstPar[2])
XLim <- max(UseData[, 1])
if(Type == "numeric") Ks.Est <- (1 - plateau)*exp(-k*c(1:XLim)) + plateau
if(Type == "binary") Ks.Est <- (0.5 - plateau)*exp(-k*c(1:XLim)) + plateau
lines(c(1:XLim), Ks.Est, lwd = 2)
text(XLim*0.7, 0.8, paste("k (slope):", round(EstPar[1], digits=2)), adj=0)
text(XLim*0.7, 0.7, paste("plateau :", round(EstPar[2], digits=2)), adj=0)
EstPar
}
- KsPlot
サンプルデータの予測データを入れるとK's plotを描く
下のKsamplePlot関数に埋め込まれているので直接は使わない
K-sampleと予測性能に対して指数モデルを当てはめて、予測性能の上がり方をチェックする
- KsamplePlot
K's plotの本体関数。
データから、サンプリング、クロスバリデーション、グラフ描画を行う。
X：説明変数行列（2変数以上）
y：結果変数ベクトル
Ksample：サンプリング数の指定
Method：回帰方法の指定、デフォルトは線形回帰、今はlmとsvmのみ実装
Medhodの拡張は簡単で、「develop model」のところの「if(Method == "lm") Model <- lm(y~., data = TrainData)」を追加するだけです。
下の行のpredictに繋げれればいいので、ランダムフォレストも追加できます。
クロスバリデーションはHalf-CV（2-fold CVの片方だけ）を3回行っています。
取り急ぎ完成しているのはここまでです。
まだデータによってはたまにエラーが出たりするので完璧じゃないのですが、使ってみて下さい。
この手順で予測モデルの構築や選択を行えば、大規模データでも素早く探索できると思います。
現在プログラムにしているのは連続値の予測だけですが、2値予測の場合はAUCを指標にすることで同じロジックを使えます。
説明変数行列と結果変数ベクトルさえあればKsamplePlot関数が使えますので、是非使ってみて下さい。
これ、実は今年度出したD論の一部なんです・・……(-｡-) ボソッ
また拡張したら記事にします。
それでは。
修正１
- KsamplePlot関数の中にミスがあったので修正。
- 連続値のオプションにニューラルネットを追加
- パッケージcaretを使えるように修正
- caret = "Yes"、Model = "使用したモデル名"で利用できる
- 2値予測も可能
修正２
- Epi::ROCをcaTools::colAUCに修正
- sampleのサンプリングによってエラーが出る
- エラーは2種類
- yが2グループ共少なくとも1つはサンプルされてないといけない
- svd(X, nu = 0)のエラーは不明
修正３
- MARSの追加（brutoは動かすとRがクラッシュするので追加できず）
- LASSOの追加（CVして最適なラムダを推定する、Ksampleは最低でも40は必要）
- ランダムフォレストの追加
動作報告があったら教えて頂けると幸いです。
サンプルプログラム
set.seed(1) x1 <- rnorm(1000000) set.seed(2) x2 <- rnorm(1000000) set.seed(3) y <- 2*x1 + x2**2 + rnorm(1000000) X1 <- data.frame(x1 = x1, x2 = x2) X2 <- data.frame(x1 = x1, x2 = x2, x3 = x2**2) y <- y KsResult1 <- KsamplePlot(X1, y) KsResult2 <- KsamplePlot(X2, y)
*1:y - yhat)**2) } EstPar <- Par[Object == min(Object), ][1, ] if(Type == "numeric") plot(UseData, xlim = c(0, max(UseData[, 1]
