walkingsdbgの日記

データ分析技術（機械学習・統計等）についてのメモを書きます。

トップ > 競馬 > 競走馬の「強さ」をベイズ統計モデリングで推定する

競走馬の「強さ」をベイズ統計モデリングで推定する

競馬

はじめに

学生時代にほんの少し競馬をかじったことがあるのですが、競馬の着順予測は考えるべき要素が多く、自分の頭の中でそれらの要素を考慮して予測するのは難しいと感じたので、１年足らずで止めてしまいました。ですが、「StanとRでベイズ統計モデリング」10.2節の棋士の「強さ」を推定する手法を知り、この手法を活用すれば、着順予測をシステマティックに行えるのではないかと考えるようになりました。そこで、以下では着順予測に考慮すべき要素の一つである競走馬の「強さ」を統計モデリングにより推定します。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

使用データ

以下のようなレース結果のシミュレーションデータを使用します（1438頭、290レース分）*1。このデータは、後述の統計モデルから生成したものです。

> d <- read.csv(file='input/sim-data.csv',header = TRUE)
> d[1:20,c("RaceID","OoA","horseID")]
   RaceID OoA horseID
1       1   1     625
2       1   2     525
3       1   3     628
4       1   4     343
5       1   5     110
6       1   6     341
7       1   7     527
8       1   8       4
9       1   9     136
10      1  10       6
11      1  11       3
12      1  12     629
13      1  13     624
14      1  14     626
15      1  15     627
16      2   1     495
17      2   2     492
18      2   3     115
19      2   4     496
20      2   5     262

ここで、列名の意味は以下の通りです。

  • RaceID：レースのインデックス
  • OoA：着順
  • horseID：競走馬のインデックス

統計モデル

競馬の着順が決まるメカニズムを以下のように仮定します。

  1. 各競走馬には固有の「強さ」 \mu[n]が存在する（ nは競走馬のインデックス）。
  2. レースで発揮できる「強さ」（＝パフォーマンス）にはバラつき \sigma[n]が存在し、このパフォーマンスの大小でレースの着順が決まる。

パフォーマンスが平均 \mu[n]・標準偏差 \sigma[n]の正規分布に従うものとし、 これらの仮定を統計モデルで表現すると以下のようになります。

{\displaystyle performance[g,i,r ] \sim {\rm Normal} \left( \mu[HorseID [g,i,r ] ] ,\sigma[HorseID [g,i,r ] ] \right) , \\ (g=1,\ldots, G[r], \ \ i=1,\ldots, r, \ \ r = 7,\ldots,18) \\ performance[g,1,r ] \lt performance[g,2,r ] \lt \ldots \lt performance[g,r,r ] \\ (g=1,\ldots, G[r], \ \ i=1,\ldots, r, \ \ r = 7,\ldots,18) \\ \mu[n] \sim {\rm Normal} \left( 0, \sigma_{\mu} \right), \ \ \sigma[n] \sim {\rm Gamma}(10,10) \\ (n = 1,\ldots,N) }

ここで、 gはレースのインデックス、 iは各レースの着順（を反転したもの）のインデックス、 rは各レースの競走馬数、 G[r]は競走馬数 rのレース数、 Nは競走馬数を表します。また、 performanceは各レース・各競走馬のパフォーマンス、 HorseID [g,i,r ]は競走馬のインデックス（競走馬数 rのレース gにおいて i着だった競走馬）を表します。最終行は、競走馬の「強さ」 \mu[n]とそのバラつき \sigma[n]を一意に推定できるように事前分布を設定しています*2。上記の統計モデルに基づき、レース結果のデータ HorseIDから、 performance[g,i,r ], \ \mu[n], \ \sigma[n], \ \sigma_{\mu}を推定します。

上記の統計モデルのstanでの実装例は以下の通りです。なお、このコードは「StanとRでベイズ統計モデリング」10.2節を参考にしています（というよりほとんど同じです）。レース毎に競走馬数が変わるため、以下では愚直に競走馬数ごとに配列を作成しています。もっとスマートな書き方をご存じの方がいらっしゃれば、ご指摘いただけるとありがたいです。。。

data {
  int N;      // num of horses
  int G[18];  // num of races
  int<lower=1, upper=N> LW7[G[7],7];  // order of arrival of each race (7th:LW7[,1],6th:LW[,2],...1st:LW7[,7])
  int<lower=1, upper=N> LW8[G[8],8];  
  int<lower=1, upper=N> LW9[G[9],9];  
  int<lower=1, upper=N> LW10[G[10],10];  
  int<lower=1, upper=N> LW11[G[11],11];  
  int<lower=1, upper=N> LW12[G[12],12];  
  int<lower=1, upper=N> LW13[G[13],13];  
  int<lower=1, upper=N> LW14[G[14],14];  
  int<lower=1, upper=N> LW15[G[15],15];  
  int<lower=1, upper=N> LW16[G[16],16];  
  int<lower=1, upper=N> LW17[G[17],17];  
  int<lower=1, upper=N> LW18[G[18],18];  
}

parameters {
  ordered[7] performance7[G[7]];
  ordered[8] performance8[G[8]];
  ordered[9] performance9[G[9]];
  ordered[10] performance10[G[10]];
  ordered[11] performance11[G[11]];
  ordered[12] performance12[G[12]];
  ordered[13] performance13[G[13]];
  ordered[14] performance14[G[14]];
  ordered[15] performance15[G[15]];
  ordered[16] performance16[G[16]];
  ordered[17] performance17[G[17]];
  ordered[18] performance18[G[18]];
  
  vector[N] mu;
  real<lower=0> s_mu;
  vector<lower=0>[N] s_pf;
}

