walkingsdbgの日記

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

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

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

競馬

はじめに

前回記事（競走馬の「強さ」をベイズ統計モデリングで推定する - walkingsdbgの日記）の続きです。 競走馬の「強さ」をレース結果から推定する方法を前回記事で説明しましたが、着順に影響を与える他の要素として、 騎手の「強さ」があります。例えば、11月24日（日）に開催されたジャパンカップ（ＧⅠ）では、オイシン・マーフィー騎手の的確な進路選択がスワーヴリチャードを1着に導きました。そこで今回は、競走馬に加えて騎手の「強さ」も同時に推定する方法を説明します。

使用データ

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

> d <- read.csv(file='input/sim-data.csv',header = TRUE)
> d[1:20,c("RaceID","OoA","horseID","jockeyID")]
   RaceID OoA horseID jockeyID
1       1   1     343       34
2       1   2     629       53
3       1   3     525        9
4       1   4     625       36
5       1   5     136       20
6       1   6       3       83
7       1   7       4        4
8       1   8     624       41
9       1   9     628       61
10      1  10     110       73
11      1  11     527       54
12      1  12     341        3
13      1  13       6       15
14      1  14     626       91
15      1  15     627       22
16      2   1     493       27
17      2   2     262       58
18      2   3     172       67
19      2   4     495       13
20      2   5     497       52

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

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

統計モデル

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

  1. 競走馬、騎手にはそれぞれ固有の「強さ」 \mu_h[n], \mu_j[m]が存在する（ n, mはそれぞれ競走馬、騎手のインデックス）。
  2. 競走馬、騎手がレースで発揮できる「強さ」（＝パフォーマンス）にはそれぞれバラつき \sigma_h[n], \sigma_j[m]が存在し、競走馬と騎手のパフォーマンスの総和が大きい順にレースの着順が決まる。

競走馬、騎手のパフォーマンスはそれぞれ平均 \mu_h[n]・標準偏差 \sigma_h[m]の正規分布、平均 \mu_j[n]・標準偏差 \sigma_j[m]の正規分布に従うものとし、これらの仮定を統計モデルで表現すると以下のようになります。

{\displaystyle \begin{align} pf_h [g,i,r ] \sim {\rm Normal} \left( \mu_h[HorseID [g,i,r ] ] ,\sigma_h[HorseID [g,i,r ] ] \right) , \tag{1} \\ (g =1,\ldots, G[r], \ \ i=1,\ldots, r, \ \ r = 7,\ldots,18) \\ pf_j [g,i,r ] \sim {\rm Normal} \left( \mu_j[JockeyID [g,i,r ] ] ,\sigma_j[JockeyID [g,i,r ] ] \right) , \tag{2} \\ (g=1,\ldots, G[r], \ \ i=1,\ldots, r, \ \ r = 7,\ldots,18) \\ pf[g,i,r ] = pf_h [g,i,r ] + pf_j [g,i,r ], \tag{3}\\ (g=1,\ldots, G[r], \ \ i=1,\ldots, r, \ \ r = 7,\ldots,18) \\ pf[g,1,r ] \lt pf[g,2,r ] \lt \ldots \lt pf[g,r,r ], \tag{4}\\ (g=1,\ldots, G[r], \ \ i=1,\ldots, r, \ \ r = 7,\ldots,18) \\ \mu_h[n] \sim {\rm Normal} \left( 0, \sigma_{\mu_h} \right), \ \ \sigma_h[n] \sim {\rm Gamma}(10,10), \tag{5}\\ (n = 1,\ldots,N) \\ \mu_j[m] \sim {\rm Normal} \left( 0, \sigma_{\mu_j} \right), \ \ \sigma_j[m] \sim {\rm Gamma}(10,10), \tag{6} \\ (m = 1,\ldots,M) \end{align} }

ここで、上式で新しく登場した記号の意味は以下の通りです。

  • g：レースのインデックス
  • i：各レースの着順のインデックス
  • r：各レースの競走馬数
  • G[r]：競走馬数 rのレース数
  • N：競走馬数
  • pf_h：競走馬のパフォーマンス
  • pf_j：騎手のパフォーマンス
  • pf：競走馬と騎手のパフォーマンスの合計
  • HorseID [g,i,r ]：競走馬のインデックス（競走馬数 rのレース gにおいて r-i+1着だった競走馬）
  • JockeyID [g,i,r ]：騎手のインデックス（競走馬数 rのレース gにおいて r-i+1着だった騎手）

式(5)(6)では、競走馬と騎手の「強さ」 \mu_h, \mu_jとそのバラつき \sigma_h,\sigma_jを一意に推定できるように事前分布を設定しています。 上記の統計モデルに基づき、レース結果のデータ HorseID,JockeyIDから、パラメータ pf, \ \mu_h, \ \mu_j, \ \sigma_h, \ \sigma_j, \ \sigma_{\mu_h}, \ \sigma_{\mu_j}を推定します。

上記の統計モデルのstanでの実装例は以下の通りです。なお、以下のコードでは式(1)~(3)をまとめて次式のように表記しています*1。

{\displaystyle pf[g,i,r ] \sim {\rm Normal} \left( \mu_h[HorseID [g,i,r ] ] + \mu_j[JockeyID [g,i,r ] ] , \sqrt{\sigma^2_h[HorseID [g,i,r ] ] + \sigma^2_j[JockeyID [g,i,r ] ]} \right) }

