ひらめの日常

日常のメモをつらつらと

トップ > 統計 > Rで2次元正規分布を描画する

Rで2次元正規分布を描画する

統計 R

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

github.com

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が多次元正規分布を描画するときには必要になります。
流れとしては、

  1. sigmaという分散共分散行列の定義
  2. rmvnormの引数に、平均と分散共分散行列を渡す
  3. 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")

今回は、dmvnormという関数を用います。

  1. x1, x2に軸用の値を代入
  2. dmvnormの第一引数に軸ラベル、そして平均と分散共分散行列を渡す
  3. dmvnormは渡された次元数を元に、多次元正規分布を生成
  4. 多次元正規分布の出力をzに代入しperspでグラフ描画

こんな感じになります。
f:id:thescript1210:20180425000954p:plain

参考

こちらのサイトなどを参考にさせていただきました。
http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html