はじめに
先日はRとXgboostのインストールおよび動作確認をしたので、本日はKaggleのチュートリアルであるタイタニックのタスクに参加する。「Data」にtrain.csvとtest.csvがあるのでダウンロードする。
関連投稿(追記)
データの確認
入っているのは以下のデータ。
属性名 | 値 | 備考 |
---|---|---|
PassengerId | 固有ID | |
Survived | 生存 | 0:No, 1:Yes |
Pclass | クラス | 1:1st, 2:2nd, 3:3rd |
Name | 名前 | (見方がわからない) |
Sex | 性別 | male,female |
Age | 年齢 | ブランクあり、推定値は小数点x.xx |
SibSp | 同乗する兄弟・夫婦の数 | 合計値 |
ParCh | 同情する両親・子供の数 | 合計値 |
Ticket | チケットナンバー | 同じチケット番号がある、表記法に一貫性がない |
Fare | 運賃 | |
Cabin | キャビン | ブランクあり、複数付いている場合がある |
Embarked | 乗船した場所 | C:Cherbourg, Q:Queenstown, S:Sourthamption |
なお、親戚や友人関係の情報はない。 名前は見方がよくわからないが、少なくとも「○○, □□」という形式にはなっていて、○○の部分がFamilyNameとして使えそうな印象。□□の部分は正直よくわからなくて、ここだけでFastName+FamilyNameの形式にも見えるし、括弧がついてる場合があってその括弧の中に別の人の名前がある。括弧の中の人名は検索してもヒットしない。MrとかMrsとかついてるので既婚か未婚かくらいの属性は追加で付けられそう。 チケットは重複があり、表記法に規則性は見つけられなかった。若干、Pclassと連動しているように思える。チケットの番号とキャビンが連動していれば、なんとなく宿泊位置関係が見えてくるかもしれない。 乗船場所からはPclassと関係がありそう。Sourthamptionからの乗船は3rdクラスが多い。そして3rdクラスは死亡率が高い。 年齢については推測するのが良いようです。年齢以外の素性を使って、生存判定が一番あがるようにEMアルゴリズム等でチューニングするのが一般的っぽいです。
素性設計1
まずは以下の素性でやってみます。
属性名 | 値 | 備考 |
---|---|---|
Pclass | クラス | 1:1st, 2:2nd, 3:3rd |
FamilyName | 名前 | (頭から最初のコンマまで) |
Sex | 性別 | male,female |
Age | 年齢 | (ブランク->-1, 四捨五入, 文字列化) |
Embarked | 乗船した場所 | C:Cherbourg, Q:Queenstown, S:Sourthamption |
とてもシンプル。識別器はXgboostを使います。パラメータはデフォルト値です。Ageが不明な場合、下手な推定しないで不明(-1)としました。また、数値のままで扱うと素性の強さを表すことになってしまうので、文字列としました。本来であればレンジを決めた方がいいですが、面倒なので四捨五入くらいにしておきました。 5-fold cross validationの結果は正解率80.4%でした。少し設計を変えてみましょう。
素性設計2
以下の素性でやってみます。
属性名 | 値 | 備考 |
---|---|---|
Pclass | クラス | 1:1st, 2:2nd, 3:3rd |
FamilyStatus | 家族の生死 | U:家族情報なし, A:生者が多い, D:死者が多い |
Sex | 性別 | male,female |
Age | 年齢 | (ブランク->-1, 四捨五入, 文字列化) |
Embarked | 乗船した場所 | C:Cherbourg, Q:Queenstown, S:Sourthamption |
一緒に居た家族の生死をFamilyStatusという素性にしました。生死に関わる一番の部分は災害時にどこに居たかだと考え、家族は近くにいるものなので追加した素性です。テスト時はトレイン時のデータに加え、テストで予測した結果も使います。もちろん、CVで不正が発生しないように気をつけてコーディングします。生者と死者が同数の場合は「D:死者が多い」とします。これは生者に比べて死者が多かったという事実を加味してのことです。家族かどうかはFamilyName素性が一致するかどうかをみます。 5-fold cross validationの結果は正解率79.9%でした。微妙に下がってますね。有意差検定していないのでなんともいえませんが、多分誤差の範囲だと思います。せっかくなのでこのままいきます。
パラメータの最適化
本来であれば、交差検定に粗密探索を組み合わせてパラメータチューニングをします。面倒なのでいくつか実験して適当に決めました。Xgboostのパラメータは以下のとおりです。
パラメータ | 値 |
---|---|
nrounds | 10 |
eta | 0.1 |
gamma | 0.1 |
max_depth | 6 |
min_child_weight | 10 |
subsumple | 1 |
colsumple_bytree | 1 |
objective | binary:logistic |
5-fold cross validationの結果は正解率81.3%でした。まじめにやればもう少し上がると思います。
Kaggleへ投稿
このアルゴリズムでどこまでいくか試してみます。Kaggleのタイタニックのページに行って「Make a Submission」から投稿します。投稿フォーマットはCSVで「PassengerId,Survived」のフォーマットになっていればOKです。 さてさて、気になる結果は・・・ 正解率79.4%で1293位(4491人中)でした。 思ったより層が厚いですね。ちょっと悔しいです。上位層はランダムフォレストを使っているそうなので、次はランダムフォレストを使ってみましょうか。もしやったら、こちらの記事に追記します。
終わりに
機械学習では「機械学習の枠組みに素性選択が入ってるから、思いつく素性を全部入れとけ」という人もいます。確かにそういう側面もありますが、機械学習の大原則は「各次元が独立であること」です。なので、似たような意味合いを持つものを入れるのはナンセンスです。私のこれまでの経験でも、素性を何でもとにかく多く突っ込めば性能が上がる、ということはなかったです。むしろノイズっぽいものを入れると性能が下がります。ランダムフォレストはアルゴリズム的にそのあたりは比較的強く出来ているので、何も考えずにポンポン投げ込むならランダムフォレストが良いという印象です。 さて、毎日ブログを書くと決意してこのブログを始めましたが、さすがにプログラム書いて実験しながらだと少し時間がかかっちゃいます。この記事の投稿のために、3日かかりました。今回はじめてRに触れましたので、使った時間の大半はRの使い方をググってたという感じです。少し慣れてきたので、次のタスクからはもう少し早くサイクルを回せそうです。 肝心のKaggleについてですが、まだチュートリアルをやっただけですが、遊び感覚で機械学習やデータマイニングができていいですね。懸賞金付きのタスクも遊び感覚でやってみようと思います。
ソースコード
library(xgboost)
library(Matrix)
set.seed(1)
# Load Data
data = read.csv("train.csv")
data$Pclass = factor(ifelse(is.na(data$Pclass), as.character(-1), as.character(data$Pclass)))
data$Age = factor(ifelse(is.na(data$Age), as.character(-1), as.character(round(data$Age))))
data$FamilyName = factor(gsub(",.+$","",data$Name))
# Create data for k-fold cross validation
cv = function(d, k) {
n = sample(nrow(d), nrow(d))
d.randomized = data[n,] # randomize data
n.residual = k-nrow(d)%%k
d.dummy = as.data.frame(matrix(NA, nrow=n.residual, ncol=ncol(d)))
names(d.dummy) = names(d)
d.randomized = rbind(d.randomized, d.dummy) # append dummy for residuals
d.splitted = split(d.randomized, 1:k)
for (i in 1:k) {
d.splitted[[i]] = na.omit(d.splitted[[i]])
}
d.splitted
}
# train data
cv.train = function(d, k) {
d.train = as.data.frame(matrix(0, nrow=0, ncol=ncol(d[[1]])))
names(d.train) = names(d[[1]])
for (i in 1:length(d)) {
if (i != k) {
d.train = rbind(d.train, d[[i]])
}
}
d.train
}
# test data
cv.test = function(d, k) {
d[[k]]
}
# 家族の生死のリスト準備
familystatus = function(d) {
rows = nrow(d)
rtn = NULL
for (i in 1:rows) {
if (d[i,"SibSp"]==0 && d[i,"Parch"]==0) {
# do nothing
} else {
name = as.character(d[i,"FamilyName"])
survive = d[i,"Survived"]
if (survive==0) survive=-1
if (is.null(rtn[name]) || is.na(rtn[name])) {
rtn[name]=survive
} else {
rtn[name]=rtn[name]+survive
}
}
}
rtn
}
# 家族の生死素性化
addfamilystatus = function(d,l) {
rows = nrow(d)
rtn = NULL
for (i in 1:rows) {
if (d[i,"SibSp"]==0 && d[i,"Parch"]==0) {
rtn = c(rtn,"U")
} else {
name = as.character(d[i,"FamilyName"])
if (is.na(l[name])) {
rtn = c(rtn,"U")
} else if (l[name]>0) {
rtn = c(rtn,"A")
} else {
rtn = c(rtn,"D")
}
}
}
rtn = factor(rtn)
rtn
}
# 家族の生死のリストアップデート
updatefamilystatus = function(d,l) {
rows = nrow(d)
rtn = l
for (i in 1:rows) {
if (d[i,"SibSp"]==0 && d[i,"Parch"]==0) {
# do nothing
} else {
name = as.character(d[i,"FamilyName"])
survive = d[i,"Survived"]
if (survive==0) survive=-1
if (is.null(rtn[name]) || is.na(rtn[name])) {
rtn[name]=survive
} else {
rtn[name]=rtn[name]+survive
}
}
}
rtn
}
## load data as sparse matrix
k = 5
data.splitted = cv(data, k)
## Xgboost
rtn.gold = NULL
rtn.pred = NULL
for (i in 1:k) {
data.train = cv.train(data.splitted, i)
fsurvivelist = familystatus(data.train)
data.train$FamilyStatus = addfamilystatus(data.train,fsurvivelist)
train = sparse.model.matrix(Survived~Pclass+FamilyStatus+Sex+Age+Embarked, data.train)
train.label = data.train$Survived
# train model
bst = xgboost(data=train,
label=train.label,
nrounds=10,
eta=0.1,
gamma=0.1,
max_depth=6,
min_child_weight=10,
subsumple=1,
colsumple_bytree=1,
objective="binary:logistic",
verbose=0)
# test
data.test = cv.test(data.splitted, i)
rows = nrow(data.test)
for (j in 1:rows) {
dtest = data.test
dtest$FamilyStatus = addfamilystatus(dtest,fsurvivelist)
test = sparse.model.matrix(Survived~Pclass+FamilyStatus+Sex+Age+Embarked, dtest[j,])
test.label = data.test[j,"Survived"]
# predict test data
pred = predict(bst, test)
prediction = as.numeric(pred > 0.5)
table(test.label, prediction)
acc = mean(prediction == test.label)
rtn.gold = c(rtn.gold,test.label)
rtn.pred = c(rtn.pred,prediction)
# update fsurvivelist
dtest = dtest[j,]
dtest$Survived = prediction
fsurvivelist = updatefamilystatus(dtest,fsurvivelist)
}
}
table(rtn.gold, rtn.pred)
acc = mean(rtn.gold == rtn.pred)
print(paste("Xgboost test-total-accuracy=", acc))