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

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

多変量正規分布とは
- 求めたいパラメータ
使用する関数
- mvnorm関数
- optim関数
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となりました。