The One with ...

研究用のメモなどを公開しています

多次元正規分布のパラメータ推定

今日はStanを使って多次元正規分布のパラメータを推定する練習をやってみましょう(授業の続き).

まず推定用のデータをRで生成します.

library(MASS) # for mvrnorm( )
N0 <- 300 #データ個数
m1 <- 0;m2 <- 0;s1 <- 1;s2 <- 1;r <- 0.6 #パラメータ設定
mu = c(m1,m2)  #平均ベクトルの定義
sigma = matrix(c(s1^2, r*s1*s2, r*s1*s2, s2^2), 2, 2)#共分散行列の定義
#r=(r*s1*s2)/(s1*s2) 共分散行列は基準化しているので,rは相関
data0 = mvrnorm(N0, mu, sigma)
#2次元正規乱数をN0個を生成  相関cor(exv[,1],exv[,2]) 

2次元の正規分布にしたがう乱数を発生させる関数はmvrnorm(N0, mu, sigma) です(MASSパッケージの関数).引数のN0が個数,mu, sigmaは平均ベクトルと共分散行列です. 相関rをダイレクトに指定するため基準化しています.

この平均ベクトルと共分散行列がStanで推定すべきパラメータです.正しく推定できれば,自分が定義した 真のパラメータに近い値が出てくるはずです.

次の図はmvrnorm( )を使って生成したデータの散布図です.相関が0.6なので,右方向に点が上がっていく様子が見て取れます.多次元正規分布を使えば, このように指定した相関を持つデータ列を簡単に作ることができます.

f:id:hamada7418:20170712163246p:plain

コード

R上の分析実行コードですStanコードは例によってアヒル本記載のコードを流用しています.

Stanコードの中で注目すべきは,共分散行列用の型であるcov_matrixという型を使っている点です. 具体的には,Stanのパラメータブロックが次のように書かれています.

parameters {
  vector[D] mn;
  cov_matrix[D] cov;
}

この型は,対称な正定値行列に対応しています.

Stanにはそのほかにも単体用のsimplexや, 順序を持ったベクトルに対するorderedなどの型も用意されています.

さて推定結果fitの中身を確認してみましょう.

            mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
mn[1]       0.03    0.00 0.06   -0.09   -0.01    0.03    0.07    0.14  2853    1
mn[2]      -0.03    0.00 0.06   -0.14   -0.07   -0.03    0.01    0.08  2660    1
cov[1,1]    1.05    0.00 0.08    0.90    0.99    1.04    1.10    1.22  2852    1
cov[1,2]    0.59    0.00 0.07    0.47    0.55    0.59    0.63    0.73  2640    1
cov[2,1]    0.59    0.00 0.07    0.47    0.55    0.59    0.63    0.73  2640    1
cov[2,2]    0.98    0.00 0.08    0.84    0.93    0.98    1.03    1.14  3000    1
lp__     -236.19    0.04 1.53 -239.99 -236.92 -235.88 -235.09 -234.19  1367    1

Samples were drawn using NUTS(diag_e) at Wed Jul 12 16:05:24 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

おおむね,うまく推定できています*1

plot(fit)で95%ベイズ信頼区間を確認しましょう.

f:id:hamada7418:20170712173334p:plain

stan_hist(fit)でパラメータのヒストグラムを描きます.

f:id:hamada7418:20170712173337p:plain

パラメータの真値( m1 <- 0;m2 <- 0;s1 <- 1;s2 <- 1;r <- 0.6 の部分)をいろいろ変えながら 自分で試してみましょう.

解析的な事後分布の導出

さてここからは事前共役分布を使った事後分布の導出について話します.まず1次元の場合,{ \displaystyle  \mu   }に対する事前分布を{ \displaystyle
\mu \sim N( \mu_{0}, \sigma_{0}^{2}) }とすると,事後分布のカーネル

{  \displaystyle p ( \mu | y ) \propto \text{exp} \left\{ - \frac{\mu - \mu^{*})}{2 \sigma^{*2}} \right\} }

となる.ここで

{  \displaystyle \mu^{*}=w \mu_{0}+(1-w) \bar{y}, \sigma^{*2}=\frac{1}{1/\sigma_{0}^{2} + n/ \sigma^{2}} } である.

またウェイトは

{ \displaystyle w = \frac{1 / \sigma_{0}^{2} }{1 / \sigma_{0}^{2} + n/ \sigma^{2}} }

である.つまり \mu の事後分布は正規分布となる.

(続く.はてぶろのTeX表記,エスケープが多すぎて安定しない…)

*1:当たり前やん,と思うかもしれませんが,「正解」が存在する推定の問題を自分で作っているのだと 考えてください.もちろん経験的には意味はありません.しかし分析手法を理解するために,これは最適な方法の一つです