スポンサーリンク

ベイズモデリングで男子プロテニスの強さを分析してみた

Python Stan データビジュアライゼーション データ分析

久しぶりにベイズモデルをやります。

今回は、Kaggleのデータセットで公開されている、男子プロテニスの試合結果データがありましたので、これを使って各プレイヤーの強さをモデル化してみます。

- Association of Tennis Professionals Matches　-ATP tournament results from 2000 to 2017-

そもそも私が趣味がテニスであり、プロテニスもよく観ていますので、スポーツアナリティクスをやる場合には、親近感のあるテニスのデータで何かしらやってみたいと思っていました。

なので、たびたびテニスの専門用語などを普通に書いてしまうかと思いますので、読みづらかったらすみません。

今回もコードはカーネルで投稿していたり。

kaggle kernel: https://www.kaggle.com/itoeiji/player-s-strength-analysis-using-bayesian

GitHubは以下。

GitHub: https://github.com/Gin04gh/datascience/tree/master/kaggle_dataset_association_of_tennis_professionals_matches

データ

データは、１試合ごとのツアー名、ツアーレベル、ツアー開催日、サーフェス、何回戦か、勝った選手の情報と得点情報、負けた選手の情報と得点情報などが格納されています。

Pandasに読み込んだコードが以下になります。

cols = [
    'tourney_id', # Id of Tournament
    'tourney_name', # Name of the Tournament
    'surface', # Surface of the Court (Hard, Clay, Grass)
    'draw_size', # Number of people in the tournament
    'tourney_level', # Level of the tournament (A=ATP Tour, D=Davis Cup, G=Grand Slam, M=Masters)
    'tourney_date', # Start date of tournament
    'match_num', # Match number
    'winner_id', # Id of winner
    'winner_seed', # Seed of winner
    'winner_entry', # How the winner entered the tournament
    'winner_name', # Name of winner
    'winner_hand', # Dominant hand of winner (L=Left, R=Right, U=Unknown?)
    'winner_ht', # Height in cm of winner
    'winner_ioc', # Country of winner
    'winner_age', # Age of winner
    'winner_rank', # Rank of winner
    'winner_rank_points', # Rank points of winner
    'loser_id',
    'loser_seed',
    'loser_entry',
    'loser_name',
    'loser_hand',
    'loser_ht',
    'loser_ioc',
    'loser_age',
    'loser_rank',
    'loser_rank_points',
    'score', # Score
    'best_of', # Best of X number of sets
    'round', # Round
    'minutes', # Match length in minutes
    'w_ace', # Number of aces for winner
    'w_df', # Number of double faults for winner
    'w_svpt', # Number of service points played by winner
    'w_1stIn', # Number of first serves in for winner
    'w_1stWon', # Number of first serve points won for winner
    'w_2ndWon', # Number of second serve points won for winner
    'w_SvGms', # Number of service games played by winner
    'w_bpSaved', # Number of break points saved by winner
    'w_bpFaced', # Number of break points faced by winner
    'l_ace',
    'l_df',
    'l_svpt',
    'l_1stIn',
    'l_1stWon',
    'l_2ndWon',
    'l_SvGms',
    'l_bpSaved',
    'l_bpFaced'
]
df_matches = pd.concat([
    pd.read_csv('./atp_matches_2000.csv', usecols=cols),
    pd.read_csv('./atp_matches_2001.csv', usecols=cols),
    pd.read_csv('./atp_matches_2002.csv', usecols=cols),
    pd.read_csv('./atp_matches_2003.csv', usecols=cols),
    pd.read_csv('./atp_matches_2004.csv', usecols=cols),
    pd.read_csv('./atp_matches_2005.csv', usecols=cols),
    pd.read_csv('./atp_matches_2006.csv', usecols=cols),
    pd.read_csv('./atp_matches_2007.csv', usecols=cols),
    pd.read_csv('./atp_matches_2008.csv', usecols=cols),
    pd.read_csv('./atp_matches_2009.csv', usecols=cols),
    pd.read_csv('./atp_matches_2010.csv', usecols=cols),
    pd.read_csv('./atp_matches_2011.csv', usecols=cols),
    pd.read_csv('./atp_matches_2012.csv', usecols=cols),
    pd.read_csv('./atp_matches_2013.csv', usecols=cols),
    pd.read_csv('./atp_matches_2014.csv', usecols=cols),
    pd.read_csv('./atp_matches_2015.csv', usecols=cols),
    pd.read_csv('./atp_matches_2016.csv', usecols=cols),
    pd.read_csv('./atp_matches_2017.csv', usecols=cols),
])
df_matches = df_matches.dropna(subset=['tourney_date'])
df_matches['year'] = df_matches['tourney_date'].apply(lambda x: int(str(x)[0:4]))
display(df_matches.head())
print(len(df_matches))

