はじめに
前回記事(競走馬の「強さ」をベイズ統計モデリングで推定する - walkingsdbgの日記)の続きです。
競走馬の「強さ」をレース結果から推定する方法を前回記事で説明しましたが、着順に影響を与える他の要素として、
騎手の「強さ」があります。例えば、11月24日(日)に開催されたジャパンカップ(GⅠ)では、オイシン・マーフィー騎手の的確な進路選択がスワーヴリチャードを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:騎手のインデックス
統計モデル
競馬の着順が決まるメカニズムを以下のように仮定します。
- 競走馬、騎手にはそれぞれ固有の「強さ」が存在する(はそれぞれ競走馬、騎手のインデックス)。
- 競走馬、騎手がレースで発揮できる「強さ」(=パフォーマンス)にはそれぞれバラつきが存在し、競走馬と騎手のパフォーマンスの総和が大きい順にレースの着順が決まる。
競走馬、騎手のパフォーマンスはそれぞれ平均・標準偏差の正規分布、平均・標準偏差の正規分布に従うものとし、これらの仮定を統計モデルで表現すると以下のようになります。
ここで、上式で新しく登場した記号の意味は以下の通りです。
- :レースのインデックス
- :各レースの着順のインデックス
- :各レースの競走馬数
- :競走馬数のレース数
- :競走馬数
- :競走馬のパフォーマンス
- :騎手のパフォーマンス
- :競走馬と騎手のパフォーマンスの合計
- :競走馬のインデックス(競走馬数のレースにおいて着だった競走馬)
- :騎手のインデックス(競走馬数のレースにおいて着だった騎手)
式(5)(6)では、競走馬と騎手の「強さ」とそのバラつきを一意に推定できるように事前分布を設定しています。
上記の統計モデルに基づき、レース結果のデータから、パラメータを推定します。
上記の統計モデルのstanでの実装例は以下の通りです。なお、以下のコードでは式(1)~(3)をまとめて次式のように表記しています*1。
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);
}
推定結果
馬・騎手の「強さ」の推定値(事後分布の中央値)の妥当性を、馬・騎手の「強さ」の真値と比較することで評価します。
馬の「強さ」の推定結果
馬の「強さ」の推定値と真値の上位5頭をそれぞれ以下に示します。
「強さ」の真値の中で1位であるhorseID:616の馬は、推定値の中でも2位にランクインしていることから、ある程度適切に推定できています。
一方で、「強さ」の真値の中で2~5位であるhorseID:1324, 1080, 955, 1246の馬は、推定値の上位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
>
> 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。
>
> 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())
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人には含まれていません。
>
> 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
>
> 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の騎手に比べて少ないです。
>
> 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())
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 に置いています。
次回の記事は、統計モデルの改善は一旦置いておいて、統計モデルを使った着順予測方法と、予測した着順に基づく最適な馬券の買い方について書こうと考えています。