まずはデータを集める

ということで、競馬の統計解析をしたいわけなのだが、解析するためのデータがなければ何も始まらない。
まずは、競馬のデータを以下のサイトからスクレイピングして取ってくることにする。
netkeiba.com - 競馬データベース
netkeiba.comでスピード指数（ある基準を元に走破タイムを数値化したもの）や馬場指数（馬場コンディションを数値化したもの）を閲覧するには有料会員に登録する必要がある。私は有料会員に登録した上でスピード指数や馬場指数まで含めてスクレイピングを行った。
以下にスクレイピング＆素性作成用のScalaコードを公開する。
github.com
ちなみにデータ解析はデータを解析できる形に持っていくまでが全工程の九割を占めると言われている。実際私もこのスクレイピング＆素性作成用スクリプトを作成するのに数週間はかけている*4。このスクリプトを無料で使える皆さんは幸運である。
作成された素性は最終的にSQLiteに格納されるようになっている。このコードを使うのにnetkeiba.comの有料会員に登録する必要はないが、その場合はスピード指数や馬場指数のカラムにはNULL値が入ることになるので気をつけて欲しい。

何を予測するのか

データが集まった所で、次に「何を」予測するのか決めよう。
私が調べた限りでは、競馬の予測には２つの方法がある。*5

1. あるレースに関する情報を入力として、そのレースが荒れるか否か（一番人気の馬が一着になるかどうか）を予測
2. ある馬の過去の勝率や騎手の勝率などを入力として、その馬がレースにおいて何着になるかを予測

例えば以下の本では前者の方法で予測を行っている。

実践データマイニング―金融・競馬予測の科学

• 作者: 月本洋
• 出版社/メーカー: オーム社
• 発売日: 1999/12
• メディア: 単行本
• 購入: 5人 クリック: 74回
• この商品を含むブログ (3件) を見る

つまり、個別の馬に関するデータを入力とし、その馬がレースで一着になるかどうかの二値を出力とする統計モデルを作成するわけである。
なお今回は、予測するのはレースの着順ではなくあくまでも「一着になるかどうか」の二値だけにする。

なぜこうするのかというと、競馬ではレースの途中で騎手が「このままでは上位になれないな」と気付いたとき、馬を無駄に疲れさせないためにあえて遅く走らせることがあるのだそうだ（競馬は着順が上位じゃないと賞金が貰えないため）。つまり、着順が上位ならばその馬には実力があると言えるが、着順が下位だからといって必ずしも実力が無いとはいえないのだ。だから「厳密な着順の数値」ではなく「一着になるかどうかの二値」だけを予測するシンプルなモデルを作成したほうがうまくいくようだ。*6（参考：Identifying winners of competitive events: A SVM-based classification model for horserace prediction）

気をつけないといけないのは、一着になった馬は少ない一方で、一着にならなかった馬はたくさんいるということだ。このままだと学習データが不均衡*7になってしまい、予測モデルを作成すると偏ったモデルが出来てしまう。不均衡データを扱う方法はいくつかあるが、今回は面倒臭いので多い方のクラス（一着にならなかった馬）のデータをサンプリングで減らしてしまうことにする。

何を統計モデルの入力とするのか

次に問題なのは、統計モデルの入力に何の変数を使うかだ。
このモデルの入力として、私は以下の素性を使うことにした。

このリストは、私が競馬関連の本とか論文とかを読んで「なんとなく良さそう」と思った変数をかき集めただけなので、これらの変数を使うことに必然性があるわけではない。他の変数を使った場合にどうなるか気になるという方は自分でコードを弄って試すべし。

randomForestを使って予測してみる

予測モデルの作成にはRのrandomForestパッケージを使うことにする。random forestとは2001年にLeo Breiman によって提案された教師あり学習のアルゴリズムである。このブログを見に来るような人には解説の必要はないかもしれないが、ざっくり言うと、decision treeはbias-variance分解で言うところのvariance（学習結果の不安定性）が高いのでbaggingと素性のsamplingを適用してみたらvarianceが下がって汎化性能アップしました、というアルゴリズムがrandom forestである。*8
それでは、実際にRのrandomForestパッケージを使って予測モデルを作成してみよう。