開催日 tourney_date が様々な形式で書かれていますので、開催年だけ別カラム year に変換しておきました。

様々な得点情報などもデータに含まれていますが、今回は、開催した年と、勝ったか負けたか、その時の勝ちプレイヤー、負けプレイヤーの情報だけで、各プレイヤーの強さを分析してみます。

テニス選手の強さモデル

下記の書籍が参考になりました。

ベイズ的な統計モデリングの書籍の中でも、具体的な問題設定とその実装例コードが多く記載されている書籍で、動かしながら学習できてオススメです。

StanとRでベイズ統計モデリング
松浦 健太郎、石田 基広

今回は上記書籍のp.189に記載されている、棋士の強さのモデリングで使われた以下の統計モデル

    \[performance[g,1]{\sim}Normal(\mu[Loser], \sigma_{pf}[Loser]), \hspace{1em} g=1{\ldots}G\]

    \[performance[g,2]{\sim}Normal(\mu[Winner], \sigma_{pf}[Winner]), \hspace{1em} g=1{\ldots}G\]

    \[performance[g,1] < performance[g,2], \hspace{1em} g=1{\ldots}G\]

    \[\mu[n]{\sim}Normal(0, \sigma_\mu), \hspace{1em} n=1{\ldots}N\]

    \[\sigma_{pf}[n]{\sim}Gamma(10, 10), \hspace{1em} n=1{\ldots}N\]

を使ってみます。

考え方としては、

  • Loser : 負けたプレイヤーのインデックス
  • Winner : 勝ったプレイヤーのインデックス

とします。

各選手の強さを \mu[n], 勝負ムラを \sigma_{pf}[n] として、1回の勝負で発揮する力（パフォーマンス）は、平均 \mu[n] , 標準偏差 \sigma_{pf}[n] の正規分布から生成されると考えます。

つまり、各選手にはそれぞれの強さの値を持っているのですが、試合ごとに見れば、調子が良くてパフォーマンスを発揮できたという時もあれば、不調で本来の力を出すことができなかった時などがあり、その変化の大きさを、力を発揮できたりできなかったりのムラとして考えます。

そして、勝負の結果はパフォーマンスの大小で決まる（大きかった方が勝つ）と考えます。

対象期間は、ひとまずは直近のみを使うとして、2015〜2017.02頃（直近期間）の戦績のデータとします。（2018年の直近までデータが欲しかったけど、今は無いようです）

対象選手は、この期間中の戦績データ数を確認しながら手動で選択します。

本当は世界ランキングTOP100位内くらいの選手すべてでやってみたいところですが、処理が膨大になりますし、結果の解釈も難しくなるので、適当に数名だけ選択することにします。

まずは分析対象の選手を選定します。

完全に個人的な興味と、データ数の兼ね合いで、以下のプレイヤーのリストを用意しました。

arr_target_player = np.array([
    'Roger Federer',
    'Rafael Nadal',
    'Novak Djokovic',
    'Andy Murray',
    'Stanislas Wawrinka',
    'Juan Martin Del Potro',
    'Milos Raonic',
    'Kei Nishikori',
    'Gael Monfils',
    'Tomas Berdych',
    'Jo Wilfried Tsonga',
    'David Ferrer',
    'Richard Gasquet',
    'Marin Cilic',
    'Grigor Dimitrov',
    'Dominic Thiem',
    'Nick Kyrgios',
    'Alexander Zverev'
])

df_tmp = df_matches[
    (df_matches['year'] >= 2015) &
    (df_matches['winner_name'].isin(arr_target_player)) &
    (df_matches['loser_name'].isin(arr_target_player))
]

arr_cnt = []
arr_rate = []
for player in arr_target_player:
    cnt_win = len(df_tmp[df_tmp['winner_name'] == player])
    cnt_lose = len(df_tmp[df_tmp['loser_name'] == player])
    arr_cnt.append(cnt_win+cnt_lose)
    arr_rate.append(cnt_win/(cnt_win+cnt_lose))

fig, axs = plt.subplots(ncols=2, figsize=(15, 5))

