RのFactoRizationMachinesパッケージを用いたFMのページビューデータへの適用

Googleタグマネージャーで集めたアクセスログデータを用いて、前回と同様に記事のレコメンドにチャレンジしてみようと思います。FactoRizationMachinesパッケージという便利そうなパッケージの存在も知れたことから、今回は以前から気になっていたFactorization Machineを扱います。

【目次】
・Factorization Machine(FM)の概要
・パッケージ紹介とインストール
・サンプルデータの構造把握と前処理
・FMの実行
・結果
・参考文献

Factorization Machine(FM)の概要

組み合わせ特徴量を扱う教師あり学習モデル。行列分解とSVMを合体させた手法。
スパースになりやすいデータの予測問題で扱う。
1ユーザーのある商品に対しての評価を、1評価1行として表して、ユーザーとアイテムの交互作用の特徴ベクトルを扱う。
相互作用項に関して、時間や文脈などを自由に入れられる。
相互作用項を次元圧縮する際の要素数を事前に決める必要がある。
Matrix Factorizationよりも精度が良いとされている。特徴量エンジニアリングなどで使われているようです。（Click-Through Rate Prediction）

パッケージ紹介とインストール

FactoRizationMachinesパッケージは線形SVMと2次のFMと高次のFMを実行することができ、引数で正則化項も加えることができます。現段階においては回帰のみで分類問題への適用は今後の開発となるようです。CRANから普通に
install.packages(‘FactoRizationMachines’)
でインストールします。libFMexeパッケージの場合は、libFMをインストールしてパスを指定しておく必要がありますが、このパッケージに関しては不要となります。

サンプルデータの構造把握と前処理

FactoRizationMachinesパッケージのサンプルコードにおいては、MovieLensのデータがサンプルデータとして載せられていました。ユーザーのID(整数)、映画のID(整数)、評価(整数、5段階)、日時(整数)からなるデータに対して、sparseMatrixに変換していました。

今回は、前回の投稿非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみると同じデータを使って、アクセスログデータに適用しようと思います。FactoRizationMachinesの形式に合わせるために、このブログのアクセスログも、クッキーのIDを整数に、記事のIDを整数に、閲覧回数を5段階(5以上を5に変換)に、日時を整数に変更しています。

library(RGA)
library(tidyverse)
library(reshape2)

authorize()
prof <-list_profiles()

start_date <- "2017-07-04"
end_date   <- "2017-10-07"

#demention until 7
accesslogdata <- get_ga(profileId = prof$id[2],
                        start.date = start_date,
                        end.date = end_date,
                        dimensions = "ga:pagePath,
                        ga:dimension1,
                        ga:dateHourMinute,
                        ga:deviceCategory,
                        ga:userType,
                        ga:fullReferrer",
                        sort = "-ga:sessions",
                        metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",
                        fetch.by = "day")

accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))
accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))
accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

#ユーザーIDの整数化
user <- unique(accesslogdata$dimension1)
user <- data.frame(dimension1=user) %>% mutate(user_number=1:n())

#記事IDの整数化
articles <- unique(accesslogdata$pagePath)
articles <- data.frame(pagePath=articles) %>% mutate(article_number=1:n())

