[データ分析] 陸上の神はサイコロを振るのか？

<img width=500 alt="decathlon_1499492047_83444.jpg" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/b902d63d-b6da-3f97-1f76-d5f99ffb994a.jpeg">
出典　https://sportsmatik.com/sports-corner/sports-know-how/decathlon/rules

# 着想など
今回のデータ分析は、スポーツの様なやってみるまで結果が分からない様な現象が、実は何らかの法則に従っていたら面白いなと思って始めました。ランキングのパフォーマンスの頻度分布を説明する尤もらしいモデルもできたので、今後何回かに分けて、このモデルを検証していこうと思います。また、この記事は、WAのデータ収集コード（ https://github.com/kyaFUK/make-athletes-toplist ）を使って欲しくて書き始めたものなので、この記事が面白かったら、LGTMとgithubの方もよろしくお願いします。

# Abstract
今回の解析で、陸上競技の記録からシーズンランキングを予測できる数理モデルを設計しました。陸上競技の記録を確率変数$X$としたとき、シーズンランキングで観測される記録は以下の分布関数に従うことが分かりました。

```math
P(X \leq x) = 1 - \frac{1}{N}\exp\left((-ax+b)^{1/n}\right)
```
ただし、
* $N$: サンプル数
* $a,b, n$: 推定するパラメータ
とします。また、Track種目の記録は$(平均速度):=(距離)/(タイム)$と変換します。

次回以降では、この分布関数を極値理論を踏まえて分析します。

# データの詳細
WA(World Athletics、 国際陸連)のホームページに掲載されている、Season/All-time Top-listを用いました。
データ収集のためのコードも自作しました。最終的には、5万件を超えるデータを扱っているので、今後さらに解析を進める上では手作業でデータを集めることは困難であり、機械的にデータを収集するツールが必須だと考えられます。ソースコードは以下で公開しています。
https://github.com/kyaFUK/make-athletes-toplist

データの中身はシーズンランキングがそのまま入っていて、これの記録の部分を使います。
![2019-100m.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/6ddea47b-783a-d814-27a0-2a45921067e5.png)

# 計算上の工夫
今後、計算手法を拡張しやすい様に、いくつか工夫を行います。
1. Track種目の記録は、(距離)/(タイム)に変換します。
2. 累積分布関数$P(X \leq x)$の代わりに、$P(X > x) := 1 - P(X \leq x)$を計算します。
3. 上位20人は外れ値として除外します。

