2024-05-10

長期投資における投資比率の最適化(2/2)

はじめに

　前回の記事では、期待対数効用のもとで最適な投資比率が満たすべき方程式を導出しました。しかし、この方程式を具体的に解くには、投資対象の金融資産を選び、さらに、それらの収益率の確率分布を求める必要があります。そこで本記事では、金融資産の選定と確率分布の推定を行い、投資比率の最適解を導出します。さらに、推定された収益率の確率分布を用いて資産運用のシミュレーションを行い、様々な投資比率における結果を比較します。なお、登場する数式・記号の定義は前回の記事に従います。

投資対象の金融資産の選定

　適当なインデックスファンド１つを投資対象として選びます。この理由は以下の通りです。

資産運用が簡単になるから。
(とある強い仮定の下で)ボラティリティの意味でリスクをほぼ最小化するから。
NISAを利用できるものが多いから。

　理由１について、投資する金融資産を1種類に絞ることで、データの取得、収益率の確率分布のモデリング、そして資産の売買の手間が減らせます。あまり投資に時間と労力を割きたくないと考える自分にとっては魅力的です。
　理由２の仮定とは、「全ての投資家が、金融資産の収益率の平均・分散・共分散について同じ情報を持ち、また、これらの情報に基づくとある共通のリスク回避的投資方法に従う」というものです。この仮定の下で、収益率のボラティリティを最小化する投資方法は、無リスク資産と市場平均ポートフォリオの組合せであることが知られています（１ファンド定理）*1。市場平均ポートフォリオを個人投資家が購入するのは手間ですが、インデックスファンドが市場平均ポートフォリオを近似するため、これに投資することがボラティリティをほぼ最小化すると考えられます。注意点として、１ファンド定理の仮定は実際の市場で厳密に成立しません。しかし、この仮定にある程度納得感があるのと、理由１で述べたように投資の負担が非常に小さくなる（インデックスファンドへの投資だけ考えればよい）ので、この仮定を受け入れて投資を行うことにします。
　理由３は、(決められた上限額まで)運用益が非課税であるため、その分だけ資産を増やせることを意味します。本来は運用益に対して20%の税金がかかるので、多くのインデックスファンドに対してNISAを利用して投資できるのは魅力的です。

金融資産の収益率の確率分布のモデリング

　以下のインデックスファンドの月次収益率のデータを用いて、1期先の収益率の確率分布 $p(R_{t+1} | R_{1:t})$ を推定します。
emaxis.am.mufg.jp
　時系列データをモデリングする際には、主にトレンド・自己相関・周期性に注意する必要がありますが、これらをグラフからざっくり判断します。

月次収益率の時系列

　顕著なトレンドは見当たりません。

月次収益率の自己相関（コレログラム）

　自己相関も小さいです。また、12期前の自己相関が小さいことから、年単位の周期性もなさそうです。

上のグラフから、月次収益率 $R_t$ に明確なトレンド・自己相関・周期性は見当たらないことがわかります。したがって、月次収益率 $R_t$ は各時刻で独立かつ同一の分布に従うと仮定します。
　月次収益率 $R_t$ の確率分布をモデリングするため、このヒストグラムを見てみます。

上図から示唆されることは以下の通りです。

単峰性である。
左に裾が長い。

これらの特徴を表現できる確率分布として、以下のGumbel分布を用います。

$R_t \sim f \left( \left. R_t \right| \mu, \beta \right) = \frac{1}{\beta}\exp\left\{-\frac{R_t - \mu}{\beta} \exp\left( -\frac{R_t - \mu}{\beta} \right) \right\}$

この分布は2つの未知パラメータ $\mu,\beta$ を持つため、これらを最尤法で推定します。ただし、Gumbel分布は「右に」裾が長いため、収益率の符号を反転させたデータに対してパラメータを推定します。その結果、以下の推定値が得られました。

$\hat{\mu} = -0.037, \hat{\beta} = 0.040$

この推定値を用いて、収益率のヒストグラムとGumbel分布（をy軸に関して反転させたもの）を重ねて書いたものが下図です。収益率の分布の性質が表現できているように見えます。

　ただし、限られたデータから推定した $\hat{\mu}, \hat{\beta}$ はばらつきを持っています。さらに、将来もこの分布に従うかは不明なので、保守的な推定値を代わりに使いたいです。そこで、収益率の最頻値を表すパラメータ $\mu$ については、最尤推定値の代わりに、その近似的な信頼区間の上限 ${\hat{\mu}}_{ub}$ を用います。

${\hat{\mu}}_{ub} = \hat{\mu} + \frac{z_{{\alpha}/{2}}}{nI_{1}\left( \mu \right)}$

ここで、 $z_{{\alpha}/{2}}$ は標準正規分布の両側 $\alpha$ 点、 $n$ はデータ数、 $I_{1}\left( \mu \right)$ はデータ1点あたりのフィッシャー情報量を表します。このフィッシャー情報量を厳密に計算することは困難なので、以下のように近次します。

${ \displaystyle \begin{eqnarray} I_{1}\left( \mu \right) &=& E \left[ \frac{{\partial}^2}{\partial {\mu}^2}\left( \log f \left( \left. R_t \right| \mu, \beta \right) \right) \right] \\ &=& \frac{1}{{\beta}^2}E \left[ \exp\left( -\frac{R_t - \mu}{\beta} \right) \right] \\ &\approx& \frac{1}{{\hat{\beta}}^2}E \left[ \exp\left( -\frac{R_t - \hat{\mu}}{\hat{\beta}} \right) \right] \\ &\approx& \frac{1}{{\hat{\beta}}^2} \frac{1}{N_\mathrm{MC}}\sum_{i=1}^{N_\mathrm{MC}} \exp\left( -\frac{R_{t,i} - \hat{\mu}}{\hat{\beta}} \right) \end{eqnarray} }$

ここで、2行目から3行目ではパラメータ $\mu, \beta$ を最尤推定値 $\hat{\mu}, \hat{\beta}$ に置き換え、また、3行目から4行目ではモンテカルロ法による期待値の近似を行っています。なお、 $N_\mathrm{MC}$ はモンテカルロ法のサンプルサイズを表します。 $N_\mathrm{MC}=100000$ として上式を計算すると、 ${\hat{\mu}}_{ub}=-0.028$ と推定されました。これを用いて、各時刻の収益率 $R_t$ は独立に以下の分布に従うと仮定します。

$R_t \sim \frac{1}{\hat{\beta}}\exp\left\{-\frac{R_t - {\hat{\mu}}_{ub}}{\hat{\beta}} \exp\left( -\frac{R_t - {\hat{\mu}}_{ub}}{\hat{\beta}} \right) \right\} \tag{1}$

各期の投資比率の最適解の導出

　前節で得られた収益率の確率分布を用いて、前回の記事で求めた以下の方程式を解くことで、インデックスファンドへの投資比率の最適解を求めます。ただし、この方程式の右辺は解析的に計算できないので、モンテカルロ法で近似します。

$E \left[ \left. \frac{R_{t+1}-R_{f,t}}{\left( 1+R_{f,t}+C_t+\omega^{*} \left( R_{t+1}-R_{f,t} \right) \right)} \right| S_{t} \right]=0$

この左辺を様々な $\omega$ に対して計算したのが下図の黒の曲線です。縦軸が左辺の値、横軸が投資比率 $\omega$ を表します。また、赤の直線が右辺を表しており、これと黒の曲線の交点が投資比率の最適解 $\omega^{*}$ を表します。

上図より、最適解は150~160％であり、100%を大きく上回っていることが分かります。これは、借金をしながらインデックスファンドに投資することを意味するので、非常にリスクが高い戦略であることが分かります。前回の記事で、期待対数効用の最大化は(期待値の最大化より)リスク回避的と述べましたが、それでもなおリスクが大きすぎることが分かりました。

資産運用のシミュレーションによる評価

　前節では、期待対数効用のもとでの最適解ですら現実的でないことを示しました。そこで本節では、投資比率を適当な範囲に制限してシミュレーションを行い、損失リスクと収益のバランスを見ることで、自分にとってベストな投資割合を探します。シミュレーションの設定は以下の通りです。