#データの結合
accesslogdata <- accesslogdata %>% mutate(date=as.numeric(as.POSIXlt(as.Date(format(substr(accesslogdata$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d"))))
accesslogdata <- accesslogdata %>% left_join(articles,by="pagePath")
accesslogdata <- accesslogdata %>% left_join(user,by="dimension1")

#最終接触日の抽出
last_date <- accesslogdata %>% select(user_number,article_number,date) %>% 
                arrange(user_number,desc(date))
last_date <- last_date[!duplicated(last_date[c("user_number","article_number")]),]

#各ユーザーの各記事に対するページビュー数の集計
accesslogdata_sum <- accesslogdata %>% group_by(article_number,user_number) %>% summarise(pageview=n()) %>% 
                        select(user_number,article_number,pageview)
accesslogdata_sum <- accesslogdata_sum %>% mutate(pageview=ifelse(pageview < 5,pageview,5))

accesslogdata_sum <- accesslogdata_sum %>% left_join(last_date, by = c("user_number" = "user_number", "article_number" = "article_number"))

library(RGA)

library(tidyverse)

library(reshape2)

authorize()

prof <-list_profiles()

start_date <- "2017-07-04"

end_date <- "2017-10-07"

#demention until 7

accesslogdata <- get_ga(profileId = prof$id[2],

start.date = start_date,

end.date = end_date,

dimensions = "ga:pagePath,

ga:dimension1,

ga:dateHourMinute,

ga:deviceCategory,

ga:userType,

ga:fullReferrer",

sort = "-ga:sessions",

metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",

fetch.by = "day")

accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))

accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))

accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

#ユーザーIDの整数化

user <- unique(accesslogdata$dimension1)

user <- data.frame(dimension1=user) %>% mutate(user_number=1:n())

#記事IDの整数化

articles <- unique(accesslogdata$pagePath)

articles <- data.frame(pagePath=articles) %>% mutate(article_number=1:n())

#データの結合

accesslogdata <- accesslogdata %>% mutate(date=as.numeric(as.POSIXlt(as.Date(format(substr(accesslogdata$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d"))))

accesslogdata <- accesslogdata %>% left_join(articles,by="pagePath")

accesslogdata <- accesslogdata %>% left_join(user,by="dimension1")

#最終接触日の抽出

last_date <- accesslogdata %>% select(user_number,article_number,date) %>%

arrange(user_number,desc(date))

last_date <- last_date[!duplicated(last_date[c("user_number","article_number")]),]

#各ユーザーの各記事に対するページビュー数の集計

accesslogdata_sum <- accesslogdata %>% group_by(article_number,user_number) %>% summarise(pageview=n()) %>%

select(user_number,article_number,pageview)

accesslogdata_sum <- accesslogdata_sum %>% mutate(pageview=ifelse(pageview < 5,pageview,5))

accesslogdata_sum <- accesslogdata_sum %>% left_join(last_date, by = c("user_number" = "user_number", "article_number" = "article_number"))

FMの実行

デフォルトの設定c(1, 10)では線形のウェイトが有効で、2次の項の要素数が10で正則化項なしのFMを実行することになります。引数に関する詳しい情報はPackage ‘FactoRizationMachines’に書かれています。今回はサンプルを参考に正則化項ありでモデルを実行します。まず、アクセスログデータに対して、ユーザーのIDからなる整数ベクトル、記事のIDからなる整数ベクトル、セッションのあった日時のデータからなる整数ベクトルを作成し、sparseMatrix関数を用いて元データを変形し、80%のデータをトレーニングに、20%のデータをテストに割り当てます。さらに、テストデータに関して、予測値との平均二乗誤差を計算します。

library(FactoRizationMachines)
library(Matrix)

user = accesslogdata_sum$user_number
items = accesslogdata_sum$article_number + max(user)
wdays = (as.POSIXlt(accesslogdata_sum$date,origin="2017-07-04")$wday+1)+max(items)

# Transform access log to feature form
data = sparseMatrix(i=rep(1:nrow(accesslogdata_sum),3),j=c(user,items,wdays),giveCsparse=F)
target = accesslogdata_sum$pageview

# Subset data to training and test data
set.seed(123)
subset = sample.int(nrow(data),nrow(data)*.8)
data.train = data[subset,]
data.test = data[-subset,]
target.train = target[subset]
target.test = target[-subset]

# Predict ratings with Support Vector Machine with linear kernel
model = SVM.train(data.train,target.train)

# RMSE resulting from test data prediction
sqrt(mean((predict(model,data.test)-target.test)^2))

# Predict ratings with second-order Factorization Machine
# with second-order 10 factors (default) and regularization
model = FM.train(data.train,target.train,regular=0.1)
# RMSE resulting from test data prediction
sqrt(mean((predict(model,data.test)-target.test)^2))

library(FactoRizationMachines)

library(Matrix)

user = accesslogdata_sum$user_number

items = accesslogdata_sum$article_number + max(user)

wdays = (as.POSIXlt(accesslogdata_sum$date,origin="2017-07-04")$wday+1)+max(items)

# Transform access log to feature form

data = sparseMatrix(i=rep(1:nrow(accesslogdata_sum),3),j=c(user,items,wdays),giveCsparse=F)

target = accesslogdata_sum$pageview

# Subset data to training and test data

set.seed(123)

subset = sample.int(nrow(data),nrow(data)*.8)

data.train = data[subset,]

data.test = data[-subset,]

target.train = target[subset]

target.test = target[-subset]

# Predict ratings with Support Vector Machine with linear kernel

model = SVM.train(data.train,target.train)

# RMSE resulting from test data prediction

sqrt(mean((predict(model,data.test)-target.test)^2))

# Predict ratings with second-order Factorization Machine

# with second-order 10 factors (default) and regularization

model = FM.train(data.train,target.train,regular=0.1)

# RMSE resulting from test data prediction

sqrt(mean((predict(model,data.test)-target.test)^2))

結果

各モデルについての平均二乗誤差を計算しています。
線形モデルや高次元モデルよりも、2次の項を持つFMが精度が高いようです。

> model = SVM.train(data.train,target.train)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.514266
> model = FM.train(data.train,target.train,regular=0.1)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.369983
> model　=　HoFM.train(data.train,target.train,c(1,3,1),regular=0.1)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.388247

> model = SVM.train(data.train,target.train)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.514266

> model = FM.train(data.train,target.train,regular=0.1)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.369983

> model　=　HoFM.train(data.train,target.train,c(1,3,1),regular=0.1)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.388247

こちらは、この中で性能の良かった2次の項を持つFMの予測結果とテストデータの結果をプロットしたものです。4点を超える値をあまり予測できていないようです。今回はサンプルを回しただけなので、本来であれば次元の数kや正則化のセッティングをいろいろいじったり、相互作用項を新しく追加するなどして精度を高めることが必要です。

結果の比較だけでは仕事で使えないので、実際に予測した結果を取り出したいと思います。
実際に運用するとなると、ページIDを所与として、ページビュー数を0とおいて（型をそろえるため。NULLだとエラーになった）、任意のタイミング(date)を想定して、モデルにデータを適用し、評価の高いものをサジェストするスタイルになるのではないでしょうか。

この結果だと、ユーザー98に記事1を見せることに対して4.02点が与えられています。

参考文献

Factorization Machinesを今更読みました
 Factorization Machines
High-order factorization machines with R #tokyor 61
Factorization Machinesのおはなし。
libFMexeを動かすまで (R Wrapper for the libFM Executable参照記事)
一歩Matrix Factorization、二歩Factorization Machines、三歩Field-aware Factorization Machines…『分解、三段突き！！』
[論文] Factorization Machines (ICDM 2010) 読んだ 22:41
Factorization machines with r
Factorization Machinesについて調べてみた

人工知能学会全国大会2017のWebマーケティングで参考になりそうな研究9選

はじめに

今更ではありますが、2017年5月に開かれた人工知能学会全国大会2017の公開資料の中から、私の本業であるWebマーケティングで役に立ちそうな研究を独断と偏見で9本ほど集めてみました。思っていたよりもWebマーケティングぽい領域が多かったので、社内での分析業務におけるアイデアに直結しそうな気がしています。

ちなみに、全ての資料はこちらにあります。

基本的に各研究ごとに
・目的と結果
・対象となるデータ
・手法の概要
・PDFのリンク
について記していきます。

紹介する研究

今回は以下の研究を紹介したいと思います。
「オンライン広告におけるスパース性と遅れコンバージョンを考慮した予測モデル」「テキストマイニングを用いた転職サイトの会員離脱予測」「インターネット広告におけるスパースなユーザー行動ベクトルからのユーザー特徴抽出モデル」「ユーザーのページビュー系列からのコンバージョン予測」「SNSにおけるユーザの行動のモチベーション要因の分析」「状況に応じた楽曲推薦に向けたソーシャルデータ分析」「テキストマイニングを用いた口コミ分析による点数評価の信頼性確認手法」「不動産仲介マーケティングのためのユーザ行動予測」「SUUMOでの不動産データ活用の取り組みと未来」

オンライン広告におけるスパース性と遅れコンバージョンを考慮した予測モデル

目的と結果

オンライン広告のCVR予測の精度向上を目的としている。
これまでの課題としては、
「質的変数による特徴量のスパース性」
「広告クリックからのCVまでの期間の長さ（期間が空いてからCVしたりするケースがあるので正例だったはずが負例と扱われる可能性がある）」
などがあった。提案手法により従来手法に比べ高い精度でCVR予測が可能であることが示された。

対象となるデータ

過去3週間分の広告配信データでテスト期間は7日間

手法の概要

次元圧縮
- Factorization Machine(FM)を用いて、従来のロジスティック回帰では考慮できていなかった交互作用をモデリングし、低次元のベクトル内積に分解して表現している。
遅れCV問題の対応
- 遅れてCVが発生する事象に関しては、指数分布を用いて遅れてCVする事象をモデリングするDelayed Feedback Modelを構築することで、本当は事後的にCVしそうなユーザーの推定を行っている。
予測
- これらの手法を合わせた、Delayed Feedback Factorization Machineを構築し、SGD(確率的勾配降下法)やAdaGradを用いてパラメータを推定している。

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/110.pdf

テキストマイニングを用いた転職サイトの会員離脱予測

目的と結果

登録情報や履歴書のテキスト情報を用いて特徴量とし、転職サイトの登録会員の離脱予測を行っている。ブラックボックスでは無い形での結果の解釈ができ、予測精度も高いことが示された。

対象となるデータ

2009年10月1日〜2016年9月28日までの履歴書、職務経歴書、自己PR書、希望条件
離脱会員の定義は「登録して1ヶ月以内に一度も応募しない会員」としている。

手法の概要

STEP1：Pythonのmojimojiというライブラリを利用してカナを全角に、数字とアルファベットは半角に変換。
STEP2：MeCabを用いて名詞のみ抽出
STEP3：「職務」「経歴」「業務」「内容」などを不用語（stop word）とした。
STEP4：Grid Searchでパラメータチューニングしたランダムフォレストを用いて変数重要度（Feature Importance）を出し、重要度の上位90%を占める単語以外を削除
STEP5：登録情報のデータ（年齢、転職回数、住所）を数値データに変換
STEP6：選定した特徴量を用いて決定木アルゴリズムを適用

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/949.pdf

インターネット広告におけるスパースなユーザー行動ベクトルからのユーザー特徴抽出モデル

目的と結果

広告のCTR予測を行う際の特徴量が非常にスパースであることから、ユーザーをクラスタリングして広告配信に効果的なユーザーを抽出することを目的としている。研究ではCVRの高いユーザー群の抽出を可能にしている。

対象となるデータ

2016年6月1日〜2016年6月21日の間のクリックやコンバージョンのデータ（株式会社アイモバイルが提供）
クリック実績のある広告数、コンバージョン実績のある広告数、クリックされた広告番号、コンバージョンが発生した広告番号、ユーザー番号など

手法の概要

高次元スパースデータを低次元に落とし込むために、オートエンコーダを用いている。
オートエンコーダにより得られた特徴量をもとにユーザーのクラスタリングを行っている。(オートエンコーダは入力された情報をなるべく欠損しないように、効率的に圧縮する構造を持たせたニューラルネットワークのこと。)
データにおけるゼロの成分を確率的に欠損とみなして計算処理を早めている。
学習モデルはAdaGrad(学習率を学習に合わせて最適化する手法)
訓練データ（特徴量の数がクリックで23339個、コンバージョンで5619個）にスパースデータ許容型オートエンコーダで500次元まで圧縮し、さらに通常のオートエンコーダで146次元まで圧縮している。
次元圧縮して獲得したデータに対してk-means法でクラスタリングを行い、CVRなどを比較している。

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/12.pdf

ユーザーのページビュー系列からのコンバージョン予測

目的と結果

Web広告において、ユーザーの閲覧履歴からコンバージョンを予測することを目的に、Recurrent Neural Network(RNN)を用いた結果、非時系列のSVMよりもわずかばかりかの精度向上が観察された。データ数が少ない場合はSVMに軍配があがる結果となっている。

対象となるデータ

2016年7月1日〜2016年10月31日の広告の接触履歴データなど（広告の業種別データやユーザーのページビュー数）を特徴量に用いて、2016年11月1日〜11月30日までの期間を予測対象とした。データは株式会社マイクロアドが提供している模様。

手法の概要

Long Short-Term Memory RNN(LSTM)
- ソフトアテンションモデル
  - 可変長の系列に対して内積で重みを算出し、足し合わせることで一部から情報を取り出すモデル。途中の層まで複数業種のデータを用いることで、複数の業種などの情報を共有することができる。（最後の隠れ層には業種別の特徴量を用いている。）
  - 勾配の計算に関してはAdamを用いて、実装はKerasを用いている。活性化関数はTanh、ドロップアウトはGaussianDropout
    としている。
- RNNで用いる特徴量としては、特定のWebページを640種類に分類して、1日のうちに各トピックにページビューがあったかなかったかを1-0の変数にするなどの前処理を行っている。

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/765.pdf

SNSにおけるユーザの行動のモチベーション要因の分析

目的と結果

SNSユーザーのモチベーションの要因（一ヶ月先に投稿するかどうか）をユーザー同士の関係性（staticな関係）と、そのユーザーの周りのインタラクション（dynamicな関係）という観点から分析し、ユーザー同士の関係性とユーザー周りのインタラクションを同時に用いたモデルにおいて、それらがモチベーションの要因につながりうる傾向が観察された。

対象となるデータ

RoomClipという部屋の家具やインテリアの様子の写真を投稿するSNSのデータ
- 気に入った写真に「いいね」を付けたり、コメントを送れる。
2015年3月の一ヶ月間を学習期間とし、2015年4月の投稿の有無を予測対象とする。

手法の概要

ネットワーク特徴量
- 各ユーザーに対して、
  - 中心性（次数中心性、媒介中心性、近接中心性）
  - クラスタ係数、PageRank、フォロー数とフォロワー数の比率、フォローしているユーザのアクティブ度
フィードバック特徴量
- 特定期間中にユーザが周りのユーザから自身の投稿に対してもらった「いいね」「コメント」に対して
  - 反応の数
  - 反応を行ったユーザの数
  - 反応を行ったユーザのうち、相互フォロー関係にあるユーザからのものの割合
予測器としてSVMを利用

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/534.pdf

状況に応じた楽曲推薦に向けたソーシャルデータ分析

目的と結果

ユーザへの楽曲推薦システムの構築を目的に、楽曲とその再生時の状況の関係を収集したソーシャルデータより分析している。取得したデータのうちの7%ほどではあるが、ツィート情報から楽曲と再生時の状況との関係を抽出可能であるということが示された。

対象となるデータ

2016年8月10日〜2016年9月3日の間に投稿されたツィート（66879件）

手法の概要

ハッシュタグ#nowplaying付きのツィートを収集して、アーティスト名や楽曲名や楽曲再生時の状況（同一ユーザの5分前後の投稿内容）を収集
取得したツィートからアーティスト名、楽曲名を除去して再生時の状況を取得
取得したデータからデータベースを作成

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/102.pdf

テキストマイニングを用いた口コミ分析による点数評価の信頼性確認手法

目的と結果

製品などの口コミの評価レビューに関する信頼性を確かめることを目的としている。レビューデータに対して感性評価を行い、さらにその結果に対して総合評価点数別に感性表現（Goodなど）の伴った表現の頻出ランキングを作成することで、点数が高い理由を把握しやすくなった。

対象となるデータ

ソニーの製品サイトにおける、ユーザの製品に対する評価とコメントのデータ（1406件）

手法の概要

テキストマイニングスタジオによる感性評価（ポジティブ、中立、ネガティブ）
総合評価点数別（5段階）の「ポジティブ、中立、ネガティブ」の割合の調査
総合評価4の口コミかつ中立な表現にフィルターをかけて、感性情報を有する単語と係り受け関係となった単語の頻度ランキングのTOP10を抽出

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/1115.pdf

不動産仲介マーケティングのためのユーザ行動予測　Prediction of user behaviors for real estate brokerage marketing

目的と結果

不動産紹介業において、顧客の属性情報と来店初期の行動から、来店者が最終的に商品を購買するかどうかの予測モデル(不均衡データにおける予測タスク)を構築することでマーケティングの広告効果測定に適用することを目的としている。結果として、顧客の来店後30分間という短い時間における行動だけでも識別力のある予測が可能であることが示されている。

対象となるデータ

オンライン不動産賃貸仲介店舗iettyの顧客データベース
- 希望物件の情報
  - 引っ越し希望日
  - 住みたいエリアの駅
  - 間取り
  - 広さ
  - 家賃
- 顧客属性
  - 在住の地域
  - 性別
  - 年齢
- 行動データ
  - 登録時
    - フォームの記入にかかった時間
    - フリーコメントの長さ
    - 利用デバイス
  - 登録後
    - 初回チャット発言までの経過時間
    - 初回物件評価までの経過時間

手法の概要

STEP1：顧客の属性および登録初期(30分以内)の行動から特徴量を作成し、2ヶ月以内に契約すれば正例、それ以外を負例とする。（正例は619件、負例は33212件）
STEP2：検証用データを無作為に20%取り置く
STEP3：残り80%のデータに対して負例に対する、ランダムアンダーサンプリングと、正例に対するSMOTE(Synthetic Minority Over-sampling Technique)を用いたオーバーサンプリングによりサンプル数を調整する。
STEP4：ランダムフォレストを用いて、無効な特徴量の取捨選択や離散化のレンジ調整する。
STEP5：全特徴量（83次元）を二値化して、RBFカーネル(Gaussian カーネル)のSVMを適用し、様々なCやγの組み合わせを試行して、二値分類を行う。
STEP6：評価を予測確率の平均とAUCの2種類で行う。

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/466.pdf

SUUMOでの不動産データ活用の取り組みと未来

目的と結果

物件検索ユーザーが比較検討を行いやすくすることを目的に様々な取り組みをしている。

対象となるデータ

不動産ポータルサイトSUUMOのデータ
- 行動履歴(Webサイト、アプリ)
- 物件情報(不動産情報)
- 街情報(独自で収集)
- アンケート

手法の概要

実活用
- 相場推定
- 単純な加算平均ではなく、外れ値や時系列変動、データ量などを考慮して推定している。
- リアルタイム物件レコメンド
- ランキング学習により推定したモデルを利用し、リアルタイムに計算を行っている。個人ごとにパーソナライズしたランキングも用いている。
- バンディットアルゴリズムによる配信施策の最適化
- メールやプッシュ通知などで件名やデザイン、配信時間など様々な選択肢があり、その最適なパターンを見つけるために、バンディットアルゴリズムを適用している。
実験段階
- 間取り画像の分析
- ユーザーの求めている間取りが何なのかを画像から抽出しようとしている。
- ユーザーの検索行動、嗜好の推定
- ユーザーの検索行動から、嗜好性を掴み取り、レコメンドする際の理由付けができるようにモデルを作ろうとしている。

PDFのリンク

https://kaigi.org/jsai/webprogram/2017/pdf/743.pdf

感想

特徴量の削減のためにランダムフォレストを適用する流れ
高次元スパースデータに対するニューラルネットワークの適用
CVの時系列予測のためにRNNの適用
SNS分析におけるネットワーク特徴量とフィードバック特徴量という概念
口コミデータ分析の作法(点数別に感性表現の係り受けを見る)
など勉強になりました。Kaggleと違って日本語なのが読みやすかったですね。

参考文献

初めてのディープラーニング –オープンソース”Caffe”による演習付き
 SMOTE で不均衡データの分類

Kaggleで使われた特徴量エンジニアリングとアルゴリズムまとめ

Kaggleにおいて、人によってはソリューションを書いてくれているものがあります。特徴量エンジニアリングやアルゴリズムなど業務での参考になるだろうと思われるので、仕事で関わりそうなデータという観点から4つほどですが事例を集めてみました。このような事例のまとまった本とかあったらいいのですが。

基本的に各事例ごとに
・目的
・どんなデータセット
・どんな特徴量エンジニアリング
・どんなアルゴリズム
・リンク
を記していきます。

Walmart Recruiting – Store Sales Forecasting

目的

ウォルマートの部門ごとの売上予測

どんなデータセット

外部のデータ
・気温
・物価指数
・スーパーボウルの時期
など

どんな特徴量エンジニアリング

を行った程度のようです。

どんなアルゴリズム

機械学習の手法
3. Random Forest
4. Linear Regression
5. K nearest regression
6. Principle Component Regression

以上の6つのモデルから平均値をとって予測し、部門ごとにモデルを作ったようです。

リンク

6 bad models make 1 good model: Power of Ensemble Learning

Algorithmic Trading Challenge

目的

大規模なトレーディングにおけるマーケットの反応を予測

どんなデータセット

どんな特徴量エンジニアリング

どんなアルゴリズム

リンク

Winning the Kaggle Algorithmic Trading Challenge with the Composition of Many Models and Feature Engineering

Predict Closed Questions on Stack Overflow

目的

質問が閉じられるかどうかを予測

どんなデータセット

どんな特徴量エンジニアリング

どんなアルゴリズム

Vowpal WabbitというMicrosoftやYahoo!が出資している機械学習ツールを使い、多クラス分類を行ったようです。カテゴリカルな変数をone-hotに変換する必要がなく、楽に分析できるようです。

リンク

Predicting closed questions on Stack Overflow

vowpal_wabbit

Basic Multiclass Classification

Outbrain Click Prediction

目的

ユーザーにレコメンドするコンテンツのクリック予測

どんなデータセット

どんな特徴量エンジニアリング

どんなアルゴリズム

リンク

2nd place solution | team brain-afk

感想

調べてみて、複数の時系列モデルの予測結果の平均値で予測する手法、特徴量の選択をstepAICのようにランダムフォレストで行う手法、Vowpal Wabbit（今回の2つのソリューションで扱われていた）などを新たに知れました。Kaggleには他にもKernelという手法をシェアする場があるので、その情報も今後キャッチアップしていきたいと思います。

参考情報

機械学習コンペティションの進展と今後の展開

SlideShareの機械学習に関するスライドを大量に集めてみた

SlideShareの検索機能に満足できなかったので、Googleのsiteコマンド検索（　「機械学習 site:slideshare.net　」）で集めてき
た、機械学習に関するスライド520件のURLに対して、タイトル・ビュー数・公開日をWebスクレイピングにより取得してみました。

以下の表はビュー数の順番で並んでいます。WordPressのTablepressというプラグインを使っているので、簡単に検索ないし100件表示なども行えます。数の多いものから見ていこうと思います。

(取得用の稚拙なPythonコードは最下部に載っています。)
(URLがリンクになるように改良しました。)

[table id=3 /]

一応、SlideShareのURLのCSVファイル（N行1列のデータ）があれば取得できるPythonコード（2系）を載せておきます。

import urllib, urllib2, csv
from bs4 import BeautifulSoup
from BeautifulSoup import BeautifulSoup
from urllib2 import Request, urlopen, URLError, HTTPError
import sys, string, codecs, time, re
from httplib import BadStatusLine, IncompleteRead

#CSVでリストを読み込む。
f = open('url_list_ss.csv', 'rb')
dataReader = csv.reader(f)

#結果の出力用のリストを作る。
data01 =[]
data02 =[]
data03 =[]
data04 =[]
data05 =[]

for row in dataReader:
       for url in row:
              try:
                     res = urllib2.urlopen(url)
                     soup = BeautifulSoup(res.read(),fromEncoding="utf-8")
                     time.sleep(5.0) #sleep(秒指定)
                     
                     for title in soup.findAll("div",{"class":"small-10 columns"}):
                            for views in soup.findAll("div",{"class":"small-2 columns text-right format-views"}):
                                for times in soup.findAll("time",{"itemprop":"datePublished"}):
                                            data01.append(url)
                                            data02.append(''.join(title.findAll(text=True)))
                                            data03.append(''.join(views.findAll(text=True)))
                                            data04.append(''.join(times.findAll(text=True)))
                                            data05.append("contentslist")
                                            
                                            #データ結合
                                            data = zip(data01,data02,data03,data04,data05)

                                            #CSV出力
                                            f= open('ss_extract_result.csv','w')
                                            writecsv = csv.writer(f,lineterminator='\n')
                                            for x in range(len(data)):
                                                writecsv.writerows([data[x]])
                                                                        

              except HTTPError, e:
                     print e.code
                     
              except BadStatusLine:
                     print "could not fetch"
                     
              except IncompleteRead:
                     print "IncompleteRead"
                     
              except IndexError:
                     print "IndexError"

import urllib, urllib2, csv

from bs4 import BeautifulSoup

from BeautifulSoup import BeautifulSoup

from urllib2 import Request, urlopen, URLError, HTTPError

import sys, string, codecs, time, re

from httplib import BadStatusLine, IncompleteRead

#CSVでリストを読み込む。

f = open('url_list_ss.csv', 'rb')

dataReader = csv.reader(f)

#結果の出力用のリストを作る。

data01 =[]

data02 =[]

data03 =[]

data04 =[]

data05 =[]

for row in dataReader:

for url in row:

try:

res = urllib2.urlopen(url)

soup = BeautifulSoup(res.read(),fromEncoding="utf-8")

time.sleep(5.0) #sleep(秒指定)

for title in soup.findAll("div",{"class":"small-10 columns"}):

for views in soup.findAll("div",{"class":"small-2 columns text-right format-views"}):

for times in soup.findAll("time",{"itemprop":"datePublished"}):

data01.append(url)

data02.append(''.join(title.findAll(text=True)))

data03.append(''.join(views.findAll(text=True)))

data04.append(''.join(times.findAll(text=True)))

data05.append("contentslist")

#データ結合

data = zip(data01,data02,data03,data04,data05)

#CSV出力

f= open('ss_extract_result.csv','w')

writecsv = csv.writer(f,lineterminator='\n')

for x in range(len(data)):

writecsv.writerows([data[x]])

except HTTPError, e:

print e.code

except BadStatusLine:

print "could not fetch"

except IncompleteRead:

print "IncompleteRead"

except IndexError:

print "IndexError"

Japan.R 2016のスライドまとめ

まだ手に入れていないスライドもあるので随時更新しますが、Japan.R 2016(connpass)のスライドをまとめています。後日、登場したパッケージなどのサンプルコードも載せていく予定です。

目次

・石田基広さんのキーノート
・ホクソエムとは何だったのか（ホクソエムさん）
・Rと探索的データ分析で、国連での日本の立ち位置を可視化する（安田洋介さん）
・マウス操作でかんたん予測分析（鈴木了太さん）
・高速・省メモリにlibsvm形式でダンプする方法を研究してみた（@hskksk）
・Rでてんしょくかつどう（@Med_KU）
・RStudio vs Emacs（@y__mattu）
・randomforestで高次元の変数重要度見る（@siero5335）
・Rで本を作りたい（前田和寛さん）
・28歳でプログラミングを始めた話（市川太祐さん）
・LDA-Visパッケージのご紹介（@doradora09）
・【e2d3R】E2D3からDot-Bar-Chartのご紹介（楠本一哲さん）
・このIRのグラフがすごい！上場企業2016（@ito_yan）
・Rでカルマンフィルタをしたい（@tetsuroito）
・PPAP（仮）（@yutannihilation）
・スライド未公開、ユーザーの状態遷移に関する分析のお話（@sanoche16）
・私とR（高栁慎一さん）
・めくってもめくってもサンプル画像（服部恵美さん）
・木と電話と選挙(causalTree)（安井翔太さん）
・スライド未公開、dplyrの話（@tomomoto）
・てかLINEやってる？(仮)（@wonder_zone）
・心理学における「再現性」の問題とBayes Factor（@NSushi）

・石田基広さんのキーノート

スライド未公開です。

・Linux使い
・ヘブライ語の意味構造を代数学でやっていた
・S/R言語の生みの親はJohn Chambers
　以下の二つは最近書かれた本だそうです。
　Software for Data Analysis: Programming with R (Statistics and Computing)
　Extending R (Chapman & Hall/CRC The R Series)
・S→S-plus→Rの順番で発展
・purrrを最近使い始めたそうです。
・XLConnectパッケージを使って、大学教員の採点活動を効率化しているそうです。

・ホクソエムとは何だったのか（ホクソエムさん）

匿名技術者集団ホクソエムの2016年の成果
・densratio（ densratio: Density Ratio Estimation ）
・githubinstall
（ githubinstall: A Helpful Way to Install R Packages Hosted on GitHub ）
・healthplanet（ Wrapper package for healthplanet api ）
・RODBCDBI
（ RODBCDBI: Provides Access to Databases Through the ODBC Interface ）
・jpmesh（ jpmesh: Utilities for Japanese Mesh Code ）

起業されたとのことです。懸命に頑張って下さい！
株式会社ホクソエム

awesomeな人材が必要とのことで、awesomeな方はアプライしてみてはいかがでしょうか。

・Rと探索的データ分析で、国連での日本の立ち位置を可視化する（安田洋介さん）

スライド未公開です。
国連のデータを使って、Exploratoryを用いた探索的データ分析の実演をされていました。

・マウス操作でかんたん予測分析（鈴木了太さん）

R AnalyticFlow
Rで実践！データサイエンス～初めの一歩から高度な応用まで～

・高速・省メモリにlibsvm形式でダンプする方法を研究してみた（@hskksk）

高速・省メモリにlibsvm形式でダンプする方法を研究してみた from Keisuke Hosaka

・Rでてんしょくかつどう（@Med_KU）

Rmd でreveal.js のhtml スライドプレゼンテーション

・RStudio vs Emacs（@y__mattu）

RStudio vs Emacs Japan.R 2016

・randomforestで高次元の変数重要度見る（@siero5335）

Randomforestで高次元の変数重要度を見る #japanr LT from Akifumi Eguchi

・Rで本を作りたい（前田和寛さん）

Rで本を作りたい

・28歳でプログラミングを始めた話（市川太祐さん）

・医療関連のアプリ開発でデータサイエンスを駆使しようとしているそうです。

スライド未公開です。
スライドがシェアされ次第載せます。

・LDA-Visパッケージのご紹介（@doradora09）

20161127 doradora09 japanr2016_lt from Nobuaki Oshiro

・【e2d3R】E2D3からDot-Bar-Chartのご紹介（楠本一哲さん）

スライドは未公開です。
E2D3をRで表示する試みのようです。
Experiments with e2d3 in R

・このIRのグラフがすごい！上場企業2016（@ito_yan）

スライド未公開です。後日シェアしていただけるようです。

今年のプレゼンはネタを集めるにしても、これが限界でした…楽しんで頂けたでしょうか？スライドは後日アップします。 #JapanR

— いとーちゃん (@ito_yan) 2016年11月27日

・Rでカルマンフィルタをしたい（@tetsuroito）

・PPAP（仮）（@yutannihilation）

Tidyverseとは from yutannihilation

・スライド未公開、ユーザーの状態遷移に関する分析のお話（@sanoche16）

スライドがシェアされ次第載せます。

・私とR（高栁慎一さん）

RjpWiki
統計・データ解析
 統計解析フリーソフト R の備忘録頁 ver.3.1
seekR(R限定の検索エンジン)
からだにいいもの
アブラタニブログってなんでしょう。油谷さんのブログ？

・めくってもめくってもサンプル画像（服部恵美さん）

Rのサンプルコードはあるけれども、どんな図ができるのかはわからない。そこで、サンプルコードとグラフを大量にまとめているサイトを作ったそうです。検索性は未知数ですが、暇なときに眺めておきたいですね。
R Graphical Manual

・木と電話と選挙(causalTree)（安井翔太さん）

木と電話と選挙(causalTree) from Shota Yasui

・スライド未公開、dplyrの話（@tomomoto）

スライドがシェアされ次第載せます。

・てかLINEやってる？(仮)（@wonder_zone）

てかLINEやってる？ (Japan.R 2016 LT) #JapanR from cancolle

・心理学における「再現性」の問題とBayes Factor（@NSushi）

スライドは後日公開とのことです。

今日はjapan.Rお疲れ様でしたー！楽しかったですありがとうございましたー！
ネット環境がないので心理学の再現性およびベイズファクターおよびにこにースライドは明日アップしますー。

— numba shushi (@NSushi) 2016年11月27日

scikit-learnのモジュールのGitHubでの利用頻度を調べてみた

『Python機械学習プログラミング』を読んで、scikit-learnのモジュールは充実しているなぁと感じたんですが、実際にWebサイトでUser Guide( http://scikit-learn.org/stable/user_guide.html )を見た所、この本に載り切らないような数多くの機械学習手法に応じたモジュールが用意されていました。そこで、世のデータサイエンティストはどのモジュールを良く使っているのだろうと気になったので、GitHubのSearchでヒットしたCodeの数を各モジュール単位で集計してみました。検索クエリは「scikit-learn + モジュール名」なので、正確なものではないのですが、相対的な利用頻度を見るぶんには使えるのではないかと思われます。

データ集計方法

・User Guideに登場するscikit-learnのモジュール名を集めています。
・教師付き学習か教師無し学習かどうかの判断は、User Guideで紹介されているモジュールかどうかで判断しています。
・GitHubのSearchで「scikit-learn + モジュール名」でヒットした件数をそのまま使っています。(2016年9月22日時点)

可視化コード

Jupyterで実行しています。

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import pylab

module_df = pd.read_csv("scikit_learn_modulelist.csv")

pylab.figure(figsize=(5, 25))
sns.set_style("whitegrid")
ax = sns.barplot(x="Hitscode", y="Module", data=module_df)

%matplotlib inline

import matplotlib.pyplot as plt

import seaborn as sns

import pandas as pd

import pylab

module_df = pd.read_csv("scikit_learn_modulelist.csv")

pylab.figure(figsize=(5, 25))

sns.set_style("whitegrid")

ax = sns.barplot(x="Hitscode", y="Module", data=module_df)

教師付き学習編

圧倒的に多いのがSVM(Support Vector Machine)を扱っているSVCモジュールで、続いて定番のロジスティック回帰やRandom Forestが使われているようです。統計解析ではメジャーなはずの線形回帰が5位なのは、初歩的なのであまりコードがアップされていないのかもしれません。GBDTのモジュールももう少し上位にくるかと思ったんですが15位でした。DMLCのXGBoostモジュールを使っているのかもしれませんね。私も実際のところXGBoostを使ってますし。

教師無し学習編

主成分分析やK-mean法など因子分解などのモジュールが上位を占めています。LDA(Latent Dirichlet Allocation)がもっと上位に来ると思ったんですが、思えばGensimの方が充実しているなぁと思うので、このランキングは妥当なのかもしれません。私もLDAなどはGensimを使っていますし。

収集を終えて

・社内だとデータサイエンティストの方がいないので、scikit-learnのモジュールの利用状況を知れてマニアックな共感をすることができた。
・SVMは実践例が豊富そうなので分析事例を探せば良い発見があるかもしれない。
・scikit-learnのUser Guideは充実していたので、時間を作って向き合ってみたいと思った。

XGBoostのパラメータチューニング実践 with Python

以前の投稿で紹介したXGBoostのパラメータチューニング方法ですが、実際のデータセットに対して実行するためのプログラムを実践してみようと思います。プログラム自体はAnalytics_Vidhya/Articles/Parameter_Tuning_XGBoost_with_Example/XGBoost models.ipynbに載っているのですが、データセットがついていません。そこで、前回の投稿(不均衡なデータの分類問題について with Python)で赤ワインのデータセットを手に入れているので、こちらのデータセットを用います。誤植なのかところどころ、うまく回らなかったところがあったので、手直しをしています。

以下の工程に従って進みます。結構長いですが、辛抱強く実践してみて下さい。
・ライブラリの読み込み
・データの読み込み
・前処理
・学習用データとテスト用データの作成
・XGBoostの予測結果をもとに、AUCの数値を返すための関数の定義
・モデルの実行
・チューニング

ライブラリの読み込み

import pandas as pd
import numpy as np
import xgboost as xgb
from xgboost.sklearn import XGBClassifier
from sklearn import cross_validation, metrics
from sklearn.grid_search import GridSearchCV

import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 12, 4

import pandas as pd

import numpy as np

import xgboost as xgb

from xgboost.sklearn import XGBClassifier

from sklearn import cross_validation, metrics

from sklearn.grid_search import GridSearchCV

import matplotlib.pylab as plt

%matplotlib inline

from matplotlib.pylab import rcParams

rcParams['figure.figsize'] = 12, 4

データの読み込み

#importing the red wine data
wine_df = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", sep=";")

1 2	#importing the red wine data wine_df = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", sep=";")

前処理

#ユニークIDを行ごとに割り当てる。
wine_df['ID'] = range(1, len(wine_df) + 1)

#ワインの質に関するデータを0-1データに置換する。qualityが6よりも小さかったら0、それ以外は1とする。
Y = wine_df.quality.values
wine_df.quality = np.asarray([1 if  i>=6 else 0 for i in Y])
wine_df.head(10)

#ユニークIDを行ごとに割り当てる。

wine_df['ID'] = range(1, len(wine_df) + 1)

#ワインの質に関するデータを0-1データに置換する。qualityが6よりも小さかったら0、それ以外は1とする。

Y = wine_df.quality.values

wine_df.quality = np.asarray([1 if i>=6 else 0 for i in Y])

wine_df.head(10)

学習用データとテスト用データの作成

#学習用データとテスト用データの作成
msk = np.random.rand(len(wine_df)) < 0.8 #乱数を発生させて0.8よりも小さいデータを選ぶ
train = wine_df[msk]
test = wine_df[~msk]

train.shape, test.shape
((1236, 13), (363, 13))

target='quality'
IDcol = 'ID'

#訓練データの目的変数の確認
train[target].value_counts()
1    659
0    577
Name: quality, dtype: int64

#学習用データとテスト用データの作成

msk = np.random.rand(len(wine_df)) < 0.8 #乱数を発生させて0.8よりも小さいデータを選ぶ

train = wine_df[msk]

test = wine_df[~msk]

train.shape, test.shape

((1236, 13), (363, 13))

target='quality'

IDcol = 'ID'

#訓練データの目的変数の確認

train[target].value_counts()

1 659

0 577

Name: quality, dtype: int64

XGBoostの予測結果をもとに、AUCの数値を返すための関数の定義

XGBoostの予測結果から、AUCの数値を返し、特徴量に応じた重要度を出力するためのプログラムです。

#テスト結果を格納するデータフレームの生成
test_results = pd.DataFrame(data=test.ID)

#関数の定義
def modelfit(alg, dtrain, dtest, predictors,useTrainCV=True, cv_folds=5, early_stopping_rounds=50):
    
    if useTrainCV:
        xgb_param = alg.get_xgb_params()
        xgtrain = xgb.DMatrix(dtrain[predictors].values, label=dtrain[target].values)
        xgtest = xgb.DMatrix(dtest[predictors].values)
        cvresult = xgb.cv(xgb_param, xgtrain, num_boost_round=alg.get_params()['n_estimators'], nfold=cv_folds,
            metrics=['auc'], early_stopping_rounds=early_stopping_rounds, show_progress=False)
        alg.set_params(n_estimators=cvresult.shape[0])
    
    #Fit the algorithm on the data
    alg.fit(dtrain[predictors], dtrain[target],eval_metric=['auc'])
        
    #Predict training set:
    dtrain_predictions = alg.predict(dtrain[predictors])
    dtrain_predprob = alg.predict_proba(dtrain[predictors])[:,1]
        
    #Print model report:
    print "\nModel Report"
    print "Accuracy : %.4g" % metrics.accuracy_score(dtrain[target].values, dtrain_predictions)
    print "AUC Score (Train): %f" % metrics.roc_auc_score(dtrain[target], dtrain_predprob)
    
    # Predict on testing data:
    dtest['predprob'] = alg.predict_proba(dtest[predictors])[:,1]
    #results = test_results.merge(dtest[['ID','predprob']], on='ID')
    print 'AUC Score (Test): %f' % metrics.roc_auc_score(dtest[target], dtest['predprob'])
                
    feat_imp = pd.Series(alg.booster().get_fscore()).sort_values(ascending=False)
    feat_imp.plot(kind='bar', title='Feature Importances')
    plt.ylabel('Feature Importance Score')

#テスト結果を格納するデータフレームの生成

test_results = pd.DataFrame(data=test.ID)

#関数の定義

def modelfit(alg, dtrain, dtest, predictors,useTrainCV=True, cv_folds=5, early_stopping_rounds=50):

if useTrainCV:

xgb_param = alg.get_xgb_params()

xgtrain = xgb.DMatrix(dtrain[predictors].values, label=dtrain[target].values)

xgtest = xgb.DMatrix(dtest[predictors].values)

cvresult = xgb.cv(xgb_param, xgtrain, num_boost_round=alg.get_params()['n_estimators'], nfold=cv_folds,

metrics=['auc'], early_stopping_rounds=early_stopping_rounds, show_progress=False)

alg.set_params(n_estimators=cvresult.shape[0])

#Fit the algorithm on the data

alg.fit(dtrain[predictors], dtrain[target],eval_metric=['auc'])

#Predict training set:

dtrain_predictions = alg.predict(dtrain[predictors])

dtrain_predprob = alg.predict_proba(dtrain[predictors])[:,1]

#Print model report:

print "\nModel Report"

print "Accuracy : %.4g" % metrics.accuracy_score(dtrain[target].values, dtrain_predictions)

print "AUC Score (Train): %f" % metrics.roc_auc_score(dtrain[target], dtrain_predprob)

# Predict on testing data:

dtest['predprob'] = alg.predict_proba(dtest[predictors])[:,1]

#results = test_results.merge(dtest[['ID','predprob']], on='ID')

print 'AUC Score (Test): %f' % metrics.roc_auc_score(dtest[target], dtest['predprob'])

feat_imp = pd.Series(alg.booster().get_fscore()).sort_values(ascending=False)

feat_imp.plot(kind='bar', title='Feature Importances')

plt.ylabel('Feature Importance Score')

モデルの実行

predictors = [x for x in train.columns if x not in [target, IDcol]]
xgb1 = XGBClassifier(
        learning_rate =0.1,
        n_estimators=1000,
        max_depth=5,
        min_child_weight=1,
        gamma=0,
        subsample=0.8,
        colsample_bytree=0.8,
        objective= 'binary:logistic',
        nthread=4,
        scale_pos_weight=1,
        seed=27)

modelfit(xgb1, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.
Stopping. Best iteration: 237

Model Report
Accuracy : 1
AUC Score (Train): 1.000000
AUC Score (Test): 0.875199

predictors = [x for x in train.columns if x not in [target, IDcol]]

xgb1 = XGBClassifier(

learning_rate =0.1,

n_estimators=1000,

max_depth=5,

min_child_weight=1,

gamma=0,

subsample=0.8,

colsample_bytree=0.8,

objective= 'binary:logistic',

nthread=4,

scale_pos_weight=1,

seed=27)

modelfit(xgb1, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.

Stopping. Best iteration: 237

Model Report

Accuracy : 1

AUC Score (Train): 1.000000

AUC Score (Test): 0.875199

チューニング

max_depthとmin_child_weightの数値をチューニングするためのプログラムです。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test1 = {
    'max_depth':range(3,10,2),
    'min_child_weight':range(1,6,2)
}
gsearch1 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=1000, max_depth=5,
                                        min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1, seed=27), 
                       param_grid = param_test1, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch1.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test1 = {

'max_depth':range(3,10,2),

'min_child_weight':range(1,6,2)

}

gsearch1 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=1000, max_depth=5,

min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1, seed=27),

param_grid = param_test1, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch1.fit(train[predictors],train[target])

gsearch1.grid_scores_, gsearch1.best_params_, gsearch1.best_score_

([mean: 0.76728, std: 0.03045, params: {'max_depth': 3, 'min_child_weight': 1},
  mean: 0.76649, std: 0.03378, params: {'max_depth': 3, 'min_child_weight': 3},
  mean: 0.76540, std: 0.03620, params: {'max_depth': 3, 'min_child_weight': 5},
  mean: 0.76509, std: 0.03183, params: {'max_depth': 5, 'min_child_weight': 1},
  mean: 0.76430, std: 0.02988, params: {'max_depth': 5, 'min_child_weight': 3},
  mean: 0.76221, std: 0.03336, params: {'max_depth': 5, 'min_child_weight': 5},
  mean: 0.77162, std: 0.03335, params: {'max_depth': 7, 'min_child_weight': 1},
  mean: 0.76575, std: 0.03585, params: {'max_depth': 7, 'min_child_weight': 3},
  mean: 0.76277, std: 0.03511, params: {'max_depth': 7, 'min_child_weight': 5},
  mean: 0.77235, std: 0.03283, params: {'max_depth': 9, 'min_child_weight': 1},
  mean: 0.76452, std: 0.03414, params: {'max_depth': 9, 'min_child_weight': 3},
  mean: 0.76114, std: 0.03561, params: {'max_depth': 9, 'min_child_weight': 5}],
 {'max_depth': 9, 'min_child_weight': 1},
 0.77235073909956886)

gsearch1.grid_scores_, gsearch1.best_params_, gsearch1.best_score_

([mean: 0.76728, std: 0.03045, params: {'max_depth': 3, 'min_child_weight': 1},

mean: 0.76649, std: 0.03378, params: {'max_depth': 3, 'min_child_weight': 3},

mean: 0.76540, std: 0.03620, params: {'max_depth': 3, 'min_child_weight': 5},

mean: 0.76509, std: 0.03183, params: {'max_depth': 5, 'min_child_weight': 1},

mean: 0.76430, std: 0.02988, params: {'max_depth': 5, 'min_child_weight': 3},

mean: 0.76221, std: 0.03336, params: {'max_depth': 5, 'min_child_weight': 5},

mean: 0.77162, std: 0.03335, params: {'max_depth': 7, 'min_child_weight': 1},

mean: 0.76575, std: 0.03585, params: {'max_depth': 7, 'min_child_weight': 3},

mean: 0.76277, std: 0.03511, params: {'max_depth': 7, 'min_child_weight': 5},

mean: 0.77235, std: 0.03283, params: {'max_depth': 9, 'min_child_weight': 1},

mean: 0.76452, std: 0.03414, params: {'max_depth': 9, 'min_child_weight': 3},

mean: 0.76114, std: 0.03561, params: {'max_depth': 9, 'min_child_weight': 5}],

{'max_depth': 9, 'min_child_weight': 1},

0.77235073909956886)

より細かい数値で再度最適なパラメータを探します。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test2 = {
    'max_depth':[4,5,6,7,8,9],
    'min_child_weight':[1,2,3,4,5,6]
}
gsearch2 = GridSearchCV(estimator = XGBClassifier( learning_rate=0.1, n_estimators=1000, max_depth=5,
                                        min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test2, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch2.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test2 = {

'max_depth':[4,5,6,7,8,9],

'min_child_weight':[1,2,3,4,5,6]

}

gsearch2 = GridSearchCV(estimator = XGBClassifier( learning_rate=0.1, n_estimators=1000, max_depth=5,

min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test2, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch2.fit(train[predictors],train[target])

gsearch2.grid_scores_, gsearch2.best_params_, gsearch2.best_score_

([mean: 0.76820, std: 0.03336, params: {'max_depth': 4, 'min_child_weight': 1},
  mean: 0.76431, std: 0.02792, params: {'max_depth': 4, 'min_child_weight': 2},
  mean: 0.76171, std: 0.03356, params: {'max_depth': 4, 'min_child_weight': 3},
  mean: 0.76257, std: 0.03277, params: {'max_depth': 4, 'min_child_weight': 4},
  mean: 0.76128, std: 0.03661, params: {'max_depth': 4, 'min_child_weight': 5},
  mean: 0.75902, std: 0.03280, params: {'max_depth': 4, 'min_child_weight': 6},
  mean: 0.76509, std: 0.03183, params: {'max_depth': 5, 'min_child_weight': 1},
  mean: 0.76426, std: 0.02974, params: {'max_depth': 5, 'min_child_weight': 2},
  mean: 0.76430, std: 0.02988, params: {'max_depth': 5, 'min_child_weight': 3},
  mean: 0.76262, std: 0.02992, params: {'max_depth': 5, 'min_child_weight': 4},
  mean: 0.76221, std: 0.03336, params: {'max_depth': 5, 'min_child_weight': 5},
  mean: 0.76655, std: 0.03397, params: {'max_depth': 5, 'min_child_weight': 6},
  mean: 0.77066, std: 0.02936, params: {'max_depth': 6, 'min_child_weight': 1},
  mean: 0.76422, std: 0.03038, params: {'max_depth': 6, 'min_child_weight': 2},
  mean: 0.76126, std: 0.03021, params: {'max_depth': 6, 'min_child_weight': 3},
  mean: 0.76334, std: 0.03176, params: {'max_depth': 6, 'min_child_weight': 4},
  mean: 0.76347, std: 0.03245, params: {'max_depth': 6, 'min_child_weight': 5},
  mean: 0.76437, std: 0.03546, params: {'max_depth': 6, 'min_child_weight': 6},
  mean: 0.77162, std: 0.03335, params: {'max_depth': 7, 'min_child_weight': 1},
  mean: 0.76140, std: 0.03245, params: {'max_depth': 7, 'min_child_weight': 2},
  mean: 0.76575, std: 0.03585, params: {'max_depth': 7, 'min_child_weight': 3},
  mean: 0.76345, std: 0.03518, params: {'max_depth': 7, 'min_child_weight': 4},
  mean: 0.76277, std: 0.03511, params: {'max_depth': 7, 'min_child_weight': 5},
  mean: 0.75858, std: 0.03375, params: {'max_depth': 7, 'min_child_weight': 6},
  mean: 0.77487, std: 0.03377, params: {'max_depth': 8, 'min_child_weight': 1},
  mean: 0.76740, std: 0.03472, params: {'max_depth': 8, 'min_child_weight': 2},
  mean: 0.76048, std: 0.03267, params: {'max_depth': 8, 'min_child_weight': 3},
  mean: 0.76288, std: 0.03773, params: {'max_depth': 8, 'min_child_weight': 4},
  mean: 0.76138, std: 0.03045, params: {'max_depth': 8, 'min_child_weight': 5},
  mean: 0.76233, std: 0.03652, params: {'max_depth': 8, 'min_child_weight': 6},
  mean: 0.77235, std: 0.03283, params: {'max_depth': 9, 'min_child_weight': 1},
  mean: 0.76929, std: 0.03267, params: {'max_depth': 9, 'min_child_weight': 2},
  mean: 0.76452, std: 0.03414, params: {'max_depth': 9, 'min_child_weight': 3},
  mean: 0.76152, std: 0.03731, params: {'max_depth': 9, 'min_child_weight': 4},
  mean: 0.76114, std: 0.03561, params: {'max_depth': 9, 'min_child_weight': 5},
  mean: 0.76551, std: 0.03394, params: {'max_depth': 9, 'min_child_weight': 6}],
 {'max_depth': 8, 'min_child_weight': 1},
 0.77486987248915451)

gsearch2.grid_scores_, gsearch2.best_params_, gsearch2.best_score_

([mean: 0.76820, std: 0.03336, params: {'max_depth': 4, 'min_child_weight': 1},

mean: 0.76431, std: 0.02792, params: {'max_depth': 4, 'min_child_weight': 2},

mean: 0.76171, std: 0.03356, params: {'max_depth': 4, 'min_child_weight': 3},

mean: 0.76257, std: 0.03277, params: {'max_depth': 4, 'min_child_weight': 4},

mean: 0.76128, std: 0.03661, params: {'max_depth': 4, 'min_child_weight': 5},

mean: 0.75902, std: 0.03280, params: {'max_depth': 4, 'min_child_weight': 6},

mean: 0.76509, std: 0.03183, params: {'max_depth': 5, 'min_child_weight': 1},

mean: 0.76426, std: 0.02974, params: {'max_depth': 5, 'min_child_weight': 2},

mean: 0.76430, std: 0.02988, params: {'max_depth': 5, 'min_child_weight': 3},

mean: 0.76262, std: 0.02992, params: {'max_depth': 5, 'min_child_weight': 4},

mean: 0.76221, std: 0.03336, params: {'max_depth': 5, 'min_child_weight': 5},

mean: 0.76655, std: 0.03397, params: {'max_depth': 5, 'min_child_weight': 6},

mean: 0.77066, std: 0.02936, params: {'max_depth': 6, 'min_child_weight': 1},

mean: 0.76422, std: 0.03038, params: {'max_depth': 6, 'min_child_weight': 2},

mean: 0.76126, std: 0.03021, params: {'max_depth': 6, 'min_child_weight': 3},

mean: 0.76334, std: 0.03176, params: {'max_depth': 6, 'min_child_weight': 4},

mean: 0.76347, std: 0.03245, params: {'max_depth': 6, 'min_child_weight': 5},

mean: 0.76437, std: 0.03546, params: {'max_depth': 6, 'min_child_weight': 6},

mean: 0.77162, std: 0.03335, params: {'max_depth': 7, 'min_child_weight': 1},

mean: 0.76140, std: 0.03245, params: {'max_depth': 7, 'min_child_weight': 2},

mean: 0.76575, std: 0.03585, params: {'max_depth': 7, 'min_child_weight': 3},

mean: 0.76345, std: 0.03518, params: {'max_depth': 7, 'min_child_weight': 4},

mean: 0.76277, std: 0.03511, params: {'max_depth': 7, 'min_child_weight': 5},

mean: 0.75858, std: 0.03375, params: {'max_depth': 7, 'min_child_weight': 6},

mean: 0.77487, std: 0.03377, params: {'max_depth': 8, 'min_child_weight': 1},

mean: 0.76740, std: 0.03472, params: {'max_depth': 8, 'min_child_weight': 2},

mean: 0.76048, std: 0.03267, params: {'max_depth': 8, 'min_child_weight': 3},

mean: 0.76288, std: 0.03773, params: {'max_depth': 8, 'min_child_weight': 4},

mean: 0.76138, std: 0.03045, params: {'max_depth': 8, 'min_child_weight': 5},

mean: 0.76233, std: 0.03652, params: {'max_depth': 8, 'min_child_weight': 6},

mean: 0.77235, std: 0.03283, params: {'max_depth': 9, 'min_child_weight': 1},

mean: 0.76929, std: 0.03267, params: {'max_depth': 9, 'min_child_weight': 2},

mean: 0.76452, std: 0.03414, params: {'max_depth': 9, 'min_child_weight': 3},

mean: 0.76152, std: 0.03731, params: {'max_depth': 9, 'min_child_weight': 4},

mean: 0.76114, std: 0.03561, params: {'max_depth': 9, 'min_child_weight': 5},

mean: 0.76551, std: 0.03394, params: {'max_depth': 9, 'min_child_weight': 6}],

{'max_depth': 8, 'min_child_weight': 1},

0.77486987248915451)

max_depthを8、min_child_weightを1として、他のパラメータチューニングに移ります。
続いて、gammaのチューニングを行います。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test3 = {
    'gamma':[i/10.0 for i in range(0,5)]
}
gsearch3 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=1000, max_depth=8,
                                        min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test3, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch3.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test3 = {

'gamma':[i/10.0 for i in range(0,5)]

}

gsearch3 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=1000, max_depth=8,

min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test3, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch3.fit(train[predictors],train[target])

gsearch3.grid_scores_, gsearch3.best_params_, gsearch3.best_score_

([mean: 0.77487, std: 0.03377, params: {'gamma': 0.0},
  mean: 0.77689, std: 0.03298, params: {'gamma': 0.1},
  mean: 0.77735, std: 0.03117, params: {'gamma': 0.2},
  mean: 0.78163, std: 0.03076, params: {'gamma': 0.3},
  mean: 0.78790, std: 0.03328, params: {'gamma': 0.4}],
 {'gamma': 0.4},
 0.78789976715320331)

gsearch3.grid_scores_, gsearch3.best_params_, gsearch3.best_score_

([mean: 0.77487, std: 0.03377, params: {'gamma': 0.0},

mean: 0.77689, std: 0.03298, params: {'gamma': 0.1},

mean: 0.77735, std: 0.03117, params: {'gamma': 0.2},

mean: 0.78163, std: 0.03076, params: {'gamma': 0.3},

mean: 0.78790, std: 0.03328, params: {'gamma': 0.4}],

{'gamma': 0.4},

0.78789976715320331)

gammaを0.4と置きます。
ここで、いままでにチューニングしたパラメータを用いて再度推定を行います。先ほどの0.875よりも高くなっています。

predictors = [x for x in train.columns if x not in [target, IDcol]]
xgb2 = XGBClassifier(
        learning_rate =0.1,
        n_estimators=1000,
        max_depth=8,
        min_child_weight=1,
        gamma=0.4,
        subsample=0.8,
        colsample_bytree=0.8,
        objective= 'binary:logistic',
        nthread=4,
        scale_pos_weight=1,
        seed=27)
modelfit(xgb2, train, test, predictors)


Will train until cv error hasn't decreased in 50 rounds.
Stopping. Best iteration: 120

Model Report
Accuracy : 1
AUC Score (Train): 1.000000
AUC Score (Test): 0.884028

predictors = [x for x in train.columns if x not in [target, IDcol]]

xgb2 = XGBClassifier(

learning_rate =0.1,

n_estimators=1000,

max_depth=8,

min_child_weight=1,

gamma=0.4,

subsample=0.8,

colsample_bytree=0.8,

objective= 'binary:logistic',

nthread=4,

scale_pos_weight=1,

seed=27)

modelfit(xgb2, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.

Stopping. Best iteration: 120

Model Report

Accuracy : 1

AUC Score (Train): 1.000000

AUC Score (Test): 0.884028

続いて、subsampleとcolsample_bytreeのチューニングを行います。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test4 = {
    'subsample':[i/10.0 for i in range(6,10)],
    'colsample_bytree':[i/10.0 for i in range(6,10)]
}
gsearch4 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,
                                        min_child_weight=1, gamma=0.4, subsample=0.8, colsample_bytree=0.8,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test4, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch4.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test4 = {

'subsample':[i/10.0 for i in range(6,10)],

'colsample_bytree':[i/10.0 for i in range(6,10)]

}

gsearch4 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,

min_child_weight=1, gamma=0.4, subsample=0.8, colsample_bytree=0.8,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test4, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch4.fit(train[predictors],train[target])

gsearch4.grid_scores_, gsearch4.best_params_, gsearch4.best_score_

([mean: 0.78994, std: 0.02779, params: {'subsample': 0.6, 'colsample_bytree': 0.6},
  mean: 0.78900, std: 0.03519, params: {'subsample': 0.7, 'colsample_bytree': 0.6},
  mean: 0.78509, std: 0.03202, params: {'subsample': 0.8, 'colsample_bytree': 0.6},
  mean: 0.78706, std: 0.02848, params: {'subsample': 0.9, 'colsample_bytree': 0.6},
  mean: 0.78511, std: 0.03140, params: {'subsample': 0.6, 'colsample_bytree': 0.7},
  mean: 0.78343, std: 0.03336, params: {'subsample': 0.7, 'colsample_bytree': 0.7},
  mean: 0.78939, std: 0.03203, params: {'subsample': 0.8, 'colsample_bytree': 0.7},
  mean: 0.78646, std: 0.04090, params: {'subsample': 0.9, 'colsample_bytree': 0.7},
  mean: 0.77809, std: 0.03452, params: {'subsample': 0.6, 'colsample_bytree': 0.8},
  mean: 0.78994, std: 0.03483, params: {'subsample': 0.7, 'colsample_bytree': 0.8},
  mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.8},
  mean: 0.79207, std: 0.03057, params: {'subsample': 0.9, 'colsample_bytree': 0.8},
  mean: 0.78466, std: 0.02672, params: {'subsample': 0.6, 'colsample_bytree': 0.9},
  mean: 0.78863, std: 0.03289, params: {'subsample': 0.7, 'colsample_bytree': 0.9},
  mean: 0.78905, std: 0.02660, params: {'subsample': 0.8, 'colsample_bytree': 0.9},
  mean: 0.78501, std: 0.03666, params: {'subsample': 0.9, 'colsample_bytree': 0.9}],
 {'colsample_bytree': 0.8, 'subsample': 0.8},
 0.79369231068019075)

gsearch4.grid_scores_, gsearch4.best_params_, gsearch4.best_score_

([mean: 0.78994, std: 0.02779, params: {'subsample': 0.6, 'colsample_bytree': 0.6},

mean: 0.78900, std: 0.03519, params: {'subsample': 0.7, 'colsample_bytree': 0.6},

mean: 0.78509, std: 0.03202, params: {'subsample': 0.8, 'colsample_bytree': 0.6},

mean: 0.78706, std: 0.02848, params: {'subsample': 0.9, 'colsample_bytree': 0.6},

mean: 0.78511, std: 0.03140, params: {'subsample': 0.6, 'colsample_bytree': 0.7},

mean: 0.78343, std: 0.03336, params: {'subsample': 0.7, 'colsample_bytree': 0.7},

mean: 0.78939, std: 0.03203, params: {'subsample': 0.8, 'colsample_bytree': 0.7},

mean: 0.78646, std: 0.04090, params: {'subsample': 0.9, 'colsample_bytree': 0.7},

mean: 0.77809, std: 0.03452, params: {'subsample': 0.6, 'colsample_bytree': 0.8},

mean: 0.78994, std: 0.03483, params: {'subsample': 0.7, 'colsample_bytree': 0.8},

mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.8},

mean: 0.79207, std: 0.03057, params: {'subsample': 0.9, 'colsample_bytree': 0.8},

mean: 0.78466, std: 0.02672, params: {'subsample': 0.6, 'colsample_bytree': 0.9},

mean: 0.78863, std: 0.03289, params: {'subsample': 0.7, 'colsample_bytree': 0.9},

mean: 0.78905, std: 0.02660, params: {'subsample': 0.8, 'colsample_bytree': 0.9},

mean: 0.78501, std: 0.03666, params: {'subsample': 0.9, 'colsample_bytree': 0.9}],

{'colsample_bytree': 0.8, 'subsample': 0.8},

0.79369231068019075)

より細かい範囲で再度パラメータをチューニングします。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test5 = {
    'subsample':[i/100.0 for i in range(75,90,5)],
    'colsample_bytree':[i/100.0 for i in range(75,90,5)]
}
gsearch5 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,
                                        min_child_weight=1, gamma=0.4, subsample=0.8, colsample_bytree=0.8,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test5, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch5.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test5 = {

'subsample':[i/100.0 for i in range(75,90,5)],

'colsample_bytree':[i/100.0 for i in range(75,90,5)]

}

gsearch5 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,

min_child_weight=1, gamma=0.4, subsample=0.8, colsample_bytree=0.8,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test5, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch5.fit(train[predictors],train[target])

gsearch5.grid_scores_, gsearch5.best_params_, gsearch5.best_score_

([mean: 0.78890, std: 0.03171, params: {'subsample': 0.75, 'colsample_bytree': 0.75},
  mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.75},
  mean: 0.79374, std: 0.03061, params: {'subsample': 0.85, 'colsample_bytree': 0.75},
  mean: 0.78890, std: 0.03171, params: {'subsample': 0.75, 'colsample_bytree': 0.8},
  mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.8},
  mean: 0.79374, std: 0.03061, params: {'subsample': 0.85, 'colsample_bytree': 0.8},
  mean: 0.78418, std: 0.03232, params: {'subsample': 0.75, 'colsample_bytree': 0.85},
  mean: 0.78905, std: 0.02660, params: {'subsample': 0.8, 'colsample_bytree': 0.85},
  mean: 0.78367, std: 0.03582, params: {'subsample': 0.85, 'colsample_bytree': 0.85}],
 {'colsample_bytree': 0.75, 'subsample': 0.85},
 0.79374219292158221)

gsearch5.grid_scores_, gsearch5.best_params_, gsearch5.best_score_

([mean: 0.78890, std: 0.03171, params: {'subsample': 0.75, 'colsample_bytree': 0.75},

mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.75},

mean: 0.79374, std: 0.03061, params: {'subsample': 0.85, 'colsample_bytree': 0.75},

mean: 0.78890, std: 0.03171, params: {'subsample': 0.75, 'colsample_bytree': 0.8},

mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.8},

mean: 0.79374, std: 0.03061, params: {'subsample': 0.85, 'colsample_bytree': 0.8},

mean: 0.78418, std: 0.03232, params: {'subsample': 0.75, 'colsample_bytree': 0.85},

mean: 0.78905, std: 0.02660, params: {'subsample': 0.8, 'colsample_bytree': 0.85},

mean: 0.78367, std: 0.03582, params: {'subsample': 0.85, 'colsample_bytree': 0.85}],

{'colsample_bytree': 0.75, 'subsample': 0.85},

0.79374219292158221)

続いて、reg_alphaをチューニングします。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test6 = {
    'reg_alpha':[1e-5, 1e-2, 0.1, 1, 100]
}
gsearch6 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,
                                        min_child_weight=1, gamma=0.4, subsample=0.85, colsample_bytree=0.75,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test6, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch6.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test6 = {

'reg_alpha':[1e-5, 1e-2, 0.1, 1, 100]

}

gsearch6 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,

min_child_weight=1, gamma=0.4, subsample=0.85, colsample_bytree=0.75,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test6, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch6.fit(train[predictors],train[target])

gsearch6.grid_scores_, gsearch6.best_params_, gsearch6.best_score_

([mean: 0.79377, std: 0.03058, params: {'reg_alpha': 1e-05},
  mean: 0.79068, std: 0.02953, params: {'reg_alpha': 0.01},
  mean: 0.79298, std: 0.03268, params: {'reg_alpha': 0.1},
  mean: 0.78731, std: 0.03270, params: {'reg_alpha': 1},
  mean: 0.72370, std: 0.03333, params: {'reg_alpha': 100}],
 {'reg_alpha': 1e-05},
 0.79376831622356758)

gsearch6.grid_scores_, gsearch6.best_params_, gsearch6.best_score_

([mean: 0.79377, std: 0.03058, params: {'reg_alpha': 1e-05},

mean: 0.79068, std: 0.02953, params: {'reg_alpha': 0.01},

mean: 0.79298, std: 0.03268, params: {'reg_alpha': 0.1},

mean: 0.78731, std: 0.03270, params: {'reg_alpha': 1},

mean: 0.72370, std: 0.03333, params: {'reg_alpha': 100}],

{'reg_alpha': 1e-05},

0.79376831622356758)

範囲が粗かったので、より細かくパラメータをチューニングします。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test7 = {
    'reg_alpha':[0, 0.001, 0.005, 0.01, 0.05]
}
gsearch7 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,
                                        min_child_weight=1, gamma=0.4, subsample=0.85, colsample_bytree=0.75,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test7, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch7.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test7 = {

'reg_alpha':[0, 0.001, 0.005, 0.01, 0.05]

}

gsearch7 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,

min_child_weight=1, gamma=0.4, subsample=0.85, colsample_bytree=0.75,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test7, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch7.fit(train[predictors],train[target])

gsearch7.grid_scores_, gsearch7.best_params_, gsearch7.best_score_

([mean: 0.79374, std: 0.03061, params: {'reg_alpha': 0},
  mean: 0.79433, std: 0.03076, params: {'reg_alpha': 0.001},
  mean: 0.79099, std: 0.02989, params: {'reg_alpha': 0.005},
  mean: 0.79068, std: 0.02953, params: {'reg_alpha': 0.01},
  mean: 0.79160, std: 0.02950, params: {'reg_alpha': 0.05}],
 {'reg_alpha': 0.001},
 0.79432567460197734)

gsearch7.grid_scores_, gsearch7.best_params_, gsearch7.best_score_

([mean: 0.79374, std: 0.03061, params: {'reg_alpha': 0},

mean: 0.79433, std: 0.03076, params: {'reg_alpha': 0.001},

mean: 0.79099, std: 0.02989, params: {'reg_alpha': 0.005},

mean: 0.79068, std: 0.02953, params: {'reg_alpha': 0.01},

mean: 0.79160, std: 0.02950, params: {'reg_alpha': 0.05}],

{'reg_alpha': 0.001},

0.79432567460197734)

これまでにチューニングしてきたパラメータを用いて再度推定を行います。

xgb3 = XGBClassifier(
        learning_rate =0.1,
        n_estimators=1000,
        max_depth=8,
        min_child_weight=1,
        gamma=0.4,
        subsample=0.85,
        colsample_bytree=0.75,
        reg_alpha=0.001,
        objective= 'binary:logistic',
        nthread=4,
        scale_pos_weight=1,
        seed=27)
modelfit(xgb3, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.
Stopping. Best iteration: 153

Model Report
Accuracy : 1
AUC Score (Train): 1.000000
AUC Score (Test): 0.880331

xgb3 = XGBClassifier(

learning_rate =0.1,

n_estimators=1000,

max_depth=8,

min_child_weight=1,

gamma=0.4,

subsample=0.85,

colsample_bytree=0.75,

reg_alpha=0.001,

objective= 'binary:logistic',

nthread=4,

scale_pos_weight=1,

seed=27)

modelfit(xgb3, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.

Stopping. Best iteration: 153

Model Report

Accuracy : 1

AUC Score (Train): 1.000000

AUC Score (Test): 0.880331

ブログであるように試行回数を1,000回から5,000回まで増やしてみます。

xgb4 = XGBClassifier(
        learning_rate =0.01,
        n_estimators=5000,
        max_depth=8,
        min_child_weight=1,
        gamma=0.4,
        subsample=0.85,
        colsample_bytree=0.75,
        reg_alpha=0.001,
        objective= 'binary:logistic',
        nthread=4,
        scale_pos_weight=1,
        seed=27)
modelfit(xgb4, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.
Stopping. Best iteration: 604

Model Report
Accuracy : 0.9951
AUC Score (Train): 0.999955
AUC Score (Test): 0.888000

xgb4 = XGBClassifier(

learning_rate =0.01,

n_estimators=5000,

max_depth=8,

min_child_weight=1,

gamma=0.4,

subsample=0.85,

colsample_bytree=0.75,

reg_alpha=0.001,

objective= 'binary:logistic',

nthread=4,

scale_pos_weight=1,

seed=27)

modelfit(xgb4, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.

Stopping. Best iteration: 604

Model Report

Accuracy : 0.9951

AUC Score (Train): 0.999955

AUC Score (Test): 0.888000

88.8%まで向上しました。色々と数値いじっても、1%高めるだけにとどまってしまうのですね。

とにかく、XGBoostをPythonで実行してパラメータチューニングするという一連の試行がこのコードでできそうなので、今後も使いまわしてみようと思います。

Tokyo.R#53で得たパッケージ情報とその実践

第53回のTokyo.Rで気になったパッケージの情報と実行例をいくつかあげました。スライドなどもろもろの発表はこちらの方のブログ「第53回R勉強会@東京で発表してきた」が非常に詳しく書かれています。

【目次】
・ggradarパッケージ
・proxyパッケージ
・因果推論（CBPSパッケージ）
・MXNetパッケージ
・missForestパッケージ
・RFinanceパッケージ

ggradarパッケージ

簡単にレーダーチャートを作れるパッケージです。こちらのブログを参考にしています。

install.packages("devtools")
devtools::install_github("ricardo-bion/ggradar")

1 2	install.packages("devtools") devtools::install_github("ricardo-bion/ggradar")

企業の職場環境に関してまとめられた某口コミサイトから4個ほどデータを拝借してきました。

> CompanyVoiceData
  company growth stability salary rewarding idea difficulty welfare education
1  google    5.0       5.0    4.9       5.0  4.3        5.0     5.0       4.6
2   yahoo    3.9       5.0    3.2       3.8  3.7        3.9     3.1       3.3
3 recruit    4.4       4.8    5.0       5.0  5.0        5.0     4.0       5.0
4  amazon    5.0       5.0    4.2       4.0  4.2        5.0     3.6       3.3

> CompanyVoiceData

company growth stability salary rewarding idea difficulty welfare education

1 google 5.0 5.0 4.9 5.0 4.3 5.0 5.0 4.6

2 yahoo 3.9 5.0 3.2 3.8 3.7 3.9 3.1 3.3

3 recruit 4.4 4.8 5.0 5.0 5.0 5.0 4.0 5.0

4 amazon 5.0 5.0 4.2 4.0 4.2 5.0 3.6 3.3

ggradarをそのまま使おうとすると、Circular Air Lightというフォントが必要だと怒られるので、参考のブログにある通り、OSXの場合はこちらをダブルクリックでインストールして再起動します。

先ほどのデータに対して、以下のコードを実行すれば非常に簡単にレーダーチャートが作れました。

library("ggradar")
CompanyVoiceData <- data.frame(read.csv(file ="company_voice.csv",header = TRUE))

ggradar(CompanyVoiceData, 
        grid.max = max(CompanyVoiceData[, 2:ncol(CompanyVoiceData)]),
        background.circle.colour = "#ffdd99", #背景色の指定
        background.circle.transparency = 1, #背景色の透明度を指定
        group.line.width = 2, #線の太さの指定
        group.point.size = 6, #シンボルの大きさの指定
        axis.label.size = 5, #軸ラベルサイズの指定
        gridline.min.colour = "#4b61ba", #最小円の線色の指定
        gridline.mid.colour = "#a87963", #中円の線色の指定
        gridline.max.colour = "#e1e6ea", #最大円の線色の指定
        grid.line.width = 1.5, #各円の線の太さの指定
        gridline.min.linetype = "longdash", #線種の指定
        gridline.mid.linetype = "longdash", #線種の指定
        gridline.max.linetype = "longdash") #線種の指定

library("ggradar")

CompanyVoiceData <- data.frame(read.csv(file ="company_voice.csv",header = TRUE))

ggradar(CompanyVoiceData,

grid.max = max(CompanyVoiceData[, 2:ncol(CompanyVoiceData)]),

background.circle.colour = "#ffdd99", #背景色の指定

background.circle.transparency = 1, #背景色の透明度を指定

group.line.width = 2, #線の太さの指定

group.point.size = 6, #シンボルの大きさの指定

axis.label.size = 5, #軸ラベルサイズの指定

gridline.min.colour = "#4b61ba", #最小円の線色の指定

gridline.mid.colour = "#a87963", #中円の線色の指定

gridline.max.colour = "#e1e6ea", #最大円の線色の指定

grid.line.width = 1.5, #各円の線の太さの指定

gridline.min.linetype = "longdash", #線種の指定

gridline.mid.linetype = "longdash", #線種の指定

gridline.max.linetype = "longdash") #線種の指定

proxyパッケージ

距離や類似度を計算するパッケージです。
先ほどのデータに対して類似度と距離を計算してみます。

library(proxy)
> simil(CompanyVoiceData[,-1])
          1         2         3
2 0.2286639                    
3 0.6373648 0.1339713          
4 0.6499133 0.5787506 0.4188571
> dist(CompanyVoiceData[,-1])
         1        2        3
2 3.522783                  
3 1.435270 3.401470         
4 2.269361 1.989975 2.393742

library(proxy)

> simil(CompanyVoiceData[,-1])

1 2 3

2 0.2286639

3 0.6373648 0.1339713

4 0.6499133 0.5787506 0.4188571

> dist(CompanyVoiceData[,-1])

1 2 3

2 3.522783

3 1.435270 3.401470

4 2.269361 1.989975 2.393742

こんな感じで、類似度や距離の計算ができます。

因果推論

こちらはパッケージとかそういうものではなく、既存の関数などで計算できるようです。
こちらのブログ、「調査観察データにおける因果推論(3) – Rによる傾向スコア，IPW推定量，二重にロバストな推定量の算出」に詳しく書かれています。
・glm関数での傾向スコアの算出
・傾向スコアを共変量としてlm関数で回帰分析
・コードを愚直に書いてIPW推定量の算出
・期待値の標準誤差を出すための関数を作成
・DR推定量の算出をするための関数を作成
などで、推定自体は実現できるようです。

ただし、CBPS(Covariate Balancing Propensity Score)というパッケージがあるらしく、このパッケージを用いれば因果推論の計算を行えるようです。

Package ‘CBPS’
以下のようなExampleコードが載っていたので、実行してみましたが、なかなか結果が返ってこなかったので不安になりました。計算が終わるまで10分以上はかかったと思います。

library(CBPS)
data(Blackwell)

form1<-"d.gone.neg ~ d.gone.neg.l1 + d.gone.neg.l2 + d.neg.frac.l3 + camp.length + camp.length +
deminc + base.poll + year.2002 + year.2004 + year.2006 + base.und + office"

##Fitting the models in Imai and Ratkovic (2014)
##Warning: may take a few mintues; setting time.vary to FALSE
##Results in a quicker fit but with poorer balance
fit1 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",
            iterations = NULL, twostep = TRUE, msm.variance = "full", time.vary = TRUE)
fit2 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",
            iterations = NULL, twostep = TRUE, msm.variance = "approx", time.vary = TRUE)

##Assessing balance
bal1 <- balance(fit1)
bal2 <- balance(fit2)

##Effect estimation: Replicating Effect Estimates in
##Table 3 of Imai and Ratkovic (2014)
lm1 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$glm.weights)
lm2 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$weights)
lm3 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit2$weights)
lm4 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$glm.weights)
lm5 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$weights)
lm6 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit2$weights)

library(CBPS)

data(Blackwell)

form1<-"d.gone.neg ~ d.gone.neg.l1 + d.gone.neg.l2 + d.neg.frac.l3 + camp.length + camp.length +

deminc + base.poll + year.2002 + year.2004 + year.2006 + base.und + office"

##Fitting the models in Imai and Ratkovic (2014)

##Warning: may take a few mintues; setting time.vary to FALSE

##Results in a quicker fit but with poorer balance

fit1 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",

iterations = NULL, twostep = TRUE, msm.variance = "full", time.vary = TRUE)

fit2 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",

iterations = NULL, twostep = TRUE, msm.variance = "approx", time.vary = TRUE)

##Assessing balance

bal1 <- balance(fit1)

bal2 <- balance(fit2)

##Effect estimation: Replicating Effect Estimates in

##Table 3 of Imai and Ratkovic (2014)

lm1 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$glm.weights)

lm2 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$weights)