axs[0].bar(x=arr_target_player, height=arr_cnt, color='b', alpha=0.5)
axs[0].set(xlabel='player', ylabel='cnt')
for tick in axs[0].get_xticklabels():
    tick.set_rotation(75)

axs[1].bar(x=arr_target_player, height=arr_rate, color='r', alpha=0.5)
axs[1].set(xlabel='player', ylabel='rate')
for tick in axs[1].get_xticklabels():
    tick.set_rotation(75)

plt.show()

（敬称略）

やはりおなじみのBIG4（フェデラー、ナダル、ジョコビッチ、マレー）は外せません。

ここに錦織やデル・ポトロ、ラオニッチなどの中堅、ベルディヒ、ツォンガ、フェレールなど、昔からTOP10辺りに位置する選手、そして最近活躍してきているティエム、キリオス、ズベレフなどの若手選手を入れてみました。

上記の棒グラフが、青が対象期間中の各選手の試合数（vs対象選手のみ）、赤が勝率になります。

ベイズモデリングといえども、データ数が少なすぎれば収束しない恐れがあり、デル・ポトロなどが若干データが少なめですが、収束しなかった時に外すとして、とりあえずこれで続行します。

上記の統計モデルに食わせるデータ形式を用意します。

dic_target_player = {}
for player in arr_target_player:
    if player not in dic_target_player:
        dic_target_player[player] = len(dic_target_player)+1

LW = []
for player_a in arr_target_player:
    for player_b in arr_target_player:
        df_tmp = df_matches[
            (df_matches['year'] >= 2015) &
            (df_matches['winner_name'] == player_a) &
            (df_matches['loser_name'] == player_b)
        ]

        for _ in range(len(df_tmp)):
            LW.append([dic_target_player[player_b], dic_target_player[player_a]])

        df_tmp = df_matches[
            (df_matches['year'] >= 2015) &
            (df_matches['winner_name'] == player_b) &
            (df_matches['loser_name'] == player_a)
        ]

        for _ in range(len(df_tmp)):
            LW.append([dic_target_player[player_a], dic_target_player[player_b]])

LW = np.array(LW, dtype=np.int32)

上記の LW は、「負け選手インデックス、勝ち選手インデックス」という順に、試合ごとに行を格納した２次元配列です。

これが後にパフォーマンスの大小を調整するためのデータになります。

さて、上記の統計モデルをStanで書いてみます。（といっても棋士の強さのモデリングにコードも記載されていますので、写経）

model = """
    data {
        int N;
        int G;
        int LW[G, 2];
    }
    parameters {
        ordered[2] performance[G];
        vector[N] mu;
        real s_mu;
        vector[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);
    }
"""
fit1 = pystan.stan(model_code=model, data={'N': len(dic_target_player), 'G': len(LW), 'LW': LW}, iter=1000, chains=4)
la1 = fit1.extract()
fit1

上記コードの ordered[2] performance[G]; が今回のモデルで重要なところで、制約付きパラメータで、１列目よりも２列目の方が大きいという制約を与えています。

よって performance[g, i] ~ normal(mu[LW[g, i]], s_pf[LW[g, i]]); で、 LW の負けプレイヤーインデックスの方のパフォーマンスよりも、勝ちプレイヤーインデックスの方のパフォーマンスは大きくなるよう調整を行っているということになります。

ちなみに事前分布のガンマ分布は、プロットしてみれば分かりますが、正規分布などのような釣り鐘型の分布で、平均が1付近かつ正の値を取るような形状をしています。

正の値をとって、正規分布のような形をしてほしいという事前情報を与えたい時に、こういった分布を用いるようです。

処理を回してみたところ、Rhat がだいたいが 1.00 付近となりました。

念のために、各チェインごとのサンプリングをプロットしてみた結果が以下になります。

plt.figure(figsize=(15,7))
colors = ['red', 'yellow', 'green', 'blue']
for i, player in enumerate(arr_target_player):
    for j in range(4):
        g = plt.violinplot(la1['mu'][j*500:(j+1)*500, i], positions=[i], showmeans=False, showextrema=False, showmedians=False)
        for pc in g['bodies']:
            pc.set_facecolor(colors[j])

plt.legend(['chain 1', 'chain 2', 'chain 3', 'chain 4'])
plt.xticks(list(range(len(arr_target_player))), arr_target_player)
plt.xticks(rotation=45)
plt.xlabel('player')
plt.ylabel('mu')
plt.show()


