random forestの感度分析+
教科書を読むと、random forestでモデルを学習させて、ある良い予測精度のものが得られたら、まずは各変数の重要度をみて、つぎにpartial dependency partial dependence plot(PDP)を眺めなさい、と書いてある。ある変数の動きに注目して予測値がどう変化するか見てみよう、という分析だ。
これを発展(?)させたforestFloorは、決定木の各ノードでの予測値の変化量に着目して、random forestにおける変数が動いたとき予測値がどちら向きにどのくらい動くか定量化しようというもの。
スライド中のコードは、ほとんどマニュアルの丸写しなので割愛する。
forestFloorでは、ある決定木でのある事例の予測値(もしくは、あるクラスへの所属確率)= ブートストラップ平均 + ノード1での変化量 + ... + 終端ノードの変化量、と解釈する。この解釈の下で、ブートストラップによって当該の観察事例が出てこない木や、観察事例が終端ノードに到達するまでに注目する変数を使わなかった、というケースは除外して、注目する変数による予測値の変化量をすべての木について足し合わせることで、ある観察事例に対するある変数の貢献度(feature contribution)を計算する。
このアイデアが面白いのは、回帰(→予測値)であろうが分類(→あるクラスへの所属確率)であろうが、弱学習器に決定木を使っているアンサンブル手法であれば貢献度が計算できる点だと思う。 弱点は、論文を眺めた限りでは、他の変数の影響を排除できているのかよく判らない点と、「それで結局予測はどうなるんだ」と訊かれそうな点と、原理上アンサンブルツリー縛りがある点。
とくに後者2つで、一般性(や有用性)という観点からPDP&ICEboxに譲ってしまうと思われる。Hastieたちの教科書では、partial dependenceは任意の教師あり学習に対して、ある変数についてのmarginal averageが計算できると説明されている。これに加えて、ICEBoxを使うと予測値がの動きを分解(disaggregate)でき、かつ、予測値が大きく変化する変数値のポイントを検出できる。一応、forestFloorでも似たようなことはできるし、論文でもICEboxについて言及してはいるのだが、若干歯切れの悪い流れであった。。
ちなみに、これらはmlRにて任意の教師あり学習に適用可能なwrapperとして実装されていて、下記のリンク先ではSVMをつかった謎のデータirisの3クラス分類での例が示されている。
Partial Prediction Plots - mlr tutorial
(追記)
上記のリンクがいろいろと切れてると思ったら、mlr 2.8→2.9で関数名の変更があったようだ。
LassoとOracle Property
「Oracle Property」。巷にこれほどかっこいいフレーズはそうそうないと思うので紹介を試みる。
さて、Lassoの罰則項はどういう作用を持っているか? L1 normの縮小推定を解説する図はいろんな教科書にも載っているのだが、自分は最初見たときに実感を持てなかった。Rで実際に条件を変えながらOLSと比較することで、一助になれば幸いである。
www.slideshare.net
今回使用したコード。 foreachのブロックは後半で結構時間がかかるので、イテレーションを減らすか並列化されたい。
# comparison among lm, lasso, adaptive lasso, and rigorous lasso require(dplyr) require(foreach) #install.packages("glmnet") require(glmnet) #install.packages("msgps") require(msgps) #install.packages("hdm") require(hdm) # functions --------------------------------------------------------------- # generate data according to model: y ~ 0 + B_i*x_i + N(0,1) Data_with <- function(n,p,s, silent=TRUE){ # data with sample size(n), number of variables(p) and non-zero vector(s) if(! silent){ cat(paste0(n, " obs. of ", p, " variables.\n", "Last ", length(s)," variables are non-zero.")) } beta <- c(rep(0,p-length(s)),s) X <- matrix(rnorm(n*p), ncol=p) res <-list( X = X, y = 0 + X%*%beta + rnorm(n), beta = beta, non.zero.coeff = which(beta>1e-8), non.zero.value = s ) if(! silent){ str(res) } invisible(res) } # effect of shrinkage on estimates (lm vs lasso) ---------------------------------- test_data <- Data_with(n=10^3, p=10, s= c(1,10,100)) str(test_data,max.level = 2) attach(test_data) system.time( lm.reg <- lm(y~X) ) summary(lm.reg) system.time( lasso_cv <- cv.glmnet(x = X, y=y, alpha = 1) ) coef(lasso_cv, s="lambda.1se") test_data <- Data_with(n=10^3, p=10, s= c(100,100,100)) str(test_data,max.level = 2) attach(test_data) system.time( lm.reg <- lm(y~X) ) summary(lm.reg) system.time( lasso_cv <- cv.glmnet(x = X, y=y, alpha = 1) ) par(mfrow=c(1,2)) plot(lasso_cv) abline(v=log(lasso_cv$lambda.1se), lty=2, col="blue") # plot solution path (Lasso) ------------------------------------------------------ lasso.reg <- glmnet(x = X, y=y, alpha = 1) plot(lasso.reg, xvar="lambda", label=TRUE) abline(v=log(lasso_cv$lambda.1se), lty=2, col="blue") abline(h= c(1,2,3), lty=2, col= c(1,2,3)) par(mfrow=c(1,1)) fit.coeff <- coef(lasso_cv, s="lambda.1se") %>% summary %>% transmute(beta=i-1, Var=x) coef(lasso_cv, s="lambda.1se") fit.coeff %>% str # using lasso with oracle property (adaptive lasso) ------------------------- system.time( alasso.reg <- msgps(X= X, y=c(y), penalty="alasso",gamma=1, lambda=0) ) par(mfrow=c(2,2)) plot(alasso.reg, stand.coef=FALSE, criterion="cp", main="Mallows' Cp") abline(h= c(1,2,3), lty=2, col= c(1,2,3)) plot(alasso.reg, stand.coef=FALSE, criterion="aicc",main="bias-corrected AIC (AICc)") abline(h= c(1,2,3), lty=2, col= c(1,2,3)) plot(alasso.reg, stand.coef=FALSE, criterion="gcv", main="generalized cross validation (GCV) ") abline(h= c(1,2,3), lty=2, col= c(1,2,3)) plot(alasso.reg, stand.coef=FALSE, criterion="bic", main="BIC") abline(h= c(1,2,3), lty=2, col= c(1,2,3)) par(mfrow=c(1,1)) summary(alasso.reg) ss=300 # measure: estimate error vs sample size ------------- sample_sizes <- c(5:20*10, 400, 600, 800, 1000,2000,3000, 4000) itr <- 1:10 conditions <- expand.grid(ss=sample_sizes, itr=itr) %>% arrange(ss) sse.sampsi <- foreach(Ci=1:NROW(conditions), .combine = rbind) %do% { cond <- conditions[Ci,] this_data <- Data_with(n=cond$ss, p=10, s= c(10,10,10), silent=TRUE) beta <- this_data$beta fit.lm <- lm(this_data$y~this_data$X) sse.lm <- sum( (beta - coef(fit.lm)[-1])^2 ) fit.lasso_cv <- cv.glmnet(x = this_data$X, y=this_data$y, alpha = 1) coeff.lasso_cv <- coef(fit.lasso_cv, s="lambda.1se") %>% summary %>% transmute(beta=i-1, Var=x) sse.lasso_cv <- sum( (beta[coeff.lasso_cv$beta[-1]] - coeff.lasso_cv$Var[-1])^2 ) fit.alasso <- msgps(X= this_data$X, y=c(this_data$y), penalty="alasso", gamma=1, lambda=0) coeff.alasso <- coef(fit.alasso) %>% data.frame %>% select(BIC) %>% unlist sse.alasso <- sum( (this_data$beta - coeff.alasso[-1])^2 ) data.frame(sample.size=cond$ss, lm=sse.lm, lasso=sse.lasso_cv, alasso=sse.alasso) } sse.summ.sampsi <- sse.sampsi %>% group_by(sample.size) %>% summarize_each(funs(mean)) sse.summ.sampsi with(sse.summ.sampsi,{ plot (lm ~sample.size, type="l", lwd=2, col="green", ylim=c(0,0.35), ylab="mean SSE", log="x") lines(lasso~sample.size, type="l", lwd=2, col="red") legend("topright", legend = c("lm","lasso"), col = c("green","red"), lty = 1, lwd=2) }) with(sse.summ.sampsi,{ plot (lm ~sample.size, type="l", lwd=2, col="green", ylim=c(0,0.35), ylab="mean SSE", log="x") lines(lasso ~sample.size, type="l", lwd=2, col="red") lines(alasso~sample.size, type="l", lwd=2, col="blue") legend("topright", legend = c("lm","lasso","a-lasso"), col = c("green","red","blue"), lty = 1, lwd=2) }) # example(rlasso) --------------------------------------------------------- # use Lasso with rigorous penalty system.time( nonpost_lasso.reg <- rlasso( formula, data, post = FALSE, # FALSEで、罰則だけrigorousなLasso intercept = TRUE, # TRUEで切片項(β0)の罰則付き推定をしない penalty = list(homoscedastic = FALSE, # FALSEで誤差の等分散性を仮定しない X.dependent.lambda = FALSE, # TRUEでデータ依存的なλの選択をする lambda.start = NULL, c = 1.1, gamma = 0.1/log(n)), control = list(numIter =15, tol = 10^-5, threshold = NULL), ...) ) # use Post-Lasso procedure with rigorous penalty system.time( post_lasso.reg <- rlasso( formula, data, post = TRUE, # TRUEでPost-Lasso+rigorousな罰則 intercept = TRUE, # TRUEで切片項(β0)の罰則付き推定をしない penalty = list(homoscedastic = FALSE, # FALSEで誤差の等分散性を仮定しない X.dependent.lambda = FALSE, # TRUEでデータ依存的なλの選択をする lambda.start = NULL, c = 1.1, gamma = 0.1/log(n)), control = list(numIter =15, tol = 10^-5, threshold = NULL), ...) ) # comparison results (set all=FALSE to summarize) summary(nonpost_lasso.reg, all=FALSE) summary(post_lasso.reg, all=FALSE) coef(post_lasso.reg) # measure: estimate error vs dimension size ------------- dim_sizes <- c(1:10*10, 1:5*200, 2:5*1000, 1e4, 2e4) itr <- 1:30 conditions <- expand.grid(ds=dim_sizes, itr=itr) %>% arrange(ds) sse.dim_sizes <- foreach(Ci=1:NROW(conditions), .combine = rbind) %do% { cond <- conditions[Ci,] this_data <- Data_with(n=1000, p=cond$ds, s= c(10,10,10), silent=TRUE) beta <- this_data$beta fit.lm <- lm(this_data$y~this_data$X) sse.lm <- sum( (beta - coef(fit.lm)[-1])^2 ) fit.lasso_cv <- cv.glmnet(x = this_data$X, y=this_data$y, alpha = 1) coeff.lasso_cv <- coef(fit.lasso_cv, s="lambda.1se") %>% summary %>% transmute(beta=i-1, Var=x) sse.lasso_cv <- sum( (beta[coeff.lasso_cv$beta[-1]] - coeff.lasso_cv$Var[-1])^2 ) fit.alasso <- msgps(X= this_data$X, y=c(this_data$y), penalty="alasso", gamma=1) coeff.alasso <- coef(fit.alasso) %>% data.frame %>% select(BIC) %>% unlist sse.alasso <- sum( (this_data$beta - coeff.alasso[-1])^2 ) plasso.reg <- rlasso(c(this_data$y) ~ this_data$X, post=TRUE, intercept=TRUE) coeff.plasso <-coef(plasso.reg) sse.plasso <-(this_data$beta - coeff.plasso)^2 %>% sum data.frame(dim.size=cond$ds, lm=sse.lm, lasso=sse.lasso_cv, alasso=sse.alasso, plasso=sse.plasso) } sse.summ.dim_sizes <- sse.dim_sizes %>% group_by(dim.size) %>% summarize_each((funs(mean))) sse.summ.dim_sizes with(sse.summ.dim_sizes,{ plot( lm ~dim.size, type="l", lwd=2, col="green", ylab="mean SSE", log="x") lines(lasso ~dim.size, type="l", lwd=2, col="red") lines(alasso~dim.size, type="l", lwd=2, col="blue") legend("topright", legend = c("lm","lasso","a-lasso"), col = c("green","red","blue"), lty = 1, lwd=2) abline(v=1000, lty=3) }) with(sse.summ.dim_sizes,{ plot(lasso~dim.size, type="l", lwd=2, col="red", ylim=c(0,0.1), ylab="mean SSE",log="x") lines(alasso~dim.size, type="l", lwd=2, col="blue") legend("topleft", legend = c("lasso","a-lasso"), col = c("red","blue"), lty = 1, lwd=2) abline(v=1000, lty=3) }) with(sse.summ.dim_sizes,{ plot(lasso~dim.size, type="l", lwd=2, col="red", ylim=c(0,0.1), ylab="mean SSE",log="x") lines(alasso~dim.size, type="l", lwd=2, col="blue") lines(plasso~dim.size, type="l", lwd=2, col="orange") legend("topleft", legend = c("lasso","a-lasso","post-lasso"), col = c("red","blue","orange"), lty = 1, lwd=2) abline(v=1000, lty=3) }) # measure: performance vs dimension size ------------- this_data <- Data_with(n=1000, p=2e4, s= rep(10,3),silent = FALSE) system.time( plasso.reg <- rlasso(c(this_data$y) ~ this_data$X, post=TRUE, intercept=TRUE) ) #coef(fit.lasso_cv) system.time( fit.alasso <- msgps(X= this_data$X, y=c(this_data$y), penalty="alasso", gamma=1) )
反省点その1、これほどかっこいいフレーズなのに当日聴衆には全く刺さらなかった(というか完全に滑った)。プレゼンスキルの改善がである。その2、圧縮センシングのお得感とか、スパース推定後の解釈性とか、もうちょっと使い道に触れた導入にすべきだった。最後の"rigorous lasso"では、計算時間は少なからず増えているわけで、実際の現場感としてここまでの精度が必要か?とか、いずれ費用対効果やロバスト化の議論はしてみたい。
random forestを予測以外の目的で使う
数年ほど前には最強と言われて一世を風靡していたrandom forestだが、予測以外にも使い道が提案されている。Rのパッケージから紹介したい。
予測全体の把握と仮説ルールの抽出
決定木分析が便利な理由の一つは「どういうルールでその予測が成り立っているのか」を極めて簡単に可視化出来る点。inTreesパッケージは予測ルールを集計し、適当に枝狩りして全体を要約することで、アンサンブルモデルにおいても決定木と同じような情報を可視化してくれる。
このパッケージではもう一つ、すべての木から取り出した1つずつの枝をトランザクションとみなしてアソシエーション分析する機能を提供する。メジャーな関連ルールを取り出す方に力点が置かれているが、変数同士の稀な関係が意外なアウトカムを導いてる方が、(研究上は)見つかったときに面白いので、欲を言えばLift値も評価してくれると良かった。
ノンパラメトリックな欠損値の補完
randomForestパッケージには、rfImputeという欠損値を補完する関数が用意されている。ところが、rfImputeは同じルールにより同じ予測をされる別の観察データを使って補完するため、目的変数が欠損しているデータには適用できない。
既知のデータにおいて欠測部分を推定するにはこれで構わないのだが、予測をしたいデータの前処理に使おうとするとこれでは困る。
missForestでは、補完対象の変数を目的変数とみなすことで、random forestにより欠損値を直接予測する。
もともと予測したい変数も補完プロセスに含めることができるので半教師学習と言える。その場合、イテレーションの回数と過学習の関係が気になるが、元の論文にはその言及はなかったように思う。
“rigorous” lassoを試してみた
High-Dimensional Metrics in R
CRAN - Package hdm
- やりたかったこと= oracle property保証のlassoで変数選択(n<p)
- 得たもの= 速い、一致性はそれっぽい。
- やりたいこと= ほかのパッケージとの比較、証明部分読む。
以下、試したかった部分だけ:
# install.packages("hdm") require(hdm) set.seed(1) ### modified from example(rlasso) ---------------- n <- 100 # sample size p <- 1500 # number of variables s <- 3 # nubmer of non-zero variables beta = c(rep(3,s), rep(0,p-s)) # s -> non-zero/zero coeff.vec X <- matrix(rnorm(n*p), ncol=p) y <- 1 + X%*%beta + rnorm(n) # use Lasso,not-Post-Lasso system.time( lasso.reg <- rlasso(y~X,post=FALSE, intercept=TRUE) ) summary(lasso.reg, all=FALSE) # use this option to summarize results # use Post-Lasso system.time( lasso.reg <- rlasso(y~X, post=TRUE, intercept=TRUE) ) summary(lasso.reg, all=FALSE) # use this option to summarize results yhat.lasso <- predict(lasso.reg) #in-sample prediction Xnew <- matrix(rnorm(n*p), ncol=p) # new X Ynew <- Xnew%*%beta + rnorm(n) #new Y yhat.lasso.new = predict(lasso.reg, newdata=Xnew) #out-of-sample prediction
オペレーションズ・リサーチ(Vol.58, No.5, 261)の解説記事にいろんなlassoとRの対応パッケージ表があるので、参考例のデータは微妙だがadaptive lassoとかとも比較してみたい。
http://www.orsj.or.jp/archive2/or58-05/or58_5_261.pdf
ちなみにこの記事、線形回帰モデルからlassoへの流れ、L1正則化が変数選択として機能する理由、変数選択のオラクルプロパティ、その先の展開まで、わずか数ページでまとまっていてオススメ。
詳しい解説
http://www.is.titech.ac.jp/~s-taiji/tmp/sparse_tutorial_2014.pdf
そのうちきちんと書き直します。
foreach+%dopar%をネストしたい
並列化できる評価関数があり、その調節パラメータを並列にグリッドサーチしたい。単純に考えると親のソケットクラスタから子クラスタを生成すると良いのかな、と思うがもっといい方法ないものか。
評価関数は例えば、randomForestのようなアンサンブル学習でも良いし、遺伝的アルゴリズムのような最適化でも良いし、LOOみたいなcross-validationを考えても良い。Rで簡単に書ける。
require(foreach) require(doParallel) eval_parallel <- function(){ # 親ソケットクラスタ cl <- makeCluster(8) registerDoParallel(cl) on.exit(stopCluster(cl)) # 子クラスタの格納用 cl_i <-NULL # grid searchの条件組み合わせを列挙しておく conditions <- expand.grid(a=1:3, b=letters[2:7], c=c(1.1,1.7,1.9)) foreach(i=1:NROW(conditions), .packages = c("foreach", "doParallel"), .combine = rbind) %dopar%{ # 子クラスタ cl_i[[i]] <- makeCluster(11) registerDoParallel(cl_i[[i]]) on.exit(stopCluster(cl_i[[i]])) # ハイパーパラメータのセット parameters <- conditions[i,] # iteration_sizeは弱学習器の数だったり、CVのfolding数だったり、GAの個体数だったり foreach(j=1:iteration_size, .combine = c) %dopar% { # do something } } } eval_parallel()
この例では親が8いて子が11づつ作られるので、8*(1+11)=96並列。親子のバランスはボトルネックを見て配分すればよい。ときどき親の割り当てぶんの確保を忘れる。
.combineがrbindとcになってるのは、このあとでrowMeansすれば平均精度を見てベストなパラメータセットを探せるcross-validationのイメージ。
on.exit()使ってるのは下記の記事から。
R の foreach で並列処理するとき on.exit() が便利 #rstatsj - Qiita
もうちょっとうまいやり方がありそうだけど、思いつかないので今はこれで凌いでいる。