lm3 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit2$weights)

lm4 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$glm.weights)

lm5 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$weights)

lm6 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit2$weights)

MXNet

XGBoostのパッケージを作ったチームが手がけているパッケージで、深層学習を実行できます。

インストール方法はここに書かれています。
Deep Learning for R

install.packages("drat", repos="https://cran.rstudio.com")
drat:::addRepo("dmlc")
install.packages("mxnet")

install.packages("drat", repos="https://cran.rstudio.com")

drat:::addRepo("dmlc")

install.packages("mxnet")

あれ、OSXではエラーが返ってきてライブラリが読み込めないですね。どうやら私のためにあるようなブログ「Installing mxnet for R on Yosemite」があったので、時間を見つけてチャレンジしてみようと思います。

ディープラーニングを用いた回帰分析については、Neural Network with MXNet in Five Minutesにコードがもろもろ載っていますので、チャレンジしてみると良いと思います。

リンク先に載っているのですが、一応コードを以下に記しておきます。

data(BostonHousing, package="mlbench")

train.ind = seq(1, 506, 3)
train.x = data.matrix(BostonHousing[train.ind, -14])
train.y = BostonHousing[train.ind, 14]
test.x = data.matrix(BostonHousing[-train.ind, -14])
test.y = BostonHousing[-train.ind, 14]