plt.figure(figsize=(15,7))
colors = ['red', 'yellow', 'green', 'blue']
for i, player in enumerate(arr_target_player):
    for j in range(4):
        g = plt.violinplot(la1['s_pf'][j*500:(j+1)*500, i], positions=[i], showmeans=False, showextrema=False, showmedians=False)
        for pc in g['bodies']:
            pc.set_facecolor(colors[j])

plt.legend(['chain 1', 'chain 2', 'chain 3', 'chain 4'])
plt.xticks(list(range(len(arr_target_player))), arr_target_player)
plt.xticks(rotation=45)
plt.xlabel('player')
plt.ylabel('s_pf')
plt.show()

良い感じに重なって、ちゃんと収束していそうなので、ひとまずは問題なさそうです。

改めて統合したサンプリング結果をプロットしてみます。

まずは各選手の強さ \mu の事後分布が以下になります。

plt.figure(figsize=(15,7))
cmap = matplotlib.cm.get_cmap('tab10')
for i, player in enumerate(arr_target_player):
    g = plt.violinplot(la1['mu'][:, i], positions=[i], showmeans=False, showextrema=True, showmedians=True)
    c = cmap(i%10)
    for pc in g['bodies']:
        pc.set_facecolor(c)
    g['cbars'].set_edgecolor(c)
    g['cmaxes'].set_edgecolor(c)
    g['cmedians'].set_edgecolor(c)
    g['cmins'].set_edgecolor(c)

plt.xticks(list(range(len(arr_target_player))), arr_target_player)
plt.xticks(rotation=45)
plt.xlabel('player')
plt.ylabel('latent strength')
plt.show()


# カラーマップの作成コード
from matplotlib.colors import LinearSegmentedColormap

def generate_cmap(colors):
    values = range(len(colors))
    vmax = np.ceil(np.max(values))
    color_list = []

    for v, c in zip(values, colors):
        color_list.append(( v/ vmax, c))

    return LinearSegmentedColormap.from_list('custom_cmap', color_list)

cm = generate_cmap(['white', 'violet'])
summary = np.zeros((len(arr_target_player), 4))
for i, player in enumerate(arr_target_player):
    samples = la1['mu'][:, i]
    median = np.median(samples, axis=0)
    std = np.std(samples, ddof=1)
    lower, upper = np.percentile(samples, q=[25.0, 75.0], axis=0)
    summary[i] = [median, std, lower, upper]

summary = pd.DataFrame(summary, index=arr_target_player, columns=['median', 'std', '25%', '75%'])
summary.style.background_gradient(cmap=cm, axis=0)

ヴァイオリンプロットが各選手の強さのプロット、表はそれぞれ中央値、不偏標準偏差、50%ベイズ信頼区間の下限、上限となっています。

対象期間がジョコビッチ1強時代であったため、ジョコビッチが1番強いと出ており、直感通りの結果となっています。

デル・ポトロもBIG4に引けを取らない強さとなっていますが、怪我などの影響で勝負回数が少ないからか、帯も大きく伸びていて不確実度が高いことが分かります。

個人的に意外だったのが、ワウリンカと錦織が同じくらいの強さとなっているところで、全然ワウリンカの方が強いと思っていました...。

対象期間中の戦績を可視化してみた結果が以下になります。

cm = generate_cmap(['white', 'violet'])
df_tmp = df_matches[
    (df_matches['year'] >= 2015) &
    (df_matches['winner_name'].isin(arr_target_player)) &
    (df_matches['loser_name'].isin(arr_target_player))
]

summary = np.zeros((len(arr_target_player), len(arr_target_player)), dtype=np.int32)
for i in range(len(arr_target_player)):
    for j in range(i+1, len(arr_target_player)):

        summary[i, j] = len(df_tmp[
            (df_tmp['winner_name'] == arr_target_player[i]) &
            (df_tmp['loser_name'] == arr_target_player[j])
        ])

        summary[j, i] = len(df_tmp[
            (df_tmp['loser_name'] == arr_target_player[i]) &
            (df_tmp['winner_name'] == arr_target_player[j])
        ])

fig, axs = plt.subplots(figsize=(9, 9))
im = axs.imshow(summary, cmap=cm, interpolation='nearest')
axs.grid(False)
axs.set_xticks(list(range(len(arr_target_player))))
axs.set_yticks(list(range(len(arr_target_player))))
axs.set_xticklabels(arr_target_player)
axs.set_yticklabels(arr_target_player)
axs.set_ylabel('winner')
axs.set_xlabel('loser')
for tick in axs.get_xticklabels():
    tick.set_rotation(75)
