カテゴリー: statistics

PyMCによる項目反応理論の実装する練習

tricycle

項目反応理論（IRT: Item Response Theory）は，TOEFLなどのような試験の各設問の難易度，および回答者の能力を推測するためのモデルである．基本的な考え方は，回答者の能力と設問の難易度をパラメータとするロジスティック関数をある設問の正解確率と見なし，試験の回答状況は定義された正解確率をパラメータとするベルヌーイ分布に従う，というものである．

パラメータの推定方法は最尤推定法が一般的であるが，バージョン4にアップデートされたPyMCの練習をかねて，MCMCで項目反応理論のパラメータ推定をやってみた．

以下は，Gistに置いたPyMCによるコード例であるである．どうぞご笑覧ください（Github Gistコードへの直リンク）．

主成分分析，特異値分解および潜在的意味解析

tricycle

何度やっても忘れてしまう主成分分析（PCA）と特異値分解（SVD）の関係性．忘れてしまうのは分かっていないから．ということで，再度お勉強．

長らくあやふやにしてきた主成分分析の中心化問題．今回の勉強で「主成分分析では事前にデータを中心化する必要がある」ことがようやく理解した．

なぜ中心化が必要か？それは，主成分分析が重心を中心に元データを回転させることで，新たな軸を見つける操作だから．中心化をしないと，元データは原点を中心に回転されてしまう（注：中心化しなくても，主成分分析っぽい導出はできてしまうが，これではダメなのだ（参考：主成分分析 (2) – 主成分の導出と意味））．これでは主成分分析のコンセプトから逸脱してしまう（この理解に至るまで時間がかかった（参考：回転の中心 – 何を基準に考えるか））．

ちなみに，scikit learnのPCAは事前にデータを中心化しなくても，内部で勝手にデータを中心化してくれるようである．一方で，scikit learnのSVDは，データは自動的に中心化しない．SVDは単なる行列分解手法だから，主成分分析を目的としない特異値分解をしたいときに勝手に中心化されても困る．scikit learnでPCAメソッドを使わずに，SVDを使ってPCAを行うときは，対象となる行列を事前に中心化する必要があることに注意．

なお，文書-単語行列の次元圧縮手法である潜在的意味解析（LSI）もSVDを用いるが，LSIでは元データの中心化は行わない．理由は参考サイトにも書かれているように，文書-単語行列はスパースであるため，中心化すると計算量が膨大になること．また，文書ベクトルの類似度は角度が重要であることから，中心化するとその情報が保存されないことが理由だそう（参考：Latent Semantic Indexing and Data Centering）．参考サイトにも書かれていたが，文書行列はもともとスパースなので，各次元の平均値はそもそもゼロに近い．これら理由より，LSIでは中心化は行わないようである．

参考にしたサイトの一覧：

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

tricycle

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

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

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

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

tricycle

確率分布の期待値や分散に関する諸公式をしょっちゅう忘れるので，記憶に定着させるために諸公式を導出してみたメモ．以下の式は，確率変数\(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\)が正規分布に従うことは，別途証明が必要（正規分布の再生性）．

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

tricycleOne Comment

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

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

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

tricycle4 Comments

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

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

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

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

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