# Define the input data
data <- mx.symbol.Variable("data")
# A fully connected hidden layer
# data: input source
# num_hidden: number of neurons in this hidden layer
fc1 <- mx.symbol.FullyConnected(data, num_hidden=1)

# Use linear regression for the output layer
lro <- mx.symbol.LinearRegressionOutput(fc1)

preds = predict(model, test.x)

## Auto detect layout of input matrix, use rowmajor..
sqrt(mean((preds-test.y)^2))

demo.metric.mae <- mx.metric.custom("mae", function(label, pred) {
  res <- mean(abs(label-pred))
  return(res)
})

mx.set.seed(0)
model <- mx.model.FeedForward.create(lro, X=train.x, y=train.y,
                                     ctx=mx.cpu(), num.round=50, array.batch.size=20,
                                     learning.rate=2e-6, momentum=0.9, eval.metric=demo.metric.mae)

data(BostonHousing, package="mlbench")

train.ind = seq(1, 506, 3)

train.x = data.matrix(BostonHousing[train.ind, -14])

train.y = BostonHousing[train.ind, 14]

test.x = data.matrix(BostonHousing[-train.ind, -14])

test.y = BostonHousing[-train.ind, 14]

# Define the input data

data <- mx.symbol.Variable("data")

# A fully connected hidden layer

