テニス選手の強さをベイズモデリングで分析する

Stan python

はじめに

stanとRでベイズ統計モデリングを読み終えました。 本の10章の将棋の強さと勝負ムラを推定するという内容が面白かったため、 自身の興味のあるテニスで同様のことを実施してみました。

探したら同じテーマでやっている記事があったのですが、 元体育会庭球部所属の身としてテニスは外せなかったため、入力するデータ等少し条件を変えて実施しました。
ベイズモデリングで男子プロテニスの強さを分析してみた – 戦略コンサルで働くデータサイエンティストのブログ

データの取得・確認

最新の動向が知りたいので2018年のATPの試合データをネットの海から拾いました。
https://github.com/JeffSackmann/tennis_atp

kaggleにもデータがあるのですが、今は2017年のデータまでしかないようです。
https://www.kaggle.com/gmadevs/atp-matches-dataset

まずデータを集計し、一定以上の試合に出場している 勝率が高い選手のtop10を確認してみます。

#データ読み込み
data=pd.read_csv("../input/atp_matches_2018.csv",encoding='latin-1')

#勝利数のカウント
win_count=data.groupby("winner_name")["winner_id"].count().reset_index()
win_count.columns=["name","win_count"]

#敗北数のカウント
lose_count=data.groupby("loser_name")["loser_id"].count().reset_index()
lose_count.columns=["name","lose_count"]

#データの結合
counts=pd.merge(win_count,lose_count,on="name",how="outer")
counts=counts.fillna(0)

#勝率の算出
counts["win_rate"]=counts["win_count"]/(counts["win_count"]+counts["lose_count"])
counts=counts.sort_values(by="win_rate",ascending=False)
counts[(counts["win_count"]+counts["lose_count"])>=20].head(10)

f:id:rmizutaa:20190115084249p:plain

ナダルの勝率が9割を超えてるのは驚異的ですね… この中には入ってませんが、錦織選手は30勝15敗で勝率は0.667でした。 また試合毎に得られるポイントは異なるため、ランキングとも順序が異なっていますね。 (ランキングは4大大会で2勝しているジョコビッチが1位)

出場試合数の少ない選手を入れると結果が収束しないことが予測されるため、 今回は出場試合数が20以上の99選手のデータのみを使用します。

モデル

今回は以下のような式に基づいたモデルを使用します。 この式は「stanとRでベイズ統計モデリングの」p189に記載されているものと同一になります。

performance [ g , 1 ] \sim \operatorname { Normal } \left( \mu [Loser] , \sigma _ { p f } [Loser] \right) , \quad g = 1 \ldots G

{ per formance } [ g , 2 ] \sim \operatorname { Normal } \left( \mu [ W i n n e r ] , \sigma _ { p f } [ {Winner } ] \right) , \quad g = 1 \ldots G

performance [ g , 1 ] \lt performance [ g , 2 ] , \quad g = 1 \ldots G

\mu [ n ] \sim \operatorname { Normal } \left( 0 , \sigma _ { \mu } \right) , \quad n = 1 \ldots N

\sigma _ { p f } [ n ] \sim \operatorname { Gamma } ( 10,10 ) , \quad n = 1 \ldots N

ここで、Nは選手の人数、nは選手のIDで、各選手の強さは平均0、標準偏差\sigma_{\mu}の正規分布に従うと仮定しています。

また、標準偏差の弱情報事前分布としてはガンマ分布を用いています。 ガンマ分布の行は設定しなくてもプログラムを実行することはでき、その場合は十分に広い一様事前分布が設定されるようなのですが、その場合はstanが収束しなくなりました。 今回のモデルだと、条件が勝者が敗者よりも値が大きいという部分しかなく絶対的な値の指標が存在しないため、事前分布の形をある程度指定することで値の取りうる範囲を限定させ、収束しやすくなるようです。

ちなみにGamma(10,10)の分布は以下のような形になります。

f:id:rmizutaa:20190115085107p:plain

モデルに入力するために以下のように前処理を行います。

#試合出場数が20以上の選手を抽出
player_attend=pd.concat([data["winner_name"],data["loser_name"]],axis=0).value_counts().to_frame()
player_attend=player_attend[player_attend[0]>=20]
players_use=list(player_attend.index)