今回考える累積分布の解釈については[Appendix](#累積分布の概要)に書きました。工夫1を導入した理由は、Track種目とField種目によらず「記録の値は大きくなるほど、達成が難しくなる」という性質を統一したいためです。工夫2は色々な試行錯誤の結果です。$P(X > x)$の意味は、$x$より大きな$X$が出る確率なので、$x$より良い記録を出した選手の割合、即ち記録$x$を出した選手のランキングの相対順位（上位何%）に相当します。

### 例. 100mの記録の分布関数
まずは2019年男子100mの$P(X>x)$の分布を見てみましょう。
2019年男子100mのTop1000の記録をデータとして扱います。

|男子100mの平均速度の分布|
|----|
|![100m-speed-dist.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/f531772d-0e25-38e0-09d2-d3e2d552233a.png)|


ここで、$x$軸がspeedに変わっていることに注意してください。100mのSpeedは数字が大きいほど優れていて、達成される確率が小さい現象であることを表しています。

# 提案する数理モデル
今回の解析では、将来の世界記録を統計的に予想した[Einmahl et al. (2008)](#先行研究)の式(1)を参考に、数理モデルを設計しました。他に考慮した分布については[Appendix](#有名な分布)に書いています。ある種目のシーズンランキング上で観測される陸上競技の記録を確率変数$X$としたときに、$X$は以下の分布関数に従います。

```math
P(X>x) = \frac{1}{N}\exp\left((-ax+b)^{1/n}\right)
```
ただし、
* $N$: サンプル数
* $a,b, n$: 推定するパラメータ
とします。$a,b$は単回帰によって推定し、その時の決定係数が最大になるように$n$を定めました。

# 結果
### 2019年男子100mの分布
提案した数理モデルにデータをfittingして推定した係数を示します。2019年男子100mのTop1000の記録では以下のように求まりました。

|定数|値|
|----|----|
|$a$|38.554|
|$b$|392.482|
|$n$|1.69|
この時、決定係数は$r2=0.998811$となりました。グラフも見てみましょう。

|2019年男子100mTop1000の実データと数理モデルの分布|
|----|
|![100m-valid.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/ba0761cb-df4f-9047-3b1f-65e759d79dce.png)|

$x=10.2$付近の1つのデータ点を除いて、推定した分布関数によりシーズンランキングをトレースできていることが分かります。

### 2019年のその他の種目の分布
同様にして2019年男子の10種目について、Top1000の記録が提案したモデルに従うか調べました。

|   | 100m | LJ | SP | HJ | 400m | 110mH | DT | PV | JT | 1500m |
|----|----|----|----|----|----|----|----|----|----|----|
| a |38.554  |57.895  |5.125  |150.289  |40.813  |89.736  |1.283  |71.469  |2.657  |66.174  |
| b | 392.482  | 485.392  | 117.557  | 353.790  | 370.318  | 752.332  | 91.435  | 425.122  | 234.002  | 472.058  |
| n | 1.690  | 2.140  | 1.860  | 1.980  | 1.690  | 2.230  | 1.740  | 2.280  | 2.170  | 1.860  |
| r2 | 0.998811  | 0.999389  | 0.999243  | 0.994289  | 0.997983  | 0.999277  | 0.999049  | 0.996095  | 0.999463  | 0.996589  |
いずれの種目も100mの場合と同様に提案した数理モデルに従い、表の様に係数が推定されました。これを用いると、競技の記録からシーズンランキングを推定することが可能になります。

# 考察
提案した数理モデルを再掲します。

```math
P(X \leq x) = 1 - \frac{1}{N}\exp\left((-ax+b)^{1/n}\right)
```
ただし、
* $N$: サンプル数
* $a,b, n$: 推定するパラメータ
とします。また、Track種目の記録は$(平均速度):=(距離)/(タイム)$と変換します。

ここでは、このモデルがどの程度ユニバーサルなものなのか見ていきます。
1. 2019年のTop2000のランキングとTop3000のランキングの分布も表現できるか。
2. 2015年~2019年のTop1000のランキングを混合した5000人のランキングの分布も表現できるか。
3. U20やU18のランキングでも分布を表現できるか。

以上の3つの場合の分布を調べていきます。目的は数理モデルが以下の性質を持つことを示すためです。

1.上位層のみで成り立つモデルではないことを示す。
2.特定の期間で成り立つモデルではないことを示す。
3.特定の世代で成り立つモデルではないことを示す。

### 1. 2019年男子100mのTop2000とTop3000の分布
|   | Top1000| Top2000 | Top3000 |
|----|----|----|----|
| $a$ | 38.554 |63.033  |80.907   |
| $b$ |  392.482  | 638.662  | 817.887 |
| $n$ |  1.690 | 1.910  | 2.020 |
| $r2$ | 0.998811|0.999161 |0.999313 |
データ数が増えると、係数が変化しています。サンプル数が増えるほどに決定係数$r2$も増えているので、この数理モデルはスケーラビリティがあると言えます。

| Top2000 | Top3000 |
|----|----|
|![100-metres-2000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/d4528d6b-e78d-bcdb-fde2-1347f7130431.png)|![100-metres-3000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/66d754bc-dbe7-dc2a-762d-a1450e00de28.png)|

Top2000とTop3000について、fittingした時の結果です。今回は、$x$軸が秒になっていることに注意してください。データ数が増えても、同じ数理モデルで分布を表現できていることが分かります。


### 2. 2015~2019年男子100mのTop1000のデータを用いた時の分布
|   | 2019年| 2015~2019年 | 
|----|----|----|----|
| $a$ | 38.554 |170.412  |
| $b$ |  392.482  | 1737.091  | 
| $n$ |  1.690 | 2.240  | 
| $r2$ | 0.998811|0.999047 |

|2015~2019年Top1000の5000人分のデータにfittingした結果|
|----|
|![100-metres-5years.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/9e20673e-95f4-8c00-50a2-d65d0f37d5cd.png)|

ランキングのデータを数理モデルがトレースできていることが分かります。しかし、10.5秒付近でモデルと実際のデータがわずかに乖離している様にも見えます。この原因として、年を追うごとに世界全体の競技レベルが高まり、各年の記録の分布が変化していることが考えられます。
以下の表は男子100mの100位ごとの記録の変遷です。

rank|2001|2015	|2016|	2017	|2018	|2019
|----|----|----|----|----|----|----|
100	|10.26	|10.16	|10.14	|10.17	|10.15	|10.17
200	|10.35		|10.28	|10.25	|10.27	|10.25	|10.26
300	|10.42		|10.35	|10.32	|10.32	|10.32	|10.32
400	|10.48	|10.39	|10.37	|10.37	|10.37	|10.36
500	|10.53		|10.42	|10.41	|10.40	|10.41	|10.40

2015年と2019年を比較すると、0.01~0.02秒ほど全体的に競技レベルが向上していることが分かります。この変化により、わずかながらこのモデルでは表現できないデータが存在する可能性があります。

### 3. U20U18の2019年男子100mのTop1000の分布
U18は公開データが少なかったため、サンプル数は600人で計算しました。

|   | Senior| U20 | U18 |
|----|----|----|----|
| $a$ | 38.554 |77.294    |68.298  |
| $b$ |  392.482  | 758.612 | 658.864 |
| $n$ |  1.690 |1.970  |1.970  |
| $r2$ | 0.998811|0.998799 | 0.998447  |

| U20| U18 |
|----|----|
|![100-metres-u20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/0617e32d-0d4c-0186-1371-508d99843c85.png)|![100-metres-u18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/b73bb617-8b55-25bf-9c39-183728f71c61.png)

U20U18の場合でも、上位の数データ点を除いてはモデルが実データをトレースできていることが分かります。


# まとめ
ランキング上の陸上競技の記録が従う分布関数として以下の数理モデルを提案しました。

```math
P(X \leq x) = 1 - \frac{1}{N}\exp\left((-ax+b)^{1/n}\right)
```
ただし、
* $N$: サンプル数
* $a,b, n$: 推定するパラメータ
とします。$a,b$は単回帰によって推定し、その時の決定係数が最大になるように$n$を定めました。

また、この分布関数を使うと陸上競技の記録からシーズンランキングを予想できます。さらに、この分布関数は、サンプル数、期間、世代によらず成り立つユニバーサルなモデルであることが示唆されました。

# 次回予告
Einmahl et al. (2008)を参考に、今回求めた分布関数を極値理論で表現できないか考察してみます（予定）。

# Appendix
### 累積分布関数の概要
確率論では、確率変数$X$の累積分布関数$F(x)$は以下の様に定義されています。

```math
   F(x):=P(X \leq x)
```
言い換えると、確率変数$X$を観測した時の値が、$x$以下になる確率ということです。1000人の100mのランキングに載っている記録が確率変数$X$である時の累積分布関数のグラフを見てみましょう。
<img width=500 alt="100m-cdf.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/06df2b7f-f8d7-5f95-87aa-957d5c68cdc7.png">
グラフから、上位20%に入っている選手は、10.25以下で走っていることが読み取れます。よって、このランキング上では上位20%の選手が$10.25$以下のタイムを出しているので、$x=10.25$の時の累積分布関数は、

```math
   F(10.5):=P(X \leq 10.25) = 0.2
```
となります。

### 有名な分布
![100m-speed-dist.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/f531772d-0e25-38e0-09d2-d3e2d552233a.png)
上のような分布を見た時、まず思いつくのはべき乗分布と指数分布です。
分布関数は以下のように表されます。$\alpha$は任意の実数で、$\lambda$は正の実数です。


|べき乗分布|指数分布|
|----|----|
|$F(x) = Cx^{\alpha}$|$F(x) = C e^{- \lambda x}$|
|![powerlow.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/3bba4b77-dfb2-2d52-99ef-0e106e948723.png)|![expdist.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/f656724f-11f2-ce78-bd36-4468a863b323.png)|

べき乗分布は両対数グラフで、指数分布は片対数グラフで分布が直線になるという性質があります。理由は、それぞれの分布関数に対して、`log`をとるとわかると思います。

|べき乗分布|指数分布|
|----|----|
|![logpower.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/f77571fc-6e94-8fd4-0f91-04dabce64e2e.png)|![logexp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/500bc914-8e7c-6a60-23cc-4b96ecbf219b.png)|

### 100mの分布に対し`log`をとる
$P(X>x)$のグラフを片対数、両対数で書いてみます。

|片対数|両対数|
|----|----|
|![100mylog.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/233c8937-1be6-bf1d-f415-83fa8fea3ea0.png)|![100mxylog.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613755/ad479ea1-c226-4315-081c-5ea9cb7e50b5.png)|
両グラフはどちらも直線とは言い難いので、100mの$P(X>x)$はべき乗分布でも指数分布でもないことが分かります。

## 先行研究
[1] John H. J. Einmahl and Jan R. Magnus (Dec., 2008),"Records in Athletics through Extreme-Value Theory" ,Journal of the American Statistical Association, 484,1382-1391
https://www.tandfonline.com/doi/abs/10.1198/016214508000000698
[1]では、極値理論を用いて記録の分布関数を導出し、将来達成される世界記録の上限と達成困難度を評価しました。