> library(randomForest)
> library(RSQLite)
>  
> randomRows <- function(df, n) {
+   df[sample(nrow(df),n),]
+ }
>  
> downSample <- function(df) {
+   c1 <- df[df$order_of_finish == "TRUE",]
+   c2 <- df[df$order_of_finish == "FALSE",]
+   size <- min(nrow(c1), nrow(c2))
+   rbind(randomRows(c1,size), randomRows(c2,size))
+ }
>  
> drv <- dbDriver('SQLite')
>  
> conn <- dbConnect(drv, dbname='race.db')
>  
> rs <- dbSendQuery(conn, 
+ 'select 
+    order_of_finish,
+    race_id,
+    horse_number,
+    grade,
+    age,
+    avgsr4,
+    avgWin4,
+    dhweight,
+    disRoc,
+    r.distance,
+    dsl,
+    enterTimes,
+    eps,
+    hweight,
+    jwinper,
+    odds,
+    owinper,
+    preSRa,
+    sex,
+    f.surface,
+    surfaceScore,
+    twinper,
+    f.weather,
+    weight,
+    winRun,
+    jEps,
+    jAvgWin4,
+    preOOF,
+    pre2OOF,
+    month,
+    runningStyle,
+    preLastPhase,
+    lateStartPer,
+    course,
+    placeCode,
+    race_number
+ from 
+    feature f 
+ inner join 
+    race_info r 
+ on 
+    f.race_id = r.id
+ where 
+    order_of_finish is not null
+ and 
+    preSRa is not null
+ limit 250000')
>  
> allData <- fetch(rs, n = -1)
> 
> dbClearResult(rs)
[1] TRUE
> dbDisconnect(conn)
[1] TRUE
> 
> #カテゴリ変数をファクターに変換しておく
> allData$placeCode <- factor(allData$placeCode)
> allData$month     <- factor(allData$month)
> allData$grade     <- factor(allData$grade)
> allData$sex       <- factor(allData$sex)
> allData$weather   <- factor(allData$weather)
> allData$surface   <- factor(allData$surface)
> allData$course    <- factor(allData$course)
> 
> #負担重量/馬体重を素性に追加
> allData$weightper <- allData$weight / allData$hweight
> 
> #オッズを支持率に変換
> allData$support <- 0.788 / (allData$odds - 0.1)
> allData$odds <- NULL
> 
> #着順をカテゴリ変数に変換
> allData$order_of_finish <- factor(allData$order_of_finish == 1)
> 
> #クラスバランスを50/50にする
> allData.s <- downSample(na.omit(allData))
> allData.s <- allData.s[order(allData.s$race_id),]
> 
> #今回の実験で使用するデータのサンプル数
> nrow(allData.s)
[1] 30428
> 
> #データを学習用25428サンプルとテスト用5000サンプルに分割する
> train <- allData.s[1:(nrow(allData.s)-5000),]
> test  <- allData.s[(nrow(allData.s)-5001):nrow(allData.s),]
> 
> #予測モデルを作成
> (rf.model1 <- randomForest(
+   order_of_finish ~ . - support - race_id, train))

Call:
 randomForest(formula = order_of_finish ~ . - support - race_id,      data = train) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 5

        OOB estimate of  error rate: 29.72%
Confusion matrix:
      FALSE TRUE class.error
FALSE  8420 4362   0.3412611
TRUE   3196 9450   0.2527281
>
> #素性の重要度を見てみる
> importance(rf.model1)
             MeanDecreaseGini
horse_number        276.57124
grade               191.53030
age                 210.04150
avgsr4              545.24005
avgWin4             526.77427
dhweight            296.32679
disRoc              443.31973
distance            232.20557
dsl                 371.28809
enterTimes          332.80342
eps                 682.54396
hweight             393.27570
jwinper             417.62300
owinper             366.49348
preSRa              536.27096
sex                  62.81792
surface              45.83353
surfaceScore        361.17891
twinper             348.52685
weather             123.82181
weight              165.54897
winRun              246.36929
jEps                603.00998
jAvgWin4            140.99460
preOOF              870.35176
pre2OOF             475.39642
month               737.97377
runningStyle        456.73422
preLastPhase        408.51575
lateStartPer        250.49252
course               42.43917
placeCode           564.23156
race_number         278.57604
weightper           430.13985
>
> #テストデータで予測力を見てみる
> pred <- predict(rf.model1, test)
> tbl <- table(pred, test$order_of_finish)
> sum(diag(tbl)) / sum(tbl)
[1] 0.7067173

OOBエラーとテストデータでの正解率が共に約70%になっている。50%を超えているので、このモデルに予測力があることは確かなようだ。
しかし本番はここからである。問題は、このモデルの予測力が他の馬券購入者達の予測力に勝てるかどうかだ。
「他の馬券購入者達の予測」を表すモデルとして、以下の素性だけを用いて学習したモデルを使用する。

単勝オッズから逆算された支持率は「他の馬券購入者達の予測」そのものである。だから、もし競馬市場が効率的であるならば、この支持率を使ったモデルを超える予測精度は生み出せないはずである。なので、このモデルの予測精度を超えられるかどうかが競馬市場の効率性を測る一つの目安となる。

> #支持率だけを用いて予測モデルを作成する
> (rf.model2 <- randomForest(
+   order_of_finish ~ support, train))

Call:
 randomForest(formula = order_of_finish ~ support, data = train) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 1

        OOB estimate of  error rate: 25.7%
Confusion matrix:
      FALSE  TRUE class.error
FALSE  8734  4048   0.3166954
TRUE   2486 10160   0.1965839
> 
> pred <- predict(rf.model2, test)
> tbl <- table(pred, test$order_of_finish)
> sum(diag(tbl)) / sum(tbl)
[1] 0.7379048

このモデルの予測精度は約74%である。
残念ながら私のモデルは70%なので予測力で負けている…。

> racewiseFeature <- 
+ c("avgsr4",
+ "avgWin4",
+ "dhweight",
+ "disRoc",
+ "dsl",
+ "enterTimes",
+ "eps",
+ "hweight",
+ "jwinper",
+ "owinper",
+ "preSRa",
+ "twinper",
+ "weight",
+ "jEps",
+ "jAvgWin4",
+ "preOOF",
+ "pre2OOF",
+ "runningStyle",
+ "preLastPhase",
+ "lateStartPer",
+ "weightper")
> 
> #レース毎にデータを分割する
> splited.allData <- split(allData, allData$race_id)
> 
> scaled.allData <- unsplit(
+   lapply(splited.allData,
+     function(rw) {
+       data.frame(
+         order_of_finish = rw$order_of_finish,
+         race_id = rw$race_id,
+         age = rw$age,
+         grade = rw$grade,
+         distance = rw$distance,
+         sex = rw$sex,
+         weather = rw$weather,
+         course = rw$course,
+         surface = rw$surface,
+         surfaceScore = rw$surfaceScore,
+         horse_number = rw$horse_number,
+         placeCode = rw$placeCode,
+         race_number = rw$race_number,
+         support = rw$support,
+         scale(rw[,racewiseFeature])) #ここで正規化している
+     }),
+   allData$race_id)
> 
> scaled.allData$order_of_finish = factor(scaled.allData$order_of_finish)
> 
> is.nan.df <- function(x) do.call(cbind, lapply(x, is.nan))
> scaled.allData[is.nan.df(scaled.allData)] <- 0
> 
> scaled.allData <- downSample(na.omit(scaled.allData))
> scaled.allData <- scaled.allData[order(scaled.allData$race_id),]
>
> #データを学習用とテスト用に分割する
> scaled.train <- scaled.allData[1:(nrow(scaled.allData)-5000),]
> scaled.test  <- scaled.allData[(nrow(scaled.allData)-5001):nrow(scaled.allData),]
> 
> #レース毎に正規化されたデータで予測モデルを作成
> (rf.model3 <- randomForest(
+   order_of_finish ~ . - support - race_id, scaled.train))

Call:
 randomForest(formula = order_of_finish ~ . - support - race_id,      data = scaled.train) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 5

        OOB estimate of  error rate: 28.63%
Confusion matrix:
      FALSE TRUE class.error
FALSE  8739 4047   0.3165181
TRUE   3234 9408   0.2558140
> 
> #素性の重要度を見てみる
> importance(rf.model3)
             MeanDecreaseGini
age                 138.15954
grade               157.86619
distance            192.87544
sex                  55.18635
weather              92.09389
course               33.38138
surface              33.58647
surfaceScore        287.95836
horse_number        222.12282
placeCode           537.07988
race_number         193.98961
avgsr4              858.85621
avgWin4             726.16178
dhweight            345.24014
disRoc              371.22814
dsl                 363.05980
enterTimes          357.92536
eps                1005.00112
hweight             366.85535
jwinper             471.85535
owinper             367.94282
preSRa              890.83216
twinper             381.33466
weight              336.16596
jEps                530.81950
jAvgWin4            352.48784
preOOF              794.77337
pre2OOF             500.63913
runningStyle        358.16418
preLastPhase        383.60317
lateStartPer        338.66961
weightper           359.51054
> 
> #テストデータで予測力を見てみる
> pred <- predict(rf.model3, scaled.test)
> tbl <- table(pred, scaled.test$order_of_finish)
> sum(diag(tbl)) / sum(tbl)
[1] 0.7221112

> #絶対的能力値モデルの素性に支持率を追加して予測モデルを作成
> (rf.model4 <- randomForest(
+   order_of_finish ~ . - race_id, train))

Call:
 randomForest(formula = order_of_finish ~ . - race_id, data = train) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 6

        OOB estimate of  error rate: 24.88%
Confusion matrix:
      FALSE  TRUE class.error
FALSE  8967  3793   0.2972571
TRUE   2534 10134   0.2000316
> 
> #テストデータで予測力を見てみる
> pred <- predict(rf.model4, test)
> tbl <- table(pred, test$order_of_finish)
> sum(diag(tbl)) / sum(tbl)
[1] 0.7491004
> 
> #相対的能力差モデルの素性に支持率を追加して予測モデルを作成
> (rf.model5 <- randomForest(
+   order_of_finish ~ . - race_id, scaled.train))

Call:
 randomForest(formula = order_of_finish ~ . - race_id, data = scaled.train) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 5

        OOB estimate of  error rate: 25.26%
Confusion matrix:
      FALSE  TRUE class.error
FALSE  8936  3850   0.3011106
TRUE   2572 10070   0.2034488
> 
> #テストデータで予測力を見てみる
> pred <- predict(rf.model5, scaled.test)
> tbl <- table(pred, scaled.test$order_of_finish)
> sum(diag(tbl)) / sum(tbl)
[1] 0.7457017

今後の予定

• JRDBのデータを使う。JRDBには、IDMや馬の体格、蹄など、他にはないデータが豊富に含まれているので、予測精度向上の余地がありそう。ちなみにJRDBのデータは外れ馬券裁判で話題になった卍氏も使用している。
• 今回は最終支持率をそのまま素性に追加したが、実際に利用可能なデータはレース開始直前の支持率であり、最終支持率とはズレがあるかもしれない。レース開始15分前の単勝オッズがJRDBに含まれているのでこれを使うようにする。
• 今回は中央競馬のデータを使っているが、地方競馬の方が儲けられる可能性が高そう。なぜなら地方競馬のほうが注目度が低く、市場の効率性が低そうなので。（その代わりにレース開始直前のオッズの変動が激しいみたいだけど…）
• 血統のデータを活用する。血統をどう数値化するかがちょっと悩ましいが、多分JRA-VANがやっている親馬をカテゴリ変数にしてしまう方法が一番簡単。
• レースが荒れるかどうかを予測するほうが簡単らしいので、後でそっちの方法も試そう…と思って今ちょっとだけやってみたけどあんまりうまくいかないぞこれ…
• そもそもの目的は競馬で儲けることであり、そのためには予測力ではなく回収率を高めなければならない。なので強化学習や遺伝的アルゴリズムを使って回収率が高くなるように学習させる方が望みがありそう。（教師あり学習とは違い強化学習や遺伝的アルゴリズムであれば回収率の高さそのものを最大化するように学習させることができる）

参考ページ

今回の記事を書くにあたって、私が最も参考にしたのはJRA-VANの予測モデル解説と卍氏の書籍（これとこれ）、そしてStefan Lessmannの競馬論文の三つである。「お前の解説は下手すぎて意味わからん」という方はこれらのページも参考にされたし。