# data: input source

# num_hidden: number of neurons in this hidden layer

fc1 <- mx.symbol.FullyConnected(data, num_hidden=1)

# Use linear regression for the output layer

lro <- mx.symbol.LinearRegressionOutput(fc1)

preds = predict(model, test.x)

## Auto detect layout of input matrix, use rowmajor..

sqrt(mean((preds-test.y)^2))

demo.metric.mae <- mx.metric.custom("mae", function(label, pred) {

res <- mean(abs(label-pred))

return(res)

})

mx.set.seed(0)

model <- mx.model.FeedForward.create(lro, X=train.x, y=train.y,

ctx=mx.cpu(), num.round=50, array.batch.size=20,

learning.rate=2e-6, momentum=0.9, eval.metric=demo.metric.mae)

missForest

ランダムフォレストを用いて、欠損値補完を行うためのパッケージです。目的変数が欠損していても適用できるようです。
詳しくは、スライドを見ていただいた方がいいですが、以下のプログラムで実行できました。ちなみにスライドはこちら、「Imputation of Missing Values using Random Forest」

library(missForest)
library(dplyr)

#ggplot2のデータセットを読み込む
data(diamonds, package = "ggplot2")
dia.sample <- sample_n(diamonds, size=2000)
dia.sample <- as.data.frame(dia.sample)

