乱数生成のためのコレスキー分解・固有値が負になる相関行列の対処法について

R

コレスキー分解を使用するときに行列を転置するかしないかとか、行列積の順番を忘れて紙で計算しなおすということをしたあげく
毎回似たようなコードを書いていたのであほくさいなーと思ったり、それで乱数作ろうとしたときに
相関行列の固有値が負になってコレスキー分解できない場合の対処とかまとめる。

まずはコレスキー分解可能なケースで乱数生成

#適当な正値対称行列を用意（相関0.5、分散1)
x <- matrix(0.5,nrow=2,ncol=2)
diag(x) <- 1
#コレスキー分解
x.chol <- chol(x)
#元に戻る
t(x.chol) %*% x.chol
#相関0.5になる乱数の生成
z <- matrix(rnorm(100000),ncol=2) %*% x.chol
cov(z)
colMeans(z)

と、変数zが得られる。

次にこいつの固有値がマイナスになるケースを対処

#負の固有値が出てくるように相関が-0.9ずつの相関行列を作成
x <- matrix(-0.9,nrow=3,ncol=3)
diag(x) <- 1
#コレスキー分解したいけど・・・
x.chol <- chol(x)
#> 以下にエラー chol.default(x) :  次数 3 の主対角行列が正定値ではありません

とエラーが出る。これは負の固有値が存在するため、相関行列が正定値行列になっていないため。
これをadhocに解決するために、相関行列をスペクトル分解した後、負の固有値を10^(-6)程度の小さい値にする。

#行列のスペクトル分解
x.eigen <- eigen(x)
Lambda  <- diag(x.eigen$values)
P       <- x.eigen$vectors
#元通り
P %*% Lambda %*% t(P)
#Lambdaの固有値がマイナスの箇所を10^-6程度の小さい値に置き換え
Lambda.modified <- ifelse(Lambda < 0, 10e-6, Lambda)
#新しい相関行列を定義 & 対角成分を1に
x.modified <- P %*% Lambda.modified %*% t(P)
normalization.factor <- matrix(diag(x.modified),nrow = nrow(x.modified),ncol=1)^0.5
x.modified <- x.modified / (normalization.factor %*% t(normalization.factor))

#これでコレスキー分解できる
x.modified.chol <- chol(x.modified)
z <- matrix(rnorm(300000),ncol=3) %*% x.modified.chol
cov(z)
colMeans(z)

長々と書いたが、やっぱりこの手のもパッケージになってて
MatrixパッケージにあるnearPD関数が正定値行列近似をしてくれる。
http://rss.acs.unt.edu/Rdoc/library/Matrix/html/nearPD.html

Rjpwikiにも関連する記事があるが、nearPD関数使った方が安全かも。
http://www.okada.jp.org/RWiki/?%B9%D4%CE%F3Tips%C2%E7%C1%B4#n3ade4a5

他にも、要素間に望ましい制約（行列として持っていてほしい制約）を適当にぶち込んだうえで、
元の行列に近くなるようにフロベニウスノルムを最小化する方法もある。
この辺の話を教えてくれたアレさんに感謝。