for i in range(len(arr_target_player)):
    for j in range(len(arr_target_player)):
        text = axs.text(j, i, summary[i, j], ha='center', va='center', color='black')
fig.tight_layout()
plt.show()

一気に確認できるようにしたかったため、あまりこのような可視化を行うことはないと思いますが、行が勝ちプレイヤー、列が負けプレイヤーであった試合が何回あったかをマス目の数値で表しました。

これを見ても確かに、この期間内では、錦織はワウリンカに3勝2敗で勝ち越しており、ワウリンカ、錦織共にBIG4にも勝ったりもしているため、改めて錦織選手がTOP10の中でもかなり強い部類に位置していることを認識できました。（今は怪我でランキングを落としてしまっていますが...）

また、対象期間中では、若手は中堅以上には若干劣っている様子で、その若手の中では、ティエムが強いとなっているようです。

ここはやはり2017年後半や2018年のデータが欲しいところ。

続いて、各選手の勝負ムラ \sigma_{pf} の事後分布が以下になります。

plt.figure(figsize=(15,7))
cmap = matplotlib.cm.get_cmap('tab10')
for i, player in enumerate(arr_target_player):
    g = plt.violinplot(la1['s_pf'][:, i], positions=[i], showmeans=False, showextrema=True, showmedians=True)
    c = cmap(i%10)
    for pc in g['bodies']:
        pc.set_facecolor(c)
    g['cbars'].set_edgecolor(c)
    g['cmaxes'].set_edgecolor(c)
    g['cmedians'].set_edgecolor(c)
    g['cmins'].set_edgecolor(c)

plt.xticks(list(range(len(arr_target_player))), arr_target_player)
plt.xticks(rotation=45)
plt.xlabel('player')
plt.ylabel('uneven performance')
plt.show()


cm = generate_cmap(['white', 'violet'])
summary = np.zeros((len(arr_target_player), 4))
for i, player in enumerate(arr_target_player):
    samples = la1['s_pf'][:, i]
    median = np.median(samples, axis=0)
    std = np.std(samples, ddof=1)
    lower, upper = np.percentile(samples, q=[25.0, 75.0], axis=0)
    summary[i] = [median, std, lower, upper]

summary = pd.DataFrame(summary, index=arr_target_player, columns=['median', 'std', '25%', '75%'])
summary.style.background_gradient(cmap=cm, axis=0)

見たところ、ワウリンカが勝負ムラが一番大きいようです。

上記の戦績を見てみても、BIG4に限らず、同じ選手に勝ったり負けたりしていることがよくわかります。

BIG4の中では、この期間では、若干フェデラー、ナダルが勝ったり負けたりが多く、勝負ムラが少し大きめに出ました。

デル・ポトロは、マレー、ワウリンカに1勝1敗で、それ以外は安定して勝っていそうでしたが、なぜか大きめに出ました。

データ数が少なめだったからでしょうか、少し分かりません。

チリッチもこの期間は、ジョコビッチ、マレーに勝ったことがあり、その他選手とも勝ったり負けたりが多く、結果、勝負ムラが大きく出ています。

面白いなーと思ったのが、若手選手が全体的に勝負ムラは少し高めに出ており、これはまだパフォーマンスを安定して発揮することに慣れていないということなのかなーと思いました。

テニス選手の強さの時系列モデル

期間をある程度限定して分析を行いましたが、やはり時系列で各プレイヤーの強さがどのように推移してきているのかを見てみたいと思い、上記のモデルを少し拡張してみました。

    \[performance[y][g,1]{\sim}Normal(\mu[Loser][y], \sigma_{pf}[Loser][y]), \hspace{1em} g=1{\ldots}G, y=1{\ldots}Y\]

    \[performance[y][g,2]{\sim}Normal(\mu[Winner][y], \sigma_{pf}[Winner][y]), \hspace{1em} g=1{\ldots}G, y=1{\dots}Y\]

    \[performance[y][g,1] < performance[y][g,2], \hspace{1em} g=1{\ldots}G, y=1{\ldots}Y\]

    \[\mu[n][1]{\sim}Normal(0, \sigma_{\mu}[n][1]), \hspace{1em} n=1{\ldots}N\]

    \[\mu[n][y]{\sim}Normal(\mu[n][y-1], \sigma_\mu[n][y-1]), \hspace{1em} n=1{\ldots}N, y=2{\ldots}Y\]

    \[\sigma_{pf}[n][y]{\sim}Gamma(10, 10), \hspace{1em} n=1{\ldots}N, y=1{\ldots}Y\]

    \[\sigma_\mu[n][y]{\sim}Normal(0, 1), \hspace{1em} n=1{\ldots}N, y=1{\ldots}Y\]