#既存データセットに5%の欠損を与える
dia.mis <- prodNA(dia.sample, noNA=0.05)

#補完の実行
dia.imp <- missForest(dia.mis, verbose=TRUE)
dia.imp %>% str(max.level=1)

#補完精度の推定
dia.imp$OOBerror
dia.imp <- missForest(dia.mis, verbose=TRUE, variablewise=TRUE)

#補完精度の検証
mixError(ximp = dia.imp$ximp,
         xmis = dia.mis,
         xtrue = dia.sample)

library(missForest)

library(dplyr)

#ggplot2のデータセットを読み込む

data(diamonds, package = "ggplot2")

dia.sample <- sample_n(diamonds, size=2000)

dia.sample <- as.data.frame(dia.sample)

#既存データセットに5%の欠損を与える

dia.mis <- prodNA(dia.sample, noNA=0.05)

#補完の実行

dia.imp <- missForest(dia.mis, verbose=TRUE)

dia.imp %>% str(max.level=1)

#補完精度の推定

dia.imp$OOBerror

dia.imp <- missForest(dia.mis, verbose=TRUE, variablewise=TRUE)

#補完精度の検証

mixError(ximp = dia.imp$ximp,

xmis = dia.mis,

xtrue = dia.sample)

RFinanceYJ

Yohei Sato, Nobuaki Oshiro, Shinichi Takayanagiさんたちが作った、Yahoo!ファイナンスの株価データを取得できるパッケージです。だいぶ前からあったようですが、使って分析している人は初めて見ました。どうやらYahoo!ファイナンスの仕様によって書き換えていかないといけないようです。「2015-01-20 Rでチャートを書いてみる(9)」のブログに実行可能なプログラムがあります。以下、実行可能なコードを転載いたします。

library(RFinanceYJ)

