多変量正規分布の密度関数を求めようと思うと、
R の mvtnorm パッケージの dmvnorm 関数が有名ですが、
少々計算に時間がかかるので、
繰り返しこの関数を適用する場合、
ボトルネックになってしまうことがあります。
そこで、解決方法を探していると、
計算を高速化したライブラリが提供されているようなので、
そちらを使ってみることにしました。
この記事では、その高速化版の関数と、
元の dmvnorm の計算時間を簡単に比較してみることにします。
高速化版の多変量正規分布のライブラリ
dmvnorm の計算時間の問題を解決してくれるのは、
その名の通り mvnfast というパッケージです。
このパッケージでは、
高速化された多変量正規分布の計算を行うための関数が、
まとめられています。
インストール
CRAN で提供されているので、
通常のパッケージインストールと同様に、
install.packages("mvnfast")
でインストールできます。
また、パッケージの読み込みは、
library(mvnfast)
でオッケーです。
密度関数
密度関数は、dmvnorm と同様に、
確率密度を計算する点xと、
多変量正規分布の平均値ベクトルmu、
分散共分散行列sigmaを与えて、
dmvn(X,mu,sigma)
で計算します。
計算時間の比較
それでは、実際に計算時間が速くなるか確かめてみましょう。
1回の計算時間は、dmvnorm でも十分に高速なので、
10000回の計算時間で比較します。
dmvnorm の計算時間
まずはお馴染みの dmvnorm の計算時間です。
system.time(sapply(1:10000,function(i){dmvnorm(1:10,1:10,diag(1:10))}))
結果は、1.37 秒でした。
dmvn の計算時間
続いて、mvnfast パッケージの dmvn による計算時間です。
system.time(sapply(1:10000,function(i){dmvn(1:10,1:10,diag(1:10))}))
結果は、0.158 秒でした。
まさに、けた違いに速くなっていることが分かります。
まとめ
多変量正規分布の密度関数の計算を高速化してくれる
パッケージを紹介しました。
プログラムの高速化を行っている時に、
ボトルネックが dmvnorm になっていたので、
どうしようかと思っていたところ、
かなり高速化できそうです。
