多変量分布に従う乱数の生成

統計・機械学習

１．多変量正規分布に従う乱数の生成
　パラメタは期待値ベクトルΛと共分散行列∑の２つ．
作り方のイメージは，一旦N(0,1)に従う独立な標準正規分布の乱数を生成して，
それを，線形変換して期待値ベクトルΛと共分散行列∑に従う乱数に変換する．
そして，そんな線形変換は？というと，共分散行列が対称行列であるから求める事が
容易にできる．具体的な手順は下記．Rの多変量正規乱数もこれで作っているみたい．
　　
　１．標準正規分布乱数N(0,I)に従う独立なN次元のベクトルY1を生成
　２．Y2 =Q・Y1でY1を変換する時にY2がN(0,∑)に従うにはどんな変換Qを適用するか？
　　　そんなQとは，∑=PAP_1として∑を対角化するときの固有値のベクトル（A)と
　　　対角行列Pとした時に，Q=PA_1/2 (A_1/2はAの各固有値の平方根をとったベクトル)
　　　とすれば良い．証明はネット探すとある．
　３．変換後のY2に対して，Y3=Y2+Λは所望の期待値ベクトルΛと共分散行列∑に従う
　　　乱数となっている．

　実際に作る場合，Rにはすでにパッケージがあるのでそれを使うのが早くて安心．例は下記．

library(MASS)
argv=commandArgs(T)
lambda.vec = scan("lambda_vec.dat")          # 期待値ベクトルΛ
corr.mat = as.matrix(read.table("corr.mat")) # 共分散行列∑
# samples個の多変量正規乱数を生成
random.numbers = mvrnorm(samples, lambda.vec, corr.mat)


２．多変量ポアソン分布に従う乱数の生成
　http://arxiv.org/pdf/0710.5670.pdfを参考．（上の多変量正規乱数の式も書いてある）
　生成方法の概要は，多変量正規乱数を生成して，それを逆関数法で変換する．