多変量正規分布の累積密度とか

以前の記事で多変量の正規分布に従う乱数を発生させるのにMASSというパッケージを紹介した。
任意の母相関を持つ多変量の生成 - 猫も杓子も構造化

多変量の正規分布を扱えるパッケージは他にもあって、mvtnormというパッケージがある。こっちのパッケージのrmvnormという関数でも乱数を発生できる。引数も平均ベクトルと分散共分散行列の指定でOKである。

library(mvtnorm)

mu <- c(0,0)
sigma <- matrix(c(1,0.9,0.9,1), nrow = 2, ncol = 2)

dat <- rmvnorm(1000, mean = mu, sigma = sigma)

このパッケージは、乱数の生成に加えて、任意の区間の累積密度を求めたり確率点を求めたりすることができる。例えば、次のような散布図で色分けした区間のそれぞれの累積密度（確率）を求めたいとする。（正確には２変量の確率密度関数の累積密度なので、散布図で考えるのは正しくないような気もしますが、３Dグラフは見辛くて好きじゃないので、、、散布図の上に正規分布の山が乗っているイメージで読み進めください）

パッケージ内のpmvnormという関数を用いれば計算してくれる。この関数は引数に、n変量の下限と上限をベクトルで指定すると、その範囲の累積密度を求めてくれる。例えば、上記の散布図の青い部分の累積密度を求めたければ次のように指定すれば良い。

blue.l <- c(-Inf,-Inf)
blue.u <- c(0,0)

pmvnorm(lower = blue.l, upper = blue.u, mean = mu, sigma = sigma)

そうすると次のような結果が返ってくる。

[1] 0.4282169
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

この一番上の値が、累積密度である。下の2つは、絶対誤差とステータスメッセージたるものを返してくれるらしいがあまり詳しく調べていないから、どういうときに"Normal Completion"じゃない値を返すのかはちょっとよく分からない。

実際に最初に乱数で発生させたデータ（n=1000）で青色部分のデータ数を数えて見ると次の通り。1000程度のデータ数だとそれなりに誤差があったりしますね。

> sum(dat[,1] <0 & dat[,2] <0)
[1] 415

ちなみに、デフォルトでは下限が-Infで上限がInfであるので、指定しないで実行すると累積密度は当然のことながら1になる。

> pmvnorm(mean = mu, sigma = sigma)
[1] 1
attr(,"error")
[1] 0
attr(,"msg")
[1] "Normal Completion"