mvrnorm関数で「'Sigma' is not positive definite」というエラーが出たときの対処法
シミュレーションなどのために，多変量正規分布からのサンプリングを行う関数として{MASS}パッケージのmvrnorm関数があります。
この関数は次のようなプログラム例を通してサンプリングを行います。
library(MASS)
#3変量正規分布の例
#平均値として0が3つのベクトルを作る
mux=numeric(3)
#3*3でそれぞれ0.5の相関係数をとる相関行列
sigma<-matrix(0.5, 3, 3)
diag(sigma)<-1
#サンプルサイズ=100，3変数で多変量正規分布からのサンプリング
data<-mvrnorm(n=100, mu=mux, Sigma=sigma, empirical=FALSE)
empiricalというのは与えた相関行列を寸分たがわず再現するかしないか，というロジカルな引数で，TRUEだとcor(data)として相関係数を求めたときに与えた相関行列と完全に一致するため母相関行列のように使えるかと思います。
今回はFALSEとしているので，標本から推定した相関係数のように使うことができます。
さて，ここから本題ということで次のような状況を考えてみましょう。
x1-x10までの10個の変数があるときに，x1-x9までは相互に0.1で相関しており，x10だけは他の変数と0.5で相関しているという状況です。
これを前提として，10変量正規分布からのサンプリングを行うプログラムを次のように書きます。
library(MASS)
mux<-numeric(10)
sigma<-matrix(0.1, 10, 10)
sigma[10,]=0.5
sigma[,10]=0.5
diag(sigma)<-1
そうすると，
Error in mvrnorm(n = 100, mu = mux, Sigma = sigma, empirical = FALSE) :
'Sigma' is not positive definite
というエラーメッセージが返されて，データが生成されません。
この原因の確認と対処法をいかに述べようと思います。
エラーの原因の特定
メッセージを読めばわかるので言うまでもないのですけれど，相関行列の設定が不適切ということになります。
ただ，それで終わっては芸がないので，もう少し細かく見ていきます。
まず，mvrnorm関数の中身は次のようになっていることを確認します。
見やすさを重視して，スクリーンショットにしてみました。
この中で今回関心があるのは，
stop("'Sigma' is not positive definite")
という条件式までの部分になります。
かつ相関行列の指定までですので，次の部分が問題となっていると考えてよいでしょう。
eS <- eigen(Sigma, symmetric = TRUE)
ev <- eS$values
if (!all(ev >= -tol * abs(ev[1L])))
stop("'Sigma' is not positive definite")
これを，日本語に翻訳してコメントをつけてみます。
#与えた相関行列／共分散行列の固有値を求める
eS <- eigen(Sigma, symmetric = TRUE)
#固有値の計算結果（固有値，固有ベクトルが含まれる）の中から固有値だけを取り出す
ev <- eS$values
#もしすべての固有値が条件式を満たさなかった場合にはエラーメッセージを返す
if (!all(ev >= -tol * abs(ev[1L])))
stop("'Sigma' is not positive definite")
これらの内容から，今回のエラーは
if (!all(ev >= -tol * abs(ev[1L])))
という固有値に関する条件式を満たせなかったことが原因であることがわかります。
この条件式についてもうすこし細かく見てみます。
tolというのは上のスクリーンショットを見ていただければわかりますが，1e-06*1が割り当てられています。
abs関数は絶対値を返すものなので，abs(ev[1L])は1つ目の固有値の絶対値という意味になります。
また，!は否定演算子です。
したがって，この条件式は(-0.000001)*一つ目の固有値の絶対値を，得られた10個の固有値のすべてが越えなければFALSEを返すif条件式ということになります。
ためしに今回の状況における固有値の計算結果を見てみましょう。
> eS <- eigen(sigma, symmetric = TRUE)
> ev <- eS$values
> ev
[1] 2.9524175 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 -0.1524175
> tol=1e-06
> -tol * abs(ev[1L])
[1] -2.952417e-06
> ev[10]>=-tol * abs(ev[1L])
[1] FALSE
ひとつでもFALSEを返す場合にはエラーを返すようになっているので，かくして'Sigma' is not positive definiteというメッセージが表示されるわけです。
対策（力技）
さて，対策なのですが，要は与えた相関行列の組み合わせによっては固有値を計算する際に不適切なものが生じる恐れがある，ということなので，与える相関行列を最初から修正してしまえばいいじゃないかということになります。
そこで使うのが，{sfsmisc}パッケージのnearcor関数です。
nearcor {sfsmisc} | inside-R | A Community Site for R
これは不適切な相関行列を，cor関数でとりうる最も近い値にスムージングする関数とのことなので，素直に先ほどの相関行列をnearcor()の中に入れてからmvrnorm関数の相関行列として指定してあげれば，とりあえずエラーは消えるという力技的な解決法となっています。
お粗末様でした。
*1:=0.000001