無理やり上記のモデルに適当にねじ込んだ書き方でスミマセン。

  • Loser : 負けたプレイヤーのインデックス
  • Winner : 勝ったプレイヤーのインデックス

とします。

各プレイヤーの y 年度の強さを \mu[n][y], 勝負ムラを \sigma_{pf}[n][y] とし、 y 年度に行われる勝負で発揮する力（パフォーマンス）は、平均 \mu[n][y] , 標準偏差 \sigma_{pf}[n][y] の正規分布から生成されると考えます。

勝負の結果は同様に、パフォーマンスの大小で決まる（大きかった方が勝つ）とします。

各選手のある年度での強さ \mu[n][y] は、その1つ前の年度の強さ \mu[n][y-1] から生成されると考えます。

また、識別可能にするため、以下の仮定を入れます。

  1. 各選手の初年度の強さは平均 0 , 標準偏差 \sigma_{\mu}[n][1] の半正規分布に従うとする
  2. 各選手の年度別の強さの変化具合 \sigma_\mu[n][y] は半正規分布に従うとする

2は、弱情報事前分布です。

無情報事前分布でも試してみましたが、収束しませんでした。

といっても、次の年に別人のようにいきなり強くなるということはあまり無いと思いますので、この程度の仮定はあっても大丈夫かと思います。

対象期間は、2005〜2017.02頃（直近期間）の戦績のデータとします。

この期間だと、上記の若手選手などはほとんどデータが足りなくなりますので、改めて対象選手も手動で選択します。

実際に対象選手を絞り、データ数を確認してみます。

arr_target_year = np.array(list(range(2005, 2017)))

arr_target_player = np.array([
    'Roger Federer',
    'Rafael Nadal',
    'Novak Djokovic',
    'Andy Murray',
    'Stanislas Wawrinka',
    'Juan Martin Del Potro',
    #'Milos Raonic',
    'Kei Nishikori',
    #'Gael Monfils',
    'Tomas Berdych',
    #'Jo Wilfried Tsonga',
    'David Ferrer',
    #'Richard Gasquet',
    #'Marin Cilic',
    #'Grigor Dimitrov',
    #'Dominic Thiem',
    #'Nick Kyrgios',
    #'Alexander Zverev'
])

df_tmp = df_matches[
    (df_matches['winner_name'].isin(arr_target_player)) &
    (df_matches['loser_name'].isin(arr_target_player))
]

matrix_cnt = np.zeros((len(arr_target_year), len(arr_target_player)), dtype=np.float32)
matrix_rate = np.zeros((len(arr_target_year), len(arr_target_player)), dtype=np.float32)

for i, year in enumerate(arr_target_year):
    for j, player in enumerate(arr_target_player):
        cnt_win = len(df_tmp[(df_tmp['winner_name'] == player) & (df_tmp['year'] == year)])
        cnt_lose = len(df_tmp[(df_tmp['loser_name'] == player) & (df_tmp['year'] == year)])
        rate = 0 if (cnt_win+cnt_lose == 0) else cnt_win/(cnt_win+cnt_lose)
        matrix_cnt[i, j] = cnt_win+cnt_lose
        matrix_rate[i, j] = rate

for j, player in enumerate(arr_target_player):
    if j % 3 == 0:
        fig, axs = plt.subplots(ncols=3, figsize=(15, 3))
    axs[j%3].plot(arr_target_year, matrix_cnt[:, j], marker='o', color='b', alpha=0.5)
    axs[j%3].set(title=player, xlabel='year', ylabel='cnt', ylim=[0, 40])
plt.show()

for j, player in enumerate(arr_target_player):
    if j % 3 == 0:
        fig, axs = plt.subplots(ncols=3, figsize=(15, 3))
    axs[j%3].plot(arr_target_year, matrix_rate[:, j], marker='o', color='r', alpha=0.5)
    axs[j%3].set(title=player, xlabel='year', ylabel='rate', ylim=[0, 1])