#API
quoteStockTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval <- substr(time.interval,1,1)
  function.stock <- function(quote.table.item){
    if( xmlSize(quote.table.item) < 5) return(NULL) 
    d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o <- as.number(xmlValue(quote.table.item[[2]]))
    h <- as.number(xmlValue(quote.table.item[[3]]))
    l <- as.number(xmlValue(quote.table.item[[4]]))
    c <- as.number(xmlValue(quote.table.item[[5]]))
    v <- ifelse(xmlSize(quote.table.item) >= 6,as.number(xmlValue(quote.table.item[[6]])),0)
    a <- ifelse(xmlSize(quote.table.item) >= 7,as.number(xmlValue(quote.table.item[[7]])),0)
    return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))
  }
  return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))
}
quoteFundTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval <- substr(time.interval,1,1)
  function.fund <- function(quote.table.item){
    d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    if(time.interval=='monthly'){
      d <- endOfMonth(d)
    }
    c <- as.number(xmlValue(quote.table.item[[2]]))
    v <- as.number(xmlValue(quote.table.item[[3]]))
    return(data.frame(date=d,constant.value=c,NAV=v))
  }
  return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))
}
quoteFXTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval <- substr(time.interval,1,1)
  function.fx <- function(quote.table.item){
    d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o <- as.number(xmlValue(quote.table.item[[2]]))
    h <- as.number(xmlValue(quote.table.item[[3]]))
    l <- as.number(xmlValue(quote.table.item[[4]]))
    c <- as.number(xmlValue(quote.table.item[[5]]))
    return(data.frame(date=d,open=o,high=h,low=l,close=c))
  }
  return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))
}
######  private functions  #####
#get time series data from Yahoo! Finance.
quoteTsData <- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){
  r <- NULL
  result.num <- 51
  financial.data <- data.frame(NULL)
  #start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&c=\\1&a=\\2&b=\\3",since))
  #end   <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&f=\\1&d=\\2&e=\\3",date.end))
  start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&sy=\\1&sm=\\2&sd=\\3",since))
  end   <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&ey=\\1&em=\\2&ed=\\3",date.end))
  
  if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")
  
  extractQuoteTable <- function(r,type){
    if(type %in% c("fund","fx")){
      tbl <- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]
    }
    else{
      tbl <- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]
    }
    return(tbl)
  }
  
  #while( result.num >= 51 ){
  while(1){
    start.num <- start.num + 1
    quote.table <- NULL
    quote.url <- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&p=',start.num,'&tm=',substr(time.interval,1,1),sep="")
    #cat(quote.url)
    #try( r <- xmlRoot(htmlTreeParse(quote.url,error=xmlErrorCumulator(immediate=F))), TRUE)  # これだと取得時にエラーが出た。。
    try(r<-htmlParse(quote.url))
    if( is.null(r) ) stop(paste("Can not access :", quote.url))
    
    #try( quote.table <- r[[2]][[1]][[1]][[16]][[1]][[1]][[1]][[4]][[1]][[1]][[1]], TRUE )
    #try( quote.table <- extractQuoteTable(r,type), TRUE )
    try( quote.table <- xpathApply(r,"//table")[[2]], TRUE )
    
    quote.size<-xmlSize(quote.table)
    #cat(paste("size:",quote.size))
    if(xmlSize(quote.table)<=1){
      return (financial.data)
    }
    if( is.null(quote.table) ){
      if( is.null(financial.data) ){
        stop(paste("Can not quote :", x))
      }else{
        financial.data <- financial.data[order(financial.data$date),]
        return(financial.data)
      }
    }
    
    size <- xmlSize(quote.table)
    for(i in 2:size){
      financial.data <- rbind(financial.data,function.financialproduct(quote.table[[i]]))
    }
    
    #result.num <- xmlSize(quote.table)
    Sys.sleep(1)
  }
  financial.data <- financial.data[order(financial.data$date),]
  return(financial.data)  
}
#convert string formart date to POSIXct object
convertToDate <- function(date.string,time.interval)
{
  #data format is different between monthly and dialy or weekly
  if(any(time.interval==c('d','w'))){
    result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)
  }else if(time.interval=='m'){
    result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)
  }
  return(as.POSIXct(result))
}
#convert string to number.
as.number <- function(string)
{
  return(as.double(as.character(gsub("[^0-9.]", "",string))))
}
#return end of month date.
endOfMonth <- function(date.obj)
{
  startOfMonth     <- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")
  startOfNextMonth <- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")
  return(startOfNextMonth-1)
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

library(RFinanceYJ)

#API

quoteStockTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.stock <- function(quote.table.item){

if( xmlSize(quote.table.item) < 5) return(NULL)

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

v <- ifelse(xmlSize(quote.table.item) >= 6,as.number(xmlValue(quote.table.item[[6]])),0)

a <- ifelse(xmlSize(quote.table.item) >= 7,as.number(xmlValue(quote.table.item[[7]])),0)

return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))

}

return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))

}

quoteFundTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fund <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

if(time.interval=='monthly'){

d <- endOfMonth(d)

}

c <- as.number(xmlValue(quote.table.item[[2]]))

v <- as.number(xmlValue(quote.table.item[[3]]))

return(data.frame(date=d,constant.value=c,NAV=v))

}

return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))

}

quoteFXTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fx <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

return(data.frame(date=d,open=o,high=h,low=l,close=c))

}

return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))

}

###### private functions #####

#get time series data from Yahoo! Finance.

quoteTsData <- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){

r <- NULL

result.num <- 51

financial.data <- data.frame(NULL)

#start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&c=\\1&a=\\2&b=\\3",since))

#end <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&f=\\1&d=\\2&e=\\3",date.end))

start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&sy=\\1&sm=\\2&sd=\\3",since))

end <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&ey=\\1&em=\\2&ed=\\3",date.end))

if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")

extractQuoteTable <- function(r,type){

if(type %in% c("fund","fx")){

tbl <- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]

}

else{

tbl <- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]

}

return(tbl)

}

#while( result.num >= 51 ){

while(1){

start.num <- start.num + 1

quote.table <- NULL

quote.url <- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&p=',start.num,'&tm=',substr(time.interval,1,1),sep="")

#cat(quote.url)

#try( r <- xmlRoot(htmlTreeParse(quote.url,error=xmlErrorCumulator(immediate=F))), TRUE) # これだと取得時にエラーが出た。。

try(r<-htmlParse(quote.url))

if( is.null(r) ) stop(paste("Can not access :", quote.url))

#try( quote.table <- r[[2]][[1]][[1]][[16]][[1]][[1]][[1]][[4]][[1]][[1]][[1]], TRUE )

#try( quote.table <- extractQuoteTable(r,type), TRUE )

try( quote.table <- xpathApply(r,"//table")[[2]], TRUE )

quote.size<-xmlSize(quote.table)

#cat(paste("size:",quote.size))

if(xmlSize(quote.table)<=1){

return (financial.data)

}

if( is.null(quote.table) ){

if( is.null(financial.data) ){

stop(paste("Can not quote :", x))

}else{

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

size <- xmlSize(quote.table)

for(i in 2:size){

financial.data <- rbind(financial.data,function.financialproduct(quote.table[[i]]))

}

#result.num <- xmlSize(quote.table)

Sys.sleep(1)

}

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

#convert string formart date to POSIXct object

convertToDate <- function(date.string,time.interval)

{

#data format is different between monthly and dialy or weekly

if(any(time.interval==c('d','w'))){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)

}else if(time.interval=='m'){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)

}

return(as.POSIXct(result))

}

#convert string to number.

as.number <- function(string)

{

return(as.double(as.character(gsub("[^0-9.]", "",string))))

}

#return end of month date.

endOfMonth <- function(date.obj)

{

startOfMonth <- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")

startOfNextMonth <- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")

return(startOfNextMonth-1)

}

このコードでYahoo!ジャパンの株価を見てみましょう。ちなみに番号は4689です。どうやら上手く取れているようです。

> quoteStockTsData("4689.t",since="2016-01-01")
         date open high low close   volume adj_close
1  2016-05-02  476  483 475   478 18498100       478
2  2016-04-28  504  508 493   496 11966300       496
3  2016-04-27  505  511 495   497 12973800       497
4  2016-04-26  507  508 495   500  7712600       500
5  2016-04-25  513  515 506   509  7350600       509
6  2016-04-22  515  517 509   514  8908900       514
7  2016-04-21  512  517 506   514 13249900       514
8  2016-04-20  511  515 493   506 14455700       506
9  2016-04-19  516  523 511   516 13345800       516
10 2016-04-18  503  509 499   503 10275900       503
11 2016-04-15  504  519 504   513 16962900       513

> quoteStockTsData("4689.t",since="2016-01-01")

date open high low close volume adj_close

1 2016-05-02 476 483 475 478 18498100 478

2 2016-04-28 504 508 493 496 11966300 496

3 2016-04-27 505 511 495 497 12973800 497

4 2016-04-26 507 508 495 500 7712600 500

5 2016-04-25 513 515 506 509 7350600 509

6 2016-04-22 515 517 509 514 8908900 514

7 2016-04-21 512 517 506 514 13249900 514

8 2016-04-20 511 515 493 506 14455700 506

9 2016-04-19 516 523 511 516 13345800 516

10 2016-04-18 503 509 499 503 10275900 503

11 2016-04-15 504 519 504 513 16962900 513

XGBoostやパラメータチューニングの仕方に関する調査

【目次】
・XGBoostとは
・XGBoostで用いるパラメータ一覧
・XGBoostのパラメータチューニング
・参考文献

XGBoostとは

XGBoost (eXtreme Gradient Boosting) は勾配ブースティングアルゴリズムの先進的な実装例で、データサイエンスのコンペであるKaggleで話題となっていた手法です。

ブースティングアルゴリズムとは、弱識別器（weak learners）の集団を直列的に結合することで、予測における正確性を高めようとするアルゴリズムです。任意のt時点において、モデルの予測は以前のt-1時点での結果に基づき重み付けがなされます。正しく予測されたデータに対しては、重みを小さくし、誤って予測されたデータに対しては重みを大きくします。後で学習する識別器ほど、誤ったデータに集中して学習を進めることになります。

以下はブースティングのイメージ図です。

STEP1では全ての学習データに対して、等しい重み付けで学習を行い、決定境界を引きます。これを弱学習器による学習と言います。このケースでは毒キノコを2つ当てており、キノコを５つ当てています。
STEP2ではSTEP1で正しく識別されたデータの重みが下げられ、誤って識別されたデータの重みが上げられています。高く重み付けがなされたデータは決定境界で正しく識別されていますが、他のデータは誤って分類されています。
STEP3においてもSTEP2と同様の傾向があります。このような弱学習器による処理を繰り返すことで識別性能を高めていきます。

最終的にはこのような決定境界を引くことができるような識別器を求めていきます。

勾配ブースティングの勾配とは、ブースティングアルゴリズムにおけるパラメータ推定での最適化手法が勾配降下法に従っているという意味での勾配です。以上が勾配ブースティングモデルの簡素な説明です。

XGBoostで用いるパラメータ一覧

XGBoostで用いるパラメータに関して、大きく分けて3つあります。

1.全体パラメータ・・・XGBoost全体を司る。

2.ブースターパラメータ・・・各ステップでツリーなどのブースティングを司る。

3.学習タスクパラメータ・・・最適化タスクを司る。

以下、3つのパラメータについて、「パラメータ名」・「デフォルトの値」・「役割」・「引数」を表にしています。

1.全体パラメータ

パラメータ名	デフォルトの値	役割	引数
booster	gbtree	実行するモデルのタイプをツリーモデルか線形モデルのどちらかを指定できる。	gbtree: ツリーモデル gblinear: 線形モデル
silent	0	モデルの実行結果を出力するかどうかを決めることができる。モデルを理解する上で、0のままにしておく方が良いとされている。	0:出力結果を表示する。 1:出力結果を表示しない。
nthread	not set	並列処理のためのコア数などを指定できる。フルコアで実行したい場合は何も指定しなければ自動的にフルコアになる。

2.ブースターパラメータ

パラメータ名	デフォルトの値	役割	引数
eta	0.3	学習率を調整できる。小さくすることで、モデルの頑健性を高めることができる。	0.01〜0.2の値になることが多いらしい。
min_child_weigh	1	子ノードにおいて観察されるデータの重み付けの合計値の最小の値で、過学習を避けるために用いられる。高い値にすることで特定のサンプルに観察されるような傾向を学習することを避けられる。ただし、高くし過ぎるとフィッティングが悪くなる。
max_depth	6	木の深さの最大値過学習を制御するために用いられる。高いと過学習しやすくなる。	3〜10の値になることが多いらしい。
max_leaf_nodes		木の終端ノードの最大値 max_depthの代わりに用いる n本を指定したら、n^2個の枝を生み出す。これが指定された場合は、max_depthは無効化される。
gamma	0	分割が、損失関数の減少に繋がる場合にのみノードの分割を行う。モデルをより保守的にする。値は損失関数に応じて大きく変わり、チューニングが必要である。
max_delta_step	0	各木のウェイトの推定に制約をかけることができる。 0の場合は制約なしで、正数値を取るとモデルがより保守的になる。通常は必要とされないが、不均衡データの分類の際に用いる。
subsample	1	各木においてランダムに抽出される標本の割合小さくすることで、過学習を避けることができるが保守的なモデルとなる。	0.5〜1の値になることが多いらしい。
colsample_bytree	1	各木においてランダムに抽出される列の割合	0.5〜1の値になることが多いらしい。
colsample_bylevel	1	各レベル単位での、分割における列のsubsample比率 subsampleとcolsample_bytreeで十分なので、あまり使わないが、探索してみるのも良いかもしれない。
lambda	1	重みに関するL2正則化項多くのデータサイエンティストは使わないが、過学習を避けるためには用いられるべき。
alpha	0	重みに関するL1正則化項高次元の場合に用いるらしい。
scale_pos_weight	1	不均衡データの際に、0以上の値を取ることで、収束を早めることができる。

3.学習タスクパラメータ

パラメータ名	デフォルトの値	役割	引数
objective	reg:linear	最小化させるべき損失関数を指定する。	binary:logistic→2項分類で確率を返す。 multi:softmax→多項分類でクラスの値を返す。（num_classでクラス数の指定が必要） multi:softprob→softmaxと同じだが、確率を返す。
eval_metric	according to objective	検証を行うためのデータの評価指標	rmse – root mean square error mae – mean absolute error logloss – negative log-likelihood error – Binary classification error rate (0.5 threshold) merror – Multiclass classification error rate mlogloss – Multiclass logloss auc: Area under the curve
seed	0	ランダムなシード番号。再現可能なデータを生み出すために、あるいはパラメータチューニングの際に用いる。

XGBoostのパラメータチューニング

複数のパラメータからなるXGBoostのチューニングは非常に複雑で、理想的なパラメータについてはケースバイケースで何とも言えないそうです。

参考文献のブログにパラメータチューニングの一般的アプローチについて触れられていたので、紹介します。

4.学習率を下げ、最適なパラメータを決定する。

具体的な実行に関するPythonスクリプトはこちらのGithubで紹介されています。（iPython）
Analytics_Vidhya/Articles/Parameter_Tuning_XGBoost_with_Example/XGBoost models.ipynb
この方法に従って、自社で抱えているモデルのチューニングにチャレンジしてみようと思います。

参考文献

Quick Introduction to Boosting Algorithms in Machine Learning
Complete Guide to Parameter Tuning in Gradient Boosting (GBM) in Python
Complete Guide to Parameter Tuning in XGBoost (with codes in Python)
勾配ブースティングについてざっくりと説明する
 xgboost のパラメータ
 OS X で XGBoost ＆ xgboost4j をビルドする手順 2016-03-07 版