最初の資産は100万円とする。
投資とは別に定期的な収入があり、毎月5万円ずつ資産が増加する。
投資対象は上記のインデックスファンドであり、その月次収益率は式(1)に従う。
月に一度、一定の投資比率となるようリバランスを行う。
運用期間は35年とする。

上記の設定のもとで、投資比率を10%、30%、50%とした場合の最終的な資産の分布を以下に示します。なお、赤の直線は、投資比率が0%である(=貯金のみを行う)場合の最終的な資産を表します。

上図より、投資割合を高めるほど、最終的な資産の取りうる金額の範囲が大きくなることが分かります。リスク回避を重視する場合は、左側の裾を見て性能を評価するのが良いと思います。例えば、老後に必要な貯金額を2000万円と仮定すると、最終的な資産の分布の最小値がこれを上回るような投資比率（約10％）を選ぶのが良さそうです。

まとめ

　本記事では、投資比率の最適解を求めるため、インデックスファンド１種類を投資対象として選び、その月次収益率の確率分布を推定しました。しかし、その分布に基づいて期待対数効用を最大化する最適解を計算した結果、借金をしながら投資を行う非常に高リスクな投資方法であることが判明しました。そこで、代わりに様々な固定比率で資産運用のシミュレーションを行い、最終的な資産の分布に基づいて、リスクの許容度に見合った投資比率を検討しました。

*1:詳細は以下の書籍をご覧ください。

金融工学入門第2版

