Ty Cobbは四割打者の夢を見るか
統計学者のCarl MorrisによるTy Cobbの打率推定の例を確認します。
これまでと同様にやってみた記事の個人的備忘録ですのでご承知を。
1. はじめに
1-1. Ty Cobbってどんな人？
約100年前に活躍したMLBのスター選手
1905年から1928年までの24年間のキャリアでシーズン打率四割超えを3度達成し、1909年には打撃全タイトル制覇を達成する等々、数々の記録を保持するMLBを代表する打者の一人
参考：Wikipediaの記事
1-2. どんな例？
Ty Cobbが真の四割打者であったかどうかを推測する分析例となっている：
Ty Cobbは選手生涯を通して打率が高く、シーズン打率四割超を3度記録
しかし本当に四割超えを記録した打者と言えるのか？
選手の成績（打率）は投手や守備、球場等も含めたランダムな要素によって左右される。以前触れた18人の打者の打率推定の例でも見た通り、観測値（最尤推定値）が真の打率（という表現は個人的には好きではないが…）を表しているとは限らない
果たして彼は真の四割打者=真の打率が四割以上となったシーズンがあったと言えるのだろうか？
Ty Cobbが上記の意味で「真の四割打者であったか」を推測したい
各シーズンにおける真の打率を推定した上でそれが一度でも四割以上となる確率を計算する
有名な統計学者であるCarl Morrisの論文「Parametric Empirical Bayes Inference: Theory and Applications」が初出（論文のリンクは貼りません。興味のある方はご自身で探してください）。また、Morris氏が2023年に逝去された際にJim Albertがブログに追悼の解説記事を掲載。同じ例を現代的に書き直している。
原論文：正規近似と線形回帰を用いて構成→（当時はMCMCのような計算機による手法が発達していなかったため）理論的・古典的
解説記事：ロジスティック回帰を用いてシミュレーションを利用→（理論値を算出するのではなく計算機による解析という意味で）解析的・現代的
1-3. 前提知識
有名な18人の打者の打率推定の例（以下記事の参考文献等を参照）と関連している（というかその一般化になっている）ため、その知識が必要。
2. 方法
2-1. 問題点と解決方法
各シーズンの真の打率を推定するときに問題となるのはTy Cobbの年齢による成績変化である。18人の打率推定では全体平均を観測値と「混ぜる」ことで推定を行っていた。今回の例で全体平均の部分に対応するものとして真っ先に思いつくのはキャリアを通しての打率や各年の打率の平均値になるのだろうが、それでは彼の年齢による成績の変化を考慮できない。
そこで打率に対する年齢曲線を利用する。各年の打率（ロジスティック回帰を行うため、正確には打率のロジット）に対して2次曲線をフィッティングさせ、その推定値を各年の「平均値」として利用する。この意味で18人の打率推定の例におけるStein推定量の一般化になっている。
2次曲線のフィッティングによる年齢曲線の構成はよくあるものだと思うが、個人的にはけっこう強い仮定な気もする。例えば「ピークを基準として左右対称な成績推移」という前提はどれくらい妥当性があるのだろう。実際GAMか何かで構成する場合もあったはず。
ただまあそこまで変な仮定でもないと思うし、それが問題になるのなら別の方法で推定しろという話だろう。今回の分析では目安としての値が欲しいだけだと思って、この方法を受け入れるものとする。
2-2. データの生成過程（概要）
大雑把には以下のようなモデルを当てはめる：
Ty Cobbの$${j}$$年目のシーズンの打数を$${AB_{j}}$$, 安打数を$${H_{j}}$$とする。
そこで$${H_{j}}$$はサンプルサイズ$${AB_{j}}$$かつ確率$${p_{j}}$$の二項分布から生成されると仮定する：
$$
\begin{array}{}
H_{j} \sim Binomial(AB_{j}, p_{j}) .
\end{array}
$$
つまり真の打率を$${p_{j}}$$と考えるのである。
さらに$${p_{j}}$$は平均$${\eta_{j}}$$、精度(precision parameter)$${K}$$のベータ分布から生成されると仮定する：
$$
\begin{array}{}
p_{j} \sim Beta(K\eta_{j}, K(1-\eta_{j})) .
\end{array}
$$
ただし$${\eta_{j}}$$は2次曲線のフィッティングによる年齢曲線の値である：
$$
\begin{array}{}
logit(\eta_{j}) = \beta_{0} + \beta_{1}j + \beta_{2}j^{2}.
\end{array}
$$
$${\beta_{0}, ~\beta_{1}, ~ \beta_{2}}$$及び$${K}$$は以下の多変量正規分布から生成・変換する：
$$
\begin{array}{l}
\theta = (\theta_{0}, \theta_{1}, \theta_{2}, \theta_{3}) \sim MultiNormal(\mu, \Sigma), \\
\beta_{i} := \theta_{i} ~(i=0,1,2), ~~~K := \exp(\theta_{3})
\end{array}
$$
以上がモデルの設定であり、$${\mu}$$と$${\Sigma}$$を推定する事になる。実際には正規分布の分散共分散行列の推定にCholesky分解を使っていたり、パラメータの一部（stanコードのsigma_vec[4]を参照）に多少の制約を入れていたりする（制約を入れないと値が暴れる君になっちゃう、指数関数で変換してるので）ので、本当にこのままの推定を行なっているわけではないことに注意。詳細は後に添付するstanコード等を参照。
このモデルをもう少し見てみると、計算によって$${p_{j}}$$の推定値$${ \hat p_{j} }$$は以下の形となる事が分かる：
$$
\begin{array}{}
\hat p_{j} &=& \frac{H_{j} + K \eta_{j}}{AB_{j} + K} \\
&=& \frac{AB_{j}}{AB_{j} + K} H_{j} + (1-\frac{AB_{j}}{AB_{j} + K})\eta_{j} \\
&=& (1 - \frac{K}{AB_{j} + K}) H_{j} + \frac{K}{AB_{j} + K} \eta_{j} .
\end{array}
$$
これより、$${\hat p_{j}}$$は観測値$${H_{j}}$$を年齢曲線による「平均値」$${\eta_{j}}$$の方向に係数$${\frac{K}{AB_{j} + K}}$$で引っ張った形のStein推定量と考える事ができる。
2-3. 計算方法
今回はStanを用いたMCMCを実施する。以下のような違いがある。
原論文：正規近似を利用した理論的な導出式に値を代入
解説記事：ラプラス近似を利用したシミュレーションによる計算
ラプラス近似を利用しているのは、おそらく生成過程の最初の分散共分散行列の推定を楽にするため？
今回：Stanを用いたMCMCによるシミュレーション
この手のシミュレーションはStanが多いと思うのでやってみた
ラプラス近似は利用しないが一部パラメータに制約あり
分散共分散行列の推定にCholesky分解を利用（したつもり）
逆ガンマ分布や逆ウィシャート分布を利用すると頑健性に問題があったり、収束が難しいケースがあるらしいため
そのためStan公式ではLKJ相関分布を事前分布としたCholesky分解の利用を推奨しているらしい
この辺はStanの有名なアヒル本の10.2あたりに書いてあったのでそのまま参考にした…が書いたコードが合っているかは不安しかない
3. 結果
3-1. 結果
まずはシーズンごとの実際の打率、年齢曲線の値、真の打率の推定値を図示する：
横軸がMLB年数を、縦軸が打率を表しており、緑が実際の打率（観測値）、青が年齢曲線（2次曲線のフィッティングによる値）、赤が真の打率の推定値（上の$${\hat p_{j} }$$）を表している。また黒の太い水平線は打率四割のラインとなっている。
縮小推定になっているため緑と青の間に赤が存在しているが、最も打率が高かった7年目においても真の打率の推定値（赤）はそのラインを超えていない。実際の値は以下の通り：
最高でも.395であるため、四割以上という推定は得られなかった。
以下、分かっている人は読まなくてOK。この辺けっこう雑というか言葉を濫用する事も多いので一応注意しておこうと思ったが、あまり上手い説明にならなかったのでメモとして残しておく。このメモ内でも言葉の濫用が酷いのでいつか書き直したい…。
この結果だけからTy Cobbは真の四割打者ではないと思えるのだろうか？答えは当然No。理由は以下の通り。
上はあくまで真の打率の推定値=推定分布の代表値であることに注意。推定値もばらつきを持つ。つまり私たちが作成したモデルによる推定値は「真の打率はこの値だ！ブレはないのだ！」と主張しているのではない。むしろ「これくらいのブレを持つ」というべき。実際、2次曲線の係数（従ってそれから算出される$${\eta_{j}}$$）やベータ分布の精度パラメータ$${K}$$は分布（の実現値）として求められているため、それを利用して算出する推定値も分布を持つ。繰り返しのようになるが、上の図はあくまでその代表値を用いたプロットである。
不正確な例え話になるが、例えば私たちが「真の打率の推定値は.395である」と言ったときに、その値(.395)が100%正しいと思っているわけではない。.380も10%くらいあり得て、.390だと20%くらいあって、.395だと30%くらいになって、.400も15%くらいあって…みたいな中で、でも一番あり得そうなのは30%の .395って所だよね、と言っている。今気にしている「四割以上か」という問題について考えた時に、前述の「.400も15%くらいあって…」を無視して良いのですか？というお話。そこは入れましょうや、という。
よって、その分布（事後分布）を算出する必要があるのだが、事後分布を求めるのは一般には難しい。そのためにStanによるシミュレーションを実施している。
もう少し細かい事を言うと、今回は2種類のシミュレーションを実施しており、1つ目は上の通りStanによるMCMCを用いた「各年の推定分布を算出するためのシミュレーション」であり、2つ目は下記の通り「最高打率が四割以上かを確認するための（=推定結果を使って問題を解くための）シミュレーション」である。2つの意味を混同しないように注意。
さて、問題となる「一度でも真の打率が四割以上となったか？」という問いに答えるため、以下のシミュレーションを実施する：
乱数を用いて彼の24年間のキャリアを擬似的に5,000回シミュレートし、各回での「推定値の最大値=24年間での推定最高打率」を取得する。そうすると5,000回分の「推定最高打率」が取得できるので、それが四割以上となった比率を計算することで「推定最高打率が四割以上となった確率」を擬似的に算出する。
結果は以下の通り：
各回での「推定最高打率」のヒストグラムは上図の通り。垂直線は打率四割を表す。多くの場合で推定最大打率は四割以上となる事が分かる。計算してみると四割以上となった回数は5,000回中4,034回、比率としては約81%となった。
以上から、Ty Cobbは真の四割打者であったと言えるだろう。
3-2. 付録：どの年で稼いでいる？
どの年度の結果がどの程度、上の結果に貢献しているかを確認する。
上で出した打率推移のグラフに各年の「四割以上率」=「5,000回のシミュレーションの中で最高打率になったかつ四割以上となった割合」を棒グラフで追加した（値の軸は右側）。これらの値を全て足すと先ほどの約81%になる。
観測値（実際の打率結果、緑）が.400以上となった7,8年目は推定値（赤）も高くなっているため、「四割以上率」の値が高い（約18%,14%）。この2年間だけで81%のうち32%分を稼いでいる。
一方でこれまた実際の打率で四割以上を記録した18年目に関してはそこまで稼げていない（約7%）。図では分かりにくいが推定値（赤）の値もやや低くなっている。これは年齢曲線（青）の影響と思われる。つまり、今回のモデルでは観測値（実際の打率）と年齢曲線の値とを混ぜ込んでいるので、片方が低いとその影響が出て来る。その意味で、今回の結果は18年目の打率四割以上の蓋然性を低く評価していると言える。
一方でその間の期間については7,8年目ほどではないがコンスタントに「四割以上率」を稼いでいる（おおよそ各年5%前後）。これも年齢曲線の影響で推定打率（赤）はそこまで落ち込まなかったことが原因と思われる。ただし、例えば16年目に打率が下がった（それでも .334だが…）際には「四割以上率」も0.2%程度にまで下がっている点には注意。
期間ごとの数値も載せておく：
実際の打率が四割以上となった3シーズンだけで合計81%の約半分（約39%分）も稼いでいる。
3-2. 付録：実行コード
今回の実行コードはこちら。
stan：
Rmd：
stanファイルを適切な場所に置いてRmdの方を叩けば動くはず。
本当はnumpyroでやろうとしたが挫折した…。
4. 感想
思い付いた順に感想をだらだら書いてみる。
この例の直接的な応用ってそんなにないと思う。引退した選手のキャリアをレビューする時なんかには似たような方法が使えるかも。ただ考え方自体は有用だと思っていて、例えば18人の打者の例では事前分布のパラメータとして全体平均の値を使っていたけれど、今回はその部分が年齢曲線の値に変わっているので、この部分を（考慮したい変数を用いたモデルの推定値等に）取り替える事で状況に応じてより良い推定が得られるはず。まあそこまで行くと単に階層モデルでベイズ推定してるだけやんと言われそうだが。
これは過去の記事2で触れるべきだったかもしれないけれど、今回の結果をなぞれば18人の打率推定のStein推定量については実は等分散の仮定は不要である事が分かる（分散既知である必要はある）。ここが技術的には一番大切かも。あの面倒な分散安定化変換なんていらんかったんや！というお話。この辺はベイズモデリングの自由度や利便性って話なんだろうけれど、他にも色々と良い性質があるはず（多分…）。
こっちでバイアス入れるってのはやっぱ面白い考え方だと思った。こういうアイデアで他にも何かできないのかな。良い感じにバイアスを入れる事でNPBからMLBに来た選手の成績推定とか。逆でも良いけど。あとは信頼できるサンプルサイズの計算とか。Stein推定量と観測値がおおよそ同じ値になるのはどのくらい？みたいな。FGのとある研究だとCronbach’s Alphaを使っていたけれどそれと違いはあったりするんだろうか。
一方で3-2で見たように年齢曲線（これは事前分布のパラメータに対応するはず）に割と引っ張られている印象もあって、事前分布の設定=バイアスの入れ方次第でそこそこ変わりそうってのはちょっと注意かなとは思った…と言いつつやっぱり明示するだけマシかもとも。この辺の想いはあるんだけれども、どう言えば良いかイマイチ分からなくて上では書かなかった。
Jim AlbertによるCarl Morrisへのインタビュー記事（Chance誌のやつ）を引用していろいろと言いたかったんだけどいつの間にか無料で閲覧できなくなってる？だとしたら非常に残念…。Morrisがいろいろと語っていて面白かった記憶があるが$64を払おうとまでは思えない。たしか彼が（スポーツ関連の応用例を多く出している事に関して）良い例を作る事は統計学の発展に大きく寄与するけれどスポーツデータってけっこう良くて…みたいな事を語っていたような。ちなみに昔の私はそれを読んでとても興奮し、知人達にいろいろと溢れる想いを語った記憶がある（思い返すとクソ恥ずかしい）。
他にも新聞の広告欄で偶然見かけたBaseball Abstractの初版（！）を買った話とか、Bill Jamesについての話とか、Stanfordかどこかでけっこうなメンツを揃えてスポアナのセミナーやっていたけれどあの内容をまとめるべきだった…みたいな話もあったような。
この辺の話をしたくて書き始めた面もあったのでやる気出なくなっちゃった…。終わり。
以上、ご覧いただきありがとうございました。