読者です 読者をやめる 読者になる 読者になる

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

統計・機械学習

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


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