GensimのHDP(Hierarchical Dirichlet Process)をクラシック音楽情報に対して試してみる

HDP(Hierarchical Dirichlet Process)いわゆる階層ディリクレ過程を実行できるモデルがPythonのGensimライブラリにあるという情報から、あまり実行例も見当たらないので、チャレンジしてみました。

HDP(Hierarchical Dirichlet Process)

HDP(Hierarchical Dirichlet Process)は文書集合全体のトピック数と文書ごとのトピック数の推定を行うことができる手法で、中華料理店フランチャイズという仕組みを用いています。通常のLDAなどでは、分析者が任意のトピック数を決める必要がありましたが、与えられたデータからその数を推定するため、その必要がないというのがHDPを使うことの利点であると思われます。

実行までの流れ

ざっくりですが、
・コーパスの準備・文書の分かち書き（名詞のみ）
・HDPの実行
という流れです。

ちなみに実行環境は
MacBook Pro
OS X Yosemite 10.10.5
2.6 GHz Intel Core i5
メモリ8GBです。

コーパスの準備

今回は、以前手に入れた某辞典サイトのクラシック音楽情報1800件のテキストデータ（1行に1件分の文字列が入っているデータで16MBくらい）があるので、それをコーパスとして使います。参考情報として挙げているブログの助けを借りて、文書単位でMeCabにより形態素解析で分かち書きした結果から、意味を持ちやすい品詞として、「名詞」に該当するもののみを結果として返す以下のPythonスクリプトを用いました。結果はtmep.txtとして出力されます。もっと良いやり方があると思いますが、目的は達成できると思います。ちなみに、MeCab Neologd(ネオログディー)という、固有名詞などに強いシステム辞書を活用してみたかったので、その利用を前提として書いています。MeCab Neologd(ネオログディー)のインストール関連の情報は参考情報にありますので、チャレンジしてみてください。（OSXかUbuntuの方が進めやすいと思います。）

#coding:utf-8
import MeCab
import sys
sys.stdout = open("tmep.txt","w")

def extractKeyword(text):
	u"""textを形態素解析して、名詞のみのリストを返す"""
	tagger = MeCab.Tagger(' -d /usr/local/lib/mecab/dic/mecab-ipadic-neologd')
	#デフォルトの辞書を用いる場合は tagger = MeCab.Tagger('')	
	encoded_text = text.encode('utf-8')
	node = tagger.parseToNode(encoded_text).next
	keywords = []
	while node:
		if node.feature.split(",")[0] == "名詞":
			keywords.append(node.surface)
		node = node.next
	return keywords

def splitDocument(documents):
	u"""文章集合を受け取り、名詞のみ空白区切りの文章にして返す"""
	splitted_documents = []
	for d in documents:
		keywords = extractKeyword(d)
		splitted_documents.append(' '.join(keywords))
	return splitted_documents


if __name__ == "__main__":
	document_text = open('music_text.txt')
	raw_documents = document_text.readlines()

	# 空白区切りの文字列を入れるリスト
	splitted_documents = splitDocument(raw_documents)

	for d in splitted_documents:
		print d
	print ''


sys.stdout = sys.__stdout__

#coding:utf-8

import MeCab

import sys

sys.stdout = open("tmep.txt","w")

def extractKeyword(text):

u"""textを形態素解析して、名詞のみのリストを返す"""

tagger = MeCab.Tagger(' -d /usr/local/lib/mecab/dic/mecab-ipadic-neologd')

#デフォルトの辞書を用いる場合は tagger = MeCab.Tagger('')

encoded_text = text.encode('utf-8')

node = tagger.parseToNode(encoded_text).next

keywords = []

while node:

if node.feature.split(",")[0] == "名詞":

keywords.append(node.surface)

node = node.next

return keywords

def splitDocument(documents):

u"""文章集合を受け取り、名詞のみ空白区切りの文章にして返す"""

splitted_documents = []

for d in documents:

keywords = extractKeyword(d)

splitted_documents.append(' '.join(keywords))

return splitted_documents

if __name__ == "__main__":

document_text = open('music_text.txt')

raw_documents = document_text.readlines()

# 空白区切りの文字列を入れるリスト

splitted_documents = splitDocument(raw_documents)

for d in splitted_documents:

print d

print ''

sys.stdout = sys.__stdout__

こちらのスクリプトをターミナルで実行します。(解析するディレクトリ下で実施しています。)

python MeCab_Norm_Extract.py music_text.txt > music_text_wakati_norm_neo.txt

1	python MeCab_Norm_Extract.py music_text.txt > music_text_wakati_norm_neo.txt

HDPの実行

以下のPythonスクリプトで実行しています。

# -*- coding: utf-8 -*-
from gensim import　models,corpora
import pandas

#ドキュメントからLDAなどの分析用コーパスを作成
corpus = gensim.corpora.TextCorpus('music_text_wakati_norm_neo')

#HDPモデルの推定
model = models.hdpmodel.HdpModel(
corpus,
id2word=corpus.dictionary,
alpha=0.1)

#各文書のトピックの重みを保存
topics = [model[c] for c in corpus]
print(topics[0])

#各トピックごとの単語の抽出（topicsの引数を-1にすることで、ありったけのトピックを結果として返してくれます。）
model.print_topics(topics=-1, topn=10)

#文書ごとに割り当てられたトピックの確率をCSVで出力
mixture = [dict(model[x]) for x in corpus]
pandas.DataFrame(mixture).to_csv("topic_for_corpus.csv")

#トピックごとの上位10語をCSVで出力
topicdata =model.print_topics(topics=-1, topn=10)
pandas.DataFrame(topicdata).to_csv("topic_detail.csv")

# -*- coding: utf-8 -*-

from gensim import　models,corpora

import pandas

#ドキュメントからLDAなどの分析用コーパスを作成

corpus = gensim.corpora.TextCorpus('music_text_wakati_norm_neo')

#HDPモデルの推定

model = models.hdpmodel.HdpModel(

corpus,

id2word=corpus.dictionary,

alpha=0.1)

#各文書のトピックの重みを保存

topics = [model[c] for c in corpus]

print(topics[0])

#各トピックごとの単語の抽出（topicsの引数を-1にすることで、ありったけのトピックを結果として返してくれます。）

model.print_topics(topics=-1, topn=10)

#文書ごとに割り当てられたトピックの確率をCSVで出力

mixture = [dict(model[x]) for x in corpus]

pandas.DataFrame(mixture).to_csv("topic_for_corpus.csv")

#トピックごとの上位10語をCSVで出力

topicdata =model.print_topics(topics=-1, topn=10)

pandas.DataFrame(topicdata).to_csv("topic_detail.csv")

HDPの結果について

topic_detail.csvの結果を見たところ、トピックの数が150個もあって、「本当にトピックの数を自動で決めれているのかなぁ」と不安に思ったのですが、実際に各文書に割り当てられているトピックの数は、先ほど出力したtopic_for_corpus.csvで見ると60個でした。そのため、今回、HDPに従って決まったトピック数は60ということになります。さらに不安に思ったので、Stack Over Flowで調べていたんですが、トピックは150個出るけど確率が割り振られていないはずと回答されていました。（ Hierarchical Dirichlet Process Gensim topic number independent of corpus size ）

出現頻度の高い上位10のトピックは以下の通りです。

['topic 0: 0.011*番 + 0.008*楽章 + 0.008*作品 + 0.007*作曲 + 0.006*演奏 + 0.006*日 + 0.006*主題 + 0.006*案内 + 0.006*月 + 0.005*的',
 'topic 1: 0.006*日 + 0.006*月 + 0.006*作品 + 0.005*オペラ + 0.005*案内 + 0.005*作曲 + 0.005*演奏 + 0.004*クラシック音楽 + 0.004*ため + 0.004*幕',
 'topic 2: 0.006*演奏 + 0.006*作品 + 0.006*番 + 0.006*案内 + 0.005*作曲 + 0.005*日 + 0.005*月 + 0.005*クラシック音楽 + 0.004*的 + 0.004*ポータル',
 'topic 3: 0.006*幕 + 0.005*の + 0.004*日 + 0.004*楽章 + 0.004*的 + 0.003*作品 + 0.003*人 + 0.003*伯爵 + 0.003*オペラ + 0.003*演奏',
 'topic 4: 0.009*ワーグナー + 0.007*幕 + 0.005*作品 + 0.005*日 + 0.005*オペラ + 0.005*的 + 0.004*月 + 0.004*作曲 + 0.004*動機 + 0.004*ジークフリート',
 'topic 5: 0.011*カルメン + 0.009*幕 + 0.007*オペラ + 0.004*ホセ + 0.004*作品 + 0.004*ため + 0.003*日 + 0.003*月 + 0.003*ビゼー + 0.003*作曲',
 'topic 6: 0.005*モーツァルト + 0.004*聖歌 + 0.004*ため + 0.004*的 + 0.003*音 + 0.003*人 + 0.003*オペラ + 0.003*の + 0.003*作品 + 0.003*作曲',
 'topic 7: 0.026*番 + 0.017*交響 + 0.006*ため + 0.005*ピアノ協奏曲 + 0.005*作品 + 0.005*案内 + 0.005*協奏曲 + 0.003*ヴァイオリン + 0.003*ヘルプ + 0.003*最近',
 'topic 8: 0.007*オペラ + 0.006*幕 + 0.005*トゥーランドット + 0.004*プッチーニ + 0.004*初演 + 0.003*作曲 + 0.003*作品 + 0.003*日 + 0.003*の + 0.003*月',
 'topic 9: 0.005*作品 + 0.004*パーヴァリ + 0.004*上演 + 0.004*音楽 + 0.004*初演 + 0.004*日 + 0.004*ため + 0.003*幕 + 0.003*月 + 0.003*中',
 'topic 10: 0.007*トスカ + 0.006*幕 + 0.005*オペラ + 0.005*作品 + 0.004*彼 + 0.004*マノン + 0.003*カヴァラドッシ + 0.003*の + 0.003*スカ + 0.003*ルピア']

['topic 0: 0.011*番 + 0.008*楽章 + 0.008*作品 + 0.007*作曲 + 0.006*演奏 + 0.006*日 + 0.006*主題 + 0.006*案内 + 0.006*月 + 0.005*的',

'topic 1: 0.006*日 + 0.006*月 + 0.006*作品 + 0.005*オペラ + 0.005*案内 + 0.005*作曲 + 0.005*演奏 + 0.004*クラシック音楽 + 0.004*ため + 0.004*幕',

'topic 2: 0.006*演奏 + 0.006*作品 + 0.006*番 + 0.006*案内 + 0.005*作曲 + 0.005*日 + 0.005*月 + 0.005*クラシック音楽 + 0.004*的 + 0.004*ポータル',

'topic 3: 0.006*幕 + 0.005*の + 0.004*日 + 0.004*楽章 + 0.004*的 + 0.003*作品 + 0.003*人 + 0.003*伯爵 + 0.003*オペラ + 0.003*演奏',

'topic 4: 0.009*ワーグナー + 0.007*幕 + 0.005*作品 + 0.005*日 + 0.005*オペラ + 0.005*的 + 0.004*月 + 0.004*作曲 + 0.004*動機 + 0.004*ジークフリート',

'topic 5: 0.011*カルメン + 0.009*幕 + 0.007*オペラ + 0.004*ホセ + 0.004*作品 + 0.004*ため + 0.003*日 + 0.003*月 + 0.003*ビゼー + 0.003*作曲',

'topic 6: 0.005*モーツァルト + 0.004*聖歌 + 0.004*ため + 0.004*的 + 0.003*音 + 0.003*人 + 0.003*オペラ + 0.003*の + 0.003*作品 + 0.003*作曲',

'topic 7: 0.026*番 + 0.017*交響 + 0.006*ため + 0.005*ピアノ協奏曲 + 0.005*作品 + 0.005*案内 + 0.005*協奏曲 + 0.003*ヴァイオリン + 0.003*ヘルプ + 0.003*最近',

'topic 8: 0.007*オペラ + 0.006*幕 + 0.005*トゥーランドット + 0.004*プッチーニ + 0.004*初演 + 0.003*作曲 + 0.003*作品 + 0.003*日 + 0.003*の + 0.003*月',

'topic 9: 0.005*作品 + 0.004*パーヴァリ + 0.004*上演 + 0.004*音楽 + 0.004*初演 + 0.004*日 + 0.004*ため + 0.003*幕 + 0.003*月 + 0.003*中',

'topic 10: 0.007*トスカ + 0.006*幕 + 0.005*オペラ + 0.005*作品 + 0.004*彼 + 0.004*マノン + 0.003*カヴァラドッシ + 0.003*の + 0.003*スカ + 0.003*ルピア']

加えて、トピックごとに文書に割り当てられた数を集計してみましたが、topic0が圧倒的に多く、コーパスの特性上、含まれやすい情報がここに集まっているのではないかと思います。幅広いテーマを抽出できるかと期待していたのですが、やたらと個別具体的な「トゥーランドット」や「ワーグナー」や「カルメン」などがトピックの上位単語に上がってきています。実行方法を間違えているかもしれないし、パラメータチューニングなどをもっと頑張れば、幅広いトピックを得ることができるかもしれないので、今後の課題としたいです。

参考情報

・トピックモデルについて
machine_learning_python/topic.md at master · poiuiop/machine_learning_python · GitHub

Factorization Machine(FM)の概要

パッケージ紹介とインストール

サンプルデータの構造把握と前処理

FMの実行

結果

参考文献

はじめに

紹介する研究

オンライン広告におけるスパース性と遅れコンバージョンを考慮した予測モデル

目的と結果

対象となるデータ

手法の概要

PDFのリンク

テキストマイニングを用いた転職サイトの会員離脱予測

目的と結果

対象となるデータ

手法の概要

PDFのリンク

インターネット広告におけるスパースなユーザー行動ベクトルからのユーザー特徴抽出モデル

目的と結果

対象となるデータ

手法の概要

PDFのリンク

ユーザーのページビュー系列からのコンバージョン予測

目的と結果

対象となるデータ

手法の概要

PDFのリンク

SNSにおけるユーザの行動のモチベーション要因の分析

目的と結果

対象となるデータ

手法の概要

PDFのリンク

状況に応じた楽曲推薦に向けたソーシャルデータ分析

目的と結果

対象となるデータ

手法の概要

PDFのリンク

テキストマイニングを用いた口コミ分析による点数評価の信頼性確認手法

目的と結果

対象となるデータ

手法の概要

PDFのリンク

不動産仲介マーケティングのためのユーザ行動予測 Prediction of user behaviors for real estate brokerage marketing

目的と結果

対象となるデータ

手法の概要

PDFのリンク

SUUMOでの不動産データ活用の取り組みと未来

目的と結果

対象となるデータ

手法の概要

PDFのリンク

感想

参考文献

不動産仲介マーケティングのためのユーザ行動予測　Prediction of user behaviors for real estate brokerage marketing