plt.show()

青が年ごとのその選手の試合数（vs対象選手）、赤が年ごとの勝率です。

全ての年で5以上くらいありそうな選手と、やはり錦織選手を入れてみました。

入力するデータを作成してみます。

dic_target_year = {}
for year in arr_target_year:
    if year not in dic_target_year:
        dic_target_year[year] = len(dic_target_year)+1

dic_target_player = {}
for player in arr_target_player:
    if player not in dic_target_player:
        dic_target_player[player] = len(dic_target_player)+1

LW = []
GY = []
for year in arr_target_year:
    for player_a in arr_target_player:
        for player_b in arr_target_player:

            df_tmp = df_matches[
                (df_matches['year'] == year) &
                (df_matches['winner_name'] == player_a) &
                (df_matches['loser_name'] == player_b)
            ]
            for _ in range(len(df_tmp)):
                LW.append([dic_target_player[player_b], dic_target_player[player_a]])
                GY.append(dic_target_year[year])

            df_tmp = df_matches[
                (df_matches['year'] == year) &
                (df_matches['winner_name'] == player_b) &
                (df_matches['loser_name'] == player_a)
            ]
            for _ in range(len(df_tmp)):
                LW.append([dic_target_player[player_a], dic_target_player[player_b]])
                GY.append(dic_target_year[year])

LW = np.array(LW, dtype=np.int32)
GY = np.array(GY, dtype=np.int32)

上記は、LW は先程と同様で「負けプレイヤーインデックス、勝ちプレイヤーインデックス」が試合数分入っており、その試合番号に対応する年度はいくつかというリストとして GY を用意しています。

うまい人はもっと分かりやすく作るのでしょうけど、Stan力が足りなくてこうなってしまいました...。

Stanでモデルを実装すると、以下のようになります。

model = """
    data {
        int N;
        int G;
        int Y;
        int GY[G];
        int LW[G, 2];
    }
    parameters {
        ordered[2] performance[G];
        matrix[N, Y] mu;
        matrix[N, Y] s_mu;
        matrix[N, Y] s_pf;
    }
    model {
        for (g in 1:G)
            for (i in 1:2)
                performance[g, i] ~ normal(mu[LW[g, i], GY[g]], s_pf[LW[g, i], GY[g]]);
        for (n in 1:N)
            mu[n, 1] ~ normal(0, s_mu[n, 1]);
        for (n in 1:N)
            for (y in 2:Y)
                mu[n, y] ~ normal(mu[n, y-1], s_mu[n, y]);
        for (n in 1:N)
            s_mu[n] ~ normal(0, 1);
        for (n in 1:N)
            s_pf[n] ~ gamma(10, 10);
    }
"""
data = {
    'N': len(dic_target_player),
    'G': len(LW),
    'Y': len(dic_target_year),
    'GY': GY,
    'LW': LW,
}
fit2 = pystan.stan(model_code=model, data=data, iter=5000, chains=4)
la2 = fit2.extract()
fit2

収束判定の確認等は省略しますが、大丈夫そうでした。

結果について可視化してみます。

まずは各プレイヤーの強さの時系列の推移を確認してみます。

plt.figure(figsize=(15,7))
cmap = matplotlib.cm.get_cmap('tab10')
for j, player in enumerate(arr_target_player):
    samples = la2['mu'][:, j, :]
    medians = np.median(samples, axis=0)
    lower, upper = np.percentile(samples, q=[25.0, 75.0], axis=0)
    c = cmap(j)
    plt.plot(arr_target_year, medians, marker='o', label=player, color=c)
    plt.fill_between(arr_target_year, lower, upper, alpha=0.2, color=c)
plt.xlabel('year')
plt.ylabel('latent strength')
plt.legend(loc='lower left', bbox_to_anchor=(1, 0.5))
plt.show()

個人的には、割りと納得の推移結果となりました。

2005年頃は、やはりフェデラー、ナダルの２強時代であった様子がよくわかります。

そして、2008年頃から、フェデラー、ナダルに迫る若手、ジョコビッチ、マレー、デル・ポトロなどの強さが、やはり上位に推移してきていることがわかります。

錦織は2014年に強さが急激に上昇しており、この年に初めて世界ランクTOP10入りしました。

デル・ポトロは2009年（初めてグランドスラムに優勝した年）に急激に上昇するも、怪我の影響で推移を落としている様子が見えています。