#playerにidを付与
id_player=pd.DataFrame()
ids =[i+1 for i in range(len(players_use))]
id_player["winner_id2"]=ids #元のファイルとかぶるので2を付与
id_player["winner_name"]=players_use
data=pd.merge(data,id_player,on="winner_name",how="left")
id_player.columns=["loser_id2","loser_name"]
data=pd.merge(data,id_player,on="loser_name",how="left")

#今回使用する選手のみに
data=data[(data["winner_name"].isin(players_use)) & (data["loser_name"].isin(players_use))]

#stan入力用
LW=data[["loser_id2","winner_id2"]].astype(int)

stanのコードは以下のようになります。

model="""
data {
  int N;  // num of players
  int G;  // num of games
  int<lower=1, upper=N> LW[G,2];  // loser and winner of each game
}

parameters {
  ordered[2] performance[G];
  vector[N] mu;
  real<lower=0> s_mu;
  vector<lower=0>[N] s_pf;
}

model {
  for (g in 1:G)
    for (i in 1:2)
      performance[g,i] ~ normal(mu[LW[g,i]], s_pf[LW[g,i]]);

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


stan_data = {'N': len(players_use), 'G': len(LW),'LW': LW}
sm = pystan.StanModel(model_code=model)
fit = sm.sampling(data=stan_data, iter=1000, chains=3)

結果

stanで推定したパラメータを確認します。

result = fit.extract()

結果閲覧用のデータフレーム作成
winners=data[["winner_id2","winner_name"]]
winners.columns=["player_id","player_name"]
losers=data[["loser_id2","loser_name"]]
losers.columns=["player_id","player_name"]

players=pd.concat([winners,losers],axis=0)
players=players.drop_duplicates()
players=players.sort_values(by="player_id")

#playerの強さの平均値と勝負ムラ
players["mu"]=np.median(result["mu"],axis=0)
players["s_pf"]=np.median(result["s_pf"],axis=0)

強さの平均値の上位10人は以下になりました。 この表で、muは強さの平均、s_pfは強さの分散を示します。

players.sort_values(by="mu",ascending=False).head(10)

f:id:rmizutaa:20190115085932p:plain

勝率の上位で出した順序とほとんど同じですね。 この中ではティエムが勝負ムラが少し大きいようですね。

次に、強さの分散（勝負ムラ）の上位10人は以下になりました。

players.sort_values(by="s_pf",ascending=False).head(10)

f:id:rmizutaa:20190115085945p:plain

ここにはあまり上位選手の名前が出てこない結果となりました。 おそらく今回抽出した99人のうち、TOP層とそれ以外だと順序による重みが異なってくることが原因かと思います。 （1位の選手が20位の選手に負けるのは稀だが、61位の選手が80位の選手に負けることはよくある。）

気になったので、ナダルと錦織の強さの分布も図示してみました。

from scipy.stats import norm
X = np.arange(-3,6,0.1)
Y = norm.pdf(X,1.580996,0.721014)
plt.plot(X,Y,color="r",label="Nadal")
Y = norm.pdf(X,0.586693,0.869935)
plt.plot(X,Y,color="g",label="Nishikori")
plt.legend()
plt.show()

f:id:rmizutaa:20190115085219p:plain

これを見る限りだと、ナダルの調子が平均以上の場合は錦織がナダルに勝利するのは難しそうです。 2018年はナダルvs錦織の対戦は一度しかなく、結果はナダルの勝利になっていますが、通算の錦織とナダルの対戦成績は2勝10敗なので、感覚的にはそんなにずれていないかな？という気がします。（対戦の相性は全く考慮していないのであくまで一つの側面からの評価になりますが。）

まとめ

2018年のATPのデータを用いてテニスの強さと、勝負ムラについての分析を行いました。このような分析が何に使えるかを考えてみましたが、例えば団体戦のオーダの組み方を考えるのに役立つかもしれません。総合力は少し高くパフォーマンスのムラが少ない選手Aと総合力はやや劣るがパフォーマンスのムラが大きい選手Bがいたとして、対戦相手毎に予測される勝率を計算してベストなオーダを組むということが感覚的ではなく数値的に行えるとかが考えられそうですね。

使用したコードは以下
https://github.com/rmizuta3/tennis_analysys