model {
  for (r in 2:18){
    for (g in 1:G[r]){
      for (i in 1:r){
        if (r==7)
          performance7[g,i] ~ normal(mu[LW7[g,i]], s_pf[LW7[g,i]]);
        else if (r==8)
          performance8[g,i] ~ normal(mu[LW8[g,i]], s_pf[LW8[g,i]]);
        else if (r==9)
          performance9[g,i] ~ normal(mu[LW9[g,i]], s_pf[LW9[g,i]]);
        else if (r==10)
          performance10[g,i] ~ normal(mu[LW10[g,i]], s_pf[LW10[g,i]]);
        else if (r==11)
          performance11[g,i] ~ normal(mu[LW11[g,i]], s_pf[LW11[g,i]]);
        else if (r==12)
          performance12[g,i] ~ normal(mu[LW12[g,i]], s_pf[LW12[g,i]]);
        else if (r==13)
          performance13[g,i] ~ normal(mu[LW13[g,i]], s_pf[LW13[g,i]]);
        else if (r==14)
          performance14[g,i] ~ normal(mu[LW14[g,i]], s_pf[LW14[g,i]]);
        else if (r==15)
          performance15[g,i] ~ normal(mu[LW15[g,i]], s_pf[LW15[g,i]]);
        else if (r==16)
          performance16[g,i] ~ normal(mu[LW16[g,i]], s_pf[LW16[g,i]]);
        else if (r==17)
          performance17[g,i] ~ normal(mu[LW17[g,i]], s_pf[LW17[g,i]]);
        else if (r==18)
          performance18[g,i] ~ normal(mu[LW18[g,i]], s_pf[LW18[g,i]]);
      }
    }
  }

  mu ~ normal(0, s_mu);
  s_pf ~ gamma(10, 10);
}

推定結果

推定された馬の「強さ」（の中央値）上位5頭と、実際の馬の「強さ」上位5頭を以下に示します。 horseIDが616と1324の馬は、実際の上位5頭と推定された上位5頭の両方に含まれていることから、ある程度適切に推定できています。 一方で、horseIDが1080, 955, 1246の馬は、実際の上位5頭には含まれているにも関わらず、推定された上位5頭には含まれていません。

> # 推定された上位5頭
> d_top5[,c("horseID","mu_p50")]                 
     horseID   mu_p50
956      956 2.352053
616      616 2.287300
1324    1324 2.197555
142      142 2.196556
118      118 2.171983
> # 実際の上位5頭
> d2 <- dplyr::distinct(d,horseID,.keep_all = TRUE)
> d2[order(-d2$mu)[1:5],c("horseID","mu")]
     horseID       mu
457      616 3.890667
1063    1324 3.179124
1042    1080 2.996332
1107     955 2.806721
1003    1246 2.773523

これは、horseIDが1080, 955, 1246の馬のデータが少ないことが原因と考えられます。実際、これらの馬の出走レース数は1レースしかありません（下表のn_race参照）。

> # 実際の上位5頭のレース数
> d %>% 
+   dplyr::select(horseID,RaceID,OoA) %>%
+   dplyr::filter(horseID == 616 | horseID == 1324 | horseID == 1080 | horseID == 955 | horseID == 1246) %>%
+   dplyr::group_by(horseID) %>%
+   dplyr::summarise(n_race = n())
# A tibble: 5 x 2
  horseID n_race
    <int>  <int>
1     616      4
2     955      1
3    1080      1
4    1246      1
5    1324      3

なお、horseIDが1080, 955, 1246の馬の「強さ」の90%ベイズ信用区間を算出すると、データ数が多い馬（horseID：616,1324）に比べて信用区間が広く、 データが少ないことに起因する推定値の不確かさを評価できていることが確認できます。

> # 実際の上位5頭の「強さ」の推定値
> d_qua %>% 
+   dplyr::filter(horseID == 616 | horseID == 1324 | horseID == 1080 | horseID == 955 | horseID == 1246) %>%
+   dplyr::mutate(BayesConf_90 = mu_p95 - mu_p05)
  horseID     mu_p05   mu_p50   mu_p95 BayesConf_90
1     616  1.2021932 2.287300 3.343293     2.141100
2     955 -0.4223376 1.045755 2.420237     2.842574
3    1080 -0.2466671 1.358909 2.843493     3.090160
4    1246 -0.3242929 1.123751 2.506408     2.830701
5    1324  1.0066156 2.197555 3.477100     2.470484

まとめ

競走馬の「強さ」によって着順が決まる統計モデルを作成し、シミュレーションデータで競走馬の「強さ」とその不確かさを推定できることを確認しました。ただし、今回作成したモデルは、騎手・脚質・枠順等、実際のレースで着順に影響を与える要素を考慮していないので、まだまだ改善の余地があると考えています。

更新履歴

  • 2019年11月25日
  • 2019年12月03日
    • 数式中の一部変数名を変更しました（ Results→ HorseID）

*1:実際にはJRA-VAN (https://jra-van.jp/) で提供されている実績データも分析していますが、この分析結果を公開するとJRA-VANの利用規約に引っかかってしまうので、シミュレーションデータで説明します。

*2:これらの事前分布を設定せずに推定すると、パラメータ \mu[n]と \sigma[n]が一意に決まらず、MCMCが収束しなくなります。