マレーがかなりジグザグしていますが、2008年頃からフェデラー、ナダルに勝つようになり、この頃に世界ランク２位まで上り詰めるものの、2014年に怪我で不調となっていますので、そのような様子がはっきりと現れています。

ちなみに詳しいところはまだよく分かっていないのですが、この時系列では、年度内での強さの比較のため、年度を跨いだ強さの比較はできないのかなと思います。

例えば、2007年のフェデラーよりも2008年のナダルの方が強そうな結果は出ていますが、2007年はフェデラーの方が多く勝った、2008年ではナダルの方が多く勝ったという結果から導いているはずなので、直接そこを比較しているわけではないと思います。

続いて、各プレイヤー、各年ごとの勝負ムラの事後分布について、可視化してみます。

cmap = matplotlib.cm.get_cmap('tab10')
for j, player in enumerate(arr_target_player):
    if j % 3 == 0:
        fig, axs = plt.subplots(ncols=3, figsize=(15, 3))
    g = axs[j%3].violinplot(la2['s_pf'][:, j, :], positions=arr_target_year, showmeans=False, showextrema=False, showmedians=False)
    c = cmap(j%10)
    for pc in g['bodies']:
        pc.set_facecolor(c)
        pc.set_alpha(0.7)
    axs[j%3].set(title=player, xlabel='year', ylabel='uneven performance')
plt.show()

マレーが、年によって勝負ムラが大きくなるような結果となりました。

ワウリンカは、直近（2015〜2017.02）の分析では勝負ムラが大きかったものの、それまでの年で言えば、比較的大きくなかったようです。

ちょっと意外だったのが、錦織選手が、直近（2015〜2017.02）の分析と同様、各年においても勝負ムラはあまり大きくない傾向にあり、勝つ相手には勝つ、負ける相手には負けるがハッキリしているようです。

確かに対BIG4においても、ナダルには割りと勝つことが多い（ただしクレーは除く笑）印象があり、ジョコビッチには苦手意識を持っていそうな印象がありますが、どうなんでしょう。

最後に、各プレイヤー、各年ごとの成長の変化の度合いの事後分布を確認してみます。

cmap = matplotlib.cm.get_cmap('tab10')
for j, player in enumerate(arr_target_player):
    if j % 3 == 0:
        fig, axs = plt.subplots(ncols=3, figsize=(15, 3))
    g = axs[j%3].violinplot(la2['s_mu'][:, j, :], positions=arr_target_year, showmeans=False, showextrema=False, showmedians=False)
    c = cmap(j%10)
    for pc in g['bodies']:
        pc.set_facecolor(c)
        pc.set_alpha(0.7)
    axs[j%3].set(title=player, xlabel='year', ylabel='change of strength')
plt.show()

ここはやはり、基本的には、強さの時系列のモデルの傾きと、大きさが連動している様子が見られます。

初年度の強さを 0 から推定させていることもあり、元から強かったフェデラー、ナダルがいきなり強くなっているように見えていますので、この辺りは無視した方が良さそうです。

マレーは、フェデラー、ナダルに勝つようになった上向きの強さの変化と、不調による下向きの強さの変化があったことが見て取れます。

他も強さの時系列と同様、デル・ポトロは、2009年で強さが大きく変化している、錦織は、2014年で強さが大きく変化していることがわかります。

感想

かなり長い投稿になってしまいましたが、楽しかったです。

今回の分析だけでもまだやり残していることは多々あるように思います。

例えば、時系列の方では、成長の変化の大きさは年齢に比例して小さくなるような仮定を入れる方が良いのではないのかとか、点数の大きさを考慮に入れたらどうなるかなどです。

ただし、点数などを考慮した場合、伝統的なモデルなどを取り上げることになるのに対して、今回の方法は、上記に紹介した書籍でも述べられているように、順位という結果だけで、それらしい強さランクをモデル化する際に強力な方法だと思いました。

今回も入力したのは「勝ったか、負けたか」という情報のみですし、書籍では他にも、アンケートの人気投票などで回答者にTOP3を回答してもらった場合などにも使えると述べられています。

その場合は、LW が「ゲームごとの行」だったのを「回答者ごとの行」にして、制約付きパラメータを ordered[3] とするイメージかと思います。

その他にも、個人的には、競馬などのレース結果で、馬の強さをモデル化してみるとかも面白そうかと思いました。

１〜5着までの結果を入力して、 ordered[5] で分析すれば良いと思います。

スポンサーリンク

コメント