バナナでもわかる話

計量経済学・統計学が専門の大学院生です。統計学・経済学・投資理論・マーケティング等々に関する勉強・解説ブログ。ときどき趣味も。極力数式は使わずイメージで説明出来るよう心掛けていますが、時々暴走します。

Rで多変量正規分布の最尤推定量を求める

今回は多変量正規分布についてやっていこうと思います。

目次

スポンサーリンク



多変量正規分布とは

n次元確率変数ベクトル Yを考えます。この時、次のような密度 f(Y;\mu,\Sigma)を持つ時、 Yは多変量正規分布に従うと言います。

 f(Y;\mu , \Sigma) = \frac{1}{(2 \pi)^\frac{n}{2}\sqrt{ |\Sigma| }}exp(-\frac{1}{2}(Y-\mu)' \Sigma^{-1} (Y-\mu))

ただし、
 \muはn次元期待値ベクトル
 \Sigma n×n分散共分散行列
 'は転置記号
 \Sigma^{-1}は分散共分散行列の逆行列
 |\Sigma|は分散共分散行列の行列式

を表すものとします。

求めたいパラメータ

ここで、最尤法で考えたいパラメータですが、次の通りになります。
期待値 \mu=(\mu_1,\mu_2,\cdots,\mu_n) n通り
分散 \sigma=(\sigma_{11},\sigma_{22},\cdots,\sigma_{nn}) n通り
共分散 cov=(\sigma_{12},\cdots,\sigma_{1n},\sigma_{23},\cdots,\sigma_{2n},\cdots,\sigma_{(n-1)n}) \begin{eqnarray*}
  && {}_n C _2 \\
\end{eqnarray*}通り

流石多変量分布!パラメータが多い!
少し工夫しながら推定を行うコードを組んでいこうと思います。

使用する関数

mvnorm関数

R既存の関数には多変量正規分布を扱う関数はありませんが、mvtnormパッケージには多変量正規分布からの乱数を発生させるrmvnorm関数や、密度を与えるdmvnorm関数がセットされています。
今回は推定に使用するデータを次のようにrmvnorm関数で発生させてみます。

library(mvtnorm)
Dim=4
mu=rep(0,Dim)
sigma=diag(Dim)
Y=rmvnorm(1000,mu,sigma)

今回は n=4(Dim)に設定し、全ての期待値を0、全ての分散を1、全ての共分散を0で試しています。
また、対数尤度を考えるにあたって、dmvnorm関数を使用しています。

optim関数

最適化はoptimで行うことにします。optimは、まずparに対してパラメータの初期値を与えて、fnに推定パラメータを引数とする関数(今回は対数尤度)を与え、control=list(fnscale=-1)と設定することで最大化問題を探索的に求めてくれます。

Rで実装

n変量正規分布からのデータ xを入れると対数尤度を返す関数norm_paraを作成します。

norm_para = function(x){
	D=dim(x)[2]
	num=dim(x)[1]
	return(function(para){
		mu=para[1:D]
		var=para[(D+1):(D+D)]
		covnum=choose(D,2)
		cov=para[(2*D+1):(2*D+covnum)]
		sig_mat=diag(var)
		cumnum=c(0,cumsum((D-1):1))
		for(i in 1:(D-1)){
			sig_mat[i,i+(1:(D-i))]=sig_mat[i+(1:(D-i)),i]=cov[(cumnum[i]+1):cumnum[i+1]]
		}
		LL=0
		for(n in 1:num){
			LL=LL+log(dmvnorm(x[n,],mean=mu,sigma=sig_mat))
		}
		
		return(LL)}
	)
}

*1更に、適当な初期値を与えてoptimを実行してやります。データは先ほどmvnorm関数の節で作ったデータを利用します。

A=optim(par=c(rep(0.5,4),rep(1,4),rep(0,choose(4,2))),fn=norm_para(Y),control=list(fnscale=-1))
optim(par=A$par, fn=norm_para(Y), method="BFGS",control=list(fnscale=-1)) 


出力結果

$par
 [1] -0.016002375  0.049480996  0.048697919  0.022775051  0.962845786
 [6]  0.967682984  0.968557138  0.988855280  0.029744115  0.005733175
[11] -0.033611495  0.017052120  0.036690593 -0.001633487

$value
[1] -5616.856

$counts
function gradient 
     122       21 

$convergence
[1] 0

$message
NULL

 parのところに推定値が並んでいます。始めの n(=4)個が期待値、次の n(=4)個が分散、その次からは共分散の推定値が cov=(\sigma_{12},\cdots,\sigma_{1n},\sigma_{23},\cdots,\sigma_{2n},\cdots,\sigma_{(n-1)n})といった形で並んでいます。

よって、今回の各パラメータの推定値は
 \hat{\mu} =(-0.016002375 , 0.049480996 , 0.048697919 , 0.022775051)
 \hat{\sigma} =(0.962845786 , 0.967682984 , 0.968557138 , 0.988855280)
 \hat{cov} = (0.029744115 , 0.005733175 , -0.033611495 , 0.017052120 , 0.036690593 , -0.001633487)

真の値は
 \mu = (0,0,0,0)
 \sigma = (1,1,1,1)
 cov = (0,0,0,0,0,0)
だったので、中々近い値が出ていますね。今回はサンプルサイズ1000で試しましたが、悪くない値だと思います。

*1:id:abrahamcow さんからのご指摘を頂き、修正しました。出力結果convergenceが0になっていないと収束したとは言えないようです。恐らく初期値問題だと思われますので、一度初期値をoptimで与えた後、再度BFGSで最適化処理を行うという方法に変更したところconvergence=0となりました。