差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
mselab:2012:stat:week3:r2 [2012/12/18 04:46] – watalumselab:2012:stat:week3:r2 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1
行 1: 行 1:
-縛り。+===== 統計工学実験 第3週 ===== 
 +==== 連絡 ====
  
-  * 学習機械にはlm()かrpart()しか使わない。 +  * [[mselab:2012:stat:week3:r2#k-means法による層別_itermax追加_nstart削除|実験室で発生したkmeans関数についての不具合へ対応]]を追記。 
-  * kmeans()を層別に使う。+==== 最終課題 ====
  
 +  * 学習用データ(tic.learn)に基づいて、V86の契約に関するモデルを構築し、訪問対象に加える条件(訪問ルール)を提案せよ。
 +    * 学習機械にはlm()かrpart()しか使わない。
 +    * 層別のモデル構築も、候補に加える。
 +    * 層別の生成には、データの集計、決定木の他に、k-means法によるクラスタリングも候補に加える。
 +  * 提案した訪問ルールの契約達成率などを、検証用データ(tic.eval)に適用して検討せよ。
 +  * 学習用データでモデルを学習し、検証用データで精度を検討する、という縛りの中で、最適な訪問ルールを定めよ。これがこの実験の最終成果物である。
 +
 +==== 準備 ====
  
 準備。 準備。
 +
 <code> <code>
-install.packages(c("mvtnorm"), dependencies=TRUE)+Sys.setenv("http_proxy"="http://130.153.8.66:8080/") 
 +install.packages(c("mvtnorm", "mvpart"), dependencies=TRUE)
 library(mvtnorm) library(mvtnorm)
 +library(mvpart)
 library(MASS) library(MASS)
 </code> </code>
 +
 +データがない人は、次のコードも実行する必要がある。
 +
 +<code>
 +tic.learn <- read.table("http://kdd.ics.uci.edu/databases/tic/ticdata2000.txt")
 +tic.eval <- read.table("http://kdd.ics.uci.edu/databases/tic/ticeval2000.txt")
 +tic.test <- read.table("http://kdd.ics.uci.edu/databases/tic/tictgts2000.txt")
 +tic.eval <- cbind(tic.eval, tic.test)
 +colnames(tic.eval)[86] <- "V86"
 +rm(tic.test)
 +</code>
 +
 +さらに先週のレポート課題と同様に、V1とV5のグループ化も済ませるとよいかもしれない。
 +あるいは今週の課題で改めて、グループ化を見直してもよい。
 +
 +==== クラスタリング ====
 +
 +データの層別を行う方法の総称。生成した層のことを、クラスタ(群)という。大きく階層クラスタリングと非階層クラスタリングに分かれる。
 +
 +階層クラスタリングは、レコード間の距離を全ての組み合わせについて算出して作成した行列を距離行列から出発する。
 +まずは最も近いレコード同士をグループ(群)にまとめる。その結果、2レコードは1つのクラスタに属するが、残りのレコードは点のままである。
 +そのクラスタと残りのレコードの間の距離行列を算出し、最も近いレコード同士、最も近いレコードとグループ、もしくは最も近いグループ同士を再びグループに\\
 +まとめる。
 +これを繰り返すことで、クラスタ(群)分けを得るのが、階層クラスタリングである。
 +階層クラスタリングには、レコード同士の距離だけでなく、レコードとグループの距離、グループ同士の距離、を定める必要がある。
 +またクラスタを生成する過程をグラフに表したものを、デンドログラムと呼ぶ。
 +
 +それに対して、クラスタ(群)ごとの群内のばらつきが最小になるように、目的関数を定め、離散最適化法を用いてクラスタを生成するのが、非階層クラスタリング\\
 +である。
 +代表的な非階層クラスタリングのひとつである$k$-means法では、各群のユークリッド距離の平均が最小になるように、クラスタ平均を付置する。
  
 === k-means法によるクラスタリング === === k-means法によるクラスタリング ===
行 99: 行 141:
     * ブースティング (Boosting)     * ブースティング (Boosting)
     * バギング (Bagging)     * バギング (Bagging)
 +
 +==== 層別 ====
  
 === k-means法による層別 === === k-means法による層別 ===
  
-k-means法は、数値変数でしか用いられないため、データは「tic.learn[, c(1:86)]」と数値変数のフィールドを指定してkmeans()関数を実行する。+k-means法は、数値変数でしか用いられないため、データは「tic.learn[, c(1:85)]」と数値変数のフィールドを指定してkmeans()関数を実行する。
 <code> <code>
-tic.learn.kmeans <- kmeans(tic.learn[, c(1:86)], centers=2, nstart=1000)+tic.learn.kmeans <- kmeans(tic.learn[, c(1:85)], centers=2, nstart=1000)
 </code> </code>
  
-k-means法の結果の「cluster」という変数を、tic.learnに追加する。+k-means法の結果を適用するには、レコードごとに「centers」の座標とのユークリッド距離を算出し、一番小さくなる層にその変数を割り当てる。 
 <code> <code>
-tic.learn$cluster <- tic.learn.kmeans$cluster+date() 
 +tic.learn$cluster <- rep(0, dim(tic.learn)[1] ) 
 +for( i in c(1:dim(tic.learn)[1]) ) { 
 +  tic.kmeans.dist <- rep(0, max(tic.learn.kmeans$cluster) ) 
 +  for( j in c(1:max(tic.learn.kmeans$cluster)) ) { 
 +    tic.kmeans.dist[j] <- sum( (tic.learn[i,c(1:85)]-tic.learn.kmeans$center[j,])^2 ) 
 +  } 
 +  tic.learn$cluster[i] <- sort.list(tic.kmeans.dist)[1] 
 +} 
 +date() 
 +</code> 
 + 
 +念のため、kmeansの結果と比較して、計算に誤りがないことを確認するためにクロス集計を行う。この結果が対角であれば、計算はあっている。 
 + 
 +<code> 
 +table(tic.learn$cluster,tic.learn.kmeans$cluster)
 </code> </code>
  
行 116: 行 176:
 === 変数の値を直接指定した層別(V1の例) === === 変数の値を直接指定した層別(V1の例) ===
  
-例えばV1で層別する場合、clusterという変数をtic.learn内に作り、まず0で埋める。 +例えばV1で層別する場合、まずその変数を集計して、何種類の値があるかを確認する。 
-そして、V1の値を参照しながら、層の番号を入力していく。+ 
 +<code> 
 +table(tic.learn$V1) 
 +</code> 
 + 
 +V1は1から41の値を取ることが確認できたら、次のように、clusterという変数をtic.learn内に作り、まず0で埋める。 
 +そして、V1の値ごとに層の値を指定するコードを書き、Rに実行させる。下記は、先週に行った決定木分析を踏まえた層別の例である。
  
 <code> <code>
行 164: 行 230:
 </code> </code>
  
-0の値が残っていないかを確認するために、集計する。+最後に、0の値が残っていないかを確認するために、集計する。 
 <code> <code>
 table(tic.learn$cluster) table(tic.learn$cluster)
 </code> </code>
 +
 +=== 変数の値を直接指定した層別(V2とV4の組み合わせの例) ===
 +
 +
 +例えばV2とV4の組み合わせで、下記のように層別することを考える。
 +
 +|V2/V4|1|2|3|4|5|6|
 +|1|1|1|1|1|3|3|
 +|2|2|2|2|2|3|3|
 +|3|3|3|3|3|3|3|
 +|4|3|3|3|3|3|3|
 +|5|3|3|3|3|3|3|
 +|6|3|3|3|3|3|3|
 +|7|3|3|3|3|3|3|
 +|8|3|3|3|3|3|3|
 +|9|3|3|3|3|3|3|
 +|10|3|3|3|3|3|3|
 +
 +この場合も、まずその変数を集計して、何種類の値があるかを確認する。
 +
 +<code>
 +table(tic.learn$V2, tic.learn$V4)
 +</code>
 +
 +すべての組み合わせを書くのは煩雑なので、最初にcluster変数は3で埋めてしまう。
 +
 +<code>
 +tic.learn$cluster <- rep(3, dim(tic.learn)[1])
 +</code>
 +
 +clusterが1になる箇所と2になる箇所だけ、書く。
 +
 +<code>
 +tic.learn$cluster[tic.learn$V2==1 & tic.learn$V4==1] <- 1
 +tic.learn$cluster[tic.learn$V2==1 & tic.learn$V4==2] <- 1
 +tic.learn$cluster[tic.learn$V2==1 & tic.learn$V4==3] <- 1
 +tic.learn$cluster[tic.learn$V2==1 & tic.learn$V4==4] <- 1
 +tic.learn$cluster[tic.learn$V2==2 & tic.learn$V4==1] <- 2
 +tic.learn$cluster[tic.learn$V2==2 & tic.learn$V4==2] <- 2
 +tic.learn$cluster[tic.learn$V2==2 & tic.learn$V4==3] <- 2
 +tic.learn$cluster[tic.learn$V2==2 & tic.learn$V4==4] <- 2
 +</code>
 +
 +
 +最後に、確認のために、集計する。
 +
 +<code>
 +table(tic.learn$cluster)
 +</code>
 +
 +=== 検証用データの層別 ===
 +
 +検証用データに層別変数 cluster を加えるには、上のコードの「tic.learn」をすべて「tic.eval」に変えればよい。
 +
 +==== 層別後の解析 ====
 +
  
 === 層別後の解析 === === 層別後の解析 ===
  
 +層ごとに回帰分析を適用したり、決定木を適用したりする。
  
 +<code>
 +for( i in sort(unique(tic.learn$cluster)) ) {
 +  tic.learn.cluster <- tic.learn[tic.learn$cluster == i,]
 +  tic.working.lm <- lm(V86~V2+V3+V4, data=tic.learn.cluster)
 +  print(summary(tic.working.lm))
 +}
 +</code>
 +
 +全体に回帰分析を行うのに較べて、残差の標準偏差 (Residual standard error) が小さくなっている層と、大きくなっている層があることを確認できるだろうか。
 +
 +<code>
 +print(summary(lm(V86~V2+V3+V4, data=tic.learn)))
 +</code>
 +
 +=== 考えたルールに基づく対象限定 (層別なし) ===
 +
 +各変数に閾値を設けてルールを生成したとする。
 +たとえば、「V47が5.5以上かつV44が1未満」または「V47が5.5以上かつV1が{1,3,6,8,12,20}のどれか」、というルールは
 +次のように記す。
 +
 +<code>
 +(tic.learn$V47>5.5 & tic.learn$V44<1) | (tic.learn$V47>5.5 & (tic.learn$V1==1 |tic.learn$V1==3 | tic.learn$V1==6 | tic.learn$V1==8 | tic.learn$V1==12 | tic.learn$V1==20) ) 
 +</code>
 +
 +「&」が「かつ(AND)」、「|」が「または(OR)」である。
 +
 +このルールを検証用データに適用するには、
 +
 +<code>
 +tic.learn$visit <- (tic.learn$V47>5.5 & tic.learn$V44<1) | (tic.learn$V47>5.5 & (tic.learn$V1==1 |tic.learn$V1==3 | tic.learn$V1==6 | tic.learn$V1==8 | tic.learn$V1==12 | tic.learn$V1==20) ) 
 +</code>
 +
 +と、訪問するか否かを二値(TRUE, FALSE)で表すオブジェクトを生成する。
 +このモデルに予測に基づいた訪問の成果を検証するには、訪問対象のリストtic.learn$visitと正解V86のクロス集計を行えばよい。
 +
 +<code>
 +table(tic.learn$visit)
 +table(tic.learn$visit, tic.learn$V86)
 +</code>
 +
 +参考までに、上の訪問ルールを評価用のデータtic.evalに適用した場合、次のようになる。
 +
 +<code>
 +table(tic.eval$visit)
 +FALSE  TRUE 
 + 3029   971 
 +
 +table(tic.eval$visit, tic.eval$V86)
 +tic.eval.visit    0    1
 +         FALSE 2878  151
 +         TRUE   884   87
 +</code>
 +
 +ここでは、訪問対象に884+87=971人を選定し、そのうちの87人が実際に契約してくれる人だったことになる。
 +契約率は87/971=8.96%。また誤判別率は
 +<code>
 +(884+151)/4000
 +</code>
 +で25.9%となる。
 +
 +=== 回帰分析に基づく対象限定 (層別なし) ===
 +
 +まず、すべての変数を用いて回帰分析を行う。
 +
 +<code>
 +tic.lm <- lm(V86~V1gr+V2+V3+V4+V5gr+V6+V7+V8+V9+V10
 +            +V11+V12+V13+V14+V15+V16+V17+V18+V19+V20
 +            +V21+V22+V23+V24+V25+V26+V27+V28+V29+V30
 +            +V31+V32+V33+V34+V35+V36+V37+V38+V39+V40
 +            +V41+V42+V43+V44+V45+V46+V47+V48+V49+V50
 +            +V51+V52+V53+V54+V55+V56+V57+V58+V59+V60
 +            +V61+V62+V63+V64+V65+V66+V67+V68+V69+V70
 +            +V71+V72+V73+V74+V75+V76+V77+V78+V79+V80
 +            +V81+V82+V83+V84+V85, data=tic.learn)
 +print(tic.lm)
 +summary(tic.lm)
 +plot(tic.lm)
 +</code>
 +
 +上で「object 'V1gr' not found」などとエラーが表示されたら、先週のデータの読み込みからやり直し。
 +この結果をtic.lmに収めたら、次はAICで変数選択をさせる。
 +
 +<code>
 +tic.lm.aic <- stepAIC(tic.lm)
 +</code>
 +
 +推定したモデルに基づいて、個別のレコードごとの契約の予測を立てる。その予測値が__0.05__以上なら訪問してみることとした場合、次のようになる。
 +
 +<code>
 +tic.learn$visit <- predict(tic.lm.aic, newdata=tic.learn)>0.05
 +</code>
 +
 +集計してみると。
 +<code>
 +table(tic.learn$visit)
 +table(tic.learn$visit, tic.learn$V86)
 +</code>
 +
 +=== 回帰分析に基づく対象限定 (層別あり) ===
 +
 +学習したモデルに基づいて、訪問対象を狭めるには、predict()という関数を用いて、訪問対象か否かというリストを作成する。
 +まずクラスタ1について、設定まで調整したモデルを、学習用データ(tic.learn)から得る。(rpart関数に指定してあるモデルは適当なので、各自の得たものを用いること)
 +
 +<code>
 +tic.lm.1 <- stepAIC(lm(V86~., data=tic.learn[tic.learn$cluster==1, ]))
 +</code>
 +
 +次に、このモデル(ここではtic.rpart.1)を、そのまま学習用データ(tic.learn)に適用して、契約してくれるか否かの予測を行う。
 +この際、下に記した0.05という閾値は調整(増やしたり減らしたり)が必要かもしれない。
 +
 +<code>
 +tic.learn$visit <- rep(0, dim(tic.learn)[1])
 +tic.learn$visit[tic.learn$cluster==1] <- predict(tic.lm.1, newdata=tic.learn[tic.learn$cluster==1,])>0.05
 +</code>
 +
 +同様のことを、クラスタごとに、すべてのクラスタについて行う。(ここではk=2)
 +
 +<code>
 +tic.lm.2 <- stepAIC(lm(V86~., data=tic.learn[tic.learn$cluster==2, ]))
 +tic.learn$visit[tic.learn$cluster==2] <- predict(tic.lm.2, newdata=tic.learn[tic.learn$cluster==2,])>0.05
 +tic.lm.3 <- stepAIC(lm(V86~., data=tic.learn[tic.learn$cluster==3, ]))
 +tic.learn$visit[tic.learn$cluster==3] <- predict(tic.lm.3, newdata=tic.learn[tic.learn$cluster==3,])>0.05
 +...
 +</code>
 +
 +この層別解析の結果をまとめる。
 +
 +<code>
 +table(1*(predict(tic.lm.1)>0.05), tic.learn$V86[tic.learn$cluster==1])
 +table(1*(predict(tic.lm.2)>0.05), tic.learn$V86[tic.learn$cluster==2])
 +table(1*(predict(tic.lm.3)>0.05), tic.learn$V86[tic.learn$cluster==3])
 +...
 +</code>
 +
 +全体の集計は、上のテーブルを足しても良いし、集計し直しても良い。
 +
 +<code>
 +table(tic.learn$visit, tic.learn$V86)
 +</code>
 +
 +=== 決定木に基づく対象限定 (層別なし) ===
 +
 +学習したモデルに基づいて、訪問対象を狭めるには、predict()という関数を用いて、訪問対象か否かというリストを作成する。
 +まず、設定まで調整したモデルを、学習用データ(tic.learn)から得る。
 +
 +<code>
 +tic.rpart <- rpart(V86~., data=tic.learn, control=c(cp=0.005))
 +</code>
 +
 +次に、このモデル(ここではtic.rpart)を、そのまま学習用データ(tic.learn)に適用して、契約してくれるか否かの予測を行う。
 +この際、下に記した__0.05__という閾値は、確率が0.05以上なら訪問する、という意味なので、調整(増やしたり減らしたり)が必要かもしれない。
 +
 +<code>
 +tic.learn$visit <- predict(tic.rpart, newdata=tic.learn)[,2]>0.05
 +</code>
 +
 +このモデルに予測に基づいた訪問の成果を検証するには、訪問対象のリストtic.visitと検証用データの正解V86のクロス集計を行えばよい。
 +
 +<code>
 +table(tic.learn$visit)
 +table(tic.learn$visit, tic.learn$V86
 +</code>
 +
 +参考までに、評価用データに適用した場合は、次のようになる。
 +
 +<code>
 +tic.eval$visit <- predict(tic.rpart, newdata=tic.eval)[,2]>0.05
 +</code>
 +
 +<code>
 +table(tic.eval$visit)
 +
 +tic.eval$visit 
 +FALSE  TRUE
 + 2389  1611
 +
 +table(tic.eval$visit, tic.eval$V86)
 +
 +tic.eval$visit    0    1
 +         FALSE 2310   79
 +         TRUE  1452  159
 +</code>
 +
 +ここでは、訪問対象に1452+159=1611人を選定し、そのうちの159人が実際に契約してくれる人だったことになる。契約率は159/1452=11.0%。
 +また誤判別率は
 +<code>
 +(79+1452)/4000
 +</code>
 +で38.275%となる。
 +
 +=== 決定木に基づく対象限定 (層別あり) ===
 +
 +学習したモデルに基づいて、訪問対象を狭めるには、predict()という関数を用いて、訪問対象か否かというリストを作成する。
 +まずクラスタ1について、設定まで調整したモデルを、学習用データ(tic.learn)から得る。(rpart関数に指定してあるモデルは適当なので、各自の得たものを用いること)
 +
 +<code>
 +tic.rpart.1 <- rpart(V86~., data=tic.learn[tic.learn$cluster==1, ], control=c(cp=0.005))
 +</code>
 +
 +あるいは
 +
 +<code>
 +tic.rpart.1 <- rpart(V86~., data=subset(tic.learn, (cluster==1) ), control=c(cp=0.005))
 +</code>
 +
 +とsubset()関数を用いて、データの一部を用いるように指定する。
 +
 +次に、このモデル(ここではtic.rpart.1)を、そのまま学習用データ(tic.learn)に適用して、契約してくれるか否かの予測を行う。
 +この際、下に記した0.05という閾値は調整(増やしたり減らしたり)が必要かもしれない。
 +
 +<code>
 +tic.learn$visit <- rep(0, dim(tic.learn)[1])
 +tic.learn$visit[tic.learn$cluster==1] <- predict(tic.rpart.1, newdata=tic.learn[tic.learn$cluster==1,])>0.05
 +</code>
 +
 +同様のことを、クラスタごとに、すべてのクラスタについて行う。(ここではk=2)
 +
 +<code>
 +tic.rpart.2 <- rpart(V86~., data=tic.learn[tic.learn$cluster==2, ], control=c(cp=0.005))
 +tic.learn$visit[tic.learn$cluster==2] <- predict(tic.rpart.2, newdata=tic.learn[tic.learn$cluster==2,])>0.05
 +tic.rpart.3 <- rpart(V86~., data=tic.learn[tic.learn$cluster==3, ], control=c(cp=0.005))
 +tic.learn$visit[tic.learn$cluster==3] <- predict(tic.rpart.3, newdata=tic.learn[tic.learn$cluster==3,])>0.05
 +...
 +</code>
 +
 +この層別解析の結果をまとめる。
 +
 +<code>
 +table(1*(predict(tic.rpart.1)>0.05), tic.learn$V86[tic.learn$cluster==1])
 +table(1*(predict(tic.rpart.2)>0.05), tic.learn$V86[tic.learn$cluster==2])
 +table(1*(predict(tic.rpart.3)>0.05), tic.learn$V86[tic.learn$cluster==3])
 +...
 +</code>
 +
 +全体の集計は、上のテーブルを足しても良いし、集計し直しても良い。
 +
 +<code>
 +table(tic.learn$visit, tic.learn$V86)
 +</code>
 +
 +ここも層別なしの解析と比較しておく。
 +
 +=== 検証 ===
 +
 +以上のように、層別も加えて分析して作り上げたモデルが、目標を達成しているかどうか、を検証用データ (tic.eval) で確認する。
 +
 +
 +==== レポート提出について ====
 +
 +レポート提出要領:下記「XXXXXXX」は各自の学籍番号(半角文字)で置き換えること
 +
 +^項目^指定^
 +|提出期限|実験実施の翌週の月曜日の午後7時0分まで (でも今回は12月25日から冬休みなので、12月25日の午後7時でいい)|
 +|提出方法|電子メールに添付 (宛先は配付資料に記載)|
 +|ファイル形式|Wordファイル (LaTeXで作成する場合は、dvipdfmxでPDFに変換すること)|
 +|メールの件名|統計実験3レポート提出(XXXXXXX)|
 +|レポートファイルの名称|統計実験3_XXXXXXX.doc あるいは 統計実験3_XXXXXXX.docx|
 +|提出部数|レポートは各自1通ずつ。{{:mselab:report-header-2012.doc|レポートの表紙}}に、共同実験者の学籍番号と氏名を記すこと。|
 +
 +
 +==== 実験当日の修正 ====
 +
 +自分のMacだと問題が発生せず、実験室のPCだとRが落ちる不具合があります。メモリが少ないためかもしれません。Rが勝手に終了する場合には、kmeansのコードは次のように nstart を削除し、iter.max を追加する修正をお願いします。
 +
 +=== k-means法による層別 (iter.max追加、nstart削除) ===
 +
 +k-means法は、数値変数でしか用いられないため、データは「tic.learn[, c(1:85)]」と数値変数のフィールドを指定してkmeans()関数を実行する。
 +<code>
 +tic.learn.kmeans <- kmeans(tic.learn[, c(1:85)], centers=2, iter.max=100)
 +</code>
 +
 +k-means法の結果を適用するには、レコードごとに「centers」の座標とのユークリッド距離を算出し、一番小さくなる層にその変数を割り当てる。
 +
 +<code>
 +date()
 +tic.learn$cluster <- rep(0, dim(tic.learn)[1] )
 +for( i in c(1:dim(tic.learn)[1]) ) {
 +  tic.kmeans.dist <- rep(0, max(tic.learn.kmeans$cluster) )
 +  for( j in c(1:max(tic.learn.kmeans$cluster)) ) {
 +    tic.kmeans.dist[j] <- sum( (tic.learn[i,c(1:85)]-tic.learn.kmeans$center[j,])^2 )
 +  }
 +  tic.learn$cluster[i] <- sort.list(tic.kmeans.dist)[1]
 +}
 +date()
 +</code>
 +
 +念のため、kmeansの結果と比較して、計算に誤りがないことを確認するためにクロス集計を行う。この結果が対角であれば、計算はあっている。
 +
 +<code>
 +table(tic.learn$cluster,tic.learn.kmeans$cluster)
 +</code>
 +
 +centers=に指定する値は、いろいろ変えてみるといい。