FKが一番上手いのは誰? 「過去5年間の成功率」で見る“ベスト・フリーキッカー”TOP20

という記事を見かけた。
FKが一番上手いのは誰? 「過去5年間の成功率」で見る“ベスト・フリーキッカー”TOP20

現在、最もフリーキックが上手い選手は誰なのか……。イギリス誌『FourFourTwo』が、「2013−14シーズンからの5年間で最もフリーキックの成功率が高かった選手TOP20」という記事を掲載した。
FKでのゴール数はよく伝えられるが、その成功率は影に隠れがちだ。
弱小チームの場合、ビッグクラブと比べて好位置でのFK自体が少なく、本当にすごいキッカーが埋もれてしまうこともある。また、優れたキッカーの多いチームに在籍している場合も、蹴るチャンスは少なくなってしまう。
一方、チームで常にキッカーを任されている選手は、蹴る機会が多い分、成功率を維持するのが難しくなってくる。
そういった条件も影響し、「過去5年間のFK成功率TOP20」には意外な選手の名前も登場した。
なお、データは1月23日時点のもので、5年間で最低20本以上のFK蹴った選手だけが対象となっている。また、『FourFourTwo』は意外にもランク外となった“名キッカー”も紹介している。

ということなので、28人分のフリーキック回数とゴール数が載っていたので、これをベイズ的に解析してフリーキック得点率を推定する。
単純に二項分布に当てはめるだけ。

特に何も考えない場合の成功率はこんな感じになる。色のついた部分は80% 信用区間であり、ベイズなので「このフリーキックで得点になる確率」と思ってもらってよい。
クリスティアーノロナウド(p[2]) やメッシ(p[3]) は蹴った回数が多いため、推定の精度があがり、幅は狭くなる。
ミラレム・ピアニッチ(p[17]) は76回蹴って11ゴールなので、他の選手と比べて少し幅が狭い。

　
さて、何も考えずにフリーキックを蹴って、入る確率が0% から100% まで一様分布、と期待するのは、サッカーをしていて非常に不自然である。
というわけで、得点確率p を適当に0% (ゴールが決まりにくい）に偏ったベータ分布からサンプリングすることにしよう。
いま、beta(1, 1) ならば一様分布のピンクだが、beta(1, 10) とすると得点確率0% から低いところに非常にかたよる。

　
この条件で同じようにすると、成功確率は過小評価気味になる。

二項分布はベータ関数 で表現することができて、これは

qbeta(alpha/2, dat$g+1, dat$n-dat$g+1)

で計算できる。ベータ関数のshape パラメータは -1 なので計算するときには +1 を忘れないこと。
このとき、事前分布が一様ならば、何も乗算しなくてよいが、事前分布をベータ分布で偏らせるならば、事前分布は であり、これを尤度関数（上のR コードのところ、つまり）にかけるだけで事後分布になる。つまり、 と のところに各々足す。
　
これ読んでおけばいいと思う（唐突な宣伝、COI なし）。

Rで楽しむベイズ統計入門[しくみから理解するベイズ推定の基礎] (Data Science Library)

　
期待としては、クリスティアーノロナウドやメッシといった5-10% くらいがせいぜい達成できるFK成功率だと思っていたが、頑張れば20% もいけるのだな、という印象。

library(rstan)
txt <- "name n g
ズラタン・イブラヒモヴィッチ 91 4
クリスティアーノ・ロナウド 112 7
リオネル・メッシ 165 11
ハカン・チャルハノール 99 10
ケヴィン・デ・ブライネ 35 4
アレクシス・サンチェス 42 5
ネイマール 41 5
ウィリアン 31 4
スソ 22 3
ガエル・カクタ 22 3
ウィサム・ベン・イェデル 22 3
ウェイン・ルーニー 29 4
ロメル・ルカク 21 3
ヤヤ・トゥーレ 28 4
ズラトコ・ユヌゾヴィッチ 35 5
ダヴィド・アラバ 42 6
ミラレム・ピアニッチ 76 11
マルクス・ズットナー 27 4
マウリシオ・レモス 20 3
アンリ・セヴェ 20 3
セヤド・サリホヴィッチ 20 3
マルヴィン・プラッテンハルト 37 6
フェデリコ・ヴィヴィアーニ 30 5
ハメス・ロドリゲス 36 6
マルコス・アロンソ 22 4
フィリペ・コウチーニョ 26 5
パウロ・ディバラ 38 8
フアン・マタ 22 5
"


dat <- read.table(text=txt, header=TRUE)
code <- "
data{
  int N_kicker;
  int N[N_kicker];
  int G[N_kicker];
  vector<lower=0>[2] s;
}
parameters{
  vector<lower=0, upper=1>[N_kicker] p;
}
model{
  p ~ beta(s[1], s[2]);
  G ~ binomial(N, p);
}
"

m <- stan_model(model_code=code)

cols <- c("hotpink", "skyblue")
input <- list(N=dat$n, G=dat$g, N_kicker=nrow(dat), s=c(1, 1))
fit <- sampling(m, data=input, seed=1234)

g <- stan_plot(fit, show_density=T, pars="p", outer_level=0.95, ci_level=0.8, fill_color=cols[1])
g + xlim(0, 0.4)

x <- seq(0, 1, length=1000)
d <- cbind(dbeta(x, 1, 1), dbeta(x, 1, 10))
matplot(x, d, type="l", lwd=4, lty=1, col=cols)

# サンプリングからの分布を考えるならば
ex <- extract(fit)
alpha <- 0.05
ci <- apply(ex$p, 2, quantile, c(alpha/2, 1-alpha/2))

# 一様な事前分布
qbeta(alpha/2, dat$g+1, dat$n-dat$g+1)
# ベータ分布
qbeta(alpha/2, dat$g+1+0, dat$n-dat$g+1+9)