作者:デービッド G.ルーエンバーガー
日経BP マーケティング(日本経済新聞出版

Amazon

2024-03-13

長期投資における投資比率の最適化(1/2)

はじめに
問題設定
再帰式の導出
動的計画法による最適解の導出
まとめ
参考文献

はじめに

　新NISAをきっかけに投資を始めようと思い立ったのですが、「どの金融商品にどれくらいの比率で投資すればよいのか」が分からないので、投資を最適化問題として定式化し、その最適解を求めることでこの比率を決めます。

問題設定

　投資比率の決定を最適化問題として表すには、まず「何を最大化するのか」を決める必要があります。これを目的関数と呼びます。投資なので「最終的な資産の期待値」を目的関数とするのが自然に感じますが、これは破産リスクの高い投資に繋がります。例えば、以下のギャンブルを考えます。

所持金10万円
コイントスをして、表が出れば賭け金が10倍に、裏が出れば賭け金が1/10になる
表が出る確率と裏が出る確率は1/2で等しい
所持金の何割を投資するのが最適か？

このギャンブルへの投資比率をαとし、資産の期待値を最大化するようにαを決めてみます。

表が出た場合、所持金は10(1-α)+100α=10+90α
裏が出た場合、所持金は10(1-α)+α=10-9α
期待値は1/2×(10+90α)+1/2×(10-9α)=10+(81/2)α

上式より、αが大きいほど期待値も大きくなるので、α=1、すなわち、所持金全額をこのギャンブルに賭けることが最適解となります。したがって、資産の期待値を最大化すると、破産確率が高い投資比率を選んでしまうことが分かります。
　そこで、最終的な資産の期待値の代わりに、最終的な資産の「対数」の期待値を最大化します。これを期待対数効用と呼びます。期待対数効用を選んだ理由は以下の2点です。

期待値をそのまま最大化するよりもリスク回避的な投資比率を選択できる。
長期投資の最適化問題の計算が簡単になる。

1つ目の理由を、上記のギャンブルの例を使って説明します。資産の対数の期待値は1/2×log(10+90α)+1/2×log (10-9α)と書けるので、(計算過程は省略しますが)最適解はα=1/2となります。期待値の最大化と比較して、期待対数効用の最大化は投資割合を控えめにすることが分かります。2つめの理由については、後述の最適化問題を解く中で説明します。
　期待対数効用を目的関数として用いると、投資比率を決めるために解くべき最適化問題は以下のように表すことができます。

$\mathrm{max}_{\omega_{t:T-1}} E _t \left[ \mathrm{log} \left(W_{T} \right) \right]$

ここで、上式の記号・変数の意味は以下の通りです。

$E _t \left[ \cdot \right]$ :時刻tまでの「情報」で条件付けられた期待値(「情報」の具体的な中身は後述します)
$\omega_{t:T-1}$ :各時刻 $k$ における金融資産への投資比率 $\omega_{k}$ を時刻tから $T-1$ までまとめた集合
$W_T$ : 終端時刻 $T$ における資産

再帰式の導出

　前節で定めた最適化問題を具体化して解きやすくするため、以下の仮定を置きます。

時刻t-1における資産 $W_{t-1}$ が与えられたとき、時刻tにおける資産 $W_{t}$ は次式で表すことができる。

${ \displaystyle \begin{eqnarray} W_{t}&=&(1+R_{p,t}+C_{t-1})W_{t-1} \\ &=&(1+\tilde{R}_{p,t})W_{t-1} \end{eqnarray} }$

- $R_{p,t}$ :時刻tにおけるポートフォリオの収益率
- $C_{t-1}$ :時刻tにおける定期的な収入等による資産の増加割合(簡単のため、一時刻前のt-1の時点で既知とする)
- 表記を簡単にするため、2行目では $\tilde{R}_{p,t}=R_{p,t}+C_{t-1}$ とした。

時刻tにおけるポートフォリオの収益率 $R_{p,t}$ は次式で表される。

$R_{p,t}=R_{f,t-1}+\sum_{i=1}^{n} \omega_{i,t}\left( R_{i,t}-R_{f,t-1} \right)$

- $R_{f,t-1}$ :時刻tにおける無リスク資産の収益率(簡単のため、一時刻前のt-1の時点で既知とする)
- $\omega_{i,t}$ :時刻tにおける金融資産iへの投資比率(i=1,...,n)。 $\omega_{i,t} \geq 0$ を満たす。
- $R_{i,t}$ :時刻tにおける金融資産iの収益率(i=1,...,n)。これは、収益率に関連する変数 $X_ {t}$ に依存する。

これまで登場した確率変数の関係は、下図のグラフィカルモデルで表される(簡単のため、図中では金融資産を指す添字iを省略している)。

上記の仮定のもとで、時刻tにおける情報 $S_{t} = \left[ X_{t}, R_{t}, W_{t}, C_{t}, R_{f,t} \right]^{\mathrm{T}}$ が与えられているとき、目的関数は以下のように書くことができます。

${ \displaystyle \begin{eqnarray} \mathrm{max}_{\omega_{t:T-1}} E \left[ \left. \mathrm{log} \left(W_{T} \right) \right| S_t \right] &=& \mathrm{max}_{\omega_{t:T-1}} E \left[ \left. \mathrm{log}\left(W_{t} \right) + \sum_{k=t}^{T-1} \mathrm{log} \left(1+\tilde{R}_{p,k+1} \right) \right| S_t \right] \\ &=& \mathrm{max}_{\omega_{t:T-1}} E \left[ \left. E \left[ \left. \mathrm{log}\left(W_{t} \right) + \sum_{k=t}^{T-1} \mathrm{log} \left(1+\tilde{R}_{p,k+1} \right) \right| S_{t+1},S_t \right] \right| S_t \right] \\ &=& \mathrm{max}_{\omega_{t:T-1}} E \left[ \left. E \left[ \left. \mathrm{log}\left(W_{t} \right) + \sum_{k=t}^{T-1} \mathrm{log} \left(1+\tilde{R}_{p,k+1} \right) \right| S_{t+1} \right] \right| S_t \right] \\ &=& \mathrm{log}\left(W_{t} \right) + \mathrm{max}_{\omega_{t}} E \left[ \left. \mathrm{log} \left(1+\tilde{R}_{p,t+1} \right) + \mathrm{max}_{\omega_{t+1:T-1}} E \left[ \left. \sum_{k=t+1}^{T-1} \mathrm{log} \left(1+\tilde{R}_{p,k+1} \right) \right| S_{t+1} \right] \right| S_t \right] \end{eqnarray} }$

ここで、

1行目の式変形では、 $W_{t}=(1+\tilde{R}_{p,t})W_{t-1}$ を代入し、さらに対数関数の性質を利用することで、目的関数を、各期の収益率の対数の和に分解しました。
2行目の式変形では、期待値の繰返し公式( $E \left[X \right] = E \left[ E \left[ \left. X \right| Y \right] \right]$ )を用いました。
3行目の式変形では、 $S_{t}$ のマルコフ性( $p\left( \left. S_{t+1} \right| S_{0},\ldots, S_{t}\right) = p\left( \left. S_{t+1} \right| S_{t}\right)$ )を用いました*1。
4行目の式変形では、 $\mathrm{log}\left(W_{t} \right)$ が定数であることと、目的関数の単調性を用いました。

上式より、 $J \left(S_t \right) = \mathrm{max}_{\omega_{t:T-1}} E \left[ \left. \mathrm{log} \left(W_{T} \right) \right| S_t \right] - \mathrm{log} \left(W_{t} \right) = \mathrm{max}_{\omega_{t:T-1}} E \left[ \left. \sum_{k=t}^{T-1} \mathrm{log} \left(1+\tilde{R}_{p,k+1} \right) \right| S_t \right]$ と定義すると、以下のような再帰式を得ることができます。

$J \left(S_t \right) = \mathrm{max}_{\omega_{t}} E \left[ \left. \mathrm{log} \left(1+\tilde{R}_{p,t+1} \right) + J \left(S_{t+1} \right) \right| S_t \right]$

動的計画法による最適解の導出

　前節で導出した再帰式を、以下のように終端時刻から後ろ向きに解くことで、全時刻の最適解を求めることができます。このような帰納的な解法を総称して動的計画法と呼びます。

$J \left(S_{T-1} \right)$ を解いて最適解 $\omega_{T-1}^{*}$ を求める。
手順1で求めた最適解 $\omega_{T-1}^{*}$ を $J \left(S_{T-2} \right)$ に代入し、一時刻前の最適解 $\omega_{T-2}^{*}$ を求める。
手順2と同様の手続きを繰り返すことで、全時刻の最適解 $\omega_{t:T-1}^{*}$ を求める。

　手順1から具体的に計算していきます。 $J \left(S_{T-1} \right)$ は次のように表されます。

$J \left(S_{T-1} \right) = \mathrm{max}_{\omega_{T-1}} E \left[ \left. \mathrm{log} \left(1+\tilde{R}_{p,T} \right) \right| S_{T-1} \right]$

上式より目的関数は凹関数なので、最適性の一次の必要条件を満たす $\omega_{T-1}$ が大域的最適解に一致します。よって、目的関数を $\omega_{i,T-1}$ について微分することで、最適解 $\omega_{T-1}^{*}$ が満たす方程式は以下のようになります。

$E \left[ \left. \frac{R_{i,T}-R_{f,T-1}}{\left( 1+\tilde{R}_{p,T}\left(\omega_{T-1}^{*}\right) \right)} \right| S_{T-1} \right]=0 \quad (i=1,\ldots,n)$

この方程式を解くには金融資産の収益率の確率分布 $p\left( \left. R_{T} \right| X_{T-1} \right)$ が必要なのですが、ここでは解が得られたとして次の手順に進みます。
　手順2では、手順1で求めた解を用いて次式を解きます。

${ \displaystyle \begin{eqnarray} J\left(S_{T-2}\right) &=& E \left[ \mathrm{log} \left(1+\tilde{R}_{p,T-1}\right) + \left. J \left({S}_{T-1} \right)\right| S_{T-2} \right] \\ &=& E \left[\left. \mathrm{log} \left(1+\tilde{R}_{p,T-1}\right) + \mathrm{log} \left( 1+\tilde{R}_{p,T}\left(\omega_{T-1}^{*}\right) \right) \right| S_{T-2} \right] \end{eqnarray} }$

第二項は $\omega_{T-2}$ に依存しないため、第一項を $\omega_{i,T-2}$ について微分することで、最適解 $\omega_{T-2}^{*}$ が満たす方程式は以下のようになります。

$E \left[ \left. \frac{R_{i,T-1}-R_{f,T-2}}{\left( 1+\tilde{R}_{p,T-1}\left(\omega_{T-2}^{*}\right) \right)} \right| S_{T-2} \right]=0 \quad (i=1,\ldots,n)$

　手順1,2を見てみると、最適な投資比率 $\omega_{t}^{*}$ を決めるには、一期先の収益率の期待対数効用 $E \left[\left. \mathrm{log} \left(1+\tilde{R}_{p,t+1}\right) \right| S_{t} \right]$ を最大化すればよいことが分かります。したがって、任意の時刻tにおける最適解は次式の解に一致します。

$E \left[ \left. \frac{R_{i,t+1}-R_{f,t}}{\left( 1+\tilde{R}_{p,t+1}\left(\omega_{t}^{*}\right) \right)} \right| S_{t} \right]=0 \quad (i=1,\ldots,n)$

まとめ

資産の期待対数効用の最大化問題として投資を定式化し、最適な投資比率が満たす条件を動的計画法により求めました。この条件から具体的な投資比率を求めるには金融資産の一期先の収益率の確率分布が必要なのですが、これについては別の記事で書きます。

参考文献

デービッド・G・ルーエンバーガー．金融工学入門．第2版，日本経済新聞出版，2015

金融工学入門第2版

作者:デービッド G.ルーエンバーガー
日経BP マーケティング(日本経済新聞出版

Amazon

　期待値を最大化する戦略が破産のリスクに直面することを、18章の回転盤のギャンブルの例で分かりやすく示しています。本記事の例題は本書を参考にしました。

梶野洸ほか. 強化学習から信頼できる意思決定へ. サイエンス社，2024．

強化学習から信頼できる意思決定へ (AI/データサイエンスライブラリ“基礎から応用へ” 5)

Amazon

　再帰方程式の導出において重要な目的関数の性質(単調性・繰返し法則等)が3章1〜3節に分かりやすく書かれており、目的関数の設定および再帰式の導出に役立ちました。

安達智彦，池田昌幸．長期投資の理論と実践: パーソナル・ファイナンスと資産運用．東京大学出版会，2019．

長期投資の理論と実践: パーソナル・ファイナンスと資産運用

作者:安達智彦,池田昌幸
東京大学出版会

Amazon

　本記事では資産の最大化のみを考慮しましたが、本書では、資産だけでなく消費の最適化も含めたより一般的な問題の解を3章1〜3節で導出しています。

C. M. ビショップ著，元田浩ほか訳．パターン認識と機械学習（下）．丸善出版，2017．

パターン認識と機械学習下 (ベイズ理論による統計的予測)

作者:C.M. ビショップ
丸善出版

Amazon

　グラフィカルモデルに基づく条件付き独立性のチェック方法が本書の8.2節に示されており、(本記事では割愛しましたが)マルコフ性が成り立つことを確認する際に参考にしました。

*1:この性質のおかげで、期待値計算に必要な確率分布が投資期間の長さに依らずシンプルになります。本記事ではマルコフ性の確認を割愛します。

2021-08-12

新型コロナウイルス陽性者数の時間変化を日本地図に描く

はじめに

　データ可視化手法の勉強のため、新型コロナウイルス陽性者数の時間変化を日本地図にヒートマップで描きます。本記事は備忘録的な意味合いが強いので、読みにくい点が多いかと思います。あらかじめご了承ください。

作図方法

　以下の手順に従ってグラフを作成します。

日本地図のデータを取得する
新型コロナウイルスの陽性者数のデータを取得する
日本地図と陽性者数のデータを結合する
日本地図に陽性者数のヒートマップを載せる
ヒートマップの時間変化をアニメーションで表現する

なお、グラフ作成に使用したパッケージを先に示しておきます。

library(jsonlite)
library(dplyr)
library(sf)
library(animation)
library(ggplot2)
library(lubridate)

日本地図のデータを取得する

　Rのパッケージ{NipponMap}のjpn.shpに日本地図のデータがあります。パッケージのインストール後、インストール先フォルダから以下のように読込みます。*1

d_map <- read_sf(system.file("shapes/jpn.shp", package = "NipponMap")[1],
                 crs = "+proj=longlat +datum=WGS84")

新型コロナウイルスの陽性者数のデータを取得する

　2020年4月末からの陽性者数のデータが以下のサイトにjson形式で公開されています。
https://opendata.corona.go.jp/api/Covid19JapanAll
ダウンロード後、以下のように読込みます。

d_covid_json <- fromJSON("data/Covid19JapanAll.json")

このデータには日次陽性者数が含まれていない（累計陽性者数（npatient）として保存されている）ため、以下のように累計陽性者数の差分を取ることで日次陽性者数を計算します*2。

d_covid <- d_covid_json$itemList %>%
  rename(pref = name_jp,
         npatients_acc = npatients) %>%
  mutate(date = as.Date(date),
         npatients_acc = as.numeric(npatients_acc)) %>%
  group_by(pref) %>%
  arrange(date) %>%
  mutate(npatients = c(min(npatients_acc),diff(npatients_acc))) %>%
  ungroup()

さらに、作図用のデータとして、新規陽性者数の週次平均を以下のように計算します。

d_covid <- d_covid %>%
  mutate(date_week = floor_date(date,"week")) %>%
  group_by(date_week,SP_ID,pref) %>%
  summarise(npatients = mean(npatients) + 10^(0)) %>% #後で対数表記を行うので、log(0)を避けるために微小な値を足しておく
  ungroup()

日本地図と陽性者数のデータを結合する

日本地図のデータと陽性者数のデータには都道府県の項目が共通しているので、これをキーにして結合します。ただし、日本地図のデータ内の都道府県はID表記である一方、陽性者数のデータ内の都道府県は漢字表記であるため、陽性者数のデータの都道府県をIDに変換してから結合します。

d_map$SP_ID <- as.numeric(d_map$SP_ID) #陽性者数のデータと結合するため、numeric型に揃えておく
d_covid_map <- d_covid %>%
    mutate(
         SP_ID = case_when(
           pref == "北海道" ~ 1,
           pref == "青森県" ~ 2,
           pref == "岩手県" ~ 3,
           pref == "宮城県" ~ 4,
           pref == "秋田県" ~ 5,
           pref == "山形県" ~ 6,
           pref == "福島県" ~ 7,
           pref == "茨城県" ~ 8,
           pref == "栃木県" ~ 9,
           pref == "群馬県" ~ 10,
           pref == "埼玉県" ~ 11,
           pref == "千葉県" ~ 12,
           pref == "東京都" ~ 13,
           pref == "神奈川県" ~ 14,
           pref == "新潟県" ~ 15,
           pref == "富山県" ~ 16,
           pref == "石川県" ~ 17,
           pref == "福井県" ~ 18,
           pref == "山梨県" ~ 19,
           pref == "長野県" ~ 20,
           pref == "岐阜県" ~ 21,
           pref == "静岡県" ~ 22,
           pref == "愛知県" ~ 23,
           pref == "三重県" ~ 24,
           pref == "滋賀県" ~ 25,
           pref == "京都府" ~ 26,
           pref == "大阪府" ~ 27,
           pref == "兵庫県" ~ 28,
           pref == "奈良県" ~ 29,
           pref == "和歌山県" ~ 30,
           pref == "鳥取県" ~ 31,
           pref == "島根県" ~ 32,
           pref == "岡山県" ~ 33,
           pref == "広島県" ~ 34,
           pref == "山口県" ~ 35,
           pref == "徳島県" ~ 36,
           pref == "香川県" ~ 37,
           pref == "愛媛県" ~ 38,
           pref == "高知県" ~ 39,
           pref == "福岡県" ~ 40,
           pref == "佐賀県" ~ 41,
           pref == "長崎県" ~ 42,
           pref == "熊本県" ~ 43,
           pref == "大分県" ~ 44,
           pref == "宮崎県" ~ 45,
           pref == "鹿児島県" ~ 46,
           pref == "沖縄県" ~ 47)) %>%
  left_join(.,d_map %>% select(SP_ID,geometry),key = "SP_ID")

日本地図に陽性者数のヒートマップを載せる

Rのパッケージ{ggplot2}のgeom_sf()を使えば、日本地図をプロットすることができます。ただし、入力データは予めdata.frameからsfオブジェクト（都道府県と属性情報がセットになった47行のオブジェクト）に変換する必要があります。そこで、パッケージ{sf}をインストールし*3、このパッケージに含まれる関数st_as_sf()を用いてsfオブジェクトに変換後、以下のように日本地図上に陽性者数のヒートマップを描きます。

#地図上の沖縄県の位置を右上に移動する
    d_covid_map$geometry[grep("沖縄県", d_covid_map$pref)] <-
        d_covid_map$geometry[grep("沖縄県", d_covid_map$pref)] + c(5, 17)
#沖縄県の移動先に補助線を引く
    layer_autoline_okinawa <- function(
                    x = c(129, 132.5, 138),
                    xend = c(132.5, 138, 138),
                    y = c(40, 40, 42),
                    yend = c(40, 42, 46),
                    size = ggplot2::.pt / 15
                  ){
                     ggplot2::annotate("segment",
                         x = x,
                         xend = xend,
                         y = y,
                         yend = yend,
                         size = .pt / 15
                  )}

#日本地図上に陽性者数のヒートマップを作成
    plotdate_s <- as.Date('2020-04-26') #適当に日付を設定
    plotdate_e <- plotdate_s + 6
    sf_plot <- st_as_sf(d_covid_map %>% filter(date_week == plotdate_s))
    p <- ggplot(sf_plot)
    p <- p + geom_sf(aes(fill = npatients))
    p <- p + layer_autoline_okinawa()
    p <- p + scale_fill_gradient(low="white", high="#E60012", trans="log10", limits=c(10^(0),10^4), name="陽性者数")
    p <- p + labs(title = paste0("都道府県別新規陽性者数(週次平均)\n",as.character(plotdate_s),"~",as.character(plotdate_e))) 
    p <- p + theme(
      panel.background = element_rect(fill = "midnightblue"),
      panel.grid  = element_blank(),
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank()
    )
    plot(p)

ここで、プログラム冒頭で、日本地図上の沖縄県の位置を南西から北西に変更しています。これは、沖縄県の位置が南西のままだと日本地図のサイズが大きくなり、本州が見づらくなるためです。なお、沖縄県の位置変更については、以下の記事・つぶやきを参考にさせていただきました。
rpubs.com

勉強になります！沖縄の位置をずらすだけであれば、
df$geometry[str_which("Okinawa", df$name)] <-
df$geometry[str_which("Okinawa", df$name)] + c(5, 17)
※ c(5, 17) は任意の数字
で書き換えることができると思います。ご参考になれば幸いです。
— Shoei (@shoei05) 2021年5月30日

ヒートマップの時間変化をアニメーションで表現する

前節で作成した日本地図上のヒートマップは、ある特定の時期における陽性者数を表すものです。このヒートマップを時期に応じて変化させるため、以下のような手順でアニメーションを作成します。

前節のヒートマップを、時期を変化させてN枚作成する。
N枚のヒートマップを、パッケージ{animation}を利用して1つのgifファイルにする

これは以下のプログラムで実現できます。

#GIFアニメーションのオプション設定
ani.options(loop = TRUE, outdir=getwd(), convert = 'program/ImageMagick-7.1.0-4-portable-Q16-x64/convert.exe')

#GIFアニメーションの作成
plotdate_list <- d_covid_map %>% arrange(date_week) %>% distinct(date_week) %>% pull(date_week)
saveGIF({
  for (n in 1:length(plotdate_list)){
    #各週の陽性者数のグラフを描く
    plotdate_s <- plotdate_list[n]
    plotdate_e <- plotdate_s + 6
    sf_plot <- st_as_sf(d_covid_map %>% filter(date_week == plotdate_s))
    p <- ggplot(sf_plot)
    p <- p + geom_sf(aes(fill = npatients))
    p <- p + layer_autoline_okinawa()
    p <- p + scale_fill_gradient(low="white", high="#E60012", trans="log10", limits=c(10^(0),10^4), name="陽性者数")
    p <- p + labs(title = paste0("都道府県別新規陽性者数(週次平均)\n",as.character(plotdate_s),"~",as.character(plotdate_e))) 
    p <- p + theme(
      panel.background = element_rect(fill = "midnightblue"),
      panel.grid  = element_blank(),
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank()
    )
    plot(p)
  }
}, movie.name="covid19_JapanPref.gif", interval=0.25)

このプログラムでは、saveGIF()というgifファイルを作成・保存する関数を使用し、引数として、各週の陽性者数（週次平均）のヒートマップをforループで作成しています。冒頭のani.option()ではgifの詳細設定（gifのループ有無・出力先フォルダ指定・gif作成に必要なexeファイルの読込）を行っています。ここで注意しないといけないのは、ani.option()で指定しているように、gif作成にはImageMagickというソフトが別途必要であることです。
imagemagick.org
このソフトには色々なバージョンがあってややこしいのですが、以下のサイトのアドバイスに従い、portable版をダウンロード・解凍しました。ImageMagickを利用するには、上記プログラムのように、ani.option()にてImageMagickのconvert.exeへのパスを指定すればOKです。
higuma.github.io

実行結果

　これまでの手順を実行すると、以下のようなgifアニメーションが作成できます。期間は2020年5月～2021年7月です。東京等から全国に感染が広がる様子と、緊急事態宣言による感染者数の増減がよく見えます。

f:id:walkingsdbg:20210812000422g:plain — 新型コロナウイルス陽性者数の時間変化を表すアニメーション

まとめ

　Rのパッケージ{NipponMap}・{sf}・{animation}とgif作成ソフトImageMagicを用いて、新型コロナウイルスの新規陽性者数（週次平均）の変化を2020年5月～2021年7月の期間で地図上に描きました。ある指標の地域別での時間変化を概観するには便利な方法だと感じたので、仕事でも機会があれば使ってみようと思います。

補足

本記事のプログラムは以下にまとめております。詳細に興味がある方はこちらをご覧ください。
github.com
　また、本記事の作成にあたって、以下の記事を参考にさせていただきました。
dinov.sakura.ne.jp
totech.hateblo.jp

*1:インストール後、library(NipponMap)のようにパッケージを読込む必要は無いことに注意してください。

*2:データの中には、なぜか累計陽性者数が時間とともに減少する（＝日次陽性者数を計算すると負になる）部分があるので、その日の陽性者数は0で置き換えています。

*3:パッケージ{sf}を読込む前に、最新版の{Rcpp}をインストールすることを推奨します。私の環境では1.0.7より古いバージョンのRcppを使うと、sfの読込時にエラーが出ました。

2021-07-31

統計的仮説検定によるウマ娘のトレーニング失敗率の妥当性検証

はじめに

　スマートフォン／PCアプリ「ウマ娘　プリティダービー」が配信されてからはや半年弱。多くのトレーナーがウマ娘の育成に力を注いでいると思います。ところで、ウマ娘のトレーニングの失敗率が、表記の数値より高いと感じたことはありませんか？（失敗率が数％でも普通に失敗したり...）この失敗率の検証について、既に他のトレーナーの方がYoutubeに動画を公開されたりしています。しかし、検証結果の信頼性を統計的に評価したものは私が見た限りではありませんでした。そこで、本記事では「トレーニング失敗率が表記の数値より高いか」を統計的仮説検定を用いて検証します。

f:id:walkingsdbg:20210728223153p:image:h480:w260 — 無残にもトレーニングで散るウマ娘

検証方法

統計的仮説検定とは

　統計的仮説検定とは、その名の通り、確率的な出来事に対する仮説をデータに基づいて検証する方法のことです。例えば、コイントスを10回行って10回とも表が出た場合、このコインはどこかおかしいと考える方が多いと思います。これは、以下のように推論した結果と考えることができます。

コインが普通のコイン（表・裏が出る確率がそれぞれ50%）と仮定する
コイントスを10回行った結果、10回とも表という結果が得られた
表が出る確率が50％と仮定すると、この結果は極めて起こりにくいため、この仮定が間違っている（コインがおかしい）と考える

このように、仮説の下で期待される結果（例：表裏が５回ずつ出る）と観測結果（表が10回出る）との差が、単に偶然によるものなのか、それとも仮説が誤っているのかを確率に基づいて判断する方法が統計的仮説検定です。

２種類の誤り

　上記の例題のように、統計的仮説検定では、観測結果に基づいて2種類の仮説のどちらかを採用します。初めに立てる仮説（例：コインが普通）を帰無仮説と呼び、帰無仮説と対立する仮説（例：コインがおかしい）を対立仮説と呼びます。統計的仮説検定には、どちらの仮説を採用するかによって２種類の誤りを犯すリスクがあります。１つは、帰無仮説が正しいのに対立仮説を採用してしまう誤りで、もう一つは、対立仮説が正しいのに帰無仮説を採用してしまう誤りです。１つ目の誤りを第１種の過誤、２つ目の誤りを第２種の過誤と呼びます。統計的仮説検定では、これらの誤りを犯す確率を閾値以下に抑えるように仮説の採用基準を定めることで、正しい仮説を高い確率で採用することが可能です。

		正しい仮説
		帰無仮説	対立仮説
採用する仮説	帰無仮説	〇	×（第２種の過誤）
採用する仮説	対立仮説	×（第１種の過誤）	〇

仮説の採用基準

　それでは、仮説の採用基準をどのように定めればよいでしょうか？ここからは、ウマ娘の失敗率検証に絞って話を進めます。コイントスの例と同様に考えると、例えば「失敗率10%のトレーニングを100回行い、20回以上トレーニングに失敗すれば、失敗率の表記がおかしい」のように定めることができます。しかし、第１種の過誤・第２種の過誤を閾値（例：10％）以下に抑えるために、トレーニングの回数・失敗数の閾値を何回にすればよいかは自明ではありません。そこで、以下では第１種の過誤・第２種の過誤の確率とトレーニング回数等との関係式を求め、その関係式から逆算して、トレーニング回数等を設定します。

※ここからは、第１種の過誤・第２種の過誤の確率計算のために技術的な話が多くなるため、結果のみに関心がある方は「データの収集方法」までスキップして差し支えありません。

第１種の過誤の確率

　本節では、第１種の過誤の確率を求めます。そのために、まずトレーニングの失敗回数が従う確率分布を設定します。トレーニング回数を $n$ 、失敗回数を $x$ 、失敗率を $p$ とするとき、失敗回数 $x$ は以下の二項分布に従うと仮定します。

$\displaystyle p(x) = nCxp^x(1-p)^{n-x} \tag{1} \label{1}$

ここで、 $p(x)$ は失敗回数が $x$ になる確率を表します。二項分布は、コイントスのように二種類の結果が得られる確率的な出来事を表す分布です。トレーニングも成功／失敗の二種類の結果が現れるため、この分布を仮定します。
次に、失敗率 $p$ に対して、帰無仮説 $H_0$ と対立仮説 $H_1$ を以下のように設定します。

帰無仮説 $H_0$ : $p = p_0$
対立仮説 $H_1$ : $p > p_0$

ここで、 $p_0$ はゲーム内で表示される失敗率（以下、表示失敗率）を表します。つまり、帰無仮説 $H_0$ は真の失敗率が表示通りであること、対立仮説 $H_1$ は真の失敗率が表示失敗率より大きい（＝不正がある）ことを表します。さらに、仮説の採用基準を以下のように設定します。
「 $u_0 = \frac{\hat{p}-p_0}{\sqrt{p_0(1-p_0)/n}} \ \ (\hat{p} = \frac{x}{n})$ とし、 $u_0 \le c$ のとき帰無仮説 $H_0$ を採用し、 $u_0 > c$ のとき対立仮説 $H_1$ を採用する」
ここで、 $c$ は分析者が設定する閾値です。採用基準には、前節のように失敗回数 $x$ を用いた方が分かりやすいのですが、 $u_0$ の方がトレーニング回数等の設定を行いやすいので、 $u_0$ を使います。
$u_0$ の直感的なイメージは次の通りです。 $u_0$ の分子を見ると分かるように、実験から得られる失敗率の推定値 $\hat{p}$ と表示失敗率 $p_0$ の差を計算していて、この差が大きいほど、真の失敗率が表示失敗率より大きい（＝不正がある）可能性が高いと考えることができます。
　準備が長くなりましたが、以上の前提から第１種の過誤の確率を求めます。帰無仮説が成り立つとき、トレーニング回数 $n$ が十分大きければ、中心極限定理より、 $u_0$ が従う確率分布は標準正規分布で近似することができます。これを用いると、第１種の過誤が起きる確率は

${ \displaystyle \begin{eqnarray} \mathrm{Pr} \left(u_0 > c \right) = 1 - \mathrm{Pr} \left(u_0 \le c \right) \approx 1 - \Phi(c) \tag{2} \label{1stError} \end{eqnarray} }$

と表すことが可能です。ここで、 $\Phi(x)$ は標準正規分布の累積分布関数を表します。上式より、第１種の過誤が生じる確率を $\alpha$ に抑えたいとき、標準正規分布の上側 $\alpha$ %点を $z_{\alpha}$ とすると、閾値 $c$ を
$c = z_{\alpha} \tag{3} \label{c}$
と設定すればよいことが分かります。

第２種の過誤の確率

　次に、第２種の過誤の確率を計算します。第２種の過誤は対立仮説が正しいのに帰無仮説を採用してしまう誤りなので、対立仮説（ $p > p_0$ ）が成り立つと仮定して次式を計算すればよいです。

${ \displaystyle \begin{eqnarray} \mathrm{Pr} \left(u_0 \le c \right) &=& \mathrm{Pr} \left(u_0 \le z_{\alpha} \right) \\ &=& \mathrm{Pr} \left(\frac{\hat{p}-p_0}{\sqrt{p_0(1-p_0)/n}} \le z_{\alpha} \right) \\ &=& \mathrm{Pr} \left( \left\{ \frac{\hat{p}-p}{\sqrt{p(1-p)/n}} + \frac{p-p_0}{\sqrt{p(1-p)/n}} \right\} \sqrt{ \frac{p(1-p)}{p_0(1-p_0)} } \le z_{\alpha} \right) \\ &=& \mathrm{Pr} \left( u \le z_{\alpha}A-B \right) \\ &\approx& \Phi(z_{\alpha}A-B) \tag{4} \label{2ndError} \end{eqnarray} }$

ここで、1行目では前節の結果（ $c = z_{\alpha}$ ）を適用しました。また、4行目の $u, A, B$ は以下のように定義しました。さらに5行目では、トレーニング回数 $n$ が十分大きければ、中心極限定理より $u$ の従う確率分布が標準正規分布で近似できることを利用しました。
$u = \frac{\hat{p}-p}{\sqrt{p(1-p)/n}}, \ A = \sqrt{ \frac{p_0(1-p_0)}{p(1-p)} }, \ B = \frac{p-p_0}{\sqrt{p(1-p)/n}} \tag{5} \label{5}$
式(\ref{2ndError})を見ると、第２種の過誤の確率はトレーニング回数 $n$ と真の失敗率 $p$ に依存することが分かります。第２種の過誤の確率と $n,p$ との関係を視覚化するため、様々な $n,p$ に対して第２種の過誤の確率を以下に示します。ここで、作図のために $\alpha=0.1, \ p_0=0.05$ としています。

f:id:walkingsdbg:20210730231955p:plain — α=0.1,p0=0.05の場合の第２種の過誤の確率。実線,破線,点線はそれぞれn=10,100,1000に対応する。

上図より、 $n$ が大きいほど第２種の過誤の確率は小さくなり、また、真の失敗率 $p$ が表示失敗率 $p_0$ に近いほど第２種の過誤の確率は大きくなることが分かります。したがって、 $p_0$ の近くで第２種の過誤の確率を小さくするためには $n$ を大きくする必要があり、その分データ取得が大変になることが分かります。一方、今は「真の失敗率 $p$ が表示失敗率 $p_0$ より不正と言える程大きいか」に関心があるため、 $p_0$ の近くで第２種の過誤の確率を小さくする必要性は低いです。そこで、本分析では、真の失敗率 $p$ が表示失敗率 $p_0$ より $\Delta_0$ 以上大きい場合に、第２種の過誤の確率を $\beta$ 以下に抑えられるように $n$ を求めます。上図より、 $p$ が $p_0$ に最も近い場合、つまり $p=p_0+\Delta_0$ の場合に第２種の過誤の確率が最も大きいので、この時の確率が $\beta$ になるように $n$ を計算します。式(\ref{2ndError})より、この条件は次式で表すことができます。
$z_{\alpha}A(p_0+\Delta_0)-B(p_0+\Delta_0,n) = z_{1-\beta} \tag{6} \label{6}$
ここで、 $A, \ B$ がそれぞれ $p,n$ に依存することを示すために $A(p), \ B(p,n)$ と表記しました。上式を $n$ について解くと、トレーニング回数 $n$ は次式で表すことができます。
$n = \frac{(p_0+\Delta_0)(1-(p_0+\Delta_0))}{{\Delta_0}^2} \left( z_{\alpha}A(p_0+\Delta_0) + z_{\beta} \right)^2 \tag{7} \label{n}$

各種条件の設定

　ここまでの検討を整理します。仮説の採用基準を以下のように設定し、第１種の過誤の確率、第２種の過誤の確率をそれぞれ $\alpha, \beta$ に抑えるような $c, n$ はそれぞれ式(\ref{c})、式(\ref{n})で表されることを示しました。

帰無仮説 $H_0$ : $p = p_0$
対立仮説 $H_1$ : $p > p_0$
仮説の採用基準:「 $u_0 = \frac{\hat{p}-p_0}{\sqrt{p_0(1-p_0)/n}} \ \ (\hat{p} = \frac{x}{n})$ とし、 $u_0 \le c$ のとき帰無仮説 $H_0$ を採用し、 $u_0 > c$ のとき対立仮説 $H_1$ を採用する」

本節では、 $\alpha, \beta, \Delta_0, p_0$ の値を設定することで $c, \ n$ の値を決めます。もちろん、 $\alpha, \beta$ を小さくするほど分析結果の信頼性は高くなり、また、 $\Delta_0$ を小さくするほど $p_0$ との僅かな差を検出しやすくなります。しかし、その分 $n$ が大きくなって実験が大変になるので、今回は $\alpha=0.1, \beta=0.1, \Delta_0=0.05$ と設定します。また、次節のデータ収集の都合上、 $p_0=0.05$ に設定します。このとき $c, \ n$ は、式(\ref{c})、式(\ref{n})より次式で与えられます。
$c = z_{0.1} = 1.296, \ \ n = 176.3... < 177 \tag{8} \label{candn}$
したがって、表示失敗率5%のもとでトレーニングを177回行えば、表示失敗率の妥当性に関する検証を行うことができます。次節ではこのトレーニングデータの収集方法について説明します。

データの収集方法

トレーニング失敗率の検証を行うにあたり、育成対象をスーパークリークに絞ります。なぜなら、スーパークリークは育成期間中に失敗率を一定値に固定することができるため、実験が行いやすいからです。通常、トレーニングの失敗率は体力によって変動するため、ある失敗率に対してデータを収集するのが大変です。しかし、スーパークリークは、クラシック級3月前半～クラシック級10月前半の間「小さなほころび」というバッドステータスが付与され、（体力が十分にあっても）常に失敗率が5%になります。そこで、この期間を利用して、失敗率を5%に固定しつつ177回分のトレーニングデータを集めました。*1

f:id:walkingsdbg:20210723113851p:image:h480:w260 — 「小さなほころび」状態のスーパークリーク。体力が十分にあっても失敗率が常に5%になる

検証結果

表示失敗率5%のもとで177回トレーニングを行った結果、失敗回数は7回、仮説の採用基準で用いる $u_0$ は $-0.638$ でした。これは閾値 $z_{0.1}=1.296$ より小さいため、帰無仮説（真の失敗率が表示失敗率と同じ）を採用します。帰無仮説を採用する際には第２種の過誤を犯すリスクがありますが、真の失敗率と表記失敗率との差が $\Delta_0=0.05$ 以上の場合にこの過誤を犯す確率を10%以下に抑えているので、この結果から「トレーニング失敗率は概ね表記通り（＋5%以上の乖離はない）」と言えます。
　

まとめ

第一種の過誤・第二種の過誤の確率を10%以下に抑えた上で統計的仮説検定を行った結果、「トレーニング失敗率は概ね表記通り（＋5%以上の乖離はない）」という結論が得られました。「失敗率が表記の数値より高い！」とこれまで感じていたのは、単に失敗したときのショックが大きく、その時のことを強く覚えているからだったのでしょう。。。

補足

本分析で使用したデータ・プログラムは以下にアップしていますので、興味がある方はご覧ください。
github.com

*1:ただし、効率的にデータを集めるため、「小さなほころび」が無い状態でも失敗率が5%のトレーニングがあれば、そのデータも分析対象に含めています

2020-04-18

統計モデルに基づく着順予測と、分枝限定法による馬券組合せ最適化についての検討

はじめに

本記事では、前回記事（競走馬と騎手の「強さ」をベイズ統計モデリングで推定する - walkingsdbgの日記）で紹介した統計モデルから予測した着順に基づき、最適な馬券の組合せを求める方法について説明します。

はじめに
統計モデルに基づく着順予測
馬券組合せ最適化
- 問題設定
- 分枝限定法による最適解の求め方
  - 分枝限定法とは
  - 上界の求め方
レースデータを使用した検証
まとめ

統計モデルに基づく着順予測

前回記事では、競馬の着順が決まるメカニズムを以下のように仮定し、統計モデルを構築しました。

競走馬と騎手にはそれぞれ固有の「強さ」が存在する。
競走馬と騎手がレースで発揮できる「強さ」（＝パフォーマンス）にはそれぞれバラつきが存在し、パフォーマンスは正規分布に従う。
競走馬と騎手のパフォーマンスの総和が大きい順にレースの着順が決まる。

このモデルに基づいて着順予測を行うには、モデルから生成した標本を利用します。例えば、ある馬が１着になる確率の予測方法は次の通りです。

競走馬と騎手のパフォーマンスを正規分布から生成し、そのパフォーマンスの大小でレース着順を決める、というシミュレーションを $N$ 回実施する。
$N$ 回のシミュレーションの内、ある馬 $i$ が１着になった回数を $n_i$ 回とすると、その馬が１着になる確率 $p_i$ は $p_i = n_i/N$ と求まる。

馬券組合せ最適化

競馬の着順が前回記事の統計モデルに従って決まる場合、どの馬が何％の確率で１着になるかを上記の方法で予測することが可能です。しかし、この予測結果だけでは、どのような馬券を購入すればよいかが分かりません（例えば、１着になる確率が最大の馬に１点張りすることもできますし、外れるリスクを避けて他の馬の馬券を合わせて購入することも可能です）。そこで以下では、着順予測結果から購入すべき馬券の組合せを求める方法について述べます。

問題設定

私は、以下の条件を満たす範囲で利益（＝払戻額－購入額）の期待値を最大化するような馬券の組合せを求めることを目指します。なお、ここでは単勝の馬券のみを買うことを前提とします*1。

予算は $Ｗ$ ×100円以内
損をする（利益が０円未満となる）確率が $\alpha$ 以下である

この利益の最大化問題を数式に直すと次のようになります。

評価関数

$\displaystyle \max_{z_{ij}} \sum_{i=1}^{N} \sum_{j=0}^{W} 100j \left(R_{i}p_{i}-1 \right)z_{ij} \tag{1} \label{1}$

制約条件

${ \displaystyle \begin{eqnarray} \sum_{i=1}^{N} \sum_{j=0}^{W} jz_{ij} \leq W \tag{2} \label{2} \\ \mathrm{Pr} \left(\sum_{i=1}^{N} \sum_{j=0}^{W} 100j(r_{i}-1)z_{ij} < 0 \right) \leq \alpha \tag{3} \label{3}\\ \sum_{j=0}^{W} z_{ij} = 1 \ \left(i = 1,2,\ldots,N \right) \tag{4} \label{4}\\ z_{ij} \in \{0,1\} \ \left(i = 1,2,\ldots,N, \ \ j = 0,1,\ldots, W \right) \tag{5} \label{5} \end{eqnarray} }$
ここで、式中に現れる記号の意味は以下の通りです。

$N$ ：馬券の総数
$R_i$ ：馬券 $i$ のオッズ
$p_i$ ：馬券 $i$ が当たる確率
$z_{ij}$ ：馬券 $i$ を $j \times 100$ 円で購入する／購入しないを指定する変数
$r_i$ ：馬券 $i$ の回収率（確率変数。当たれば $R_i$ 、外れれば $0$ ）

また、上式の意味は以下の通りです。

式\eqref{1}：利益の期待値[円]

（払戻額の期待値 $\sum_{i=1}^{N} \sum_{j=0}^{W} 100jz_{ij}R_{i}p_{i}$ [円]－購入額 $\sum_{i=1}^{N} \sum_{j=0}^{W} 100jz_{ij}$ [円]）

式\eqref{2}：購入額が予算を超えないようにする制約（条件1に対応）
式\eqref{3}：損をする確率がαを超えないようにする制約（条件2に対応）
式\eqref{4}：同じ馬券 $i$ を重複購入しないようにする制約
式\eqref{5}：購入する場合は1、購入しない場合は0。ただし $z_{i0} = 1$ の場合、馬券 $i$ を一切購入しないことを表す

分枝限定法による最適解の求め方

式\eqref{1}\eqref{2}\eqref{3}\eqref{4}\eqref{5}を満たす最適な馬券を求めるための一番シンプルな方法は、全ての馬券の組合せを探索することです。しかし、これは計算コストの観点で現実的ではありません。例えば、予算1万円で10頭立てのレースの馬券を購入する場合 $2^{100\times10} \approx 10^{301}$ 通りの組合せが存在するため*2、全ての馬券の組合せを探索することは困難です。そこで、私は分枝限定法という手法を用いて馬券の組合せを最適化しました。

分枝限定法とは

分枝限定法は、以下の分枝操作・限定操作という２種類の操作を繰り返すことで効率的に最適解を探索する方法です。

分枝操作：ある商品（ここでは馬券）を買う場合／買わない場合で場合分けすることにより、元の問題を分割する。
限定操作：分枝操作により分割された問題（部分問題）を評価して以下のことが判明した場合、部分問題の探索を打ち切る。
- 制約条件（ここでは式\eqref{2}\eqref{3}\eqref{4}\eqref{5}）を満たさない
- 部分問題の最適解の上界が暫定解（探索済の組合せの中での最適解）より小さい

例えば、予算1500円で、馬券Ａ・Ｂ・Ｃの中から最適な馬券の組合せを探索することを考えます。なお、簡単のためにそれぞれの馬券の値段は１０００円であるとします。また、分枝操作によりそれぞれの馬券を買う場合／買わない場合で場合分けし、下図のように①→②→③と探索を進めたとします。このとき、③では馬券Ａ・Ｂを購入して2000円を支払っているため、既に予算をオーバーしています。したがって、③以降のノードの探索を限定操作により打ち切ることが可能です。

f:id:walkingsdbg:20200416135923p:plain — 限定操作の一例（予算オーバーの場合）。番号のついたノードは探索順を表し、赤いノードは限定操作によって探索を打ち切られたノードを表す

もう一つ限定操作の例を挙げます。予算が3000円である点を除いて、問題設定は上記例と同様です。まず、下図のように①→②→③→④→⑤と探索を進め、暫定解としてノード⑤の3000円を得たとします。次に、ノード⑥において、部分問題の最適解の上界が後述の方法で2000円と判明したとします。このとき、上界（2000円）が暫定解（3000円）を下回っているため、⑥以降のノードを探索してもより良い解を得ることができないと分かります。したがって、限定操作により⑥以降のノードの探索を打ち切ることが可能です。

f:id:walkingsdbg:20200416143244p:plain — 限定操作の一例（上界が暫定解を下回る場合）。青のノードは暫定解を表す

このように、分枝限定法では限定操作を利用することで探索回数を減らすことが可能です。
なお、分枝限定法については以下の本・講義資料を参考にしました。より詳しい解説はこれらの資料をご覧ください。

組合せ最適化とアルゴリズム (インターネット時代の数学シリーズ 8)
- 作者:幹雄, 久保
- 発売日: 2000/12/10
- メディア: 単行本
http://www.dais.is.tohoku.ac.jp/~shioura/teaching/mp11/mp11-06.pdf

上界の求め方

分枝限定法では、元の問題の制約条件を緩めた緩和問題を解くことで最適解の上界を求めます。ここでは、非線形制約\eqref{3}\eqref{5}を線形制約に緩和し、線形計画法によって上界を計算します。まず、式\eqref{5}で定めた変数 $z_{ij}$ の定義域を{0,1}から実数に緩和します。次に、式\eqref{3}は「損をする（利益が０円未満となる）確率が $\alpha$ 以下である」という条件を表していますが、この条件を「全ての購入馬券が外れる確率が $\alpha$ 以下である」という条件に緩和します。緩和した条件は以下の式\eqref{6}のように $z_{i0}$ について線形であるため、式\eqref{1}\eqref{2}\eqref{4}\eqref{6}を線形計画法で解くことにより最適解の上界を求めることが可能になります。
${ \displaystyle \begin{eqnarray} (全ての購入馬券が外れる確率) &=& 1 - (少なくとも1つの購入馬券が当たる確率) \\ &=& 1 - \sum_{i=1}^N p_i\left(1-z_{i0}\right) \leq \alpha \tag{6} \label{6} \end{eqnarray} }$

レースデータを使用した検証

上述の着順予測と馬券組合せ最適化を用いて、過去のレースデータを用いたシミュレーションで回収率が100%を超えるか検証します。*3

使用データ

訓練データ：2017年1月～2018年11月のG1・G2・G3・オープンレース（計536レース）
テストデータ：2018年12月のG1レース（計5レース）

シミュレーション条件

$N = 2000$
$W = 100$
$\alpha = 0.2$
分枝限定法の探索回数上限*4： $10^5$ 回
馬券１口あたりの購入額*5： $500$ 円

結果・考察

テストデータ（計5レース）の購入額・払戻額・利益を以下に示します。なお、購入額が0円のレースは、上記の探索回数において暫定解が得られなかったレースを表します。

レース名	購入額[円]	払戻額[円]	利益[円]
チャンピオンズカップ(G1)	10,000	0	-10,000
阪神ジュベナイルフィリーズ(G1)	0	0	0
朝日杯フューチュリティステークス(G1)	10,000	11,500	1,500
有馬記念(G1)	0	0	0
ホープフルステークス	9,000	0	-9,000
合計	29,000	11,500	-17,500

上表より、総利益は-17,500円（回収率は約40%）という残念な結果になりました。。。この原因を調べるため、全ての購入馬券が外れたレース（チャンピオンズカップとホープフルステークス）について購入馬券と実際の着順を見ます。

チャンピオンズカップ

                 馬名 OoA   prob  odds weight
1    ルヴァンスレーヴ   1 0.1485   1.9     NA
2    ウェスタールンド   2 0.0630  32.1      5
3      サンライズソア   3 0.0980   8.2     15
4    アンジュデジール   4 0.0415  40.1      5
5    オメガパフューム   5 0.1350   9.2     15
6    サンライズノヴァ   6 0.1140  10.0     10
7        ノンコノユメ   7 0.0705  13.4     10
8              ミツバ   8 0.0500 139.5      5
9  ヒラボクラターシュ   9 0.0615  95.3      5
10     アスカノロマン  10 0.0010 296.0     NA
11   ケイティブレイブ  11 0.0710   5.0     20
12     センチュリオン  12 0.0590 251.8      5
13 インカンテーション  13 0.0665  75.5      5
14 アポロケンタッキー  14 0.0060 204.6     NA
15           パヴェル  15 0.0145  27.3     NA

ホープフルステークス

                 馬名 OoA   prob  odds weight
1  サートゥルナーリア   1 0.1965   1.8     NA
2  アドマイヤジャスタ   2 0.1075   6.2     15
3      ニシノデイジー   3 0.1990   6.6     15
4  コスモカレンドゥラ   4 0.0785  41.6      5
5  ブレイキングドーン   5 0.0975   9.9     10
6    ヴァンドギャルド   6 0.1215  10.4     10
7      ヒルノダカール   7 0.0285  68.7      5
8    キングリスティア   8 0.0175  24.4      5
9    ミッキーブラック   9 0.0625  35.8      5
10   マードレヴォイス  10 0.0165 436.3      5
11   ジャストアジゴロ  11 0.0255  70.3      5
12   ハクサンタイヨウ  12 0.0245 430.9      5
13       タニノドラマ  13 0.0245 171.0      5

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

OoA：着順
prob：統計モデルから推定した勝率
odds：単勝オッズ
weight：購入額（weight×100円）。NAは購入していないことを表す

購入馬券を見ると、1番人気の馬が含まれていないことが分かります。これは、1番人気の馬について、統計モデルから算出した勝率が低い割にオッズが高くないため、利益最大化の観点で良くない馬券と判断されたことが原因と考えられます。この戦略は、統計モデルから算出した勝率が妥当であれば良い戦略です。しかし、今回のシミュレーションでは馬券を購入した3レース中2レースで全滅という結果だったので、統計モデルによる勝率の推定精度が低いと考えています。参考情報ですが、以下ブログによるとG1レースで1番人気かつ単勝1倍台の馬の勝率は44.8%（本統計モデルでは15~19%程度）なので、現在使用している統計モデルは、1番人気の馬の勝率を過小評価していると考えています。
jra-van.jp

まとめ

統計モデルに基づく着順予測と分枝限定法を組み合せて、馬券組合せ最適化の方法を検討しました。過去のレースデータを用いた検証では回収率が約40％となったため、実戦投入には更なる改善が必要と分かりました。回収率が低い原因として統計モデルによる勝率の推定精度が低いことが考えられるため、今後は統計モデルの推定精度向上を行います。

*1:これは、制約条件\eqref{3}の計算がややこしくなる・問題の規模が大きくなることを避けるためです。利益最大化の観点では、単勝馬券に絞ることは選択肢を狭めるのであまり良くないです。

*2:1頭あたり、100円の単勝馬券を買う／買わない、200円の単勝馬券を買う／買わない、…10000円の単勝馬券を買う／買わない、という $2^{100}$ 通りの組合せがあるため、10頭立てのレースでは $2^{100\times10}$ 通りになります。

*3:本節の解析に使用したプログラムはhttps://github.com/walkingsdbg/horseracingにUPしています。ただし、使用データはJRA-VAN (https://jra-van.jp/)の有料サービスから取得したものであるため、非公開とさせていただきます。

*4:上限回数に達した場合、探索を途中で打ち切り、これまでの探索で得られた暫定解を最適解の近似解として出力します。

*5:1口100円の場合、問題の規模が大きすぎるためか探索が全然進まなかったため、1口500円にすることで問題の規模を小さくしています。

2019-12-15

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

競馬

はじめに

前回記事（競走馬の「強さ」をベイズ統計モデリングで推定する - 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：騎手のインデックス

統計モデル

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

競走馬、騎手にはそれぞれ固有の「強さ」 $\mu_h[n], \mu_j[m]$ が存在する（ $n, m$ はそれぞれ競走馬、騎手のインデックス）。
競走馬、騎手がレースで発揮できる「強さ」（＝パフォーマンス）にはそれぞれバラつき $\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位に位置しています。

2019-11-20

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

競馬

はじめに

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

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

作者:松浦健太郎
出版社/メーカー: 共立出版
発売日: 2016/10/25
メディア: 単行本

使用データ

以下のようなレース結果のシミュレーションデータを使用します（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：競走馬のインデックス

統計モデル

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

各競走馬には固有の「強さ」 $\mu[n]$ が存在する（ $n$ は競走馬のインデックス）。
レースで発揮できる「強さ」（＝パフォーマンス）にはバラつき $\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日
- 馬を表す英単語のスペルが誤っていたため修正しました（誤：horce→正：horse）
- コードをGitHubにUPしましたhttps://github.com/walkingsdbg/horseracing
2019年12月03日
- 数式中の一部変数名を変更しました（ $Results$ → $HorseID$ ）

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

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