data {
  int N;      // num of horses
  int M;      // num of jockeys
  int G[18];  // num of races
  int HorseID[sum(G),18,18];
  int JockeyID[sum(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_h;
  vector[M] mu_j;
  real<lower=0> s_mu_h;
  real<lower=0> s_mu_j;
  vector<lower=0>[N] s_pf_h;
  vector<lower=0>[M] s_pf_j;
}

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_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==8)
          performance8[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==9)
          performance9[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==10)
          performance10[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==11)
          performance11[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==12)
          performance12[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==13)
          performance13[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==14)
          performance14[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==15)
          performance15[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==16)
          performance16[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==17)
          performance17[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
        else if (r==18)
          performance18[g,i] ~ normal(mu_h[HorseID[g,i,r]]+mu_j[JockeyID[g,i,r]], sqrt( s_pf_h[HorseID[g,i,r]]^2 + s_pf_j[JockeyID[g,i,r]]^2 ) );
      }
    }
  }

  mu_h ~ normal(0, s_mu_h);
  mu_j ~ normal(0, s_mu_j);
  s_pf_h ~ gamma(10, 10);
  s_pf_j ~ gamma(10, 10);
}

推定結果

馬・騎手の「強さ」の推定値（事後分布の中央値）の妥当性を、馬・騎手の「強さ」の真値\mu_h, \ \mu_jと比較することで評価します。

馬の「強さ」の推定結果

馬の「強さ」の推定値と真値の上位5頭をそれぞれ以下に示します。 「強さ」の真値の中で1位であるhorseID：616の馬は、推定値の中でも2位にランクインしていることから、ある程度適切に推定できています。 一方で、「強さ」の真値の中で2～5位であるhorseID：1324, 1080, 955, 1246の馬は、推定値の上位5頭には含まれていません。

> # 推定された上位5頭
> d_top5_h[,c("horseID","mu_p50")]                 
    horseID   mu_p50
455     455 2.271581
616     616 2.266670
343     343 2.183429
321     321 2.137451
83       83 2.047325
> # 実際の上位5頭
> d2 <- dplyr::distinct(d,horseID,.keep_all = TRUE)
> d2[order(-d2$mu_horse)[1:5],c("horseID","mu_horse")]
     horseID mu_horse
457      616 3.890667
1063    1324 3.179124
1042    1080 2.996332
1107     955 2.806721
1006    1246 2.773523

これは、前回記事と同様に、horseID：1324, 1080, 955, 1246の馬のデータが少ないことが原因と考えられます。 実際、これらの馬の出走レース数（下表のn_race）はhorseID：616の馬に比べて少ないです*2。

> # 実際の上位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

騎手の「強さ」の推定結果

騎手の「強さ」の推定値と真値の上位5頭をそれぞれ以下に示します。 「強さ」の真値の中で1～3位であるjockeyID：94, 48, 29の騎手は、推定値の中でも上位3位にランクインしていることから、ある程度適切に推定できています。 一方で、「強さ」の真値の中で4～5位であるjockeyID：136, 44の騎手は、推定値の上位5人には含まれていません。

> # 推定された上位5人
> d_top5_j[,c("jockeyID","mu_p50")]                 
    jockeyID   mu_p50
29        29 2.214286
48        48 1.808931
94        94 1.546535
103      103 1.269752
1          1 1.255741
> # 実際の上位5人
> d2 <- dplyr::distinct(d,jockeyID,.keep_all = TRUE)
> d2[order(-d2$mu_jockey)[1:5],c("jockeyID","mu_jockey")]
    jockeyID mu_jockey
30        94  2.779812
67        48  2.362584
23        29  2.319888
130      136  2.051941
74        44  2.044647

これも、馬の「強さ」の推定結果の節で述べたように、出走レース数の少なさが原因と考えています。 実際、これらの騎手の出走レース数（下表のn_race）はhorseID：94, 48, 29の騎手に比べて少ないです。

> # 実際の上位5人のレース数
> d %>% 
+   dplyr::select(jockeyID,RaceID,OoA) %>%
+   dplyr::filter(jockeyID == 94 | jockeyID == 48 | jockeyID == 29 | jockeyID == 136 | jockeyID == 44) %>%
+   dplyr::group_by(jockeyID) %>%
+   dplyr::summarise(n_race = n())
# A tibble: 5 x 2
  jockeyID n_race
     <int>  <int>
1       29    106
2       44      2
3       48     23
4       94     11
5      136      4

まとめ

競走馬と騎手の「強さ」によって着順が決まる統計モデルを作成し、シミュレーションデータからそれらの「強さ」をある程度推定できることを確認しました。 なお、本記事で使用したデータ・コードは https://github.com/walkingsdbg/horseracing に置いています。

次回の記事は、統計モデルの改善は一旦置いておいて、統計モデルを使った着順予測方法と、予測した着順に基づく最適な馬券の買い方について書こうと考えています。

*1:パラメータ pfは、正規分布に従う確率変数 pf_h, pf_jの和で表されるので、正規分布に従います。また、 pf_h, pf_jは独立であるため、 pfの平均と分散は、 pf_h, pf_jの平均の和 \mu_h + \mu_jと分散の和 \sigma^2_h + \sigma^2_jでそれぞれ表されます。

*2:ただし、horseID：1324の馬はレース数3とそこそこデータがあるので、「強さ」の推定値の中では7位に位置しています。