NN基礎からCNNと画像解析入門 - Chainer
大分前にちょこっと画像解析と畳込みニューラルネット（CNN）を勉強してて、周りに話そうかなと思ってスライド作ったけど、忙しくてそれどころじゃなくなったので、とりあえず保存。
Chainerのコードとかそのうちまとめて保存する予定...。
p.2 -11までがCNNに進むために理解が必要なニューラルネットワーク（NN）基礎。
p.12-23が画像データとCNNとChainerでの実装例。
参考書は、おなじみの以下書籍とマニュアル。
- 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
- 出版社/メーカー: 丸善出版
- 発売日: 2012/04/05
- メディア: 単行本（ソフトカバー）
- 購入: 6人 クリック: 33回
- この商品を含むブログ (18件) を見る
深層学習 Deep Learning (監修:人工知能学会)
- 作者: 麻生英樹,安田宗樹,前田新一,岡野原大輔,岡谷貴之,久保陽太郎,ボレガラダヌシカ,人工知能学会,神嶌敏弘
- 出版社/メーカー: 近代科学社
- 発売日: 2015/11/05
- メディア: 単行本
- この商品を含むブログ (1件) を見る
A Small Study on Segmentation by Clustering Analysis
業務で、Webログや契約者データを用いて、クラスター分析を使ったセグメンテーション（顧客分類）に関して考える必要があったので、それに関するメモ。
１）セグメンテーションに関して
以下の文献で、参考になることが多かったので、引用とメモ。
中村 2008 “マーケット・セグメンテーション - 購買履歴データを用いた販売機会の発見”
マーケット・セグメンテーション―購買履歴データを用いた販売機会の発見 (専修大学商学研究所叢書)
- 作者: 中村博
- 出版社/メーカー: 白桃書房
- 発売日: 2008/03
- メディア: 単行本
- クリック: 5回
- この商品を含むブログを見る
セグメンテーションアプローチ（顧客分類し、グループ別に異なったコミュニケーション）の際の注意点。使用する変数に関して。
スライドp.2
セグメンテーションアプローチにおいて、これらの評価指標があいまいな場合、実行不可能、もしくは結果に繋がらないと考えられる。
スライドp.3
実際に分析する場合、データ取得が容易だからと言って、性別や年齢といった変数を用いても効果は低い。
コンバージョンに直接繋がる影響の大きい変数は、一般的に抽象的かつ取得困難。
２）クラスター分析に関しての実験
セグメンテーションアプローチにおいて、これらの評価指標があいまいな場合、実行不可能、もしくは結果に繋がらないと考えられる。
スライドp.3
実際に分析する場合、データ取得が容易だからと言って、性別や年齢といった変数を用いても効果は低い。
コンバージョンに直接繋がる影響の大きい変数は、一般的に抽象的かつ取得困難。
２）クラスター分析に関しての実験
「何でこのように分けたのだ」、「どのような変数がビジネス上重要なのか」といった質問に、クラスター分析は非常に弱い。クラスター分析は教師なし学習。データの一つの見方を提供する分析であるり、使用する変数によって安定性を失ってしまう分析であり、使い方を気をつけないと、間違った結論を導いてしまう可能性があるので注意が必要だと思う。
スライドp.3-7では、あるビジネス目標を達成するためのセグメンテーションにクラスター分析を用いる場合、なんでも変数を投入してクラスター分析してしまっていいのか、に関して実験してみた。
スライドp.3-7では、あるビジネス目標を達成するためのセグメンテーションにクラスター分析を用いる場合、なんでも変数を投入してクラスター分析してしまっていいのか、に関して実験してみた。
結論としては、セグメンテーションアプローチをするといった何らかのビジネス目標がある場合、データの取得が容易な変数だけを集めてクラスター分析にかけると言ったアプローチは危険で、クラスター分析に投入する変数は吟味すべきであると考える。
スライドp.8では、簡単であるが、クラスター分析によるセグメンテーションアプローチの提案。
Stackingに関して（２）
前回の続き。
highschoolstudent.hatenablog.com
サンプルデータの作成。あまり細かいことは考えず適当なパラメータ設定で適当に作成。
library(MASS) # 乱数の作成（多変量正規分布） library(glmnet) # Ridge回帰 library(lme4) # 混合モデル library(lattice) # lme4による分析結果の可視化
n_obs <- 2000 # サンプル数 n_r <- 20 # rの水準数 n_x <- 30 # Xの列数 n_r_each <- n_obs/n_r # 各rの水準のサンプル数 ### r(random effect)の作成 ### s <- 5 # random effectの標準偏差 # rのmeanの生成（20水準） set.seed(123) mean_r <- floor(rnorm(n=n_r, mean=0, sd=15)) mean_r # mean_rのそれぞれの要素を20回繰り返す mean_r_vec <- rep(mean_r, times=1, each=n_r_each) # n_obs行のベクトル # 乱数(sd=sの正規分布)を加える set.seed(123) r_vec <- mean_r_vec + rnorm(n=n_obs, mean=0, sd=s) # 各rのID列の作成 r_id=rep(1:n_r, times=1, each=n_r_each) ### X（2000×30）の作成 ### mu <- rep(0, n_x) # 各変数の平均は0 sig <- diag(n_x) # 相関係数行列（30×30） sig[sig==0] <- 0.1 # 30の変数間に0.1の相関を与える sig_var <- sig*4 # 分散共分散行列。分散を4とする set.seed(123) x_mtrx <- mvrnorm(n_obs, mu=mu, Sigma=sig_var) # 1000*30の行列 X <- data.frame(x_mtrx) ### 目的変数yの作成 ### # すべての係数(β)は1とする。重み付け無しで30列を足し合わせる logit_p <- apply(X, MARGIN=1, sum) # randon effectを加える logit_r_p <- logit_p + r_vec # scaleしておく p <- 1 / (1+exp(-1*scale(logit_r_p))) # yの生成 set.seed(123) y <- rbinom(n=n_obs, size=1, prob=p) ### data.frameにまとめる ### data1 <- data.frame(r_id, mean_r_vec, X, p, y) # 列IDを付けておく data1 <- data.frame(OBS_ID=1:nrow(data1), data1) # データ型の変更 data1[, "r_id"] <- as.factor(data1[, "r_id"]) data1[, "y"] <- as.factor(data1[, "y"]) data1[, "OBS_ID"] <- as.character(data1[, "OBS_ID"])
> str(data1) 'data.frame': 2000 obs. of 35 variables: $ OBS_ID : chr "1" "2" "3" "4" ... $ r_id : Factor w/ 20 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ... $ mean_r_vec: num -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ... $ X1 : num -0.267 2.431 -0.354 -1.632 -2.469 ... $ X2 : num 0.0797 -0.6252 -3.987 -0.219 0.9425 ... $ X3 : num -2.678 1.573 0.518 -3.111 0.12 ... $ X4 : num 3.29 -0.87 -7.05 1.1 -1.74 ... $ X5 : num 3.235 -0.832 2.129 -1.876 2.369 ... $ X6 : num -1.846 -0.7481 0.0862 0.16 -2.193 ... $ X7 : num 2.383 -2.76 2.649 0.34 -0.901 ... $ X8 : num -1.725 -1.1587 -1.2911 0.325 0.0011 ... $ X9 : num 0.0625 1.6856 0.11 -0.6957 3.046 ... $ X10 : num -0.0845 2.1115 -2.3013 3.282 1.1797 ... $ X11 : num 0.6461 1.1761 -3.7507 1.4679 0.0739 ... $ X12 : num 1.024 2.317 -1.685 -0.367 -2.538 ... $ X13 : num -1.27 -3.483 -0.253 0.399 2.68 ... $ X14 : num -2.569 -0.559 -3.942 -3.517 -0.387 ... $ X15 : num 2.5588 -2.0188 0.0141 0.8221 0.2959 ... $ X16 : num 2.227 2.08 -0.565 0.555 0.261 ... $ X17 : num -3.413 1.206 1.578 0.226 -1.649 ... $ X18 : num 3.189 -0.286 -1.177 1.409 -1.41 ... $ X19 : num 4.053 0.355 0.173 -2.395 0.607 ... $ X20 : num -0.356 -0.632 -1.924 0.464 1.393 ... $ X21 : num 2.0457 -0.4395 0.7572 -2.1578 -0.0979 ... $ X22 : num 2.877 1.097 -1.977 1.85 0.151 ... $ X23 : num -0.447 -0.177 -1.276 -1.001 -0.904 ... $ X24 : num -0.506 0.924 -0.718 0.292 1.549 ... $ X25 : num 2.343 0.305 -3.97 -3.443 -0.222 ... $ X26 : num -1.731 -0.769 2.971 2.385 -0.441 ... $ X27 : num 1.6675 -0.0147 -2.1689 1.4975 -2.7424 ... $ X28 : num 0.0978 3.6201 -2.9536 1.036 1.4886 ... $ X29 : num -2.901 0.138 -3.443 -1.862 -1.06 ... $ X30 : num 0.1366 -0.6662 0.0819 3.1395 -0.1976 ... $ p : num 0.491 0.429 0.162 0.374 0.364 ... $ y : Factor w/ 2 levels "0","1": 1 2 1 2 2 1 1 1 1 1 ...
plot(data1[,"r_id"],data1[,"p"], xlab="r_id", ylab="p")
rによって、適当に散らばったデータ。
Train(Stacking用に分割), Testセットの作成。
Train:Test=1600:400。Trainは8分割し、各Stacking用データは200とする。
ID列（dataSet_ID）を作成し、管理。
Stacking用に8分割したTrainセットは1,2,3,4,5,6,7,8とIDを付け、Testセットを99とする。
# Train(Stacking), Test IDの作成 # Train -> 1,2,3,4,5,6,7,8 # Test -> 99 each_ID <- c(rep(1,10),rep(2,10),rep(3,10),rep(4,10),rep(5,10),rep(6,10),rep(7,10),rep(8,10),rep(99,20)) dataSet_ID <- rep(each_ID, n_r) data2 <- data.frame(data1, dataSet_ID) # Train, Test dataの分割 trainData <- data2[data2[,"dataSet_ID"]!=99, ] testData <- data2[data2[,"dataSet_ID"]==99, ]
では、Ridge回帰で、1stステージモデルを作成する。
glmnetパッケージを使用。
highschoolstudent.hatenablog.com
# Trainデータの説明変数 exp_vars <- as.matrix(trainData[,c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10", "X11","X12","X13","X14","X15","X16","X17","X18","X19","X20", "X21","X22","X23","X24","X25","X26","X27","X28","X29","X30")] ) # Trainデータの目的変数 #target_var <- trainData[,"y"] target_var <- as.vector(trainData[,"y"]) # Testデータの説明変数 exp_vars_test <- as.matrix(testData[,c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10", "X11","X12","X13","X14","X15","X16","X17","X18","X19","X20", "X21","X22","X23","X24","X25","X26","X27","X28","X29","X30")] )
まず、Cross Validationで最適なλを見つける。
cv.glmnetは、k+1回数glmnetを走らせる。
k回でλを決定し、残り1回でパラメータを推定する。
前回述べた通り、この分割はStackingの分割と同とする必要はない。
しかし、foldid引数にIDが指定できる、8分割で十分だろうという事で、Cross ValidatioとStackingの分割は同じとして処理する。
set.seed(123) fit1st_ridgeCV <- cv.glmnet( x=exp_vars, y=target_var, family="binomial", #nfold=8, foldid=trainData[,"dataSet_ID"], # CVのIDを指定 keep=TRUE, # CVのIDを保存 alpha=0 ) plot(fit1st_ridgeCV) # CV(k=8)による、最適なλ lambdaMin <- fit1st_ridgeCV$lambda.min
> lambdaMin [1] 0.1540261
見つけたλ（lambdaMin）によって、2ndモデル用の変数（Ridge回帰による予測値）を作成。
Trainデータの処理。
8個に分けた各Stackingセットを一つ抜いて、残りのデータで固定したλ（lambdaMin）でglmnetを実行し、パラメータを推定、抜いておいたセットにこの推定したモデルをあてはめて、予測値を計算する。
これをStackingに分けた数（8回）実行する。
# Stackingセット（= CVセット）のID cv_id <- as.integer(levels(as.factor(trainData[,"dataSet_ID"]))) # 1 2 3 4 5 6 7 8 # 8回のあてはめ結果オブジェクトの格納用リスト list_res <- list() # 予測値格納用のdf（col1:OBS_ID、col2:予測obsのID、col3:予測値） pred_y <- NULL for(id in cv_id){ #print(paste("cv_id:",id)) # パラメータ推定用データ train <- trainData[trainData[,"dataSet_ID"]!=id, ] #print("train - cv_id") #print(table(train[,"dataSet_ID"])) x <- as.matrix(train[,c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10", "X11","X12","X13","X14","X15","X16","X17","X18","X19","X20", "X21","X22","X23","X24","X25","X26","X27","X28","X29","X30")]) y <- as.vector(train[,"y"]) # あてはめ用データ target <- trainData[trainData[,"dataSet_ID"]==id, ] #print("target - cv_id") #print(table(target[,"dataSet_ID"])) x_target <- as.matrix(target[,c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10", "X11","X12","X13","X14","X15","X16","X17","X18","X19","X20", "X21","X22","X23","X24","X25","X26","X27","X28","X29","X30")]) fitRidge <- glmnet( # 指定したλ（lambdaMin）でリッジ回帰 x=x, y=y, family="binomial", lambda=lambdaMin, alpha=0 ) list_res[[id]] <- fitRidge # あてはめ結果オブジェクトの保存 # あてはめように抜いてあったデータに、パラメータ推定用データで推定したモデルをあてはめて、予測値を推定。あてはめデータも保存。 pred_y <- rbind(pred_y, data.frame(target[,"OBS_ID"], mean(target[,"dataSet_ID"]), predict(fitRidge, newx=x_target))) #print("---------------------") } names(pred_y) <- c("OBS_ID", "target_id", "pred_input") # リッジ回帰による予測で作成した列を含めたdata.frame（Trainデータ）の作成 trainData2 <- merge(x=trainData, y=pred_y, by.x="OBS_ID", by.y="OBS_ID", sort=FALSE, all=TRUE) # 列の絞り込みと順序の入れ替え trainData2 <- trainData2[, c("OBS_ID","r_id","mean_r_vec", "X1","X2","X3","X4","X5","X6","X7","X8","X9","X10", "X11","X12","X13","X14","X15","X16","X17","X18","X19","X20", "X21","X22","X23","X24","X25","X26","X27","X28","X29","X30", "p","y","dataSet_ID","pred_input")]
Trainデータの処理。
Trainデータ全体で、固定したλ（lambdaMin）でglmnetを実行し、パラメータを推定、それにTestデータをあてはめ予測値を作成。
fitRidge_test <- glmnet( # 指定したλ（lambdaMin）でリッジ回帰 x=exp_vars, y=target_var, family="binomial", lambda=lambdaMin, alpha=0 ) # Testデータでの予測値 ypred <- predict(fitRidge_test, newx=exp_vars_test) # リッジ回帰による予測で作成した列を含めたdata.frame（Testデータ）の作成 testData2 <- data.frame( testData, pred_input=as.vector(ypred) )
TrainとTestデータを一つにまとめる。
data3 <- rbind(trainData2, testData2) # pred_input（線形結合）をpに直す ypredp <- 1 / (1 + exp(-1*data3[,"pred_input"])) data3 <- data.frame(data3, pred_input_p=ypredp)
以上、1stステージモデル。
では、2ndステージモデルに入る。
# Train, Test dataの分割 Sec_trainData <- data3[data3[,"dataSet_ID"]!=99, ] Sec_testData <- data3[data3[,"dataSet_ID"]==99, ]
以前、混合モデルを使用した際は、glmmMRを使ったが、実際はlme4パッケージの方がメジャーなようなので、こちらを使う。
highschoolstudent.hatenablog.com
predict関数もサポートしており、変量効果の各水準に対する予測値(BLUP(Best Linior Unbiased Predictor))を考慮した予測が可能。
詳細に関しては、lme4: Mixed-effects modeling with R 参照。
ランダム切片モデルとする。
# Intercept model fit_mixed <- glmer( y ~ pred_input + (1|r_id), # 1は切片の意味 data=Sec_trainData, family=binomial, verbose=1 )
> summary(fit_mixed) Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [ glmerMod] Family: binomial ( logit ) Formula: y ~ pred_input + (1 | r_id) Data: Sec_trainData AIC BIC logLik deviance df.resid 1978.0 1994.2 -986.0 1972.0 1597 Scaled residuals: Min 1Q Median 3Q Max -2.8686 -0.8164 -0.3042 0.8195 3.2091 Random effects: Groups Name Variance Std.Dev. r_id (Intercept) 0.3299 0.5743 Number of obs: 1600, groups: r_id, 20 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.003755 0.139643 -0.027 0.979 pred_input 1.168448 0.094078 12.420 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) pred_input 0.017
ranef関数に、結果オブジェクトを引数で渡すと、変量効果の各水準に対する予測値(BLUP(Best Linior Unbiased Predictor))を出力できる。
condVar=TRUEとすると、各変量効果の水準の予測値の95%予測区間も出力する。
dotplot(ranef(fit_mixed, condVar=TRUE))
では、2ndステージモデルの予測値を出力する。
# Trainデータの予測値 predTrain_fit_mixed <- 1 / ( 1 + exp(-1*predict(fit_mixed, newdata=Sec_trainData)) ) # Testデータの予測値 predTest_fit_mixed <- 1 / ( 1 + exp(-1*predict(fit_mixed, newdata=Sec_testData)) ) # data.frameにまとめる Sec_trainData2 <- data.frame(Sec_trainData, pred_y=predTrain_fit_mixed ) Sec_testData2 <- data.frame(Sec_testData, pred_y=predTest_fit_mixed) data4 <- rbind(Sec_trainData2, Sec_testData2)
上で使用したprediction関数であるが、引数re.form=NULL（デフォルト）としておくと、変量効果の予測値（BLUP）を含めた予測となる。
allow.new.levelsという引数があり、allow.new.levels=TRUE（デフォルトはFALSE）とすると、新しいデータに対する予測の際、モデル推定のデータになかった変量効果の水準が含まれた場合、変量効果の予測を考慮しない値が返される。
デフォルト（FALSE）だと、モデル推定のデータになかった変量効果の水準が含まれた場合、Errorが返るとのこと。
> head(data4) OBS_ID r_id mean_r_vec X1 X2 X3 X4 X5 1 1 1 -9 -0.2670289 0.07965359 -2.6781329 3.2920218 3.2347046 2 2 1 -9 2.4309568 -0.62516429 1.5726213 -0.8702877 -0.8316664 3 3 1 -9 -0.3538252 -3.98699188 0.5181232 -7.0503276 2.1287953 4 4 1 -9 -1.6324691 -0.21901088 -3.1107837 1.1008922 -1.8763152 5 5 1 -9 -2.4687447 0.94246812 0.1196235 -1.7423328 2.3693425 6 6 1 -9 1.2053775 -2.38738177 0.7406409 -2.9720772 -1.2437253 ### 省略 ### X27 X28 X29 X30 p y dataSet_ID pred_input 1 1.6675496 0.09780506 -2.9014770 0.13657163 0.4909926 0 1 0.49382763 2 -0.0146797 3.62010783 0.1383445 -0.66618073 0.4288412 1 1 0.05697817 3 -2.1689243 -2.95355951 -3.4428897 0.08190659 0.1619845 0 1 -0.89058592 4 1.4974864 1.03595553 -1.8618701 3.13950267 0.3741005 1 1 0.02697507 5 -2.7424436 1.48861771 -1.0603866 -0.19760381 0.3637205 1 1 -0.26863659 6 0.1702958 1.57345489 -2.8409584 -1.65512597 0.1465225 0 1 -1.05824906 pred_input_p pred_y 1 0.6210077 0.5308346 2 0.5142407 0.4044539 3 0.2909889 0.1833040 4 0.5067434 0.3960387 5 0.4332418 0.3170407 6 0.2576442 0.1557724
pred_yが最終予測値である。
なお、パフォーマンスに関してはAUCで確認している。
結果だけだが、Trainデータで0.7404149、Testデータで0.7289979となる。
Stackingを行わなかった場合のモデルとの比較等に関しては省略。
以上。
Stackingに関して （１）
予測モデルのコンペで使われているStackingという手法に関して。
モデルをEnsembleしまくって、マルチステージモデルにして予測精度を高めるときに使われる手法とのこと。
実務家な私も、業務でマルチステージモデル(2ステージ)を組む必要があったので、そのためのメモ。
kaggleなどで発明され、まだ教科書にはなっていない方法論のようで、以下を参照。
具体的には以下のような課題（データ）があった。
目的変数y(0,1のBinomial)、説明変数rとXでモデルを作りたい。
rは個々のセールスマンを示すようなカテゴリカルデータで、変量効果として扱いたい。
Xは営業日誌のようなテキストデータから単語抽出して、TF-IDF変換した横長なデータ。
混合モデル（Logistic GLMM）を作成したいが、横長なXをそのままロジスティックモデルの変数としてそのまま加えるのは不適切、といった理由で、まずXをRidge回帰によって処理をしてから、その予測値を2段階モデルの説明変数として、ｒが変量効果となる混合モデルを作成することで対処する。
Ridge回帰に関しては、以下を参照。
highschoolstudent.hatenablog.com
混合モデルに関しては、以下を参照。
highschoolstudent.hatenablog.com
マルチステージモデルを組む際、ちゃんとしたStackingの手順を踏まないと予測値が上がらない、と、こういったのに詳しい同僚に言われたので、上の参考文献（KAGGLE ENSEMBLING GUIDE）を参考にしつつ、具体的手順を記してみる。
コアとなる考え方だが、最終モデル（例では2nd Model）以前では、ペアとなっている目的変数をモデルパラメータ推定の際使ってはいけない、と言うことだと思う。
注：文献を読んでの私なりの理解なので、間違ってたり、他にもやり方がある可能性有り
データをTrain（学習）データとTest（テスト）データに分ける。
Trainデータをさらに、パラメータ推定のためにいくつかに分ける。図では説明を簡単にするために2つに分けたが、もう少し細かく分割しても良いと思う。
なお、この分割は、モデルの複雑さを決定するためのk分割交差検証のためではない。
まず、Trainデータ全体でk分割交差検証法（kは適当な値）などを用い、モデルの複雑さを決定する。今回の例では、Ridge回帰のλ(Complexity Parameter)を決定する。
決定したモデルの複雑さ（λ）を用い、X1とy1でモデル（Mod1）を作成（パラメータを推定）。そのモデルをX2にあてはめ、予測値x2を作成。
決定したモデルの複雑さ（λ）を用い、X2とy2でモデル（Mod2）を作成（パラメータを推定）。そのモデルをX1にあてはめ、予測値x1を作成。
決定したモデルの複雑さ（λ）を用い、X1,X2とy1,y2でモデル（Mod'）を作成（パラメータを推定）。そのモデルをTestデータのX'にあてはめ、予測値x'を作成。
予測により作成する各変数（x1、x2、x'）は、ペアとなる目的変数（x1はy1、x2はy2、x'はy'）を使用してない、と言うことになる。
Trainデータにおいて、作成した予測値（x1,x2）と、取っておいたその他の変数（変量効果とするr1,r2）を説明変数、Trainデータの目的変数（y1,y2）を使って、最終モデル（Mod''）を作成する。
これをTestデータの説明変数（r',x'）にあてはめ、最終的な予測値（y*）を計算する。
では次回、シミュレーションによって作成したサンプルデータでやってみる。
単なるSQLメモ
このようなSampleデータがあったとする。
Group別に、時間（Min）が最大となる行を取り出したい。
クエリー実行後。
SELECT T1.*
FROM Sample T1
INNER JOIN(
SELECT Group, MAX(Min) AS MAX_Min
FROM Sample
GROUP BY Group
) T2 ON T1.Group = T2.Group
AND T1.Min = T2.MAX_Min
この場合、時間（Min）に重複がないことを仮定している。
またこういった、時間が分かれている（Min、Sec）Sampleデータの場合。
クエリー実行後。
SELECT T1.*
FROM Sample T1
INNER JOIN(
SELECT Group, MAX( Min||Sec ) AS MAX_MinSec
FROM Sample
GROUP BY Group
) T2 ON T1.Group = T2.Group
AND T1.Min || T1.Sec = T2.MAX_MinSec
以上、単なるSQLメモ。
モデル評価基準 追加 - リフトチャートに関して２ - R{ROCR}
highschoolstudent.hatenablog.com
highschoolstudent.hatenablog.com
上2つの記事に関連した追加。
そもそもリフトチャートと呼ばれるものは定義があいまいなようで、人によって使っているものが細かく違ったりする。
基本的に似たようなものであるのだが、”モデル評価基準 追加 - リフトチャートに関して２ - R{ROCR}”で作成したのとまたちょっと違ったリフトチャートを描いてみる。
業務上、このチャートを標準とすることとなった。確かに、シンプルで理解もしやすいのでお勧め。
以下の手順で作成できる。
１．予測値を小さいものから順に並べる。
２．いくつかのグループに分ける。（今回の19行のデータ数ではグループ数を大きくできないが、10とか20グループに分けるやり方が良く使われるらしい。）各グループ内でのデータ数（行数）は等しい。（必ずしも割り切れる訳ではないので、適当な調整を加える。）
３．各グループでの観測値の平均値を計算する。（観測値は0/1なので、各グループでの1の割合が計算される。）
４．横軸をグループ番号、縦軸を３．で計算した観測値の平均値をプロットする。すべての観測値での平均値もプロットしておくとグラフが理解しやすい。
データはこれまでと同じものを使用。
d_roc1 <- read.delim("clipboard", header=TRUE)
d_roc1
pred1 pred2 observed 1 0.95 0.90 1 2 0.90 0.80 1 3 0.80 0.75 1 4 0.60 0.70 1 5 0.50 0.60 1 6 0.85 0.95 0 7 0.75 0.85 0 8 0.70 0.65 0 9 0.65 0.55 0 10 0.55 0.50 0 11 0.45 0.45 0 12 0.40 0.40 0 13 0.35 0.35 0 14 0.30 0.30 0 15 0.25 0.25 0 16 0.20 0.20 0 17 0.15 0.15 0 18 0.10 0.10 0 19 0.05 0.05 0
関数を定義する。
fn_plot_Lift20bins <- function(pred, obs, aveline=TRUE, nbins=20, color="black"){
df0 <- data.frame(pred=pred, obs=obs)
df1 <- df0[order(df0$pred),]
n <- nrow(df1)
if(nbins*2>n) {stop(">>> Set smaller nbin. nbins*2 <= nrow")}
meanobs <- mean(df1$obs)
nrep <- floor(n/nbins) # 各binに入る行数
ndup <- n-nbins*nrep # nbinでnが割り切れない場合の対処
bin <- c()
for(i in 1:nbins){
if(i<=ndup){ bin <- c(bin, rep(i, nrep+1)) } # nbinでnが割り切れない場合,最初のbinに1つづつ加えていく
else{ bin <- c(bin, rep(i, nrep)) }
}
df2 <- data.frame(df1, bin=bin)
#print( df2 )
ratio <- tapply(df2$obs, df2$bin, mean)
df3 <- data.frame(bin=1:nbins, ratio)
print( df3 )
plot( df3$bin, df3$ratio, type="b" , xlab="Bin", ylab="Ratio", col=color)
abline(h=meanobs, col="yellow")
abline( v=1:nbins, lty=3 )
}
pred: モデル予測値のベクトル
obs: 観測値のベクトル
aveline: 全観測値の平均線をプロットするかどうか
nbins: グループ数（データ数に対し大きいとエラーになる）
color: 線の色
pred1に対するリフトチャート。
fn_plot_Lift20bins(pred=d_roc1$pred1, obs=d_roc1$observed, aveline=TRUE, nbins=5)
bin ratio 1 1 0.0000000 2 2 0.0000000 3 3 0.5000000 4 4 0.2500000 5 5 0.6666667
5グループに分けたので、各グループは3か4のデータ数となっている。（4*4+3*1=19）
各グループでの観測値の平均のデータフレームを表示し、それのプロットも出力される。
黄色の横線は、全観測値での平均値。
低い方から並べていっているので、基本的に右上がりのグラフとなる。
良いモデルでは、このリフトチャートは単調増加（monotone increasing）でカーブが急なものとなる。
今回の5グループのプロットだと、途中で平均値が落ちてしまっているが、3グループくらいにして、グループ内でのデータ数を確保すると単調増加のカーブになる。
定義した関数は重ね描きには対応してないので、今後修正の必要あり。
以上。
Gradient Boosting Treeを使ってみる - R{gbm} (Part.2)
前回に続き、いくつかの変数変換を試してみる。
highschoolstudent.hatenablog.com
Partial Dependence Plotで気づいたことを元に、以下の説明変数に対して変換を実行。
"age" - 数値変数だとノイズに敏感に見えるため、カテゴリー化してみる。
"fnlwgt" - 外れ値にキャップをかけてみる。
"capital.gain" - 外れ値にキャップをかけてみる。
"capital.loss" - 元の変数の意味が分からなくて解釈に困るが、カテゴリー化して様子を見てみる。
"hours.per.week" - 外れ値にキャップをかけてみる。
"native.country" - カテゴリカル変数の多水準問題に2つの方法（両方とも数値化する方法）で取り組んでみる。
キャップと言うのは、連続変数に上限下限を設けること。
例えば、"身長"という変数があって180cmで上限のキャップをかけた場合、それ以上高い観測値はすべて180cmに変換する。
Tree系のモデルは説明変数の外れ値にロバストであり、キャップをかけることによってパフォーマンスが大幅に改善したりするとは考えられない気がするが、それでも念のために外れ値の処理はしておいたほうが良いとのこと（同僚より）なので、実行しておく。
カテゴリー化の理由だが、Plotを見た感じで学習データにオーバーフィッティングしている感じがしたので、カテゴリー化することによりスムーズにしようというのが一つで、もう一つが解釈しやすくなる（Plotで見やすくなる）のではないかという期待から。
いずれにせよそこまで科学的な根拠があると言うわけではないが、良くやられているので、とりあえずやってみる、といった要素が強い。
キャップ関数を定義。変数（データフレームの列）を引数とし、ベクトルを返す。
# Capをかける関数、カラムと最小/最大値を指定
fn_cap <- function(num_col, low=min(num_col), high=max(num_col)){
num_col[num_col<=low] <- low
num_col[num_col>=high] <- high
return(num_col)
}
カテゴリー化の関数を定義。変数（データフレームの列）を引数とし、ベクトルを返す。
分け方だが、上記理由を考慮した完全な主観。
# age(numeric)をカテゴリー化する関数
fn_age_category <- function(age_col){
age_col[age_col<20] <- "-19"
age_col[age_col>=20 & age_col<25] <- "20-24"
age_col[age_col>=25 & age_col<30] <- "25-29"
age_col[age_col>=30 & age_col<35] <- "30-34"
age_col[age_col>=35 & age_col<40] <- "35-39"
age_col[age_col>=40 & age_col<45] <- "40-44"
age_col[age_col>=45 & age_col<50] <- "45-49"
age_col[age_col>=50 & age_col<55] <- "50-54"
age_col[age_col>=55 & age_col<60] <- "55-59"
age_col[age_col>=60 & age_col<65] <- "60-65"
age_col[age_col>=65] <- "65-"
return(age_col)
}
# capital.lossをカテゴリー化する関数
fn_capLoss_category <- function(closs_col){
closs_col[closs_col==0] <- "0"
closs_col[closs_col>=0 & closs_col<1000] <- "<1000"
closs_col[closs_col>=1000 & closs_col<1250] <- "<1250"
closs_col[closs_col>=1250 & closs_col<1500] <- "<1500"
closs_col[closs_col>=1500 & closs_col<1750] <- "<1750"
closs_col[closs_col>=1750 & closs_col<2000] <- "<2000"
closs_col[closs_col>=2000] <- ">=2000"
return(closs_col)
}
データ加工は、adult2_2から出発。
後で実行する多水準カテゴリカル変数の変換のときに必要であるからだが、学習データ（学習＋検証）とテストデータを示すフラグ変数を作成しておく。
前回同様、24420までが学習＋検証（0とする）、24421から最後の32561までがテスト（1とする）。
test_col <- c(rep(0,24420),rep(1,(nrow(adult2_2)-24420)))
adult2_3 <- data.frame(adult2_2, Test=test_col)
head(adult2_3)
age workclass fnlwgt education education.num marital.status
1 34 Private 191856 Some-college 10 Married
2 40 Private 205047 HS-grad 9 Married
3 35 Private 230054 7th-8th 4 Married
4 27 Private 279608 5th-6th 3 Married
5 34 Private 226872 Bachelors 13 Divorced
6 61 State-gov 151459 10th 6 Never-married
occupation relationship race sex capital.gain capital.loss
1 Exec-managerial Wife White Female 7298 0
2 Craft-repair Not-in-family White Male 0 0
3 Machine-op-inspct Husband White Male 0 0
4 Farming-fishing Husband White Male 0 0
5 Sales Not-in-family White Female 0 0
6 Other-service Not-in-family Black Female 0 0
hours.per.week native.country income_class income Test
1 40 United-States >50K 1 0
2 40 United-States >50K 1 0
3 40 Mexico <=50K 0 0
4 40 Mexico <=50K 0 0
5 40 United-States <=50K 0 0
6 38 United-States <=50K 0 0
最終列に0/1データの"Test"カラムが作成される。
ageのカテゴリー化。
age2 <- fn_age_category(adult2_3[,"age"])
fnlwgtに500000でキャップ。
fnlwgt2 <- fn_cap(adult2_3[,"fnlwgt"], high=500000)
capital.gainに10000でキャップ。
capital.gain2 <- fn_cap(adult2_3[,"capital.gain"], high=10000)
capital.lossのカテゴリー化。
capital.loss2 <- fn_capLoss_category(adult2_3[,"capital.loss"])
hours.per.weekに80でキャップ。
hours.par.week2 <- fn_cap(adult2_3$hours.per.week, high=80)
データフレームにまとめる。
adult2_4 <- data.frame(adult2_3, age2, fnlwgt2, capital.gain2, capital.loss2, hours.par.week2)
では、カテゴリカル変数の多水準問題（High Cardinality Problem)への対処を考えてみる。
今回は逆に、数値化することによって、オーバーフィッティングへの対処を試みる。
一つ目の方法が、単純にランキング化する方法。Simple Rankingと呼ぶことにする。
学習データにおいて（これがTest列を作成しておいた理由）水準別に目的変数を集計し、値が高いものから順にランキングを付ける。
"native.country"列は42水準。この42水準別に目的変数である"income"が１となる割合を計算する。高いものから順に並べ、1から42のランクを付ける。
関数群を定義しておく。
############################
###### Simple Ranking ######
# 水準別ランキングをdfとして返す
fn_ranking_table <- function(d_f, y_name, x_name, x_new_name){
# 水準別平均値のテーブル
tbl <- tapply(d_f[,y_name], d_f[,x_name], mean)
# ランキング
rnk <- rank(-tbl, na.last="keep", ties.method="average")
# dfとして返す
tbl_df <- data.frame(names(rnk), rnk)
names(tbl_df) <- c(x_name, x_new_name)
return( tbl_df )
}
# (df, "目的変数(numeric)名", "説明変数(カテゴリー)名", "Rankingによる新しい変数名", "テスト変数名")
# テスト変数 - 0:train&validation 1:test
fn_SimpleRanking2 <- function(d_f, y_name, x_name, x_new_name, test_name=NULL){
#テストデータ(水準1)を除く
if(is.null(test_name)){ # テストデータを指定しない場合
d_f_train <- d_f
}
else if(sum(d_f[,test_name]==1)>0) { # テストデータ用の水準(1)を指定した場合
d_f_train <- d_f[d_f[,test_name]!=1,]
}
else { stop("--- Test data set??? ---") } # 不明な変数を指定した場合
#Rankingテーブルの作成
tbl_df <- fn_ranking_table(d_f_train, y_name, x_name, x_new_name)
return( merge(d_f, tbl_df, by=x_name, all.x=TRUE) )
}
############################
############################
fn_SimpleRanking2(データフレーム, "目的変数(numeric)名", "説明変数(カテゴリー)名", "Rankingによる新しい変数名", "テスト変数名(0:学習, 1:テスト)")として実行する。
temp_df <- fn_SimpleRanking2(adult2_4, y_name="income", x_name="native.country", x_new_name="native.country2_SR", test_name="Test")
では、もう一つの方法を実行。
学習データにおいて、説明変数のある水準を除いたときの目的変数の平均値をその水準の値とする。
例えば、"native.country"変数でJapan水準を除いたときの"Income"変数の平均値をJapanの値とする。
さらに、オーバーフィッティングを押さえるため、適当なノイズを加える。
Mean of Response Variable as Explanetory Variable by Leave One Out、略してLOOと呼ぶことにする。
この方法は、同僚に教えてもらった。
###########################
########### LOO ###########
# 説明変数の水準を指定すると、その水準を除いた目的変数の平均値を返す関数
# (df, "目的変数(numeric)名", "説明変数(カテゴリー)名", 説明変数の水準)
fn_mean_oneLevel <- function(d_f, y_name, x_name, x_level){
return( mean(d_f[d_f[,x_name]!=x_level,y_name]) )
}
# fn_mean_oneLevel関数を、指定した説明変数の水準分イテレーション実施して、dfを返す関数
# (df, "目的変数(numeric)名", "説明変数(カテゴリー)名", "LOOによる新しい変数名")
fn_mean_oneLevel_table <- function(d_f, y_name, x_name, x_new_name){
Num <- c()
for(lvl in levels(d_f[,x_name])){
Num <- c(Num, fn_mean_oneLevel(d_f, y_name, x_name, lvl))
}
loo_df <- data.frame(levels(d_f[,x_name]), Num)
names(loo_df) <- c(x_name, x_new_name)
return( loo_df )
}
# 元のdfにLOOを付け、返す関数
# (df, "目的変数(numeric)名", "説明変数(カテゴリー)名", "LOOによる新しい変数名", ノイズの大きさ, "テスト変数名")
# テスト変数 - 0:train&validation 1:test
# 学習データのみにある水準が含まれ、テストデータにその水準が含まれない場合の処理はしていない
fn_LOO2 <- function(d_f, y_name, x_name, x_new_name, noise_ratio=0.1, test_name=NULL){
#テストデータ(水準1)を除く
if(is.null(test_name)){ # テストデータを指定しない場合
d_f_train <- d_f
}
else if(sum(d_f[,test_name]==1)>0) { # テストデータ用の水準(1)を指定した場合
d_f_train <- d_f[d_f[,test_name]!=1,]
}
else { stop("--- Test data set??? ---") } # 不明な変数を指定した場合
#LOOテーブルの作成
loo_table <- fn_mean_oneLevel_table(d_f_train, y_name, x_name, x_new_name)
# merge
d_f_merge <- merge(d_f, loo_table, by=x_name, all.x=TRUE, sort=FALSE)
# ノイズのSD
sdiv <- sd(d_f_merge[,x_new_name])
# sdivをベースに正規分布からノイズを作成、noise_ratioで大きさを調整
noise <- rnorm(nrow(d_f), mean=0, sd=sdiv*noise_ratio)
# ノイズを加える
d_f_merge[,x_new_name] <- d_f_merge[,x_new_name] + noise
return( d_f_merge )
}
###########################
###########################
fn_LOO2(データフレーム, "目的変数(numeric)名", "説明変数(カテゴリー)名", "LOOによる新しい変数名", ノイズの大きさ, "テスト変数名(0:学習, 1:テスト)") として実行。
'ノイズの大きさ'とは、0から1の間の値を取り、0だと全くノイズを与えないという意味。
LOOによって各水準の値を計算した後、それらのばらつき（SD）を基準として、そのSDの何割かをノイズとして加えておく。
rnorm(n, mean=0, SD*'ノイズの大きさ')といったノイズを発生させ、LOOで計算した変数に加えておく。
adult2_5 <- fn_LOO2(temp_df, y_name="income", x_name="native.country", x_new_name="native.country2_LOO", noise_ratio=0.1, test_name="Test")
ノイズは0.1にしてみた。
変数が増えてややこしくなってきたので、順序を入れ替えて整理しておく。
adult3 <- adult2_5[,c("Test", "age", "age2", "workclass", "fnlwgt", "fnlwgt2", "education", "education.num", "marital.status", "occupation", "relationship", "race", "sex", "capital.gain", "capital.gain2", "capital.loss", "capital.loss2", "hours.per.week", "hours.par.week2", "native.country", "native.country2_SR",
"native.country2_LOO", "income_class","income")]
str(adult3)
'data.frame': 32561 obs. of 24 variables: $ Test : int 0 1 0 0 0 1 0 0 0 0 ... $ age : int 26 41 35 20 27 35 38 37 37 36 ... $ age2 : Factor w/ 11 levels "-19","20-24",..: 3 6 5 2 3 5 5 5 5 5 ... $ workclass : Factor w/ 7 levels "?","Federal-gov",..: 4 4 2 1 4 4 4 3 5 4 ... $ fnlwgt : int 56929 207685 35309 200061 269354 187415 220237 244803 188774 31438 ... $ fnlwgt2 : num 56929 207685 35309 200061 269354 ... $ education : Factor w/ 16 levels "10th","11th",..: 10 12 10 16 10 5 9 12 2 12 ... $ education.num : int 13 9 13 10 13 3 11 9 7 9 ... $ marital.status : Factor w/ 5 levels "Divorced","Married",..: 3 2 3 3 2 2 2 2 2 1 ... $ occupation : Factor w/ 13 levels "?","Adm-clerical",..: 4 13 12 1 11 7 3 3 11 13 ... $ relationship : Factor w/ 6 levels "Husband","Not-in-family",..: 2 1 2 4 1 1 1 1 2 5 ... $ race : Factor w/ 4 levels "Asian-Pac-Islander",..: 2 4 1 4 4 1 4 4 4 4 ... $ sex : Factor w/ 2 levels "Female","Male": 2 2 2 1 2 2 2 2 2 2 ... $ capital.gain : int 0 0 0 0 0 0 0 0 0 0 ... $ capital.gain2 : int 0 0 0 0 0 0 0 0 0 0 ... $ capital.loss : int 0 0 0 0 0 0 0 0 0 0 ... $ capital.loss2 : Factor w/ 6 levels "<1000","<1250",..: 1 1 1 1 1 1 1 1 1 1 ... $ hours.per.week : int 50 45 28 40 25 50 40 40 60 43 ... $ hours.par.week2 : int 50 45 28 40 25 50 40 40 60 43 ... $ native.country : Factor w/ 42 levels "?","Cambodia",..: 1 1 1 1 1 1 1 1 1 1 ... $ native.country2_SR : num 19 19 19 19 19 19 19 19 19 19 ... $ native.country2_LOO: num 0.241 0.24 0.241 0.242 0.241 ... $ income_class : Factor w/ 2 levels "<=50K",">50K": 1 1 1 1 1 1 2 1 1 1 ... $ income : int 0 0 0 0 0 0 1 0 0 0 ...
では、再びモデルを作成してみる。
学習(+検証)、テストデータの分割。
train3 <- adult3[adult3[,"Test"]==0,]
test3 <- adult3[adult3[,"Test"]==1,]
"native.country2_SR"と"native.country2_LOO"で、別々にモデルを作ってみる。
（元々の変数重要度が、最初に作ったモデルでも4.893%しかないので、大きく違ってくるとは考えにくい）
まず、Simple Ranking,"native.country2_SR"。
キャップをかけたりカテゴリー化した変数ももちろん用いる。
x_col_names <- c("age2", "workclass", "fnlwgt2", "education", "education.num", "marital.status", "occupation", "relationship", "race", "sex", "capital.gain2", "capital.loss2", "hours.par.week2", "native.country2_SR")
x_cols <- train3[, x_col_names]
y_col <- train3[,"income"]
train3データセットのうち、最初の16264行を学習、残りを検証。これは前回と同じ。
for_train <- 16264
前回のPatial Dependence Plotの結果を見て、単調増加(monotone increasing)に見えた変数には、その制約を課してみる。
"education.num"と"capital.gain"をmonotone increasingとする。ベクトルに1を入れる。
単調減少(monotone decreasing)の場合は-1を入れる。
"native.country2_SR"は、学習(+検証)データにおいて目的変数の値が高いものから順に並べれれているので、monotone decreasingの制約をかける。
単調増加/減少にすることにより、オーバーフィッティングの回避が行われるのだと思われる。
monotone <- c(0,0,0,0,1,0,0,0,0,0,1,0,0,-1)
0はこの制約をかけないという意味。
引数var.monotoneに、このベクトルを渡す。
mod2_1 <- gbm.fit(
x = x_cols,
y = y_col,
distribution = "bernoulli",
var.monotone = monotone,
bag.fraction = 0.5,
n.trees = 5000,
interaction.depth = 7,
n.minobsinnode = 100,
shrinkage = 0.005,
nTrain = for_train
)
変数重要度。
summary(mod2_1)
var rel.inf relationship relationship 30.1044827 capital.gain2 capital.gain2 15.2426028 occupation occupation 14.1707533 education education 11.5529362 age2 age2 8.6704758 education.num education.num 4.4322248 capital.loss2 capital.loss2 3.9274653 hours.par.week2 hours.par.week2 3.6641404 fnlwgt2 fnlwgt2 3.1119212 workclass workclass 2.3148916 marital.status marital.status 1.4242388 native.country2_SR native.country2_SR 0.8392975 sex sex 0.3368926 race race 0.2076771
native.countryをランク化したことによって、重要度が落ちた。他はあまり変化がない。
最適なTreeの数。
gbm.perf(mod2_1, plot.it=TRUE)
3132
では、予測を行う。
学習データ。
x_cols_train <- x_cols[1:for_train,]
y_col_train <- y_col[1:for_train]
検証データ。
x_cols_varidation <- x_cols[(for_train+1):nrow(x_cols),]
y_col_varidation <- y_col[(for_train+1):nrow(x_cols)]
学習データに対する予測値。
pred_mod2_1_train <- predict(mod2_1, newdata=x_cols_train, n.trees=gbm.perf(mod2_1, plot.it=FALSE), type="response")
検証データに対する予測値。
pred_mod2_1_varidation <- predict(mod2_1, newdata=x_cols_varidation, n.trees=gbm.perf(mod2_1, plot.it=FALSE), type="response")
テストデータに対する予測値。
test_x <- test3[,x_col_names]
test_y <- test3[,"income"]
pred_mod2_1_test <- predict(mod2_1, newdata=test_x, n.trees=gbm.perf(mod2_1, plot.it=FALSE), type="response")
学習データでのAUC。
fn_AUC(pred=pred_mod2_1_train, obs=y_col_train)
0.9332073
検証データでのAUC。
fn_AUC(pred=pred_mod2_1_varidation, obs=y_col_varidation)
0.9198591
テストデータでのAUC。
fn_AUC(pred=pred_mod2_1_test, obs=test_y)
0.9192887
fn_AUCに関しては、前回参照。
個人的な経験だが、今回示したようなオーバーフィッティングの回避を試みると、学習とテストデータのパフォーマンス(AUCなど)の幅が狭くなり、検証とテストデータのパフォーマンスもかなり近くなるといった傾向がある。
前回の結果と見比べてみてほしい。
一応、LOOも試しておく。
x_col_names <- c("age2", "workclass", "fnlwgt2", "education", "education.num", "marital.status", "occupation", "relationship", "race", "sex", "capital.gain2", "capital.loss2", "hours.par.week2", "native.country2_LOO")
x_cols <- train3[, x_col_names]
y_col <- train3[,"income"]
for_train <- 16264
monotone <- c(0,0,0,0,1,0,0,0,0,0,1,0,0,1)
"native.country2_LOO"はmonotone increasingとしておく。
同様にgbmを実行し、結果をmod2_2に格納。
変数重要度。
summary(mod2_2)
var rel.inf relationship relationship 30.0720371 capital.gain2 capital.gain2 15.2900749 occupation occupation 14.2201974 education education 11.1203147 age2 age2 8.7793286 education.num education.num 5.1393013 capital.loss2 capital.loss2 3.9263051 hours.par.week2 hours.par.week2 3.6922352 fnlwgt2 fnlwgt2 3.0737037 workclass workclass 2.3722413 marital.status marital.status 1.3842674 sex sex 0.3466526 native.country2_LOO native.country2_LOO 0.3382367 race race 0.2451039
最適なTreeの数。
gbm.perf(mod2_2, plot.it=TRUE)
3328
予測を行い、AUCを計算する。
学習データでのAUC。
0.9331306
検証データでのAUC。
0.9196043
テストデータでのAUC。
0.9192822
以上、今回はSimple RankingやLOO関数のメモを個人的理由とする。
