自分の実装はgithubにあります。

D次元正規分布

1次元正規分布はよく統計学とかでも出て来ますが、今回は2次元正規分布の確率密度関数を描画したいと思います。

定義

${ \displaystyle f(x) = \frac{1}{(2\pi)^{\frac{D}{2}}|\Sigma|^{\frac{1}{2}}}\exp(-\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu)) }$

${D}$ ：次元
${x}$ ：D次元ベクトル
${\Sigma}$ ：xの分散共分散行列
${\mu}$ ：D次元の分散ベクトル

以後D=2の場合を扱っていきます。

散布図

とりあえず確率密度関数を描画する前に、2次元正規分布に従う散布図をプロットしてみます。

library(mvtnorm)
sigma = matrix(c(2, 0.8, 0.8, 2), ncol = 2)
rand = rmvnorm(n = 1000, mean = c(0, 0), sigma)
x1 = rand[, 1]
x2 = rand[, 2]
plot(x1, x2)

mvtnormが多次元正規分布を描画するときには必要になります。
流れとしては、

sigmaという分散共分散行列の定義
rmvnormの引数に、平均と分散共分散行列を渡す
rmvnormは渡された次元数を元に、多次元正規分布に従った乱数を生成

f:id:thescript1210:20180425000159p:plain

確率密度関数

それでは、確率密度関数 ${f(x)}$ を描画していきます。

x1 = seq(-3, 3, length = 50)
x2 = x1
f = function(x1, x2) {
  dmvnorm(matrix(c(x1, x2), ncol = 2), mean = c(0, 0), sigma = sigma)
}
z = outer(x1, x2, f)
z[is.na(z)] = 1
op = par(bg = "white")
persp(x1, x2, z, theta = 30, phi = 30, expand = 0.5, col = "orange")