日本プロ野球(NPB)のチーム戦闘力の統計モデリング

ここ1ヶ月くらいでNumPyroの練習をやってきた(記事1記事2).よく分からなかったnumpyro.plateの意味も,ドキュメントを読んでようやくちゃんと分かった(気がする).

そろそろNumPyroの練習に切りをつけるために,何かよい課題はないかということで,2次元データのソフトクラスタリング(混合ガウス分布)を扱おうと思ったのだが,離散値を潜在変数に設定したときにNumPyroではNUTSが使えないことが分かり断念.

このまま終わるのは後味が悪いので,NumPyro短期練習の最終回として「日本プロ野球のチーム戦闘力のモデリング」を題材にしてみることにした(GitHub Gistの元記事).同じ題材を扱った事例がSlideShareにあるが,モデリングの方法が異なるので,比べてみると面白いと思う.

確率分布の平均と分散に関する諸公式の導出メモ

確率分布の期待値や分散に関する諸公式をしょっちゅう忘れるので,記憶に定着させるために諸公式を導出してみたメモ.以下の式は,確率変数\(X\)と\(Y\)がどんな確率分布に従っているかに依らず導き出される性質である.

期待値(平均)と分散

確率変数\(X\)の確率密度関数は\(f(x)\)とする.

$$ E(X) =  \int_{-\infty}^{\infty} x f(x) dx $$

\begin{align}
V(X) &= \int_{-\infty}^{\infty} (x – \mu)^2 f(x) dx \\
&= \int_{-\infty}^{\infty} x^2 f(x) dx -2 \int_{-\infty}^{\infty} \mu_x x f(x) dx + \mu^2 \int_{-\infty}^{\infty} f(x) dx \\
&= E(X^2) -2 (E(X))^2 + (E(X))^2 \\
&= E(X^2) – (E(X))^2
\end{align}

諸公式の導出

確率変数\(X\)の確率密度関数を\(f(x)\),\(E(X) = \mu_x\) ,\(V(X) = \sigma_x^2\)とする.また,確率変数\(Y\)の確率密度関数を\(g(x)\),\(E(Y) = \mu_y\) ,\(V(Y) = \sigma_y^2\)とする.また,確率変数\(X\)と\(Y\)の同時確率密度関数を\(p(x, y)\),共分散\(Cov(X,Y)\)をとする.このとき,

\begin{align}
E(X+Y)&=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} (x+y) p(x, y) dx dy \\
&= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} x p(x, y) dx dy + \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} y p(x, y) dx dy \\
&= \int_{-\infty}^{\infty}  x p(x) dx + \int_{-\infty}^{\infty}  y p(y) dy \\
&= E(X) + E(Y)
\end{align}
\begin{align}
E(aX+b)&=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} (ax+b) f(x) dx \\
&= a \int_{-\infty}^{\infty} xf(x)dx + b \int_{-\infty}^{\infty} f(x)dx\\
&= aE(X) + b
\end{align}
\begin{align}
V(aX+b)&= E((aX+b)^2) – (E(aX+b))^2 \\
&= E(a^2X^2 + 2abX + b^2) – (aE(X)+b)^2 \\
&= a^2E(X^2) + 2abE(X) + b^2 – a^2E(X)^2 – 2abE(X) – b^2 \\
&= a^2(E(X^2) – E(X)^2) \\
& = a^2 V(X)
\end{align}
\begin{align}
Cov(X,Y)&= E((X-\mu_x)(Y-\mu_y)) \\
&= E(XY-\mu_xY-X\mu_y+\mu_x\mu_y) \\
&= E(XY)-E(\mu_x)E(Y)-E(X)E(\mu_y) +E(\mu_x)E(\mu_y)\\
& = E(XY) – E(X)E(Y)
\end{align}
\begin{align}
V(X+Y)&= E((X+Y)^2) – (E(X+Y))^2 \\
&= E(X^2 + 2XY + Y^2) – (E(X) + E(Y))^2 \\
&= E(X^2)+2E(XY)+E(Y^2) – (E(X)^2 + 2E(X)E(Y) + E(Y)^2) \\
& = (E(X^2)-E(X)^2) + (E(Y)^2 – E(Y)^2) + 2(E(XY) – E(X)E(Y)) \\
& = V(X) + V(Y) + 2Cov(X,Y)
\end{align}

※ 確率変数\(X\)が正規分布\(N(\mu, \sigma^2)\)に従うとき,確率変数\(aX+b\)は平均\(a\mu +b)\),分散\(a^2\sigma^2\)の正規分布に従う.しかし,\(E(aX+b)=aE(X)+b\)をもって\(aX+b\)が正規分布に従うとは言えない.確率変数\(aX+b\)が正規分布に従うことは,別途証明が必要(正規分布の再生性).

受信メッセージ数のモデリング

先日投稿した記事「NumPyroによる本塁打率のモデリング」に続き,NumPyroの試し書きを続けている.今回の題材は「Pythonで体験するベイズ推論」の一節から.

この書籍はPyMC2を用いて確率的プログラミングを行っているのだが,個人的に題材が面白くて気になっていた.PyMC2をやろうと思って2,3年が経過して,その間,TheanoベースのPyMC3がリリースされ,その後Theanoの開発終了を受け,TensorFlowベースのPyMC4の開発が開始されることになり… NumPyroに心移りし,PyMCを学ぶのはやめることになったが,確率的プログラミングの題材としては十分に役知立ちそうだ.

NumPyroによる本塁打率のモデリング

これまで,検索閲覧時間やページ閲覧回数などの行動データの統計モデリングを行うときは,メインのプログラミング言語であるPythonではなくRを使っていた.具体的にはbrmsというパッケージを利用していた.Rとbrmsを使えば,確率モデルを数式チックに記述できる.また結果を可視化し分析するパッケージ群が豊富にあるということで,Pythonを手放してまでRを使っていた.

ところが,Rだと計算が遅い.Rでは並列化やGPU利用に少々難がある.何より,使い慣れたPythonで前処理から分析を行いたい.これら問題を解決するためのツールとして,Pythonでは確率的プログラミングのライブラリとしてTensorFlow ProbabilityPyroが公開されている.しかし,これらはTensorの扱いに慣れていない軟弱な僕にはつらい.

諦めていたところ,最近NumPyroというライブラリがあることを知った.なんとTensorではなくnumpyのarray形式でデータを扱えるというではないか.開発もPyroと同じUber AI Labsによって行われているので,メンテナンスもしっかりしてそう.ということで,NumPyroの使い方を勉強してみることにした.

題材として,以前から気になっていた@muijpさんの「野球選手が本塁打を一番打てるのは何歳のときなのかPythonStanで求める」をNumPyroで実装してみることにした.

以下は,NumPyroでの統計モデリングの実装方法に関する記事である.どうぞご笑覧ください(Github Gistコードへの直リンク).